Christian T. Seubert — Jun 23, 2013, 12:57 PM
##
# Multilevel Modeling in OpenMx - Demo Code.
##
# Based on the paper and the Mx code from the online appendix of:
# Mehta, P. D., & Neale, M. C. (2005). People are variables too: Multilevel structural equations modeling.
# Psychological Methods, 10(3), 259â284.
##
# Author: Christian T. Seubert <christian.seubert@uibk.ac.at>
# Date: 2013-06-22
##
require(OpenMx)
Loading required package: OpenMx
require(nlme)
Loading required package: nlme
#### Bivariate Multilevel Random-Intercept Modeling ####
#### Model 3: SEM/CFA (parallel-test) formulation ####
# In the SEM specification, the model is defined for each cluster j.
# Grouping of the data: Individuals in columns, clusters in rows. Within each cluster, the order of individuals is irrelvant.
# 3 clusters (rows) with a maximum of 3 observations (columns) per cluster * 2 outcomes (P & Q) = 6 observations in each cluster.
Y3<-matrix(c(11.1,20.6,13.3,16.0,14.3,15.1,
18.6,8.9,19.2,7.3,21.5,3.6,
10.4,2.7,11.1,3.2,12.8,5.0),nc=6,nr=3,byrow=T,dimnames=list(NULL,c("YP1_j","YQ1_j","YP2_j","YQ2_j","YP3_j","YQ3_j")))
# Bivariate Multilevel Model with Random Intercept 3: Estimating gamma00, rij, and u0j for two DVs.
OpenMx.mod3 <- mxModel("Bivariate random intercept model 3 SEM CFA formulation",
mxData(
observed=Y3,
type="raw"
),
# Gj: between-cluster variances & covariances for P, Q, PQ (u 0j).
# Matrix specification: diagonal with variances, offdiagonal with covariances.
# Number of rows & columns = number of DVs.
mxMatrix(
type="Symm",
nrow=2,
ncol=2,
values=c(2.1,
1,3.1),
free=c(T,
T,T),
labels=c("u0j_P",
"u0j_PQ","u0j_Q"),
byrow=TRUE,
name="G"
),
# Wij: within-cluster residual variances & covariances for P, Q, PQ (r ij).
# Matrix specification: diagonal with variances, offdiagonal with covariances.
# Number of rows & columns = number of DVs.
mxMatrix(
type="Symm",
nrow=2,
ncol=2,
values=c(2,
0,2),
free=c(T,
T,T),
labels=c("rij_P",
"rij_PQ","rij_P"),
byrow=TRUE,
name="W"
),
# Zj: random-effects design matrix.
# Number of rows = number of observations; number of columns = number of DVs.
mxMatrix(
type="Full",
nrow=6,
ncol=2,
values=c(1,0,
0,1,
1,0,
0,1,
1,0,
0,1),
free=FALSE,
byrow=TRUE,
name="Z"
),
# Xj: fixed-effects design matrix.
# Number of rows = number of observations; number of columns = number of DVs.
mxMatrix(
type="Full",
nrow=6,
ncol=2,
values=c(1,0,
0,1,
1,0,
0,1,
1,0,
0,1),
free=FALSE,
byrow=TRUE,
name="X"
),
# B: vector of fixed-effects coefficients.
# Fixed-effects parameters for P & Q (grand-means or intercepts, gamma00_P & gamma00_Q).
mxMatrix(
type="Full",
nrow=2,
ncol=1,
values=c(6,10),
free=TRUE,
labels=c("gamma00_P","gamma00_Q"),
name="B"
),
# I: Identity matrix with dimension = number of clusters.
# Used in the calculation of Rj from Wij.
mxMatrix(
type="Iden",
nrow=3,
ncol=3,
name="I"
),
# Rj: Wij is the expected within-cluster covariance for a single individual i in cluster j.
# Rj is the covariance matrix for the whole cluster j.
mxAlgebra(
expression = I %x% W, # %x% -> Kronecker product
name="R"
),
# Vj: computing model-implied covariance matrix (for cluster j).
mxAlgebra(
expression = Z %*% G %*% t(Z) + R,
name="V"
),
# Mj: computing model-implied means (for cluster j).
# other than in the Mx script, X %*% B must be transposed, otherwise mxFIMLObjective will produce an error.
# (means="M" expects a 1xn matrix, not a nx1 matrix)
mxAlgebra(
expression = t(X %*% B),
name="M"
),
# The objective function expects the means model, the covariance model, and the name of the DVs.
mxFIMLObjective(
covariance="V", # matrix of model-implied covariances
means="M", # matrix of model-implied means
dimnames=dimnames(Y3)[[2]]
)
)
OpenMx.fit3 <- mxRun(OpenMx.mod3)
Running Bivariate random intercept model 3 SEM CFA formulation
summary(OpenMx.fit3)
data:
$`Bivariate random intercept model 3 SEM CFA formulation.data`
YP1_j YQ1_j YP2_j YQ2_j
Min. :10.4 Min. : 2.7 Min. :11.1 Min. : 3.20
1st Qu.:10.8 1st Qu.: 5.8 1st Qu.:12.2 1st Qu.: 5.25
Median :11.1 Median : 8.9 Median :13.3 Median : 7.30
Mean :13.4 Mean :10.7 Mean :14.5 Mean : 8.83
3rd Qu.:14.8 3rd Qu.:14.8 3rd Qu.:16.2 3rd Qu.:11.65
Max. :18.6 Max. :20.6 Max. :19.2 Max. :16.00
YP3_j YQ3_j
Min. :12.8 Min. : 3.6
1st Qu.:13.6 1st Qu.: 4.3
Median :14.3 Median : 5.0
Mean :16.2 Mean : 7.9
3rd Qu.:17.9 3rd Qu.:10.1
Max. :21.5 Max. :15.1
free parameters:
name matrix row col Estimate Std.Error lbound ubound
1 u0j_P G 1 1 11.855 10.791
2 u0j_PQ G 1 2 -2.325 12.375
3 u0j_Q G 2 2 32.753 27.842
4 rij_P W 1 1 4.017 1.926
5 rij_PQ W 1 2 -2.474 1.926
6 gamma00_P B 1 1 14.700 2.097
7 gamma00_Q B 2 1 9.156 3.371
observed statistics: 18
estimated parameters: 7
degrees of freedom: 11
-2 log likelihood: 89.76
saturated -2 log likelihood: NA
number of observations: 3
chi-square: NA
p: NA
Information Criteria:
df Penalty Parameters Penalty Sample-Size Adjusted
AIC: 67.76 103.76 NA
BIC: 77.67 97.45 78.78
CFI: NA
TLI: NA
RMSEA: NA
timestamp: 2013-06-23 12:57:46
frontend time: 0.185 secs
backend time: 0.006 secs
independent submodels time: 0 secs
wall clock time: 0.191 secs
cpu time: 0.191 secs
openmx version number: 1.3.2-2301
# This model with two DVs cannot be computed with lme.