# 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

1. response
2. sample
3. 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

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

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}$

(insert diagram)

# Code Listing

“xxM”“SAS:“R:

## 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

### 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: 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. ## 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: 1. 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)$$. 2. 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
##  248
##
## $fit$nParameters
##  4
##
## $fit$nObservations
##  60
##
## $fit$aic
##  256
##
## $fit$bic
##  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