Multinomial Mixed Effects Models

Description

The mme package implements three multinomial area level mixed effects models for small area estimation. The first model (Model 1) is based on the area level multinomial mixed model with independent random effects for each category of the response variable (Lopez-Vizcaino et al, 2013). The second model (Model 2) takes advantage from the availability of survey data from different time periods and uses a multinomial mixed model with independent random effects for each category of the response variable and with independent time and domain random effects. The third model (Model 3) is similar to the second one, but with correlated time random effects. To fit the models, we combine the penalized quasi-likelihood (PQL) method, introduced by Breslow and Clayton (1993) for estimating and predicting th fixed and random effects, with the residual maximum likelihood (REML) method for estimating the variance components. In all models the package use two approaches to estimate the mean square error (MSE), first through an analytical expression and second by bootstrap techniques.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

Breslow, N, Clayton, D (1993). Aproximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88, 9-25.

Inverse of the Fisher information matrix of the fixed and random effects in Model 1

Description

This function calculates the inverse of the Fisher information matrix of the fixed effects (beta) and the random effects (u) and the score vectors S.beta and S.u, for the model with one independent random effect in each category of the response variable (Model 1). modelfit1 uses the output of this function to estimate the fixed and random effects by the PQL method.

Usage

Fbetaf(sigmap, X, Z, phi, y, mu, u)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

See Also

data.mme, initial.values, wmatrix, phi.mult, prmu, phi.direct, sPhikf, ci, modelfit1, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata) #data
mod=1 #type of model
datar=data.mme(simdata,k,pp,mod)
initial=datar$initial
mean=prmu(datar$n,datar$Xk,initial$beta.0,initial$u.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

#Inverse of the Fisher information matrix
Fisher=Fbetaf(sigmap,datar$X,datar$Z,initial$phi.0,datar$y[,1:(k-1)],
       mean$mean,initial$u.0)

Inverse of the Fisher information matrix of fixed and random effects in Model 3

Description

This function calculates the score vector S and the inverse of the Fisher information matrix for the fixed (beta) and the random effects (u1, u2) in Model 3. This model has two independet sets of random effects. The first one contains independent random effects u1dk associated to each category and domain. The second set contains random effects u2dkt associated to each category, domain and time period. Model 3 assumes that the u2dk are AR(1) correlated across time. modelfit3 uses the output of this function to estimate the fixed and random effect by the PQL method.

Usage

Fbetaf.ct(sigmap, X, Z, phi1, phi2, y, mu, u1, u2, rho)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix, phi.mult.ct, prmu.time, phi.direct.ct, sPhikf.ct, ci, modelfit3, msef.ct, omega, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3)
datar=data.mme(simdata3,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities) #variance-covariance

##The inverse of the Fisher information matrix and the score matrix
Fisher.beta=Fbetaf.ct(sigmap,datar$X,datar$Z,initial$phi1.0,initial$phi2.0,
         datar$y[,1:(k-1)],mean$mean,initial$u1.0,initial$u2.0,initial$rho.0)

The inverse of the Fisher information matrix of the fixed and random effects for Model 2

Description

This function calculates the score vector S and the inverse of the Fisher information matrix for the fixed (beta) and the random effects (u1, u2) in Model 2. This model has two independet sets of random effects. The first one contains independent random effects u1dk associated to each category and domain. The second set contains random effects u2dkt associated to each category, domain and time period. Model 2 assumes that the u2dk are independent across time. modelfit2 uses the output of this function to estimate the fixed and random effect by the PQL method.

Usage

