Table of Contents
- 1 Pastes: Three Level Random-Intercepts Model
- 2 xxM Model Matrices
- 3 Path Diagram
- 4 Code Listing
Pastes: Three Level Random-Intercepts Model
Pastes example included with lme4 involves a random-intercepts model of strength of pastes in different samples nested within batches. The three levels are
- response
- sample
- batch
\[ y_{i}^{1} =\nu^1 + \Lambda_{i,j}^{1,2} \times \eta_j^{2} + \Lambda_{i,j}^{1,3} \times \eta_j^{3}+ e_{i}, \]
\[ e^{1} \sim N(0, \theta^{1,1}_{1,1}) \] \[ \eta^{2} \sim N(0, \psi^{2,2}_{1,1}) \] \[ \eta^{3} \sim N(0, \psi^{3,3}_{1,1}) \]
xxM Model Matrices
Level-1 (Response): Within matrices
With a single dependent variable there are two parameters level-1: (a) residual variance \( (\theta_{1,1}) \), and (b) intercept \( \nu_{1} \)
Residual Covariance Matrix
\Theta^{1,1} =
\begin{bmatrix}
\theta^{1,1}_{1,1}
\end{bmatrix}
\]
Observed interept
\nu^{1} =
\begin{bmatrix}
\nu^{1}_{1}
\end{bmatrix}
\]
Level-2 (Sample): Within Matrices
At level-2, we have a single latent variable: random-intercept of strength, with a single parameter latent variance.
Latent Factor Covariance Matrix
The latent covariance matrix is a (1×1) matrix with single variance of the intercept:
\Psi^{2,2}=
\begin{bmatrix}
\psi_{1,1}^{2,2}
\end{bmatrix}
\]
Level-3 (Batch): Within Matrices
Likewise, at level-3, we have a single latent variable: random-intercept of strength, with a single parameter latent variance.
Latent Factor Covariance Matrix
The latent covariance matrix is a (1×1) matrix with single variance of the intercept:
\Psi^{3,3}=
\begin{bmatrix}
\psi_{1,1}^{3,3}
\end{bmatrix}
\]
Across level matrices: Sample to Response
Across-level factor-loading matrix
The factor-loading matrix \( (\Lambda^{1,2}) \) has a single element fixed to 1.0.
\Lambda^{1,2} =
\begin{bmatrix}
1.0_{i,j}^{1,2}
\end{bmatrix}
\]
Across level matrices: Batch to Response
Across-level factor-loading matrix
The factor-loading matrix \( (\Lambda^{1,3}) \) has a single element fixed to 1.0.
\Lambda^{1,3} =
\begin{bmatrix}
1.0_{i,j}^{1,3}
\end{bmatrix}
\]
Path Diagram
(insert diagram)
Code Listing
xxM
Load xxM
xxM library needs to be loaded first.
1 2 |
library(xxm) |
Prepare R datasets
For this analysis, we use data packaged with lme4.
For analysis with xxM separate response, batch and sample datasets are created as follows:
1 2 |
library(lme4) |
1 2 3 |
## Loading required package: lattice ## Loading required package: Matrix |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
### Creating datasets for xxM Lowest level is called response response <- Pastes ## First column must include unique IDs for the level Next p columns must ## include IDs of the parent levels (batch and strength) response$response <- 1:nrow(response) response <- response[c("response", "sample", "batch", "strength")] convert.magic1 <- function(obj, types) { out <- lapply(1:length(obj), FUN = function(i) { FUN1 <- switch(types[i], character = as.character, numeric = as.numeric, factor = as.factor, integer = as.integer) FUN1(obj[, i]) }) names(out) <- colnames(obj) as.data.frame(out) } ## ID columns must be integer, ys and xs must be numeric type <- c(rep("integer", 3), "numeric") response <- convert.magic1(response, type) ## Each level must have its dataset with a column of unique IDs. batch <- data.frame(batch = unique(response$batch)) sample <- data.frame(sample = unique(response$sample)) # examine datasets str(response) |
1 2 3 4 5 6 |
## 'data.frame': 60 obs. of 4 variables: ## $ response: int 1 2 3 4 5 6 7 8 9 10 ... ## $ sample : int 1 1 2 2 3 3 4 4 5 5 ... ## $ batch : int 1 1 1 1 1 1 2 2 2 2 ... ## $ strength: num 62.8 62.6 60.1 62.3 62.7 63.1 60 61.4 57.5 56.9 ... |
1 2 |
str(sample) |
1 2 3 |
## 'data.frame': 30 obs. of 1 variable: ## $ sample: int 1 2 3 4 5 6 7 8 9 10 ... |
1 2 |
str(batch) |
1 2 3 |
## 'data.frame': 10 obs. of 1 variable: ## $ batch: int 1 2 3 4 5 6 7 8 9 10 |
response dataset includes the following columns:
- response Unique response IDs (integer)
- sample IDs for the parent level (integer)
- batch IDs for the parent level (integer)
- strength the dependent varaible (numeric)
sample dataset includes the following columns:
- sample Unique batch IDs (integer)
batch dataset includes the following columns:
- batch Unique batch IDs (integer)
Construct R-matrices
For each parameter matrix, construct three related matrices:
- pattern matrix: A matrix indicating free or fixed parameters.
- value matrix: with start or fixed values for corresponding parameters.
- label matrix: with user friendly label for each parameter. label matrix is optional.
1 2 3 4 5 6 7 8 9 10 11 12 13 |
## Level-1 residual variance theta_pattern <- matrix(1, nrow = 1, ncol = 1) theta_value <- matrix(2.5, nrow = 1, ncol = 1) ## Level-1 intercept nu_pattern <- matrix(1, nrow = 1, ncol = 1) nu_value <- matrix(60.05, nrow = 1, ncol = 1) ## Level-2 and 3 latent variance psi_pattern <- matrix(1, nrow = 1, ncol = 1) psi_value <- matrix(1, nrow = 1, ncol = 1) ## Across-level factor loading matrix (fixed to 1.0) lambda_pattern <- matrix(0, nrow = 1, ncol = 1) lambda_value <- matrix(1, nrow = 1, ncol = 1) |
Construct model
xxmModel()
is used to declare level names. The function returns a model object that is passed as a parameter to subsequent statements. Variable name for the return value can be anything.
1 2 |
pastes <- xxmModel(levels = c("response", "sample", "batch")) |
1 2 |
## A new model was created. |
Add submodels
For each declared level xxmSubmodel()
is invoked to add corresponding submodel to the model object. The function adds three types of information to the model object:
- parents declares a list of all parents of the current level.
- Level with the independent variable is the parent level.
- Level with the dependent variable is the child level.
- variables declares names of observed dependent (ys), observed independent (xs) and latent variables (etas) for the level.
- data R data object for the current level.
1 2 3 4 |
## Level-1 response pastes <- xxmSubmodel(model = pastes, level = "response", parents = c("sample", "batch"), ys = "strength", xs = , etas = , data = response) |
1 2 |
## Submodel for level `response` was created. |
1 2 3 4 |
## Level-2 sample pastes <- xxmSubmodel(model = pastes, level = "sample", parents = , ys = , xs = , etas = "inter_samp", data = sample) |
1 2 |
## Submodel for level `sample` was created. |
1 2 3 4 |
## Level-3 sample pastes <- xxmSubmodel(model = pastes, level = "batch", parents = , ys = , xs = , etas = "inter_batch", data = batch) |
1 2 |
## Submodel for level `batch` was created. |
Add within-level matrices
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:
- pattern
- value
- label (optional)
1 2 3 4 |
# Level-1 residual variance pastes <- xxmWithinMatrix(model = pastes, level = "response", type = "theta", pattern = theta_pattern, value = theta_value, label = ) |
1 2 3 4 |
## ## 'theta' matrix does not exist and will be added. ## Added `theta` matrix. |
1 2 3 4 |
# Level-1 inetrcept pastes <- xxmWithinMatrix(model = pastes, level = "response", type = "nu", pattern = nu_pattern, value = nu_value, label = ) |
1 2 3 4 |
## ## 'nu' matrix does not exist and will be added. ## Added `nu` matrix. |
1 2 3 4 5 |
# Level-2 latent variance pastes <- xxmWithinMatrix(model = pastes, level = "sample", type = "psi", pattern = psi_pattern, value = psi_value, label = ) |
1 2 3 4 |
## ## 'psi' matrix does not exist and will be added. ## Added `psi` matrix. |
1 2 3 4 5 |
# Level-3 latent variance pastes <- xxmWithinMatrix(model = pastes, level = "batch", type = "psi", pattern = psi_pattern, value = psi_value, label = ) |
1 2 3 4 |
## ## 'psi' matrix does not exist and will be added. ## Added `psi` matrix. |
Add across-level matrices
Pairs of levels that share parent-child relationship have regression relationships. xxmBetweenMatrix()
is used to add corresponding rergession matrices connecting the two levels.
- Level with the independent variable is the parent level.
- Level with the dependent variable is the child level.
For each parameter matrix, the function adds the three matrices constructed earlier:
- pattern
- value
- label (optional)
1 2 3 4 |
# Level-2 to Level-1 factor-loading (fixed to 1.0) pastes <- xxmBetweenMatrix(model = pastes, parent = "sample", child = "response", type = "lambda", pattern = lambda_pattern, value = lambda_value, label = ) |
1 2 3 4 |
## ## 'lambda' matrix does not exist and will be added. ## Added `lambda` matrix. |
1 2 3 4 |
# Level-3 to Level-1 factor-loading (fixed to 1.0) pastes <- xxmBetweenMatrix(model = pastes, parent = "batch", child = "response", type = "lambda", pattern = lambda_pattern, value = lambda_value, label = ) |
1 2 3 4 |
## ## 'lambda' matrix does not exist and will be added. ## Added `lambda` matrix. |
Estimate model parameters
Estimation process is initiated by xxmRun()
. If all goes well, a q&d summary of the results is printed.
1 2 |
pastes <- xxmRun(pastes) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
## ------------------------------------------------------------------------------ ## Estimating model parameters ## ------------------------------------------------------------------------------ ## 290.9778594217 ## 285.7873333143 ## 280.7767590084 ## 276.0149840654 ## 271.5546853088 ## 267.4575220947 ## 263.8016391893 ## 261.2740756488 ## 258.5290865871 ## 257.0471463734 ## 255.9187866656 ## 255.9125449731 ## 255.7986675263 ## 255.7447647880 ## 255.4424804249 ## 254.5251352801 ## 253.1520719699 ## 250.8132821825 ## 248.3509689557 ## 247.9987917603 ## 247.9957451719 ## 247.9945071529 ## 247.9944682815 ## 247.9944658952 ## 247.9944658624 ## Model converged normally ## nParms: 4 ## ------------------------------------------------------------------------------ ## * ## 1: response_theta_1_1 :: 0.678 [ 0.000, 0.000] ## ## 2: response_nu_1_1 :: 60.053 [ 0.000, 0.000] ## ## 3: sample_psi_1_1 :: 8.434 [ 0.000, 0.000] ## ## 4: batch_psi_1_1 :: 1.199 [ 0.000, 0.000] ## ## ------------------------------------------------------------------------------ |
Estimate profile-likelihood confidence intervals
Once parameters are estimated, confidence intervals are estimated by invoking xxmCI()
. Depending on the the number of observations and the complexity of the model, xxmCI()
may take a long time to compute. xxmCI()
also prints a summary of parameter estimates and CIS.
1 2 |
pastes <- xxmCI(pastes) |
View results
A summary of results may be retrived as an R list by a call to xxmSummary()
. The returned list has two elements:
- fit is a list with five elements:
- deviance is \( -2 Log Likelihood \) for the maximum likelihood fit function.
- nParameters is the total number of unique parameters.
- nObservations is the total number of observations across all levels.
- aic is Akaike’s Information Criterion or AIC computed as \( -2ll + 2*p \).
- bic is Bayesian Information Criterion or BIC computed as \( -2ll + p*\log(n) \).
- estimates is a single table of free parameter estimates
All xxM parameters have superscripts {child, parent} and subscripts {to, from}. xxM adds a descriptive parameter label if one is not already provided by the user.
1 2 |
xxmSummary(pastes) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
## $fit ## $fit$deviance ## [1] 248 ## ## $fit$nParameters ## [1] 4 ## ## $fit$nObservations ## [1] 60 ## ## $fit$aic ## [1] 256 ## ## $fit$bic ## [1] 264.4 ## ## ## $estimates ## child parent to from label estimate ## 1 response response strength strength response_theta_1_1 0.678 ## 2 response response strength One response_nu_1_1 60.053 ## 4 sample sample inter_samp inter_samp sample_psi_1_1 8.434 ## 6 batch batch inter_batch inter_batch batch_psi_1_1 1.199 ## low high ## 1 0.4251 1.178 ## 2 58.6636 61.443 ## 4 4.6567 16.433 ## 6 0.0010 8.682 |
Free model object
xxM
model object may hog a significant amount of RAM outside of R’s workspace. This memory will automatically be released, when the workspace is cleared by a call to rm(list=ls())
or at the end of the R session. Alternatively, it is recommended that xxmFree()
may be called to release the memory.
1 2 3 |
pastes <- xxmFree(pastes) rm(list = ls()) |
Proc Mixed
1 2 3 4 5 6 7 |
PROC MIXED data=pastes method=ML covtest; CLASS sample batch; MODEL strength = /s; RANDOM intercept /subject=sample ; RANDOM intercept /subject=batch ; RUN; |
R: lmer
1 2 3 |
library(lme4) (fm01 <- lmer(strength ~ 1 + (1 | sample) + (1 | batch), REML = FALSE, Pastes)) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
## Linear mixed model fit by maximum likelihood ['lmerMod'] ## Formula: strength ~ 1 + (1 | sample) + (1 | batch) ## Data: Pastes ## AIC BIC logLik deviance ## 256.0 264.4 -124.0 248.0 ## Random effects: ## Groups Name Std.Dev. ## sample (Intercept) 2.904 ## batch (Intercept) 1.095 ## Residual 0.823 ## Number of obs: 60, groups: sample, 30; batch, 10 ## Fixed Effects: ## (Intercept) ## 60.1 |