Load xxM and data

library(xxm)
data(lranslp.xxm, package = "xxm")

Define local variables for conveinence

levels <- c("l1", "l2")
id1 <- levels
ys1 <- c("y1", "y2", "y3", "y4")
xs1 <- c("x")
l1vars <- c(levels, ys1, xs1)
es1 <- c("fw")
id2 <- "l2"
l2vars <- c("l2")
es2 <- c("int", "slp")

Construct R-matrices

For each parameter matrix, construct three related matrices:

  1. pattern matrix: A matrix indicating free or fixed parameters.
  2. value matrix: with start or fixed values for corresponding parameters.
  3. label matrix: with user friendly label for each parameter. label matrix is optional.
########################## Student Model Matrices ################ Student: factor-loading matrix:
########################## Pattern and value
ly1_pat <- matrix(1, 4, 1)
ly1_val <- matrix(0.1, 4, 1)
# Student: factor-covariance matrix: Pattern and value
ps1_pat <- diag(0, 1)
ps1_val <- diag(1, 1)
# Student: observed residual-covariance matrix: Pattern and value
th1_pat <- diag(1, 4)
th1_val <- diag(1, 4)
###################### School Model Matrices #################### School: factor-covariance
###################### matrix: Pattern and value
ps2_pat <- matrix(1, 2, 2)
ps2_val <- diag(0.2, 2)
# School: means: Pattern and value
al2_pat <- matrix(1, 2, 1)
al2_val <- matrix(c(1, 0.2), 2, 1)
###################### School - > Student matrices ############ School->Student: Beta
###################### (latent-on-latent regression) matrix: Pattern, value AND label
be12_pat <- matrix(c(0, 0), 1, 2)
be12_val <- matrix(c(1, 0), 1, 2)
be12_label <- matrix(c("one", "l1.x"), 1, 2)

Construct main model object

xxmModel() is used to declare level names. The function returns a model object that is passed as a parameter to subsequent statements.

xm <- xxmModel(levels = levels)
## A new model was created.

Add submodels to the model objects

For each declared level xxmSubmodel() is invoked to add corresponding submodel to the model object. The function adds three pieces of information: 1. parents declares a list of parents of the current level. 2. variables declares names of observed dependent (ys), observed independent (xs) and latent variables (etas) for the level. 3. data R data object for the current level.

### Submdeol: Student
xm <- xxmSubmodel(model = xm, level = "l1", parents = "l2", ys = ys1, xs = xs1, 
    etas = es1, data = l1)
## Submodel for level `l1` was created.

### Submdeol: School
xm <- xxmSubmodel(model = xm, level = "l2", parents = , ys = , xs = , etas = es2, 
    data = l2)
## Submodel for level `l2` was created.

Add Within-level parameter matrices for each submodel

For each declared level xxmWithinMatrix() is used to add within-level parameter matrices. For each parameter matrix, the function adds the three matrices constructed earlier:

# l1 within matrices (lambda, psi, and theta)
xm <- xxmWithinMatrix(xm, "l1", "lambda", ly1_pat, ly1_val, , )
## 
## 'lambda' matrix does not exist and will be added.
##  Added `lambda` matrix.
xm <- xxmWithinMatrix(xm, "l1", "psi", ps1_pat, ps1_val, , )
## 
## 'psi' matrix does not exist and will be added.
##  Added `psi` matrix.
xm <- xxmWithinMatrix(xm, "l1", "theta", th1_pat, th1_val, , )
## 
## 'theta' matrix does not exist and will be added.
##  Added `theta` matrix.
# l2 within matrices (psi and alpha)
xm <- xxmWithinMatrix(xm, "l2", "psi", ps2_pat, ps2_val, , )
## 
## 'psi' matrix does not exist and will be added.
##  Added `psi` matrix.
xm <- xxmWithinMatrix(xm, "l2", "alpha", al2_pat, al2_val, , )
## 
## 'alpha' matrix does not exist and will be added.
##  Added `alpha` matrix.

Add Across-level parameter matrices to the model

Pairs of levels that share parent-child relationship have regression relationships. xxmBetweenMatrix() is used to add corresponding parameter matrices connecting the two levels.

For each parameter matrix, the function adds the three matrices constructed earlier:

xm <- xxmBetweenMatrix(xm, "l2", "l1", "beta", be12_pat, be12_val, be12_label, 
    )
## 
## 'beta' matrix does not exist and will be added.
##  Added `beta` matrix.

Estimate model parameters

Estimation process is initiated by xxmRun(). If all goes well, a quick printed summary of results is produced.

