# Penicillin: Random-Intercept, Cross-Classified model

The penicillin example included with lme4 is designed to examine variability between multiple samples of penicillin. In this example the outcome variable,diameter, refers to the diameter of the disc in a petri dish which is successfully treated by Penicillin. There are two cross-classified parent levels: the plate, or Petri dish, used and the sample of Penicillin. Level-1 is response, level-2 is called Plate, and level-3 is Sample.

$y_{i(u,v)} = \nu + 1\times \eta_u^{p}+ 1\times\eta_v^{s} + e_{i(u,v)},$

$e \sim N(0, \Theta^{1,1})$ $\eta^{p} \sim N(0, \Psi^{p})$ $\eta^{s} \sim N(0, \Psi^{s})$

# Path Diagram # xxM representation

$y^{1}_i = \nu^1 + 1_{i,j}^{1,2} \times \eta_j^{2}+ 1_{i,k}^{1,3} \times\eta_k^{3} + e_{i}^1$

$e^1 \sim N(0, \Theta_{1,1})$ $\eta^{2} \sim N(0, \Psi^{2,2})$ $\eta^{3} \sim N(0, \Psi^{3.3})$

Note:

• Subscripts ${i, j, k}$ are not really necessary.
• Coefficients $1_{i,j}^{1,2}$ and $1_{i,k}^{1,3}$ are included for clarity. These are just fixed factor-loadings in the model.

A general multivariate cross-classified model may be simply represented as:

$y^{1} = \nu^1 + \Lambda^{1,2} \times \eta^{2}+ \Lambda^{1,3} \times\eta^{3} + e^1$

# xxM Model Matrices

## Level-1 (Response): Within matrices

### Residual Covariance Matrix

We have a single dependent variable and hence a single parameter level-1 residual variance $$(\Theta_{1,1})$$ at level-1. The residual covariance or theta matrix is a (1×1) matrix:

$\Theta^{1,1} = \begin{bmatrix} \theta^{1,1}_{1,1} \end{bmatrix}$

### Observed Variable Intercept

In this model, plate and sample both reflect deviations from, but do not capture, the grand mean (intercept) of the dependent variable diameter. As a result, we must estimate a fixed effect for the observed intercept. The estimate is a 1×1 matrix with a single element:

$\nu^{1} = \begin{bmatrix} \nu_{1}^{1} \end{bmatrix}$

## Level-2 (Plate): Within Matrices

At the Plate level, we have a single latent variable: the random-intercept of plate. Since this latent variable has a mean of 0, it is unneccesary to estimate the parameter. As a result, the only estimated parameter within the Plate level is the latent factor covariance matrix.

### 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 (Sample): Within Matrices

At the Sample level, we have a single latent variable: the random-intercept of sample. As with the Plate level, this latent variable has a mean of 0 and it is unneccesary to estimate the parameter. The only estimated parameter within the Sample level is the latent factor covariance matrix.

### 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: Plate 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,u}^{1,2} \end{bmatrix}$

## Across level matrices: Sample 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,v}^{1,3} \end{bmatrix}$

# 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 diameter, plate, and sample datasets are created as follows:

library(lme4)

## Loading required package: lattice

### Creating datasets for xxM
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)
}
pen <- Penicillin

# assign proper names
names(pen) <- c("diameter", "plate", "sample")

# create a three-level dataset with response as level-1, plate as level-2,
# and sample as level-2. create an ID column for the 'response' level with
# the same column name
pen$response <- 1:nrow(pen) # Ensure that the first column is the ID column (for)response) response <- pen[c("response", "plate", "sample", "diameter")] # ID columns (response, plate, and sample) must be integer, Outcome must be # numeric. response <- convert.magic1(response, c("integer", "integer", "integer", "numeric")) # create level-2 dataset. sample <- data.frame(sample = unique(response$sample))
plate <- data.frame(plate = unique(response$plate)) str(response)  ## 'data.frame': 144 obs. of 4 variables: ##$ response: int  1 2 3 4 5 6 7 8 9 10 ...
##  $plate : int 1 1 1 1 1 1 2 2 2 2 ... ##$ sample  : int  1 2 3 4 5 6 1 2 3 4 ...
##  $diameter: num 27 23 26 23 23 21 27 23 26 23 ...  str(sample)  ## 'data.frame': 6 obs. of 1 variable: ##$ sample: int  1 2 3 4 5 6