Fbetaf.it(sigmap, X, Z, phi1, phi2, y, mu, u1, u2)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, prmu.time, phi.direct.it, sPhikf.it, ci, modelfit2, msef.it, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2  #Type of model
data(simdata2) #data
datar=data.mme(simdata2,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

##The inverse of the Fisher information matrix of the fixed effects
Fisher=Fbetaf.it(sigmap,datar$X,datar$Z,initial$phi1.0,initial$phi2.0,
       datar$y[,1:(k-1)],mean$mean,initial$u1.0,initial$u2.0)

Add items from a list

Description

This function adds items from a list of dimension d*t, where d is the number of areas and t is the number of times periods.

Usage

addtolist(B_d, t, d)

Arguments

Value

B_d a list of dimension d.

See Also

Fbetaf.it, Fbetaf.ct, modelfit2, modelfit3

Examples

k=3  #number of categories for the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata2) # data
mod=2
datar=data.mme(simdata2,k,pp,mod)

##Add the time periods
l=addtolist(datar$X,datar$t,datar$d)

Add rows from a matrix

Description

This function adds rows from a matrix of dimension d*t*(k-1) times d*(k-1).

Usage

addtomatrix(C2, d, t, k)

Arguments

Value

C22 a matrix of dimension d*(k-1) times d*(k-1).

See Also

Fbetaf.it, Fbetaf.ct, modelfit2,modelfit3

Examples

k=3 #number of categories of the response variable
d=15 # number of areas
t=2  # number of time periods
mat=matrix(1,d*t*(k-1),d*(k-1)) # a matrix

##Add items in the matrix
mat2=addtomatrix(mat,d,t,k)

Standard deviation and p-values of the estimated model parameters

Description

This function calculates the standard deviations and the p-values of the estimated model parameters. The standard deviations are obtained from the asymptotic Fisher information matrix in the fitting algorithms modelfit1, modelfit2, modelfit3, depending of the current multinomial mixed model.

Usage

ci(a, F)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

modelfit1, modelfit2, modelfit3.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata) #data
mod=1  #Type of model
datar=data.mme(simdata,k,pp,mod)

#Model fit
result=modelfit1(pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N)
beta=result[[8]][,1] #fixed effects
Fisher=result[[3]] #Fisher information matrix

##Standard deviation and p-values
res=ci(beta,Fisher)

Function to generate matrices and the initial values

Description

Based on the input data, this function generates some matrices that are required in subsequent calculations and the initial values obtained from the function initial.values.

Usage

data.mme(fi, k, pp, mod)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

initial.values, wmatrix, phi.mult, prmu, Fbetaf, phi.direct, sPhikf, ci, modelfit1, msef, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata2) #Data
mod=2

##Needed matrix and initial values
datar=data.mme(simdata2,k,pp,mod)

Initial values for fitting algorithm to estimate the fixed and random effects and the variance components

Description

This function sets the initial values. An iterative algorithm fits the multinomial mixed models that requires initial values for the fixed effects, the random effects and the variance components. This initial values are used in modelfit1, modelfit2 and modelfit3.

Usage

initial.values(d, pp, datar, mod)

Arguments

Value

A list containing the following components, depending on the chosen model.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, wmatrix, phi.mult.it, prmu.time, Fbetaf.it, phi.direct.it, sPhikf.it, ci, modelfit2, msef.it, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)
D=nrow(simdata)
mod=1 #Type of model
datar=data.mme(simdata,k,pp,mod)

## Initial values for fixed, random effects and variance components
initial=datar$initial

Create objects of class mmedata

Description

This function creates objects of class mmedata.

Usage

mmedata(fi, k, pp)

Arguments

See Also

modelfit1, modelfit2, modelfit3

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)
r=mmedata(simdata,k,pp)

Choose between the three models

Description

This function chooses one of the three models.

Usage

model(d, t, pp, Xk, X, Z, initial, y, M, MM, mod)

Arguments

Value

the output of the function modelfit1, modelfit2 or modelfit3.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)  #data
mod=1 #Model 1
datar=data.mme(simdata,k,pp,mod)
result=model(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
datar$n,datar$N, mod)

Function used to fit Model 1

Description

