Table of Contents
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.
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:
library(lme4)
## Loading required package: lattice
## Loading required package: Matrix
### 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)
## '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 ...
str(sample)
## 'data.frame': 30 obs. of 1 variable:
## $ sample: int 1 2 3 4 5 6 7 8 9 10 ...
str(batch)
## '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.
## 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.
pastes <- xxmModel(levels = c("response", "sample", "batch"))
## 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.
## Level-1 response
pastes <- xxmSubmodel(model = pastes, level = "response", parents = c("sample",
"batch"), ys = "strength", xs = , etas = , data = response)
## Submodel for level `response` was created.
## Level-2 sample
pastes <- xxmSubmodel(model = pastes, level = "sample", parents = , ys = , xs = ,
etas = "inter_samp", data = sample)
## Submodel for level `sample` was created.
## Level-3 sample
pastes <- xxmSubmodel(model = pastes, level = "batch", parents = , ys = , xs = ,
etas = "inter_batch", data = batch)
## 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)
# Level-1 residual variance
pastes <- xxmWithinMatrix(model = pastes, level = "response", type = "theta",
pattern = theta_pattern, value = theta_value, label = )
##
## 'theta' matrix does not exist and will be added.
## Added `theta` matrix.
# Level-1 inetrcept
pastes <- xxmWithinMatrix(model = pastes, level = "response", type = "nu", pattern = nu_pattern,
value = nu_value, label = )
##
## 'nu' matrix does not exist and will be added.
## Added `nu` matrix.
# Level-2 latent variance
pastes <- xxmWithinMatrix(model = pastes, level = "sample", type = "psi", pattern = psi_pattern,
value = psi_value, label = )
##
## 'psi' matrix does not exist and will be added.
## Added `psi` matrix.
# Level-3 latent variance
pastes <- xxmWithinMatrix(model = pastes, level = "batch", type = "psi", pattern = psi_pattern,
value = psi_value, label = )
##
## '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)
# 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 = )
##
## 'lambda' matrix does not exist and will be added.
## Added `lambda` matrix.
# 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 = )
##
## '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.
pastes <- xxmRun(pastes)
## ------------------------------------------------------------------------------
## 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.
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.
xxmSummary(pastes)
## $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.
pastes <- xxmFree(pastes)
rm(list = ls())
Proc Mixed
PROC MIXED data=pastes method=ML covtest;
CLASS sample batch;
MODEL strength = /s;
RANDOM intercept /subject=sample ;
RANDOM intercept /subject=batch ;
RUN;
R: lmer
library(lme4)
(fm01 <- lmer(strength ~ 1 + (1 | sample) + (1 | batch), REML = FALSE, Pastes))
## 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