str(plate)

## 'data.frame':    24 obs. of  1 variable:
##  $plate: int 1 2 3 4 5 6 7 8 9 10 ...  response dataset includes the following columns: • response Unique response IDs (integer) • plate IDs for the parent level (integer) • sample IDs for the parent level (integer) • diameter the dependent varaible (numeric) plate dataset includes the following columns: • plate Unique plate IDs (integer) sample dataset includes the following columns: • sample Unique sample 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 (response) residual variance theta_pattern <- matrix(1, nrow = 1, ncol = 1) theta_value <- matrix(2.5, nrow = 1, ncol = 1) # level-2 (plate) latent variance psi_pattern <- matrix(1, nrow = 1, ncol = 1) psi_value <- matrix(1, nrow = 1, ncol = 1) # level-3 (sample) latent mean nu_pattern <- matrix(1, nrow = 1, ncol = 1) nu_value <- matrix(22.97, nrow = 1, ncol = 1) # plate/sample -> response: Factor-loading 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. pen <- xxmModel(levels = c("response", "plate", "sample"))  ## 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. • Levels with the independent variables are the parent levels. • 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. ### Submodel: Response pen <- xxmSubmodel(model = pen, level = "response", parents = c("plate", "sample"), ys = "diameter", xs = , etas = , data = response)  ## Submodel for level response was created.   ### Submodel: Plate pen <- xxmSubmodel(model = pen, level = "plate", parents = , ys = , xs = , etas = "plate", data = plate)  ## Submodel for level plate was created.   ### Submodel: Sample pen <- xxmSubmodel(model = pen, level = "sample", parents = , ys = , xs = , etas = "sample", data = sample)  ## Submodel for level sample 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) ### Within-Matrix: Level :: Response, Type :: observed residual-covariance pen <- xxmWithinMatrix(model = pen, level = "response", type = "theta", pattern = theta_pattern, value = theta_value, label = )  ## ## 'theta' matrix does not exist and will be added. ## Added theta matrix.   ### Within-Matrix: Level :: Response, Type :: observed variable intercept pen <- xxmWithinMatrix(model = pen, level = "response", type = "nu", pattern = nu_pattern, value = nu_value, label = )  ## ## 'nu' matrix does not exist and will be added. ## Added nu matrix.   ### Within-Matrix: Level :: Plate, Type :: latent factor variance pen <- xxmWithinMatrix(model = pen, level = "plate", type = "psi", pattern = psi_pattern, value = psi_value, label = )  ## ## 'psi' matrix does not exist and will be added. ## Added psi matrix.   ### Within-Matrix: Level :: Sample, Type ::latent factor variance pen <- xxmWithinMatrix(model = pen, level = "sample", 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) ### Between-Matrix: Child :: Response, Parent :: Plate, Type :: factor-loading pen <- xxmBetweenMatrix(model = pen, parent = "plate", child = "response", type = "lambda", pattern = lambda_pattern, value = lambda_value, label = )  ## ## 'lambda' matrix does not exist and will be added. ## Added lambda matrix.   ### Between-Matrix: Child :: Response, Parent :: Sample, Type :: ### factor-loading pen <- xxmBetweenMatrix(model = pen, 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.  ### Estimate model parameters Estimation process is initiated by xxmRun(). If all goes well, a q&d summary of the results is printed. pen <- xxmRun(pen)  ## ------------------------------------------------------------------------------ ## Estimating model parameters ## ------------------------------------------------------------------------------ ## 381.7896401031 ## 383.0593366241 ## 372.0392133728 ## 363.0400133360 ## 357.6642273898 ## 353.9375212264 ## 351.5101869798 ## 347.1445555090 ## 345.7723126374 ## 345.1403633189 ## 340.2338642579 ## 339.0994461776 ## 338.0024813307 ## 336.9734549987 ## 336.0272283759 ## 335.1795891216 ## 333.8532609313 ## 333.2730156349 ## 332.2795893391 ## 332.2264465732 ## 332.2025221517 ## 332.1884981905 ## 332.1883935085 ## 332.1883581231 ## 332.1883488233 ## 332.1883486692 ## Model converged normally ## nParms: 4 ## ------------------------------------------------------------------------------ ## * ## 1: response_theta_1_1 :: 0.302 [ 0.000, 0.000] ## ## 2: response_nu_1_1 :: 22.972 [ 0.000, 0.000] ## ## 3: plate_psi_1_1 :: 0.715 [ 0.000, 0.000] ## ## 4: sample_psi_1_1 :: 3.135 [ 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. pen <- xxmCI(pen)  ### 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(pen)  ##$fit
## $fit$deviance
##  332.2
##
## $fit$nParameters
##  4
##
## $fit$nObservations
##  144
##
## $fit$aic
##  340.2
##
## $fit$bic
##  352.1
##
##
## \$estimates
##      child   parent       to     from              label estimate     low
## 1 response response diameter diameter response_theta_1_1   0.3024  0.2360
## 2 response response diameter      One    response_nu_1_1  22.9722 21.2666
## 4    plate    plate    plate    plate      plate_psi_1_1   0.7150  0.4014
## 6   sample   sample   sample   sample     sample_psi_1_1   3.1352  1.2008
##      high
## 1  0.3962
## 2 24.6778
## 4  1.3974
## 6 12.6474


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