This function fits the multinomial mixed model with one independent random effect per category of the response variable (Model 1), like in the formulation described in Lopez-Vizcaino et al. (2013). The fitting algorithm combines the penalized quasi-likelihood method (PQL) for estimating and predicting the fixed and random effects with the residual maximum likelihood method (REML) for estimating the variance components. This function uses as initial values the output of the function initial.values

Usage

modelfit1(pp, Xk, X, Z, initial, y, M, MM)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

data.mme, initial.values, wmatrix, phi.mult, prmu, phi.direct, sPhikf, ci, Fbetaf, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)  #data
mod=1 #type of model
datar=data.mme(simdata,k,pp,mod)

#Model fit
result=modelfit1(pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
      datar$n,datar$N)

Function to fit Model 2

Description

This function fits the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another independent time and domain random effect (Model 2). The formulation is described in Lopez-Vizcaino et al. (2013). The fitting algorithm combines the penalized quasi-likelihood method (PQL) for estimating and predicting the fixed and random effects, respectively, with the residual maximum likelihood method (REML) for estimating the variance components. This function uses as initial values the output of the function initial.values.

Usage

modelfit2(d, t, pp, Xk, X, Z, initial, y, M, MM)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, prmu.time, phi.direct.it, sPhikf.it, ci, Fbetaf.it, msef.it, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2 #type of model
data(simdata2)  #data
datar=data.mme(simdata2,k,pp,mod)

##Model fit
result=modelfit2(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N)

Function used to fit Model 3

Description

This function fits the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another correlated time and domain random effect (Model 3). The formulation is described in Lopez-Vizcaino et al. (2013). The fitting algorithm combine the penalized quasi-likelihood method (PQL) for estimating and predicting the fixed and random effects, respectively, with the residual maximun likelihood method (REML) for estimating the variance components. This function uses as initial values the output of the function initial.values.

Usage

modelfit3(d, t, pp, Xk, X, Z, initial, y, M, MM, b)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix, phi.mult.ct, prmu.time, phi.direct.ct, sPhikf.ct, omega, ci, Fbetaf.ct, msef.ct, mseb

Examples

## Not run: 

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3)   #data
datar=data.mme(simdata3,k,pp,mod)

##Model fit
result=modelfit3(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N,0)
## End(Not run)

Bias and MSE using parametric bootstrap

Description

This function calculates the bias and the mse for the multinomial mixed effects models using parametric bootstrap. Three types of multinomial mixed models are considered, with one independent domain random effect in each category of the response variable (Model 1), with two random effects: the first, with a domain random effect and with independent time and domain random effect (Model 2) and the second, with a domain random effect and with correlated time and domain random effect (Model 3). See details of the parametric bootstrap procedure in Gonzalez-Manteiga et al. (2008) and in Lopez-Vizcaino et al. (2013) for the adaptation to these three models. This function uses the output of modelfit1, modelfit2 or modelfit3, depending of the current multinomial mixed model.

Usage

mseb(pp, Xk, X, Z, M, MM, resul, B, mod)

Arguments

Value

a list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

Gonzalez-Manteiga, W, Lombardia, MJ, Molina, I, Morales, D, Santamaria, L (2008). Estimation of the mean squared error of predictors of small area linear parameters under a logistic mixed model, Computational Statistics and Data Analysis, 51, 2720-2733.

See Also

data.mme, initial.values, wmatrix, phi.mult, phi.mult.it, phi.mult.ct, prmu, prmu.time, phi.direct, phi.direct.it, phi.direct.ct, sPhikf, sPhikf.it, sPhikf.ct, modelfit1, modelfit2, modelfit3, omega, Fbetaf, Fbetaf.it, Fbetaf.ct, ci.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)
mod=1  # Type of model
datar=data.mme(simdata,k,pp,mod)
##Model fit
result=modelfit1(pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],datar$n,datar$N)

B=1    #Bootstrap iterations
ss=12345 #SEED
set.seed(ss)