xx <- xxmRun(xm)
## Model specification seems reasonable.
## Validity: Checked
## Matrix Constructed
## ------------------------------------------------------------------------------
## Estimating model parameters
## ------------------------------------------------------------------------------
##                3296.1681048786 
##                3262.2204916026 
##                3224.2886019777 
##                3177.4909945365 
##                3162.6131868125 
##                3123.3414244986 
##                3059.2262312280 
##                2988.2178101527 
##                2972.9862362597 
##                2916.7072744438 
##                2878.7854717227 
##                2857.4802951707 
##                2817.5124833065 
##                2796.4374947877 
##                2747.0723309565 
##                2719.1414760471 
##                2665.0802042334 
##                2636.8943854765 
##                2611.7553331240 
##                2588.9023333240 
##                2514.0473698463 
##                2503.4578768566 
##                2500.7205999010 
##                2492.3772278300 
##                2480.6071213656 
##                2478.3038537411 
##                2477.8352425263 
##                2477.7957432990 
##                2477.7406744829 
##                2477.7391546315 
##                2477.7385137178 
##                2477.7384218124 
##                2477.7383557121 
##                2477.7383436330 
##                2477.7383429881 
##                2477.7383428371 
##                2477.7383428224 
## Model converged normally
## ------------------------------------------------------------------------------
## *
##  1:                                  l1_theta_1_1 ::      1.126 [     0.000,      0.000]
## 
##  2:                                  l1_theta_2_2 ::      1.025 [     0.000,      0.000]
## 
##  3:                                  l1_theta_3_3 ::      1.197 [     0.000,      0.000]
## 
##  4:                                  l1_theta_4_4 ::      1.431 [     0.000,      0.000]
## 
##  5:                                 l1_lambda_1_1 ::      0.917 [     0.000,      0.000]
## 
##  6:                                 l1_lambda_2_1 ::      0.839 [     0.000,      0.000]
## 
##  7:                                 l1_lambda_3_1 ::      0.690 [     0.000,      0.000]
## 
##  8:                                 l1_lambda_4_1 ::      0.592 [     0.000,      0.000]
## 
##  9:                                    l2_psi_1_1 ::      0.836 [     0.000,      0.000]
## 
## 10:                                    l2_psi_1_2 ::     -0.293 [     0.000,      0.000]
## 
## 11:                                    l2_psi_2_2 ::      0.411 [     0.000,      0.000]
## 
## 12:                                  l2_alpha_1_1 ::      0.869 [     0.000,      0.000]
## 
## 13:                                  l2_alpha_2_1 ::      0.303 [     0.000,      0.000]
## 
## ------------------------------------------------------------------------------

Estimate profile-likelihood confidence intervals

Once parameters are estimated, confidence inetrvals are estimated by invoking xxmCI() . Depending on the the number of observations and the complexity of the dependence structure xxmCI() may take very long. xxMCI() displays a table of parameter estimates and CIS.

View results

A summary of results may be retrived as an R list by a call to xxmSummary()

s <- xxmSummary(xm)
s
## $deviance
## [1] 2478
## 
## $estimates
##    child parent  to from         label estimate     low     high
## 1     l1     l1  y1   y1  l1_theta_1_1   1.1265  0.7536  1.59450
## 3     l1     l1  y2   y2  l1_theta_2_2   1.0250  0.7113  1.41107
## 6     l1     l1  y3   y3  l1_theta_3_3   1.1973  0.9047  1.56805
## 10    l1     l1  y4   y4  l1_theta_4_4   1.4312  1.1165  1.83732
## 12    l1     l1  y1   fw l1_lambda_1_1   0.9172  0.7165  1.13597
## 13    l1     l1  y2   fw l1_lambda_2_1   0.8393  0.6662  1.02581
## 14    l1     l1  y3   fw l1_lambda_3_1   0.6896  0.5317  0.86462
## 15    l1     l1  y4   fw l1_lambda_4_1   0.5921  0.4390  0.76361
## 18    l2     l2 int  int    l2_psi_1_1   0.8358  0.3884  1.71909
## 19    l2     l2 int  slp    l2_psi_1_2  -0.2934 -0.7446 -0.03198
## 20    l2     l2 slp  slp    l2_psi_2_2   0.4109  0.1331  1.00192
## 21    l2     l2 int  One  l2_alpha_1_1   0.8686  0.5422  1.25838
## 22    l2     l2 slp  One  l2_alpha_2_1   0.3032  0.0453  0.57851

Free moodel object

xxM model object may hog a large amount of RAM outside of R's memory. This memory will automatically be released, when R's workspace is cleared by a call to rm(list=ls()) or at the end of the R session. Alternatively, xxmFree() may be called to release memory.

xxmFree(xm)
##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 260118  7.0     407500 10.9   350000  9.4
## Vcells 232212  1.8     786432  6.0   786429  6.0

# rm(list=ls()) detach('package:xxm', unload=TRUE)