pen <- xxmFree(pen)
rm(list = ls()) ## Proc Mixed

PROC MIXED data=rdata method=ML covtest;
CLASS plate sample;
MODEL diameter = /s;
RANDOM intercept /subject = plate ;
RANDOM intercept /subject = sample;
RUN;


## R: lmer

library(lme4)
# LME4 model
(fm03 <- lmer(diameter ~ 1 + (1 | plate) + (1 | sample), Penicillin))

## Linear mixed model fit by REML ['lmerMod']
## Formula: diameter ~ 1 + (1 | plate) + (1 | sample)
##    Data: Penicillin
## REML criterion at convergence: 330.9
## Random effects:
##  Groups   Name        Std.Dev.
##  plate    (Intercept) 0.847
##  sample   (Intercept) 1.932
##  Residual             0.550
## Number of obs: 144, groups: plate, 24; sample, 6
## Fixed Effects:
## (Intercept)
##          23

(confint(fm03))

## Computing profile confidence intervals ...

##               2.5 %  97.5 %
## .sig01       0.6336  1.1821
## .sig02       1.0958  3.5563
## .sigma       0.4858  0.6295
## (Intercept) 21.2666 24.6778


# Obtaining the samen deviance as in xxM
(fm03.1 <- lmer(diameter ~ 1 + (1 | plate) + (1 | sample), REML = FALSE, Penicillin))

## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: diameter ~ 1 + (1 | plate) + (1 | sample)
##    Data: Penicillin
##      AIC      BIC   logLik deviance
##    340.2    352.1   -166.1    332.2
## Random effects:
##  Groups   Name        Std.Dev.
##  plate    (Intercept) 0.846
##  sample   (Intercept) 1.771
##  Residual             0.550
## Number of obs: 144, groups: plate, 24; sample, 6
## Fixed Effects:
## (Intercept)
##          23

(confint(fm03.1))

## Computing profile confidence intervals ...

##               2.5 %  97.5 %
## .sig01       0.6336  1.1821
## .sig02       1.0958  3.5563
## .sigma       0.4858  0.6295
## (Intercept) 21.2666 24.6778