##Bootstrap parametric BIAS and MSE
mse.pboot=mseb(pp,datar$Xk,datar$X,datar$Z,datar$n,datar$N,result,B,mod)

Analytic MSE for Model 1

Description

This function calculates the analytic MSE for the multinomial mixed model with one independent random effect per category of the response variable (Model 1). See Lopez-Vizcaino et al. (2013), section 4, for details. The formulas of Prasad and Rao (1990) are adapted to Model 1. This function uses the output of modelfit1.

Usage

msef(pp, X, Z, resul, MM, M)

Arguments

Value

mse is a matrix with the MSE estimator calculated by adapting the explicit formulas of Prasad and Rao (1990).

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

Prasad, NGN, Rao, JNK (1990).The estimation of the mean squared error of small area estimators. Journal of the American Statistical Association, 85, 163-171.

See Also

data.mme, initial.values, wmatrix, phi.mult, prmu, phi.direct, sPhikf, modelfit1, Fbetaf, ci, mseb.

Examples

require(Matrix)

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)  #data
mod=1 # type of model
datar=data.mme(simdata,k,pp,mod)
# Model fit
result=modelfit1(pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N)

#Analytic MSE
mse=msef(pp,datar$X,datar$Z,result,datar$N,datar$n)

Analytic MSE for Model 3

Description

This function calculates the analytic MSE for the multinomial mixed model with two independent random effects for each category of the response variable: one random effect associated with the domain and another correlated random effect associated with time and domain (Model 3). See details of the model and the expresion of mse in Lopez-Vizcaino et al. (2013). The formulas of Prasad and Rao (1990) are adapted to Model 3. This function uses the output of modelfit3.

Usage

msef.ct(p, X, result, M, MM)

Arguments

Value

mse.analitic is a matrix with the MSE estimator calculated by adapting the explicit formulas of Prasad and Rao (1990).

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

Prasad, NGN, Rao, JNK (1990).The estimation of the mean squared error of small area estimators. Journal of the American Statistical Association, 85, 163-171.

See Also

data.mme, initial.values, wmatrix, phi.mult.ct, prmu.time, phi.direct.ct, sPhikf.ct, modelfit3, Fbetaf.ct, ci, omega, mseb.

Examples

## Not run: 

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3) #data
datar=data.mme(simdata3,k,pp,mod)
##Model fit
result=modelfit3(d,t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N,0)

##Analytic MSE
msef=msef.ct(pp,datar$X,result,datar$n,datar$N)
## End(Not run)

Analytic MSE for Model 2

Description

This function calculates the analytic MSE for the multinomial mixed model with two independent random effects for each category of the response variable: one random effect associated with the domain and another independent random effect associated with time and domain (Model 2). See details of the model and the expresion of mse in Lopez-Vizcaino et al. (2013). The formulas of Prasad and Rao (1990) are adapted to Model 2. This function uses the output of modelfit2.

Usage

msef.it(p, X, result, M, MM)

Arguments

Value

mse.analitic is a matrix with the MSE estimator calculated by adapting the explicit formulas of Prasad and Rao (1990). The matrix dimension is the number of domains multiplied by the number of categories minus 1.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

Prasad, NGN, Rao, JNK (1990).The estimation of the mean squared error of small area estimators. Journal of the American Statistical Association, 85, 163-171.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, prmu.time, phi.direct.it, sPhikf.it, modelfit2, Fbetaf.it, ci, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2 #type of model
data(simdata2)
datar=data.mme(simdata2,k,pp,mod)
##Model fit
result=modelfit2(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N)

##Analytic MSE
msef=msef.it(pp,datar$X,result,datar$n,datar$N)

Model correlation matrix for Model 3

Description

This function calculates the model correlation matrix and the first derivative of the model correlation matrix for Model 3. Model 3 is the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another correlated time and domain random effect.

Usage

omega(t, k, rho, phi2)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix ,phi.mult.ct, prmu.time, phi.direct.ct, Fbetaf.ct, sPhikf.ct, ci, modelfit3, msef.ct, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3)   #data
datar=data.mme(simdata3,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)

