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
#### Univariate Multilevel Random-Intercept & Random-Slope Modeling ####
#### Model 4: 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.
# Outcome Y: 5 clusters (rows) with a maximum of 4 observations (columns) per cluster.
# Predictor X: 5 clusters (rows) with a maximum of 4 observations (columns) per cluster, appended to each line.
Y4<-read.csv("OpenMx-ML-SEM_3_Y4.csv") # 80x8=640 data values
names(Y4)
[1] "Y1_j" "Y2_j" "Y3_j" "Y4_j" "X1_j" "X2_j" "X3_j" "X4_j"
head(Y4)
Y1_j Y2_j Y3_j Y4_j X1_j X2_j X3_j X4_j
1 1.667 3.000 3.222 3.667 2.0 2.8 3.2 3.0
2 3.222 4.000 3.444 3.778 3.4 4.0 2.8 3.4
3 3.667 4.000 4.667 4.000 3.0 3.2 4.4 3.2
4 3.556 3.889 3.889 3.778 2.6 3.4 3.6 3.4
5 4.444 3.556 4.000 3.889 4.0 3.2 3.4 4.0
6 3.667 3.444 4.333 3.778 3.6 3.6 4.4 3.8
# Univariate Multilevel Model with Random Intercept & Random Slope 4: Estimating gamma00, gamma01, rij, and u0j.
OpenMx.mod4 <- mxModel("Univariate random intercept random slope model 4 SEM CFA formulation",
mxData(
observed=Y4,
type="raw"
),
# Lambda: factor-loading matrix. Corresponds to Xj and Zj matrices.
# Four observations (max. obs. in clusters = 4) by 2 factors (intercept & slope).
# First factor (random intercept): All factor loadings are fixed to value 1 (not estimated).
# Second-factor (random slope): All factor loadings are fixed to the values of the corresponding definition variable.
# -> The factor-loading for person 1 is fixed to X1; factor-loading for person 2 is fixed to X2, so on and so forth.
mxMatrix(
type="Full",
nrow=4,
ncol=2,
values=c(1,NA,
1,NA,
1,NA,
1,NA),
free=c(F,F,
F,F,
F,F,
F,F),
labels=c(NA,"data.X1_j",
NA,"data.X2_j",
NA,"data.X3_j",
NA,"data.X4_j"),
byrow=TRUE,
name="L"
),
# Psi: latent factor covariance matrix. Corresponds to the Gj matrix.
# Matrix specification: diagonal with variances, offdiagonal with covariances.
# Number of rows & columns = number of factors (intercept & slope).
mxMatrix(
type="Symm",
nrow=2,
ncol=2,
values=c(2,
0,1.5),
free=TRUE,
labels=c("u0j",
"u0j_u1j","u1j"),
byrow=TRUE,
name="P"
),
# Theta: residual covariance matrix. Corresponds to the Rj matrix.
# Diagonal matrix containing residual variances.
# Number of rows & columns = max. obs. in clusters.
mxMatrix(
type="Symm",
nrow=4,
ncol=4,
values=c(2,
0,2,
0,0,2,
0,0,0,2),
free=c(T,
F,T,
F,F,T,
F,F,F,T),
labels=c("rij",
NA,"rij",
NA,NA,"rij",
NA,NA,NA,"rij"),
byrow=TRUE,
name="T"
),
# Alpha: latent factor mean vector. Corresponds to the fixed-effects vector B.
# Number of rows = number of factors (intercept & slope).
mxMatrix(
type="Full",
nrow=2,
ncol=1,
values=c(1,
1),
labels=c("gamma00","gamma01"),
free=TRUE,
byrow=TRUE,
name="A"
),
# V: computing model-implied covariance matrix (for cluster j).
mxAlgebra(
expression = L %*% P %*% t(L) + T,
name="S"
),
# M: 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(L %*% A),
name="M"
),
# The objective function expects the means model, the covariance model, and the name of the DVs.
mxFIMLObjective(
covariance="S", # matrix of model-implied covariances
means="M", # matrix of model-implied means
dimnames=dimnames(Y4)[[2]][1:4]
)
)
OpenMx.fit4 <- mxRun(OpenMx.mod4)
Running Univariate random intercept random slope model 4 SEM CFA formulation
summary(OpenMx.fit4)
data:
$`Univariate random intercept random slope model 4 SEM CFA formulation.data`
Y1_j Y2_j Y3_j Y4_j
Min. :-5.19 Min. :-1.44 Min. :-1.40 Min. :-3.16
1st Qu.: 2.65 1st Qu.: 3.53 1st Qu.: 3.22 1st Qu.: 2.55
Median : 3.61 Median : 4.00 Median : 3.89 Median : 3.78
Mean : 3.50 Mean : 4.09 Mean : 3.84 Mean : 3.54
3rd Qu.: 4.43 3rd Qu.: 4.82 3rd Qu.: 4.67 3rd Qu.: 4.29
Max. : 8.00 Max. : 9.78 Max. : 8.02 Max. : 9.17
X1_j X2_j X3_j X4_j
Min. : 0.233 Min. :-0.391 Min. :0.576 Min. :-5.88
1st Qu.: 2.841 1st Qu.: 2.800 1st Qu.:2.987 1st Qu.: 2.60
Median : 3.479 Median : 3.600 Median :3.800 Median : 3.56
Mean : 3.626 Mean : 3.692 Mean :3.991 Mean : 3.50
3rd Qu.: 4.260 3rd Qu.: 4.317 3rd Qu.:4.610 3rd Qu.: 4.44
Max. :11.038 Max. : 9.841 Max. :9.271 Max. : 9.38
free parameters:
name matrix row col Estimate Std.Error lbound ubound
1 u0j P 1 1 1.007651 0.99167
2 u0j_u1j P 1 2 -0.311987 0.25248
3 u1j P 2 2 0.100876 0.06722
4 rij T 1 1 3.250770 0.31965
5 gamma00 A 1 1 3.758587 0.30553
6 gamma01 A 2 1 -0.003978 0.08098
observed statistics: 320
estimated parameters: 6
degrees of freedom: 314
-2 log likelihood: 1312
saturated -2 log likelihood: NA
number of observations: 80
chi-square: NA
p: NA
Information Criteria:
df Penalty Parameters Penalty Sample-Size Adjusted
AIC: 683.59 1324 NA
BIC: -64.37 1338 1319
CFI: NA
TLI: NA
RMSEA: NA
timestamp: 2013-06-23 12:57:51
frontend time: 0.18 secs
backend time: 0.107 secs
independent submodels time: 0 secs
wall clock time: 0.287 secs
cpu time: 0.287 secs
openmx version number: 1.3.2-2301
# This model is equivalent to:
Y4.lme<-read.csv("OpenMx-ML-SEM_3_Y4.lme.csv") # 320x2=640 data values + clustering variable -> equivalent to Y4
names(Y4.lme)
[1] "cluster" "Y" "X"
head(Y4.lme)
cluster Y X
1 1 1.667 2.0
2 1 3.000 2.8
3 1 3.222 3.2
4 1 3.667 3.0
5 2 3.222 3.4
6 2 4.000 4.0
lme.mod4<-lme(Y~X,random=~X|cluster,data=Y4.lme,method="ML")
summary(lme.mod4)
Linear mixed-effects model fit by maximum likelihood
Data: Y4.lme
AIC BIC logLik
1324 1346 -655.8
Random effects:
Formula: ~X | cluster
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.0038 (Intr)
X 0.3176 -0.979
Residual 1.8030
Fixed effects: Y ~ X
Value Std.Error DF t-value p-value
(Intercept) 3.759 0.3056 239 12.298 0.0000
X -0.004 0.0810 239 -0.049 0.9609
Correlation:
(Intr)
X -0.936
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.77279 -0.34387 0.04509 0.41728 3.23024
Number of Observations: 320
Number of Groups: 80
VarCorr(lme.mod4)
cluster = pdLogChol(X)
Variance StdDev Corr
(Intercept) 1.0076 1.0038 (Intr)
X 0.1009 0.3176 -0.979
Residual 3.2508 1.8030