| Title: | Multinomial Mixed Effects Models |
|---|---|
| Description: | Fit Gaussian Multinomial mixed-effects models for small area estimation: Model 1, with one random effect in each category of the response variable (Lopez-Vizcaino,E. et al., 2013) <doi:10.1177/1471082X13478873>; Model 2, introducing independent time effect; Model 3, introducing correlated time effect. mme calculates direct and parametric bootstrap MSE estimators (Lopez-Vizcaino,E et al., 2014) <doi:10.1111/rssa.12085>. |
| Authors: | E. Lopez-Vizcaino, M.J. Lombardia and D. Morales |
| Maintainer: | E. Lopez-Vizcaino <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1-6 |
| Built: | 2025-02-11 03:26:02 UTC |
| Source: | https://github.com/cran/mme |
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.
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.
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.
addtolist(B_d, t, d)addtolist(B_d, t, d)
B_d |
a list in each area. |
t |
number of time periods. |
d |
number of areas. |
B_d a list of dimension d.
Fbetaf.it, Fbetaf.ct,
modelfit2, modelfit3
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)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)
This function adds rows from a matrix of dimension d*t*(k-1) times d*(k-1).
addtomatrix(C2, d, t, k)addtomatrix(C2, d, t, k)
C2 |
a matrix of dimension d*t*(k-1) times d*(k-1). |
d |
number of areas. |
t |
number of time periods. |
k |
number of categories of the response variable. |
C22 a matrix of dimension d*(k-1) times d*(k-1).
Fbetaf.it, Fbetaf.ct,
modelfit2,modelfit3
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)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)
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.
ci(a, F)ci(a, F)
a |
vector with the estimated parameters obtained from |
F |
inverse of the Fisher Information Matrix obtained from |
A list containing the following components.
Std.dev |
vector with the standard deviations of the parameters. The parameters are sorted per category. |
p.value |
vector with the p-values of the parameters for testing H0:a=0. |
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
modelfit1, modelfit2,
modelfit3.
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)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)
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.
data.mme(fi, k, pp, mod)data.mme(fi, k, pp, mod)
fi |
input data set with (d x t) rows and 4+k+sum(pp) columns. The first four columns of the data set are, in this order: the area indicator (integer), the time indicator (integer), the sample size of each domain and the population size of each domain. The following k columns are the categories of the response variable. Then, the auxiliary variables: first, the auxiliary variables of the first category, second, the auxiliary variables of the second category, and so on. Examples of input data sets are in |
k |
number of categories of the response variable. |
pp |
vector with the number of auxiliary variables per category. |
mod |
a number specifying the type of models: 1=multinomial mixed model with one independent random effect in each category of the response variable (Model 1), 2=multinomial mixed 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) and 3= multinomial 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). |
A list containing the following components.
n |
vector with the area sample sizes. |
N |
vector with the area population sizes. |
Z |
design matrix of random effects. |
Xk |
list of matrices with the auxiliary variables per category. The dimension of the list is the number of domains |
X |
list of matrices with the auxiliary variables. The dimension of the list is the number of categories of the response variable minus one. |
y |
matrix with the response variable. The rows are the domains and the columns are the categories of the response variable. |
initial |
a list with the initial values for the fixed and the random effects obtained from |
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
initial.values, wmatrix,
phi.mult, prmu,
Fbetaf, phi.direct,
sPhikf, ci,
modelfit1, msef,
mseb
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)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)
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.
Fbetaf(sigmap, X, Z, phi, y, mu, u)Fbetaf(sigmap, X, Z, phi, y, mu, u)
sigmap |
a list with the model variance-covariance matrices for each domain. |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects. |
phi |
vector with the values of the variance components obtained from |
y |
matrix with the response variable except the reference category. The rows are the domains and the columns are the categories of the response variable minus one. |
mu |
matrix with the estimated mean of the response variable obtained from |
u |
matrix with the values of random effects obtained from |
A list containing the following components.
F.beta.beta |
the first diagonal element of the inverse of the Fisher information matrix. |
F.beta.u |
the element out of the diagonal of the inverse of the Fisher information matrix. |
F.u.u |
the second diagonal element of the inverse of the Fisher information matrix. |
S.beta |
beta scores. |
S.u |
u scores. |
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13 ,153-178.
data.mme, initial.values,
wmatrix, phi.mult,
prmu, phi.direct,
sPhikf, ci,
modelfit1, msef,
mseb.
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)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)
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.
Fbetaf.ct(sigmap, X, Z, phi1, phi2, y, mu, u1, u2, rho)Fbetaf.ct(sigmap, X, Z, phi1, phi2, y, mu, u1, u2, rho)
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects. |
phi1 |
vector with the values of the variance components for the first random effects obtained from |
phi2 |
vector with the values of the variance components for the second random effects obtained from |
y |
matrix with the response variable, except the reference category. The rows are the domains and the columns are the categories of the response variable minus one. |
mu |
matrix with the estimated mean of the response variable. |
u1 |
matrix with the values of the first random effect obtained from |
u2 |
matrix with the values of the second random effect obtained from |
rho |
vector with the values of the correlation parameter obtained from |
A list containing the following components.
F |
the inverse of the Fisher information matrix of (beta, u1, u2). |
S |
(beta, u1, u2) score vectors |
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.
data.mme, initial.values,
wmatrix, phi.mult.ct,
prmu.time, phi.direct.ct,
sPhikf.ct, ci,
modelfit3, msef.ct, omega,
mseb
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)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)
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.
Fbetaf.it(sigmap, X, Z, phi1, phi2, y, mu, u1, u2)Fbetaf.it(sigmap, X, Z, phi1, phi2, y, mu, u1, u2)
sigmap |
a list with the model variance-covariance matrices for each domain. |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
phi1 |
vector with the first variance component obtained from |
phi2 |
vector with the second variance component obtained from |
y |
matrix with the response variable, except the reference category obtained from |
mu |
matrix with the estimated mean of the response variable obtained from |
u1 |
matrix with the values of the first random effect obtained from |
u2 |
matrix with the values of the second random effect obtained from |
A list containing the following components.
F |
the inverse of the Fisher information matrix. |
S |
(beta, u1, u2) scores |
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
prmu.time, phi.direct.it,
sPhikf.it, ci,
modelfit2, msef.it,
mseb.
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)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)
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.
initial.values(d, pp, datar, mod)initial.values(d, pp, datar, mod)
d |
number of areas. |
pp |
vector with the number of auxiliary variables per category. |
datar |
output of function |
mod |
a number specifying the type of model: 1=multinomial mixed model with one independent random effect for each category of the response variable (Model 1), 2=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) and 3= 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). |
A list containing the following components, depending on the chosen model.
beta.0 |
a list with the initial values for the fixed effects beta per category. |
phi.0 |
vector with the initial values for the variance components phi of Model 1. |
phi1.0 |
vector with the initial values for the variance components phi1 of Model 2 or 3. |
phi2.0 |
vector with the initial values for the variance components phi2 of Model 2 or 3. |
u |
matrix with the initial values for the random effect for Model 1. |
u1.0 |
matrix with the initial values for the first random effect for Model 2 or 3. |
u2.0 |
matrix with the initial values for the second random effect for Model 2 or 3. |
rho.0 |
vector with the initial values for the correlation parameter for Model 3. |
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.
data.mme, wmatrix,
phi.mult.it, prmu.time,
Fbetaf.it, phi.direct.it,
sPhikf.it, ci,
modelfit2, msef.it,
mseb
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$initialk=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
This function creates objects of class mmedata.
mmedata(fi, k, pp)mmedata(fi, k, pp)
fi |
input data set with (d X t) rows and 4+k+sum(pp) columns. The first four columns of the data set are, in this order: the area indicator (integer), the time indicator (integer), the sample size of each domain and the population size of each domain. The following k columns are the categories of the response variable. Then, the auxiliary variables: first, the auxiliary variables of the first category, second, the auxiliary variables of the second category, and so on. Examples of input data set are in |
k |
number of categories of the response variable. |
pp |
vector with the number of auxiliary variables per category. |
modelfit1, modelfit2,
modelfit3
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)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)
This function chooses one of the three models.
model(d, t, pp, Xk, X, Z, initial, y, M, MM, mod)model(d, t, pp, Xk, X, Z, initial, y, M, MM, mod)
d |
number of areas. |
t |
number of time periods. |
pp |
vector with the number of the auxiliary variables per category. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
initial |
output of the function |
y |
matrix with the response variable obtained from |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
mod |
a number specifying the type of models: 1=multinomial mixed model with one independent random effect in each category of the response variable (Model 1), 2=multinomial mixed 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) and 3= multinomial 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). |
the output of the function modelfit1, modelfit2 or modelfit3.
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.
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)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)
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
modelfit1(pp, Xk, X, Z, initial, y, M, MM)modelfit1(pp, Xk, X, Z, initial, y, M, MM)
pp |
vector with the number of the auxiliary variables per category. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
initial |
output of the function |
y |
matrix with the response variable except the reference category obtained from |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
A list containing the following components.
Estimated.probabilities |
matrix with the estimated probabilities for the categories of response variable. |
Fisher.information.matrix.phi |
Fisher information matrix of the random effect. |
Fisher.information.matrix.beta |
Fisher information matrix of the fixed effect. |
u |
matrix with the estimated random effects. |
mean |
matrix with the estimated mean of the response variable. |
warning1 |
0=OK,1=The model could not be fitted. |
warning2 |
0=OK,1=The value of the variance component is negative: the initial value is taken. |
beta.Stddev.p.value |
matrix with the estimated fixed effects, its standard deviations and its p-values. |
phi.Stddev.p.value |
matrix with the estimated variance components, its standard deviations and its p-values. |
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
data.mme, initial.values,
wmatrix, phi.mult,
prmu, phi.direct,
sPhikf, ci,
Fbetaf, msef,
mseb.
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)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)
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.
modelfit2(d, t, pp, Xk, X, Z, initial, y, M, MM)modelfit2(d, t, pp, Xk, X, Z, initial, y, M, MM)
d |
number of areas. |
t |
number of time periods. |
pp |
vector with the number of the auxiliary variables per category. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects |
initial |
output of the function |
y |
matrix with the response variable obtained from |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
A list containing the following components.
Estimated.probabilities |
matrix with the estimated probabilities for the categories of response variable. |
Fisher.information.matrix.phi |
Fisher information matrix of the variance components. |
Fisher.information.matrix.beta |
Fisher information matrix of the fixed effects. |
u1 |
matrix with the estimated first random effect. |
u2 |
matrix with the estimated second random effect. |
mean |
matrix with the estimated mean of response variable. |
warning1 |
0=OK,1=The model could not be fitted. |
warning2 |
0=OK,1=The value of the variance component is negative: the initial value is taken. |
beta.Stddev.p.value |
matrix with the estimated fixed effects, its standard deviations and its p-values. |
phi.Stddev.p.value |
matrix with the estimated variance components, its standard deviations and its p-values. |
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
prmu.time, phi.direct.it,
sPhikf.it, ci,
Fbetaf.it, msef.it,
mseb
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)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)
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.
modelfit3(d, t, pp, Xk, X, Z, initial, y, M, MM, b)modelfit3(d, t, pp, Xk, X, Z, initial, y, M, MM, b)
d |
number of areas. |
t |
number of time periods. |
pp |
vector with the number of the auxiliary variables per category. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
initial |
output of the function |
y |
matrix with the response variable obtained from |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
b |
parameter that indicates the bootstrap. |
A list containing the following components.
Estimated.probabilities |
matrix with the estimated probabilities for the categories of response variable. |
Fisher.information.matrix.phi |
Fisher information matrix of phi. |
Fisher.information.matrix.beta |
Fisher information matrix of beta. |
u1 |
matrix with the estimated first random effect. |
u2 |
matrix with the estimated second random effect. |
mean |
matrix with the estimated mean of the response variable. |
warning1 |
0=OK,1=The model could not be fitted. |
warning2 |
0=OK,1=The value of the variance component is negative: the initial value is taken. |
beta.Stddev.p.value |
matrix with the estimated fixed effects, its standard deviations and its p-values. |
phi.Stddev.p.value |
matrix with the estimated variance components, its standard deviations and its p-values. |
rho |
estimated correlation parameter. |
rho.Stddev.p.value |
matrix with the estimated correlation parameter, its standard deviations and its p-values. |
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.
data.mme, initial.values,
wmatrix, phi.mult.ct,
prmu.time, phi.direct.ct,
sPhikf.ct, omega, ci,
Fbetaf.ct, msef.ct,
mseb
## 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)## 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)
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.
mseb(pp, Xk, X, Z, M, MM, resul, B, mod)mseb(pp, Xk, X, Z, M, MM, resul, B, mod)
pp |
vector with the number of the auxiliary variables per category. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
resul |
|
B |
number of bootstrap replications. |
mod |
a number specifying the type of models: 1=multinomial mixed model with one independent random effect in each category of the response variable (Model 1), 2=multinomial mixed 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) and 3= multinomial 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). |
a list containing the following components.
bias.pboot |
BIAS of the parametric bootstrap estimator of the mean of the response variable |
mse.pboot |
MSE of the parametric bootstrap estimator of the mean of the response variable |
rmse.pboot |
RMSE of the parametric bootstrap estimator of the mean of the response variable |
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.
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.
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)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)
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.
msef(pp, X, Z, resul, MM, M)msef(pp, X, Z, resul, MM, M)
resul |
the output of the function |
X |
list of matrices with the auxiliary variables obtained from |
Z |
design matrix of random effects obtained from |
pp |
vector with the number of the auxiliary variables per category. |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
mse is a matrix with the MSE estimator calculated by adapting the explicit formulas of Prasad and Rao (1990).
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.
data.mme, initial.values,
wmatrix, phi.mult,
prmu, phi.direct,
sPhikf, modelfit1,
Fbetaf, ci,
mseb.
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)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)
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.
msef.ct(p, X, result, M, MM)msef.ct(p, X, result, M, MM)
p |
vector with the number of the auxiliary variables per category. |
X |
list of matrices with the auxiliary variables obtained from |
result |
the output of the function |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
mse.analitic is a matrix with the MSE estimator calculated by adapting the explicit formulas of Prasad and Rao (1990).
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.
data.mme, initial.values,
wmatrix, phi.mult.ct,
prmu.time, phi.direct.ct,
sPhikf.ct, modelfit3,
Fbetaf.ct, ci, omega,
mseb.
## 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)## 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)
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.
msef.it(p, X, result, M, MM)msef.it(p, X, result, M, MM)
p |
vector with the number of the auxiliary variables per category. |
X |
list of matrices with the auxiliary variables obtained from |
result |
the output of the function |
M |
vector with the area sample sizes. |
MM |
vector with the population sample sizes. |
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.
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
prmu.time, phi.direct.it,
sPhikf.it, modelfit2,
Fbetaf.it, ci,
mseb
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)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)
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.
omega(t, k, rho, phi2)omega(t, k, rho, phi2)
t |
number of time periods. |
k |
number of categories of the response variable. |
rho |
vector with the correlation parameter obtained from |
phi2 |
vector with the values of the second variance component obtained from |
A list containing the following components.
Omega.d |
correlation matrix. |
First.derivative.Omegad |
Fisher derivative of the model correlation matrix. |
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.
data.mme, initial.values,
wmatrix ,phi.mult.ct,
prmu.time, phi.direct.ct,
Fbetaf.ct, sPhikf.ct, ci,
modelfit3, msef.ct,
mseb
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)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)
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.
phi.direct(sigmap, phi, X, u)phi.direct(sigmap, phi, X, u)
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
phi |
vector with the initial values of the variance components obtained from |
u |
matrix with the values of the random effects obtained from |
a list containing the following components.
phi.new |
vector with the variance components. |
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.
data.mme, initial.values,
wmatrix, phi.mult,
prmu, Fbetaf,
sPhikf, ci,
modelfit1, msef,
mseb.
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)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)
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.
phi.direct.ct(p, sigmap, X, theta, phi1, phi2, u1, u2, rho)phi.direct.ct(p, sigmap, X, theta, phi1, phi2, u1, u2, rho)
p |
vector with the number of auxiliary variables per category. |
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
theta |
matrix with the estimated log-probabilites of each category in front of the reference category obtained from |
phi1 |
vector with the initial values of the first variance component obtained from |
phi2 |
vector with the initial values of the second variance component obtained from |
u1 |
matrix with the values of the first random effect obtained from |
u2 |
matrix with the values of the second random effect obtained from |
rho |
vector with the initial values of the correlation parameter obtained from |
a list containing the following components.
phi1.new |
vector with the values of the variance component for the first random effect. |
phi2.new |
vector with the values of the variance component for the second random effect. |
rho.new |
vector with the correlation parameter. |
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.
data.mme, initial.values,
wmatrix, phi.mult.ct,
prmu.time, Fbetaf.ct
sPhikf.ct, ci,
modelfit3, msef.ct,
mseb, omega
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)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)
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.
phi.direct.it(pp, sigmap, X, phi1, phi2, u1, u2)phi.direct.it(pp, sigmap, X, phi1, phi2, u1, u2)
pp |
vector with the number of auxiliary variables per category. |
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
phi1 |
vector with the initial values of the first variance component obtained from |
phi2 |
vector with the initial values of the second variance component obtained from |
u1 |
matrix with the values of the first random effect obtained from |
u2 |
matrix with the values of the second random effect obtained from |
a list containing the following components.
phi1.new |
vector with the values of the variance component for the first random effect. |
phi2.new |
vector with the values of the variance component for the second random effect. |
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
prmu.time, Fbetaf.it
sPhikf.it, ci,
modelfit2, msef.it,
mseb
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)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)
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).
phi.mult(beta.0, y, Xk, M)phi.mult(beta.0, y, Xk, M)
beta.0 |
initial values for the fixed effects obtained in |
y |
matrix with the response variable obtained from |
Xk |
list of matrices with the auxiliary variables per category obtained from |
M |
vector with the sample size of the areas. |
phi.0 vector of inicial values for the variance components
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
data.mme, initial.values,
wmatrix, prmu,
Fbetaf, phi.direct,
sPhikf, ci,
modelfit1, msef,
mseb.
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)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)
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).
phi.mult.ct(beta.0, y, Xk, M, u1, u2)phi.mult.ct(beta.0, y, Xk, M, u1, u2)
beta.0 |
a list with the initial values for the fixed effects per category obtained from |
y |
matrix with the response variable obtained from |
Xk |
list of matrices with the auxiliary variables per category obtained from |
M |
vector with the sample size of the areas. |
u1 |
matrix with the values for the first random effect obtained from |
u2 |
matrix with the values for the second random effect obtained from |
A list containing the following components.
phi.0 |
vector of the initial values for the variance components. |
rho.0 |
vector of the initial values for the correlation parameter. |
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.
data.mme, initial.values,
wmatrix, prmu.time,
Fbetaf.ct, phi.direct.ct,
sPhikf.ct, ci,
modelfit3, msef.ct,omega,
mseb.
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)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)
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).
phi.mult.it(beta.0, y, Xk, M, u1, u2)phi.mult.it(beta.0, y, Xk, M, u1, u2)
beta.0 |
initial values for the fixed effects obtained from |
y |
matrix with the response variable obtained from |
Xk |
list of matrices with the auxiliary variables per category obtained from |
M |
vector with the sample size of the areas. |
u1 |
vector with the initial values for the first random effect obtained from |
u2 |
vector with the initial values for the second random effect obtained from |
phi.0 vector of the initial values for the variance components.
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
data.mme, initial.values,
wmatrix, prmu.time,
Fbetaf.it, phi.direct.it,
sPhikf.it, ci,
modelfit2, msef.it,
mseb.
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)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)
This function prints objects of class mme.
## S3 method for class 'mme' print(x, ...)## S3 method for class 'mme' print(x, ...)
x |
a list with the output of |
... |
further information. |
modelfit1, modelfit2,
modelfit3
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) resultk=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
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).
prmu(M, Xk, beta, u)prmu(M, Xk, beta, u)
M |
vector with the area sample sizes. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
beta |
fixed effects obtained from |
u |
values of random effects obtained from |
A list containing the following components:
Estimated.probabilities |
matrix with the estimated probabilities for the categories of response variable. |
mean |
matrix with the estimated mean of the response variable. |
eta |
matrix with the estimated log-rates of the probabilities of each category over the reference category. |
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling, 13, 153-178.
data.mme, initial.values,
wmatrix, phi.mult,
Fbetaf, phi.direct,
sPhikf, ci,
modelfit1, msef,
mseb.
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)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)
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).
prmu.time(M, Xk, beta, u1, u2)prmu.time(M, Xk, beta, u1, u2)
M |
vector with the area sample sizes. |
Xk |
list of matrices with the auxiliary variables per category obtained from |
beta |
a list with the values for the fixed effects beta per category obtained from |
u1 |
a vector with the values of the first random effect obtained from |
u2 |
a vector with the values of the second random effect obtained from |
A list containing the following components:
Estimated.probabilities |
matrix with the estimated probabilities for the categories of response variable. |
mean |
matrix with the estimated mean of the response variable. |
eta |
matrix with the estimated log-rates of the probabilities of each category over the reference category. |
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
Fbetaf.it, phi.direct.it,
sPhikf.it, ci,
modelfit2, msef.it,
mseb
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)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)
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:
simdatasimdata
A data frame with 15 rows and 9 variables in columns
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.
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)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 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:
simdata2simdata2
A data frame with 30 rows and 9 variables in columns
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.
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)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 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:
simdata3simdata3
A data frame with 40 rows and 9 variables in columns
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.
## 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)## 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)
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.
sPhikf(pp, sigmap, X, eta, phi)sPhikf(pp, sigmap, X, eta, phi)
pp |
vector with the number of the auxiliary variables per category. |
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
eta |
matrix with the estimated log-rates of probabilities of each category over the reference category obtained from
|
phi |
vector with the values of the variance components obtained from |
A list containing the following components.
S.k |
phi score vector. |
F |
Fisher information matrix of the variance component phi. |
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.
data.mme, initial.values,
wmatrix, phi.mult,
prmu, phi.direct,
Fbetaf, ci,
modelfit1, msef,
mseb.
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)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)
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.
sPhikf.ct(d, t, pp, sigmap, X, eta, phi1, phi2, rho, pr, M)sPhikf.ct(d, t, pp, sigmap, X, eta, phi1, phi2, rho, pr, M)
d |
number of areas. |
t |
number of time periods. |
pp |
vector with the number of the auxiliary variables per category. |
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
eta |
matrix with the estimated log-rates of probabilites of each category over the reference category obtained from |
phi1 |
vector with the values of the first variance component obtained from |
phi2 |
vector with the values of the second variance component obtained from |
rho |
vector with the correlation parameter obtained from |
pr |
matrix with the estimated probabilities of the response variable obtained from |
M |
vector with the area sample sizes. |
A list containing the following components.
S |
(phi1, phi2, rho) score vector. |
F |
Fisher information matrix of the variance components (phi1, phi2, rho). |
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.
data.mme, initial.values,
wmatrix, phi.mult.ct,
prmu.time, phi.direct.ct,
Fbetaf.ct, omega, ci,
modelfit3, msef.ct,
mseb
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)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)
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.
sPhikf.it(d, t, pp, sigmap, X, eta, phi1, phi2)sPhikf.it(d, t, pp, sigmap, X, eta, phi1, phi2)
d |
number of areas. |
t |
number of time periods. |
pp |
vector with the number of the auxiliary variables per category. |
sigmap |
a list with the model variance-covariance matrices for each domain obtained from |
X |
list of matrices with the auxiliary variables obtained from |
eta |
matrix with the estimated log-rates of probabilities of each category over the reference category obtained from |
phi1 |
vector with the values of the first variance component obtained from |
phi2 |
vector with the values of the second variance component obtained from |
A list containing the following components.
S |
phi score vector. |
F |
Fisher information matrix of the variance components phi1 and phi2. |
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.
data.mme, initial.values,
wmatrix, phi.mult.it,
prmu.time, phi.direct.it,
Fbetaf.it, ci,
modelfit2, msef.it,
mseb
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)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)
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.
wmatrix(M, pr)wmatrix(M, pr)
M |
vector with area sample sizes. |
pr |
matrix with the estimated probabilities for the categories of the
response variable obtained from |
W a list with the model variance-covariance matrices for each domain.
Lopez-Vizcaino, ME, Lombardia, MJ and Morales, D (2013). Multinomial-based small area estimation of labour force indicators. Statistical Modelling,13,153-178.
data.mme, initial.values,
phi.mult, prmu, prmu.time
Fbetaf, phi.direct,
sPhikf, ci,
modelfit1, msef,
mseb
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)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)