##The model correlation matrix
matrix.corr=omega(datar$t,k,initial$rho.0,initial$phi2.0)

Variance components for Model 1

Description

This function calculates the variance components for the multinomial mixed model with one independent random effect in each category of the response variable (Model 1). These values are used in the second part of the fitting algorithm implemented in modelfit1. The algorithm adapts the ideas of Schall (1991) to a multivariate model and the variance components are estimated by the REML method.

Usage

phi.direct(sigmap, phi, X, u)

Arguments

Value

a list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult, prmu, Fbetaf, sPhikf, ci, modelfit1, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)   #data
mod=1 #type of model
datar=data.mme(simdata,k,pp,mod)
initial=datar$initial
mean=prmu(datar$n,datar$Xk,initial$beta.0,initial$u.0)
#model variance-covariance matrix
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

##Variance components
phi=phi.direct(sigmap,initial$phi.0,datar$X,initial$u.0)

Variance components for Model 3

Description

This function calculates the variance components for the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another correlated time and domain random effect (Model 3). This variance components are used in the second part of the fitting algorithm implemented in modelfit3. The algorithm adapts the ideas of Schall (1991) to a multivariate model. The variance components are estimated by the REML method.

Usage

phi.direct.ct(p, sigmap, X, theta, phi1, phi2, u1, u2, rho)

Arguments

Value

a list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submitted for review.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult.ct, prmu.time, Fbetaf.ct sPhikf.ct, ci, modelfit3, msef.ct, mseb, omega

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3) #data
datar=data.mme(simdata3,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

##The variance components
phi.ct=phi.direct.ct(pp,sigmap,datar$X,mean$eta,initial$phi1.0,
       initial$phi2.0,initial$u1.0,initial$u2.0,initial$rho.0)

Variance components for Model 2

Description

This function calculates the variance components for the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another independent time and domain random effect (Model 2). This variance components are used in the second part of the fitting algorithm implemented in modelfit2. The algorithm adapts the ideas of Schall (1991) to a multivariate model. The variance components are estimated by the REML method.

Usage

phi.direct.it(pp, sigmap, X, phi1, phi2, u1, u2)

Arguments

Value

a list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, prmu.time, Fbetaf.it sPhikf.it, ci, modelfit2, msef.it, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
d=10 #number of areas
mod=2  #Type of model
data(simdata2)   #data
datar=data.mme(simdata2,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities) #variance-covariance

## The variance components
phi.it=phi.direct.it(pp,sigmap,datar$X,initial$phi1.0,initial$phi2.0,
       initial$u1.0,initial$u2.0)

Initial values for the variance components for Model 1

Description

This function is used in initial.values to calculate the initial values for the variance components in the multinomial mixed model with one independent random effect in each category of the response variable (Model 1).

Usage

phi.mult(beta.0, y, Xk, M)

Arguments

Value

phi.0 vector of inicial values for the variance components

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

data.mme, initial.values, wmatrix, prmu, Fbetaf, phi.direct, sPhikf, ci, modelfit1, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata) #data
mod=1 #type of model
datar=data.mme(simdata,k,pp,mod)
###beta values
beta.new=list()
beta.new[[1]]=matrix(c( 1.3,-1),2,1)
beta.new[[2]]=matrix(c( -1.6,1),2,1)

##Initial variance components
phi=phi.mult(beta.new,datar$y,datar$Xk,datar$n)

Initial values for the variance components in Model 3

Description

This function is used in initial.values to calculate the initial values for the variance components in the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another correlated time and domain random effect (Model 3).

Usage

phi.mult.ct(beta.0, y, Xk, M, u1, u2)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

See Also

