##
# 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 Modeling ####
#### Model 1a: mixed-effects model formulation ####
# In the MEM specification, the model is defined for all clusters and all individuals at once.
# Grouping of the data: Individuals of each cluster in consecutive columns. Within each cluster, the order of individuals is irrelvant.
# 9 variables (people) in a single row; "Yi_j" reads: Y value of individual i in cluster j.
Y1a<-matrix(c(9.03,9.54,10.80,1.22,3.12,2.18,15.75,13.69,14.13),nc=9,nr=1,
dimnames=list(NULL,c("Y1_1","Y2_1","Y3_1","Y1_2","Y2_2","Y3_2","Y1_3","Y2_3","Y3_3")))
# Univariate Multilevel Model with Random Intercept 1a: Estimating gamma00, rij, and u0j.
OpenMx.mod1a <- mxModel("Univariate random intercept model 1a MEM formulation",
mxData(
observed=Y1a,
type="raw"
),
# G: covariance matrix among cluster-effects (u 0j).
# G is a symmetric matrix with 3 rows and 3 columns (could also be declared as "Diag").
mxMatrix(
type="Symm",
nrow=3,
ncol=3,
values=c(1,
0,1,
0,0,1),
free=c(T,
F,T,
F,F,T),
labels=c("u0j",
NA,"u0j",
NA,NA,"u0j"),
byrow=TRUE,
name="G"
),
# R: covariance matrix among individual residuals (r ij).
# R is symmetric with 9 rows and 9 columns (could also be declared as "Diag").
mxMatrix(
type="Symm",
nrow=9,
ncol=9,
values=c(2,
0,2,
0,0,2,
0,0,0,2,
0,0,0,0,2,
0,0,0,0,0,2,
0,0,0,0,0,0,2,
0,0,0,0,0,0,0,2,
0,0,0,0,0,0,0,0,2),
free=c(T,
F,T,
F,F,T,
F,F,F,T,
F,F,F,F,T,
F,F,F,F,F,T,
F,F,F,F,F,F,T,
F,F,F,F,F,F,F,T,
F,F,F,F,F,F,F,F,T),
labels=c("rij",
NA,"rij",
NA,NA,"rij",
NA,NA,NA,"rij",
NA,NA,NA,NA,"rij",
NA,NA,NA,NA,NA,"rij",
NA,NA,NA,NA,NA,NA,"rij",
NA,NA,NA,NA,NA,NA,NA,"rij",
NA,NA,NA,NA,NA,NA,NA,NA,"rij"),
byrow=TRUE,
name="R"
),
# Z: random-effects design matrix: dummy variables indicating cluster membership.
# Z is a full matrix with 9 rows and 3 columns.
mxMatrix(
type="Full",
nrow=9,
ncol=3,
values=c(1,0,0,
1,0,0,
1,0,0,
0,1,0,
0,1,0,
0,1,0,
0,0,1,
0,0,1,
0,0,1),
free=FALSE,
byrow=TRUE,
name="Z"
),
# X: design matrix for fixed-effects (vector of "1").
# X is a unity matrix with 9 rows and 1 column: all elements are fixed to 1.
mxMatrix(
type="Unit",
nrow=9,
ncol=1,
free=FALSE,
name="X"
),
# B: fixed-effects coefficients vector.
# B is a single fixed-effects parameter, i.e. the mean intercept.
mxMatrix(
type="Full",
nrow=1,
ncol=1,
values=6,
free=TRUE,
labels="gamma00",
name="B"
),
# V: model-implied covariance matrix for the entire sample.
mxAlgebra(
expression = Z %*% G %*% t(Z) + R,
name="V"
),
# M: model-implied means for the entire sample.
# 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 covariance.
means="M", # matrix of model-implied mean.
dimnames=dimnames(Y1a)[[2]]
)
)
OpenMx.fit1a <- mxRun(OpenMx.mod1a)
## Running Univariate random intercept model 1a MEM formulation
summary(OpenMx.fit1a)
## data:
## $`Univariate random intercept model 1a MEM formulation.data`
## Y1_1 Y2_1 Y3_1 Y1_2
## Min. :9.03 Min. :9.54 Min. :10.8 Min. :1.22
## 1st Qu.:9.03 1st Qu.:9.54 1st Qu.:10.8 1st Qu.:1.22
## Median :9.03 Median :9.54 Median :10.8 Median :1.22
## Mean :9.03 Mean :9.54 Mean :10.8 Mean :1.22
## 3rd Qu.:9.03 3rd Qu.:9.54 3rd Qu.:10.8 3rd Qu.:1.22
## Max. :9.03 Max. :9.54 Max. :10.8 Max. :1.22
## Y2_2 Y3_2 Y1_3 Y2_3
## Min. :3.12 Min. :2.18 Min. :15.8 Min. :13.7
## 1st Qu.:3.12 1st Qu.:2.18 1st Qu.:15.8 1st Qu.:13.7
## Median :3.12 Median :2.18 Median :15.8 Median :13.7
## Mean :3.12 Mean :2.18 Mean :15.8 Mean :13.7
## 3rd Qu.:3.12 3rd Qu.:2.18 3rd Qu.:15.8 3rd Qu.:13.7
## Max. :3.12 Max. :2.18 Max. :15.8 Max. :13.7
## Y3_3
## Min. :14.1
## 1st Qu.:14.1
## Median :14.1
## Mean :14.1
## 3rd Qu.:14.1
## Max. :14.1
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 u0j G 1 1 25.5590 21.1332
## 2 rij R 1 1 0.9699 0.5599
## 3 gamma00 B 1 1 8.8289 2.9372
##
## observed statistics: 9
## estimated parameters: 3
## degrees of freedom: 6
## -2 log likelihood: 38.41
## saturated -2 log likelihood: NA
## number of observations: 1
## chi-square: NA
## p: NA
## Information Criteria:
## df Penalty Parameters Penalty Sample-Size Adjusted
## AIC: 26.41 44.41 NA
## BIC: 38.41 38.41 32.18
## CFI: NA
## TLI: NA
## RMSEA: NA
## timestamp: 2013-06-23 12:57:41
## frontend time: 0.229 secs
## backend time: 0.002 secs
## independent submodels time: 0 secs
## wall clock time: 0.231 secs
## cpu time: 0.231 secs
## openmx version number: 1.3.2-2301
# This model is equivalent to:
Y1.lme<-as.data.frame(matrix(c(9.03,9.54,10.80,
1.22,3.12,2.18,
15.75,13.69,14.13,
rep(1:3,each=3)),
nc=2,nr=9,dimnames=list(c("Y1_1","Y2_1","Y3_1","Y1_2","Y2_2","Y3_2","Y1_3","Y2_3","Y3_3"),c("Y","cluster"))))
lme.mod1<-lme(Y~1,random=~1|cluster,data=Y1.lme,method="ML")
summary(lme.mod1)
## Linear mixed-effects model fit by maximum likelihood
## Data: Y1.lme
## AIC BIC logLik
## 44.41 45.01 -19.21
##
## Random effects:
## Formula: ~1 | cluster
## (Intercept) Residual
## StdDev: 5.056 0.9848
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 8.829 3.115 6 2.834 0.0298
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.0524 -0.7595 -0.2417 0.8769 1.3178
##
## Number of Observations: 9
## Number of Groups: 3
VarCorr(lme.mod1)
## cluster = pdLogChol(1)
## Variance StdDev
## (Intercept) 25.5590 5.0556
## Residual 0.9699 0.9848
#### Model 1b: 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.
Y1b<-matrix(c(9.03,9.54,10.80,
1.22,3.12,2.18,
15.75,13.69,14.13),
nc=3,nr=3,byrow=T,dimnames=list(NULL,c("Y1_j","Y2_j","Y3_j")))
# Univariate Multilevel Model with Random Intercept 1b: Estimating gamma00, rij, and u0j.
OpenMx.mod1b <- mxModel("Univariate random intercept model 1b SEM CFA formulation",
mxData(
observed=Y1b,
type="raw"
),
# G: between-cluster variance (u 0j).
# Matrix specification: a single free element.
# Number of rows & columns = number of DVs.
mxMatrix(
type="Symm",
nrow=1,
ncol=1,
values=1.1,
free=TRUE,
labels="u0j",
name="G"
),
# R: within-cluster residual variance (r ij).
# Matrix specification: diagonal with equal variances.
# Number of rows & columns = number of observations.
mxMatrix(
type="Symm",
nrow=3,
ncol=3,
values=c(2,
0,2,
0,0,2),
free=c(T,
F,T,
F,F,T),
labels=c("r0j",
NA,"r0j",
NA,NA,"r0j"),
byrow=TRUE,
name="R"
),
# Z: random-effects design matrix (vector of 1s).
# Number of rows = number of observations.
mxMatrix(
type="Full",
nrow=3,
ncol=1,
values=c(1,
1,
1),
free=FALSE,
byrow=TRUE,
name="Z"
),
# X: fixed-effects design matrix (vector of 1s).
# Number of rows = number of observations.
mxMatrix(
type="Full",
nrow=3,
ncol=1,
values=c(1,
1,
1),
free=FALSE,
byrow=TRUE,
name="X"
),
# B: vector of fixed-effects coefficient.
# Single fixed-effects parameter (grand-mean or intercept, gamma00).
mxMatrix(
type="Full",
nrow=1,
ncol=1,
values=6,
free=TRUE,
labels="gamma00",
name="B"
),
# V: computing model-implied covariance matrix (for cluster j).
mxAlgebra(
expression = Z %*% G %*% t(Z) + R,
name="V"
),
# 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(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 covariance
means="M", # matrix of model-implied mean
dimnames=dimnames(Y1b)[[2]]
)
)
OpenMx.fit1b <- mxRun(OpenMx.mod1b)
## Running Univariate random intercept model 1b SEM CFA formulation
summary(OpenMx.fit1b)
## data:
## $`Univariate random intercept model 1b SEM CFA formulation.data`
## Y1_j Y2_j Y3_j
## Min. : 1.22 Min. : 3.12 Min. : 2.18
## 1st Qu.: 5.12 1st Qu.: 6.33 1st Qu.: 6.49
## Median : 9.03 Median : 9.54 Median :10.80
## Mean : 8.67 Mean : 8.78 Mean : 9.04
## 3rd Qu.:12.39 3rd Qu.:11.62 3rd Qu.:12.46
## Max. :15.75 Max. :13.69 Max. :14.13
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 u0j G 1 1 25.5590 21.1349
## 2 r0j R 1 1 0.9699 0.5599
## 3 gamma00 B 1 1 8.8289 2.9373
##
## observed statistics: 9
## estimated parameters: 3
## degrees of freedom: 6
## -2 log likelihood: 38.41
## saturated -2 log likelihood: NA
## number of observations: 3
## chi-square: NA
## p: NA
## Information Criteria:
## df Penalty Parameters Penalty Sample-Size Adjusted
## AIC: 26.41 44.41 NA
## BIC: 31.82 41.71 33.71
## CFI: NA
## TLI: NA
## RMSEA: NA
## timestamp: 2013-06-23 12:57:41
## frontend time: 0.088 secs
## backend time: 0.001 secs
## independent submodels time: 0 secs
## wall clock time: 0.089 secs
## cpu time: 0.089 secs
## openmx version number: 1.3.2-2301
# This model is equivalent to (see above for model specification):
summary(lme.mod1)
## Linear mixed-effects model fit by maximum likelihood
## Data: Y1.lme
## AIC BIC logLik
## 44.41 45.01 -19.21
##
## Random effects:
## Formula: ~1 | cluster
## (Intercept) Residual
## StdDev: 5.056 0.9848
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 8.829 3.115 6 2.834 0.0298
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.0524 -0.7595 -0.2417 0.8769 1.3178
##
## Number of Observations: 9
## Number of Groups: 3
VarCorr(lme.mod1)
## cluster = pdLogChol(1)
## Variance StdDev
## (Intercept) 25.5590 5.0556
## Residual 0.9699 0.9848
#### Model 2: SEM/CFA (parallel-test) formulation with additional clusters & observations ####
# 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 4 observations (columns) per cluster -> adding observations requires adjustment of matrices R, Z, X.
Y2a<-matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87),
nc=4,nr=3,byrow=T,dimnames=list(NULL,c("Y1_j","Y2_j","Y3_j","Y4_j")))
# 4 clusters (rows) with a maximum of 4 observations (columns) per cluster -> adding clusters requires no adjustments.
# The new cluster strongly affects within-cluster variance.
Y2b<-matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87,
1.02,16.23,3.43,18.65),
nc=4,nr=4,byrow=T,dimnames=list(NULL,c("Y1_j","Y2_j","Y3_j","Y4_j")))
# 4 clusters (rows) with a maximum of 4 observations (columns) per cluster -> adding clusters requires no adjustments.
# The new cluster strongly affects between-cluster variance.
Y2c<-matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87,
50.03,48.10,52.12,51.20),
nc=4,nr=4,byrow=T,dimnames=list(NULL,c("Y1_j","Y2_j","Y3_j","Y4_j")))
# Univariate Multilevel Model with Random Intercept 2: Estimating gamma00, rij, and u0j.
# This is the container model without the data argument; data Y2a, Y2b, Y2c are added below in three "super" models.
OpenMx.mod2sub <- mxModel("Univariate random intercept model 2 SEM CFA formulation",
# G: between-cluster variance (u 0j).
# Matrix specification: a single free element.
# Number of rows & columns = number of DVs.
mxMatrix(
type="Symm",
nrow=1,
ncol=1,
values=1.1,
free=TRUE,
labels="u0j",
name="G"
),
# R: within-cluster residual variance (r ij).
# Matrix specification: diagonal with equal variances.
# Number of rows & columns = number of observations.
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("r0j",
NA,"r0j",
NA,NA,"r0j",
NA,NA,NA,"r0j"),
byrow=TRUE,
name="R"
),
# Z: random-effects design matrix (vector of 1s).
# Number of rows = number of observations.
mxMatrix(
type="Full",
nrow=4,
ncol=1,
values=c(1,
1,
1,
1),
free=FALSE,
byrow=TRUE,
name="Z"
),
# X: fixed-effects design matrix (vector of 1s).
# Number of rows = number of observations.
mxMatrix(
type="Full",
nrow=4,
ncol=1,
values=c(1,
1,
1,
1),
free=FALSE,
byrow=TRUE,
name="X"
),
# B: vector of fixed-effects coefficient.
# Single fixed-effects parameter (grand-mean or intercept, gamma00).
mxMatrix(
type="Full",
nrow=1,
ncol=1,
values=6,
free=TRUE,
labels="gamma00",
name="B"
),
# V: computing model-implied covariance matrix (for cluster j).
mxAlgebra(
expression = Z %*% G %*% t(Z) + R,
name="V"
),
# 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(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 covariance
means="M", # matrix of model-implied mean
dimnames=dimnames(Y2a)[[2]]
)
)
# "Super" model for Y2a data.
OpenMx.mod2a <- mxModel(OpenMx.mod2sub,
mxData(
observed=Y2a,
type="raw"
)
)
OpenMx.fit2a <- mxRun(OpenMx.mod2a)
## Running Univariate random intercept model 2 SEM CFA formulation
summary(OpenMx.fit2a)
## data:
## $`Univariate random intercept model 2 SEM CFA formulation.data`
## Y1_j Y2_j Y3_j Y4_j
## Min. : 1.22 Min. : 3.12 Min. : 2.18 Min. : 2.99
## 1st Qu.: 5.12 1st Qu.: 6.33 1st Qu.: 6.49 1st Qu.: 5.51
## Median : 9.03 Median : 9.54 Median :10.80 Median : 8.04
## Mean : 8.67 Mean : 8.78 Mean : 9.04 Mean : 9.30
## 3rd Qu.:12.39 3rd Qu.:11.62 3rd Qu.:12.46 3rd Qu.:12.46
## Max. :15.75 Max. :13.69 Max. :14.13 Max. :16.87
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 u0j G 1 1 26.748 22.1282
## 2 r0j R 1 1 1.416 0.6676
## 3 gamma00 B 1 1 8.947 3.0056
##
## observed statistics: 12
## estimated parameters: 3
## degrees of freedom: 9
## -2 log likelihood: 51.24
## saturated -2 log likelihood: NA
## number of observations: 3
## chi-square: NA
## p: NA
## Information Criteria:
## df Penalty Parameters Penalty Sample-Size Adjusted
## AIC: 33.24 57.24 NA
## BIC: 41.36 54.54 46.54
## CFI: NA
## TLI: NA
## RMSEA: NA
## timestamp: 2013-06-23 12:57:41
## frontend time: 0.05 secs
## backend time: 0.001 secs
## independent submodels time: 0.0009999 secs
## wall clock time: 0.052 secs
## cpu time: 0.052 secs
## openmx version number: 1.3.2-2301
# This model is equivalent to:
Y2a.lme<-as.data.frame(matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87,
rep(1:3,each=4)),
nc=2,nr=12,
dimnames=list(c("Y1_1","Y2_1","Y3_1","Y4_1","Y1_2","Y2_2","Y3_2","Y4_2","Y1_3","Y2_3","Y3_3","Y4_3"),c("Y","cluster"))))
lme.mod2a<-lme(Y~1,random=~1|cluster,data=Y2a.lme,method="ML")
summary(lme.mod2a)
## Linear mixed-effects model fit by maximum likelihood
## Data: Y2a.lme
## AIC BIC logLik
## 57.24 58.7 -25.62
##
## Random effects:
## Formula: ~1 | cluster
## (Intercept) Residual
## StdDev: 5.172 1.19
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 8.947 3.139 9 2.85 0.0191
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.12555 -0.82805 -0.03803 0.56521 1.54657
##
## Number of Observations: 12
## Number of Groups: 3
VarCorr(lme.mod2a)
## cluster = pdLogChol(1)
## Variance StdDev
## (Intercept) 26.748 5.172
## Residual 1.416 1.190
# "Super" model for Y2b data.
OpenMx.mod2b <- mxModel(OpenMx.mod2sub,
mxData(
observed=Y2b,
type="raw"
)
)
OpenMx.fit2b <- mxRun(OpenMx.mod2b)
## Running Univariate random intercept model 2 SEM CFA formulation
summary(OpenMx.fit2b)
## data:
## $`Univariate random intercept model 2 SEM CFA formulation.data`
## Y1_j Y2_j Y3_j Y4_j
## Min. : 1.02 Min. : 3.12 Min. : 2.18 Min. : 2.99
## 1st Qu.: 1.17 1st Qu.: 7.93 1st Qu.: 3.12 1st Qu.: 6.78
## Median : 5.12 Median :11.62 Median : 7.12 Median :12.46
## Mean : 6.75 Mean :10.64 Mean : 7.63 Mean :11.64
## 3rd Qu.:10.71 3rd Qu.:14.32 3rd Qu.:11.63 3rd Qu.:17.32
## Max. :15.75 Max. :16.23 Max. :14.13 Max. :18.65
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 u0j G 1 1 15.264 14.631
## 2 r0j R 1 1 20.840 8.508
## 3 gamma00 B 1 1 9.168 2.262
##
## observed statistics: 16
## estimated parameters: 3
## degrees of freedom: 13
## -2 log likelihood: 99.47
## saturated -2 log likelihood: NA
## number of observations: 4
## chi-square: NA
## p: NA
## Information Criteria:
## df Penalty Parameters Penalty Sample-Size Adjusted
## AIC: 73.47 105.5 NA
## BIC: 81.45 103.6 95.31
## CFI: NA
## TLI: NA
## RMSEA: NA
## timestamp: 2013-06-23 12:57:41
## frontend time: 0.052 secs
## backend time: 0.001 secs
## independent submodels time: 0 secs
## wall clock time: 0.053 secs
## cpu time: 0.053 secs
## openmx version number: 1.3.2-2301
# This model is equivalent to:
Y2b.lme<-as.data.frame(matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87,
1.02,16.23,3.43,18.65,
rep(1:4,each=4)),
nc=2,nr=16,
dimnames=list(c("Y1_1","Y2_1","Y3_1","Y4_1","Y1_2","Y2_2","Y3_2","Y4_2","Y1_3","Y2_3","Y3_3","Y4_3","Y1_4","Y2_4","Y3_4","Y4_4"),c("Y","cluster"))))
lme.mod2b<-lme(Y~1,random=~1|cluster,data=Y2b.lme,method="ML")
summary(lme.mod2b)
## Linear mixed-effects model fit by maximum likelihood
## Data: Y2b.lme
## AIC BIC logLik
## 105.5 107.8 -49.73
##
## Random effects:
## Formula: ~1 | cluster
## (Intercept) Residual
## StdDev: 3.907 4.565
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 9.168 2.337 12 3.924 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.8934 -0.3134 -0.0201 0.3634 1.9686
##
## Number of Observations: 16
## Number of Groups: 4
VarCorr(lme.mod2b)
## cluster = pdLogChol(1)
## Variance StdDev
## (Intercept) 15.26 3.907
## Residual 20.84 4.565
# "Super" model for Y2c data.
OpenMx.mod2c <- mxModel(OpenMx.mod2sub,
mxData(
observed=Y2c,
type="raw"
)
)
OpenMx.fit2c <- mxRun(OpenMx.mod2c)
## Running Univariate random intercept model 2 SEM CFA formulation
## Warning: In model 'Univariate random intercept model 2 SEM CFA
## formulation' NPSOL returned a non-zero status code 1. The final iterate
## satisfies the optimality conditions to the accuracy requested, but the
## sequence of iterates has not yet converged. NPSOL was terminated because
## no further improvement could be made in the merit function (Mx status
## GREEN).
summary(OpenMx.fit2c)
## data:
## $`Univariate random intercept model 2 SEM CFA formulation.data`
## Y1_j Y2_j Y3_j Y4_j
## Min. : 1.22 Min. : 3.12 Min. : 2.18 Min. : 2.99
## 1st Qu.: 7.08 1st Qu.: 7.93 1st Qu.: 8.64 1st Qu.: 6.78
## Median :12.39 Median :11.62 Median :12.46 Median :12.46
## Mean :19.01 Mean :18.61 Mean :19.81 Mean :19.77
## 3rd Qu.:24.32 3rd Qu.:22.29 3rd Qu.:23.63 3rd Qu.:25.45
## Max. :50.03 Max. :48.10 Max. :52.12 Max. :51.20
##
## The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN).
##
## free parameters:
## name matrix row col Estimate Std.Error lbound ubound
## 1 u0j G 1 1 341.486 241.7417
## 2 r0j R 1 1 1.814 0.7404
## 3 gamma00 B 1 1 19.301 9.2453
##
## observed statistics: 16
## estimated parameters: 3
## degrees of freedom: 13
## -2 log likelihood: 81.44
## saturated -2 log likelihood: NA
## number of observations: 4
## chi-square: NA
## p: NA
## Information Criteria:
## df Penalty Parameters Penalty Sample-Size Adjusted
## AIC: 55.44 87.44 NA
## BIC: 63.41 85.59 77.28
## CFI: NA
## TLI: NA
## RMSEA: NA
## timestamp: 2013-06-23 12:57:41
## frontend time: 0.052 secs
## backend time: 0.002 secs
## independent submodels time: 0 secs
## wall clock time: 0.054 secs
## cpu time: 0.054 secs
## openmx version number: 1.3.2-2301
# This model is equivalent to:
Y2c.lme<-as.data.frame(matrix(c(9.03,9.54,10.80,8.04,
1.22,3.12,2.18,2.99,
15.75,13.69,14.13,16.87,
50.03,48.10,52.12,51.20,
rep(1:4,each=4)),
nc=2,nr=16,
dimnames=list(c("Y1_1","Y2_1","Y3_1","Y4_1","Y1_2","Y2_2","Y3_2","Y4_2","Y1_3","Y2_3","Y3_3","Y4_3","Y1_4","Y2_4","Y3_4","Y4_4"),c("Y","cluster"))))
lme.mod2c<-lme(Y~1,random=~1|cluster,data=Y2c.lme,method="ML")
summary(lme.mod2c)
## Linear mixed-effects model fit by maximum likelihood
## Data: Y2c.lme
## AIC BIC logLik
## 87.44 89.75 -40.72
##
## Random effects:
## Formula: ~1 | cluster
## (Intercept) Residual
## StdDev: 18.48 1.347
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 19.3 9.549 12 2.021 0.0661
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.64934 -0.76787 -0.01694 0.56410 1.33554
##
## Number of Observations: 16
## Number of Groups: 4
VarCorr(lme.mod2c)
## cluster = pdLogChol(1)
## Variance StdDev
## (Intercept) 341.486 18.479
## Residual 1.814 1.347