data.mme, initial.values, wmatrix, prmu.time, Fbetaf.ct, phi.direct.ct, sPhikf.ct, ci, modelfit3, msef.ct,omega, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3  #type of model
data(simdata3) #data
D=nrow(simdata3)
datar=data.mme(simdata3,k,pp,mod)
###Fixed effects values
beta.new=list()
beta.new[[1]]=matrix(c( 1.3,-1),2,1)
beta.new[[2]]=matrix(c( -1.6,1),2,1)
## Random effects values
u1.new=rep(0.01,((k-1)*datar$d))
dim(u1.new)=c(datar$d,k-1)
u2.new=rep(0.01,((k-1)*D))
dim(u2.new)=c(D,k-1)

## Initial variance components
phi=phi.mult.ct(beta.new,datar$y,datar$Xk,datar$n,u1.new,u2.new)

Initial values for the variance components in Model 2

Description

This function is used in initial.values to calculate the initial values for the variance components in the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect (u1) and another independent time and domain random effect (u2) (Model 2).

Usage

phi.mult.it(beta.0, y, Xk, M, u1, u2)

Arguments

Value

phi.0 vector of the initial values for the variance components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

data.mme, initial.values, wmatrix, prmu.time, Fbetaf.it, phi.direct.it, sPhikf.it, ci, modelfit2, msef.it, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata2)  #data
mod=2 #Type of model
datar=data.mme(simdata2,k,pp,mod)
D=nrow(simdata2)
###fixed effects values
beta.new=list()
beta.new[[1]]=matrix(c( 1.3,-1),2,1)
beta.new[[2]]=matrix(c( -1.6,1),2,1)
## random effects values
u1.new=rep(0.01,((k-1)*datar$d))
dim(u1.new)=c(datar$d,k-1)
u2.new=rep(0.01,((k-1)*D))
dim(u2.new)=c(D,k-1)

##Initial variance components
phi=phi.mult.it(beta.new,datar$y,datar$Xk,datar$n,u1.new,u2.new)

Print objects of class mme

Description

This function prints objects of class mme.

Usage

## S3 method for class 'mme'
print(x, ...)

Arguments

See Also

modelfit1, modelfit2, modelfit3

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=1  # Type of model
data(simdata)
datar=data.mme(simdata,k,pp,mod)
##Model fit
result=modelfit1(pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],datar$n,datar$N)
result

Estimated mean and probabilities for Model 1

Description

This function calculates the estimated probabilities and the estimated mean of the response variable, in the multinomial mixed model with one independent random effect in each category of the response variable (Model 1).

Usage

prmu(M, Xk, beta, u)

Arguments

Value

A list containing the following components:

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.

See Also

data.mme, initial.values, wmatrix, phi.mult, Fbetaf, phi.direct, sPhikf, ci, modelfit1, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata) #data
mod=1 #type of model
D=nrow(simdata)
datar=data.mme(simdata,k,pp,mod)
initial=datar$initial

##Estimated mean and probabilities
mean=prmu(datar$n,datar$Xk,initial$beta.0,initial$u.0)

Estimated mean and probabilities for Model 2 and 3

Description

This function calculates the estimated probabilities and the estimated mean of the response variable, in the multinomial mixed models with two independent random effects, one random effect associated with the area and the other associated with the time, for each category of the response variable. The first model assumes independent time and domain random effect (Model 2) and the second model assumes correlated time and domain random effect (Model 3).

Usage

prmu.time(M, Xk, beta, u1, u2)

Arguments

Value

A list containing the following components:

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicator under a multinomial mixed model with correlated time and area effects. Submited for review.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, Fbetaf.it, phi.direct.it, sPhikf.it, ci, modelfit2, msef.it, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2 #Type of model
data(simdata2) # data
datar=data.mme(simdata2,k,pp,mod)
initial=datar$initial

## Estimated mean and estimated probabilities
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)

Fisher information matrix and score vectors of the variance components for Model 1

Description

This function computes the Fisher information matrix and the score vectors of the variance components, for the multinomial mixed model with one independent random effect in each category of the response variable (Model 1). These values are used in the fitting algorithm implemented in modelfit1 to estimate the random effects. The algorithm adatps the ideas of Schall (1991) to a multivariate model. The variance components are estimated by the REML method.

Usage

sPhikf(pp, sigmap, X, eta, phi)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult, prmu, phi.direct, Fbetaf, ci, modelfit1, msef, mseb.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata) #data
mod=1 #type of model
datar=data.mme(simdata,k,pp, mod)
initial=datar$initial
mean=prmu(datar$n,datar$Xk,initial$beta.0,initial$u.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

##Fisher information matrix and score vectors
Fisher.phi=sPhikf(pp,sigmap,datar$X,mean$eta,initial$phi.0)

Fisher information matrix and score vectors of the variance components for Model 3

Description

This function computes the Fisher information matrix and the score vectors of the variance components, for the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another correlated time and domain random effect (Model 3). These values are used in the fitting algorithm implemented in modelfit3 to estimate the random effects. The algorithm adatps the ideas of Schall (1991) to a multivariate model. The variance components are estimated by the REML method.

Usage

sPhikf.ct(d, t, pp, sigmap, X, eta, phi1, phi2, rho, pr, M)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult.ct, prmu.time, phi.direct.ct, Fbetaf.ct, omega, ci, modelfit3, msef.ct, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3) #data
datar=data.mme(simdata3,k,pp, mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

## Fisher information matrix and the score vectors
Fisher.phi.ct=sPhikf.ct(datar$d,datar$t,pp,sigmap,datar$X,mean$eta,initial$phi1.0,
           initial$phi2.0,initial$rho.0,mean$estimated.probabilities,datar$n)

Fisher information matrix and score vectors of the variance components for Model 2

Description

This function computes the Fisher information matrix and the score vectors of the variance components, for the multinomial mixed model with two independent random effects for each category of the response variable: one domain random effect and another independent time and domain random effect (Model 2). These values are used in the fitting algorithm implemented in modelfit2 to estimate the random effects. The algorithm adatps the ideas of Schall (1991) to a multivariate model. The variance components are estimated by the REML method.

Usage

sPhikf.it(d, t, pp, sigmap, X, eta, phi1, phi2)

Arguments

Value

A list containing the following components.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Small area estimation of labour force indicators under a multinomial mixed model with correlated time and area effects. Submitted for review.

Schall, R (1991). Estimation in generalized linear models with random effects. Biometrika, 78,719-727.

See Also

data.mme, initial.values, wmatrix, phi.mult.it, prmu.time, phi.direct.it, Fbetaf.it, ci, modelfit2, msef.it, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2 #Type of model
data(simdata2) #data
datar=data.mme(simdata2,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
sigmap=wmatrix(datar$n,mean$estimated.probabilities)

##Fisher information matrix and score vectors
Fisher.phi=sPhikf.it(datar$d,datar$t,pp,sigmap,datar$X,mean$eta,initial$phi1.0,
           initial$phi2.0)

Dataset for Model 1

Description

Dataset used by the multinomial mixed effects model with one independent random effect in each category of the response variable (Model 1). This dataset contains 15 small areas. The response variable has three categories. The last is the reference category. The variables are as follows:

Usage

simdata

Format

A data frame with 15 rows and 9 variables in columns

Details

• area: area indicator.

• Time: time indicator.

• sample: the sample size of each domain.

• Population: the population size of each domain.

• Y1: the first category of the response variable.

• Y2: the second category of the response variable.

• Y3: the third category of the response variable.

• X1: the covariate for the first category of the response variable.

• X2: the covariate for the second category of the response variable.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata)  #data
mod=1 # type of model
datar=data.mme(simdata,k,pp,mod)
# Model fit
result=model(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N,mod)

#Analytic MSE
mse=msef(pp,datar$X,datar$Z,result,datar$N,datar$n)

B=1    #Bootstrap iterations
ss=12345 #SEED
set.seed(ss)

##Bootstrap parametric BIAS and MSE
mse.pboot=mseb(pp,datar$Xk,datar$X,datar$Z,datar$n,datar$N,result,B,mod)

Dataset for Model 2

Description

Dataset used by the multonomial mixed effects model with two independent random effects in each category of the response variable: one domain random effect and another independent time and domain random effect (Model 2). This dataset contains 10 small areas and two periods. The response variable has three categories. The last is the reference category. The variables are as follows:

Usage

simdata2

Format

A data frame with 30 rows and 9 variables in columns

Details

• area: area indicator.

• Time: time indicator.

• sample: the sample size of each domain.

• Population: the population size of each domain.

• Y1: the first category of the response variable.

• Y2: the second category of the response variable.

• Y3: the third category of the response variable.

• X1: the covariate for the first category of the response variable.

• X2: the covariate for the second category of the response variable.

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
data(simdata2)
mod=2 #type of model
datar=data.mme(simdata2,k,pp,mod)

##Model fit
result=model(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N,mod)

##Analytic MSE
msef=msef.it(pp,datar$X,result,datar$n,datar$N)

B=1    #Bootstrap iterations
ss=12345 #SEED
set.seed(ss)

##Bootstrap parametric BIAS and MSE
mse.pboot=mseb(pp,datar$Xk,datar$X,datar$Z,datar$n,datar$N,result,B,mod)

Dataset for Model 3

Description

Dataset used by the multonomial mixed effects model with two independent random effects in each category of the response variable: one domain random effect and another correlated time and domain random effect (Model 3). This dataset contains ten small areas and four periods. The response variable has three categories. The last is the reference category. The variables are as follows:

Usage

simdata3

Format

A data frame with 40 rows and 9 variables in columns

Details

• area: area indicator.

• Time: time indicator.

• sample: the sample size of each domain.

• Population: the population size of each domain.

• Y1: the first category of the response variable.

• Y2: the second category of the response variable.

• Y3: the third category of the response variable.

• X1: the covariate for the first category of the response variable.

• X2: the covariate for the second category of the response variable.

Examples

## Not run: 

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=3 #type of model
data(simdata3) #data
datar=data.mme(simdata3,k,pp,mod)
##Model fit
result=model(datar$d,datar$t,pp,datar$Xk,datar$X,datar$Z,datar$initial,datar$y[,1:(k-1)],
       datar$n,datar$N,mod)

##Analytic MSE
msef=msef.ct(pp,datar$X,result,datar$n,datar$N)

B=1    #Bootstrap iterations
ss=12345 #SEED
set.seed(ss)

##Bootstrap parametric BIAS and MSE
mse.pboot=mseb(pp,datar$Xk,datar$X,datar$Z,datar$n,datar$N,result,B,mod)
## End(Not run)

Model variance-covariance matrix of the multinomial mixed models

Description

This function calculates the variance-covariance matrix of the multinomial mixed models. Three types of multinomial mixed model are considered. The first model (Model 1), with one random effect in each category of the response variable; Model 2, introducing independent time effect; Model 3, introducing correlated time effect.

Usage

wmatrix(M, pr)

Arguments

Value

W a list with the model variance-covariance matrices for each domain.

References

Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling,13,153-178.

See Also

data.mme, initial.values, phi.mult, prmu, prmu.time Fbetaf, phi.direct, sPhikf, ci, modelfit1, msef, mseb

Examples

k=3 #number of categories of the response variable
pp=c(1,1) #vector with the number of auxiliary variables in each category
mod=2 #type of model
data(simdata2)
datar=data.mme(simdata2,k,pp,mod)
initial=datar$initial
mean=prmu.time(datar$n,datar$Xk,initial$beta.0,initial$u1.0,initial$u2.0)
##The model variance-covariance matrix
varcov=wmatrix(datar$n,mean$estimated.probabilities)