Linear mixed model covariance decomposition with random slopes- lme4

Load data

data(Orthodont, package="nlme")
Data <- Orthodont

Using lmm

library(lme4)

## Loading required package: Matrix

fit.lmer.slope <- lmer(distance ~ age + Sex  +(1+ age | Subject), data = Data) 
summary(fit.lmer.slope)

## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + Sex + (1 + age | Subject)
##    Data: Data
## 
## REML criterion at convergence: 435.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0815 -0.4568  0.0155  0.4471  3.8944 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Subject  (Intercept) 7.82249  2.7969        
##           age         0.05126  0.2264   -0.77
##  Residual             1.71621  1.3100        
## Number of obs: 108, groups:  Subject, 27
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 17.63518    0.88623  19.899
## age          0.66019    0.07125   9.265
## SexFemale   -2.14544    0.75746  -2.832
## 
## Correlation of Fixed Effects:
##           (Intr) age   
## age       -0.838       
## SexFemale -0.348  0.000

default is compound symmetry

Get x,y,z

X=(getME(fit.lmer.slope, "X"))
head(X)

##   (Intercept) age SexFemale
## 1           1   8         0
## 2           1  10         0
## 3           1  12         0
## 4           1  14         0
## 5           1   8         0
## 6           1  10         0

y=getME(fit.lmer.slope, "y")
head(y)

## [1] 26.0 25.0 29.0 31.0 21.5 22.5

Z <- getME(fit.lmer.slope, "Z")
head(Z)

## 6 x 54 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 54 column names 'M16', 'M16', 'M05' ... ]]

##                                                                                
## 1 . . . . .  . . . . . . . . . . . . . . . . . . . . . . . 1  8 . . . . . . . .
## 2 . . . . .  . . . . . . . . . . . . . . . . . . . . . . . 1 10 . . . . . . . .
## 3 . . . . .  . . . . . . . . . . . . . . . . . . . . . . . 1 12 . . . . . . . .
## 4 . . . . .  . . . . . . . . . . . . . . . . . . . . . . . 1 14 . . . . . . . .
## 5 . . . . 1  8 . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .
## 6 . . . . 1 10 . . . . . . . . . . . . . . . . . . . . . . .  . . . . . . . . .
##                                  
## 1 . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . .

dim(Z)

## [1] 108  54

Get fixed effect and random effect coefficients

bhat <- getME(fit.lmer.slope, "fixef") #getME(fit.lmer, "beta")  # fixef(fit.lmer)
bhat

## (Intercept)         age   SexFemale 
##  17.6351805   0.6601852  -2.1454431

uhat <- ranef(fit.lmer.slope) 
uhat

## $Subject
##     (Intercept)         age
## M16 -0.96698456 -0.06544230
## M05 -2.15621081  0.04447130
## M02 -1.58143839  0.02194874
## M11  0.38737823 -0.13961777
## M07 -1.40307472  0.03606404
## M08  0.37415109 -0.11799495
## M03 -0.83491587  0.02435289
## M12 -1.82593774  0.11594755
## M13 -5.59182087  0.46400729
## M14  0.51944692 -0.04982259
## M09 -1.07941522  0.11835170
## M15 -1.11909663  0.18322017
## M06  1.03469722  0.02495756
## M04  3.20171821 -0.15492789
## M01  0.96194797  0.14388309
## M10  3.04960613  0.09373458
## F10 -2.31273350 -0.13319695
## F09  0.32324282 -0.16262237
## F06 -0.07316594 -0.12598450
## F01  0.11181130 -0.12268061
## F05  1.43310624 -0.14279903
## F07  0.62044804 -0.03708906
## F02 -0.37057384  0.05450561
## F08  2.38444671 -0.16952522
## F03 -0.01384650  0.08273621
## F04  2.30508388 -0.03978828
## F11  2.62212981  0.05331079
## 
## with conditional variances for "Subject"

Get random effect covariance structure (covariance matrix expandation)

notice z matrix structure

vc <- VarCorr(fit.lmer.slope)
Lambda_new <- bdiag(replicate(length(levels(Data$Subject)), vc[["Subject"]], simplify = FALSE))
head(Lambda_new)

## 6 x 54 sparse Matrix of class "dgCMatrix"
##                                                                               
## [1,]  7.822487 -0.48495997  .         .           .         .          . . . .
## [2,] -0.484960  0.05126495  .         .           .         .          . . . .
## [3,]  .         .           7.822487 -0.48495997  .         .          . . . .
## [4,]  .         .          -0.484960  0.05126495  .         .          . . . .
## [5,]  .         .           .         .           7.822487 -0.48495997 . . . .
## [6,]  .         .           .         .          -0.484960  0.05126495 . . . .
##                                                                               
## [1,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [2,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [3,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [4,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [5,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [6,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##                   
## [1,] . . . . . . .
## [2,] . . . . . . .
## [3,] . . . . . . .
## [4,] . . . . . . .
## [5,] . . . . . . .
## [6,] . . . . . . .

dim(Lambda_new)

## [1] 54 54

“theta”: random-effects parameter estimates: these are parameterized as the relative Cholesky factors of each random effect term
use other cov structure by nlme
lme4 does not allow to specify it, and it can only fit LMMs with independent residual errors. However, the range of available variance-covariance matrices for the random effects are restricted to diagonal or general matrices.

Get residual variance and its identity matrix

sigma <- getME(fit.lmer.slope, "sigma")**2 
sigma

## [1] 1.71621

head(sigma*diag(nrow(Data)))

##         [,1]    [,2]    [,3]    [,4]    [,5]    [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 1.71621 0.00000 0.00000 0.00000 0.00000 0.00000    0    0    0     0     0
## [2,] 0.00000 1.71621 0.00000 0.00000 0.00000 0.00000    0    0    0     0     0
## [3,] 0.00000 0.00000 1.71621 0.00000 0.00000 0.00000    0    0    0     0     0
## [4,] 0.00000 0.00000 0.00000 1.71621 0.00000 0.00000    0    0    0     0     0
## [5,] 0.00000 0.00000 0.00000 0.00000 1.71621 0.00000    0    0    0     0     0
## [6,] 0.00000 0.00000 0.00000 0.00000 0.00000 1.71621    0    0    0     0     0
##      [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94] [,95]
## [1,]     0     0     0     0     0     0     0     0     0     0     0     0
## [2,]     0     0     0     0     0     0     0     0     0     0     0     0
## [3,]     0     0     0     0     0     0     0     0     0     0     0     0
## [4,]     0     0     0     0     0     0     0     0     0     0     0     0
## [5,]     0     0     0     0     0     0     0     0     0     0     0     0
## [6,]     0     0     0     0     0     0     0     0     0     0     0     0
##      [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104] [,105] [,106]
## [1,]     0     0     0     0      0      0      0      0      0      0      0
## [2,]     0     0     0     0      0      0      0      0      0      0      0
## [3,]     0     0     0     0      0      0      0      0      0      0      0
## [4,]     0     0     0     0      0      0      0      0      0      0      0
## [5,]     0     0     0     0      0      0      0      0      0      0      0
## [6,]     0     0     0     0      0      0      0      0      0      0      0
##      [,107] [,108]
## [1,]      0      0
## [2,]      0      0
## [3,]      0      0
## [4,]      0      0
## [5,]      0      0
## [6,]      0      0

Get y covariance matrix

VM <- Z%*%Lambda_new%*%t(Z)+sigma*diag(nrow(Data)) 
head(VM)

## 6 x 108 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 108 column names '1', '2', '3' ... ]]

##                                                                                
## 1 5.060294 3.194403 3.044722 2.895042 .        .        .        .        . . .
## 2 3.194403 4.965992 3.305161 3.360540 .        .        .        .        . . .
## 3 3.044722 3.305161 5.281810 3.826039 .        .        .        .        . . .
## 4 2.895042 3.360540 3.826039 6.007747 .        .        .        .        . . .
## 5 .        .        .        .        5.060294 3.194403 3.044722 2.895042 . . .
## 6 .        .        .        .        3.194403 4.965992 3.305161 3.360540 . . .
##                                                                                
## 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##                                                                                
## 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##                                        
## 1 . . . . . . . . . . . . . . . . . . .
## 2 . . . . . . . . . . . . . . . . . . .
## 3 . . . . . . . . . . . . . . . . . . .
## 4 . . . . . . . . . . . . . . . . . . .
## 5 . . . . . . . . . . . . . . . . . . .
## 6 . . . . . . . . . . . . . . . . . . .

dim(VM)

## [1] 108 108

Get fixed effect coefficients covariance matrix

vcov <- vcov(fit.lmer.slope) #fixed cov
vcov

## 3 x 3 Matrix of class "dpoMatrix"
##             (Intercept)           age     SexFemale
## (Intercept)  0.78540227 -5.292131e-02 -2.337501e-01
## age         -0.05292131  5.076869e-03 -5.668086e-16
## SexFemale   -0.23375009 -5.668086e-16  5.737502e-01

# computer correlation coefficients 
vcov@x[2]/prod(sqrt( diag(vcov(fit.lmer.slope))[-3] ))

## [1] -0.8380822

Get random effect coefficients covariance matrix (standard deviation)

vc <- VarCorr(fit.lmer.slope) 
vc ## default print method: standard dev and corr

##  Groups   Name        Std.Dev. Corr  
##  Subject  (Intercept) 2.79687        
##           age         0.22642  -0.766
##  Residual             1.31004

as.matrix(Matrix::bdiag(vc)) #random effect covariance matrix

##             (Intercept)         age
## (Intercept)    7.822487 -0.48495997
## age           -0.484960  0.05126495

# https://stackoverflow.com/questions/47307340/extracting-the-i-estimated-variance-covariance-matrix-of-random-effects-and-or
# https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect022.htm

Computer fixed effect coefficients

library(matlib)
inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X))%*%t(as.matrix(X))%*%inv(as.matrix(VM))%*%y

##           [,1]
## [1,] 17.635190
## [2,]  0.660189
## [3,] -2.145444

bhat

## (Intercept)         age   SexFemale 
##  17.6351805   0.6601852  -2.1454431

Computer covariance of fixed efffect coefficients

inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X))

##                                        
## [1,]  0.78540228 -0.05292131 -0.2337501
## [2,] -0.05292131  0.00507687  0.0000000
## [3,] -0.23375010  0.00000000  0.5737502

# standard error
sqrt(inv(t(as.matrix(X))%*%inv(as.matrix(VM))%*%as.matrix(X)))

## Warning in sqrt(inv(t(as.matrix(X)) %*% inv(as.matrix(VM)) %*% as.matrix(X))):
## NaNs produced

##                                   
## [1,] 0.8862292        NaN      NaN
## [2,]       NaN 0.07125216 0.000000
## [3,]       NaN 0.00000000 0.757463

# the following equals lmm summary
vcov(fit.lmer.slope)

## 3 x 3 Matrix of class "dpoMatrix"
##             (Intercept)           age     SexFemale
## (Intercept)  0.78540227 -5.292131e-02 -2.337501e-01
## age         -0.05292131  5.076869e-03 -5.668086e-16
## SexFemale   -0.23375009 -5.668086e-16  5.737502e-01

sqrt(vcov(fit.lmer.slope))

## Warning in sqrt(x@x): NaNs produced

## 3 x 3 Matrix of class "dsyMatrix"
##             (Intercept)        age SexFemale
## (Intercept)   0.8862292        NaN       NaN
## age                 NaN 0.07125215       NaN
## SexFemale           NaN        NaN  0.757463

default is compound symmetry

Computer random effect coefficients

comput_uhat <- (as.matrix(Lambda_new))%*%t(Z)%*%inv(as.matrix(VM))%*%(y-as.matrix(X)%*%(bhat))
cbind((comput_uhat@x),(uhat[["Subject"]]))

## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded

##    (comput_uhat@x) (Intercept)         age
## 1      -0.96698450 -0.96698456 -0.06544230
## 2      -0.06544230 -2.15621081  0.04447130
## 3      -2.15621077 -1.58143839  0.02194874
## 4       0.04447130  0.38737823 -0.13961777
## 5      -1.58143833 -1.40307472  0.03606404
## 6       0.02194874  0.37415109 -0.11799495
## 7       0.38737824 -0.83491587  0.02435289
## 8      -0.13961777 -1.82593774  0.11594755
## 9      -1.40307466 -5.59182087  0.46400729
## 10      0.03606404  0.51944692 -0.04982259
## 11      0.37415114 -1.07941522  0.11835170
## 12     -0.11799495 -1.11909663  0.18322017
## 13     -0.83491581  1.03469722  0.02495756
## 14      0.02435289  3.20171821 -0.15492789
## 15     -1.82593769  0.96194797  0.14388309
## 16      0.11594755  3.04960613  0.09373458
## 17     -5.59182082 -2.31273350 -0.13319695
## 18      0.46400729  0.32324282 -0.16262237
## 19      0.51944688 -0.07316594 -0.12598450
## 20     -0.04982259  0.11181130 -0.12268061
## 21     -1.07941517  1.43310624 -0.14279903
## 22      0.11835171  0.62044804 -0.03708906
## 23     -1.11909660 -0.37057384  0.05450561
## 24      0.18322018  2.38444671 -0.16952522
## 25      1.03469720 -0.01384650  0.08273621
## 26      0.02495756  2.30508388 -0.03978828
## 27      3.20171811  2.62212981  0.05331079
## 28     -0.15492788 -0.96698456 -0.06544230
## 29      0.96194797 -2.15621081  0.04447130
## 30      0.14388309 -1.58143839  0.02194874
## 31      3.04960605  0.38737823 -0.13961777
## 32      0.09373459 -1.40307472  0.03606404
## 33     -2.31273345  0.37415109 -0.11799495
## 34     -0.13319696 -0.83491587  0.02435289
## 35      0.32324281 -1.82593774  0.11594755
## 36     -0.16262238 -5.59182087  0.46400729
## 37     -0.07316592  0.51944692 -0.04982259
## 38     -0.12598451 -1.07941522  0.11835170
## 39      0.11181134 -1.11909663  0.18322017
## 40     -0.12268062  1.03469722  0.02495756
## 41      1.43310621  3.20171821 -0.15492789
## 42     -0.14279903  0.96194797  0.14388309
## 43      0.62044803  3.04960613  0.09373458
## 44     -0.03708906 -2.31273350 -0.13319695
## 45     -0.37057382  0.32324282 -0.16262237
## 46      0.05450561 -0.07316594 -0.12598450
## 47      2.38444667  0.11181130 -0.12268061
## 48     -0.16952523  1.43310624 -0.14279903
## 49     -0.01384653  0.62044804 -0.03708906
## 50      0.08273621 -0.37057384  0.05450561
## 51      2.30508383  2.38444671 -0.16952522
## 52     -0.03978828 -0.01384650  0.08273621
## 53      2.62212974  2.30508388 -0.03978828
## 54      0.05331080  2.62212981  0.05331079

covariance matrix of random effect coefficients

head(Lambda_new)

## 6 x 54 sparse Matrix of class "dgCMatrix"
##                                                                               
## [1,]  7.822487 -0.48495997  .         .           .         .          . . . .
## [2,] -0.484960  0.05126495  .         .           .         .          . . . .
## [3,]  .         .           7.822487 -0.48495997  .         .          . . . .
## [4,]  .         .          -0.484960  0.05126495  .         .          . . . .
## [5,]  .         .           .         .           7.822487 -0.48495997 . . . .
## [6,]  .         .           .         .          -0.484960  0.05126495 . . . .
##                                                                               
## [1,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [2,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [3,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [4,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [5,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
## [6,] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##                   
## [1,] . . . . . . .
## [2,] . . . . . . .
## [3,] . . . . . . .
## [4,] . . . . . . .
## [5,] . . . . . . .
## [6,] . . . . . . .

Computer predicted values

yhat <- X%*%(bhat)+Z%*%comput_uhat
# comput_uhat= uhat
head(yhat)

## 6 x 1 Matrix of class "dgeMatrix"
##       [,1]
## 1 25.02967
## 2 26.63781
## 3 28.24595
## 4 29.85408
## 5 21.51081
## 6 22.87508

head(fitted(fit.lmer.slope))

##        1        2        3        4        5        6 
## 25.02967 26.63781 28.24595 29.85408 21.51081 22.87508

Using lmm with two grouping factors

library(lme4) 
fit.lmer.slope2 <- lmer(distance ~ age + Sex  +(1+ age | Subject)+ (1 | age), data = Data)

## boundary (singular) fit: see help('isSingular')

summary(fit.lmer.slope2)

## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + Sex + (1 + age | Subject) + (1 | age)
##    Data: Data
## 
## REML criterion at convergence: 435.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0814 -0.4568  0.0155  0.4471  3.8944 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.  Corr 
##  Subject  (Intercept) 7.824e+00 2.797e+00      
##           age         5.127e-02 2.264e-01 -0.77
##  age      (Intercept) 3.470e-09 5.891e-05      
##  Residual             1.716e+00 1.310e+00      
## Number of obs: 108, groups:  Subject, 27; age, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 17.63523    0.88626  19.899
## age          0.66019    0.07125   9.265
## SexFemale   -2.14557    0.75746  -2.833
## 
## Correlation of Fixed Effects:
##           (Intr) age   
## age       -0.838       
## SexFemale -0.348  0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Get z2

58= 54+4 ages

Z2 <- getME(fit.lmer.slope2, "Z")
head(Z2,10)

## 10 x 58 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 58 column names 'M16', 'M16', 'M05' ... ]]

##                                                                                
## 1  . . . . .  . . . . . . . .  . . . . . . . . . . . . . . . 1  8 . . . . . . .
## 2  . . . . .  . . . . . . . .  . . . . . . . . . . . . . . . 1 10 . . . . . . .
## 3  . . . . .  . . . . . . . .  . . . . . . . . . . . . . . . 1 12 . . . . . . .
## 4  . . . . .  . . . . . . . .  . . . . . . . . . . . . . . . 1 14 . . . . . . .
## 5  . . . . 1  8 . . . . . . .  . . . . . . . . . . . . . . . .  . . . . . . . .
## 6  . . . . 1 10 . . . . . . .  . . . . . . . . . . . . . . . .  . . . . . . . .
## 7  . . . . 1 12 . . . . . . .  . . . . . . . . . . . . . . . .  . . . . . . . .
## 8  . . . . 1 14 . . . . . . .  . . . . . . . . . . . . . . . .  . . . . . . . .
## 9  . . . . .  . . . . . . . 1  8 . . . . . . . . . . . . . . .  . . . . . . . .
## 10 . . . . .  . . . . . . . 1 10 . . . . . . . . . . . . . . .  . . . . . . . .
##                                             
## 1  . . . . . . . . . . . . . . . . . 1 . . .
## 2  . . . . . . . . . . . . . . . . . . 1 . .
## 3  . . . . . . . . . . . . . . . . . . . 1 .
## 4  . . . . . . . . . . . . . . . . . . . . 1
## 5  . . . . . . . . . . . . . . . . . 1 . . .
## 6  . . . . . . . . . . . . . . . . . . 1 . .
## 7  . . . . . . . . . . . . . . . . . . . 1 .
## 8  . . . . . . . . . . . . . . . . . . . . 1
## 9  . . . . . . . . . . . . . . . . . 1 . . .
## 10 . . . . . . . . . . . . . . . . . . 1 . .

dim(Z2)

## [1] 108  58

Using lmm with nested or crossed random effects

library(lme4) 
fit.lmer.slope3 <- lmer(distance ~ age + Sex  +(1 | Sex/Subject), data = Data)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

summary(fit.lmer.slope3)

## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + Sex + (1 | Sex/Subject)
##    Data: Data
## 
## REML criterion at convergence: 437.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7489 -0.5503 -0.0252  0.4534  3.6575 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  Subject:Sex (Intercept) 3.2668   1.8074  
##  Sex         (Intercept) 0.9392   0.9691  
##  Residual                2.0495   1.4316  
## Number of obs: 108, groups:  Subject:Sex, 27; Sex, 2
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 17.70671    1.27852   13.85
## age          0.66019    0.06161   10.72
## SexFemale   -2.32102    1.56785   -1.48
## 
## Correlation of Fixed Effects:
##           (Intr) age   
## age       -0.530       
## SexFemale -0.586  0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

Get z3

29= 27+2 ages

Z3 <- getME(fit.lmer.slope3, "Z")
head(Z3,10)

## 10 x 29 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 29 column names 'M16:Male', 'M05:Male', 'M02:Male' ... ]]

##                                                             
## 1  . . . . . . . . . . . . . . 1 . . . . . . . . . . . . 1 .
## 2  . . . . . . . . . . . . . . 1 . . . . . . . . . . . . 1 .
## 3  . . . . . . . . . . . . . . 1 . . . . . . . . . . . . 1 .
## 4  . . . . . . . . . . . . . . 1 . . . . . . . . . . . . 1 .
## 5  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
## 6  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
## 7  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
## 8  . . 1 . . . . . . . . . . . . . . . . . . . . . . . . 1 .
## 9  . . . . . . 1 . . . . . . . . . . . . . . . . . . . . 1 .
## 10 . . . . . . 1 . . . . . . . . . . . . . . . . . . . . 1 .

dim(Z3)

## [1] 108  29

Remark

nested random effects occur when one factor (grouping variable) appears only within a particular level of another factor (grouping variable).
lmer (1|group1/group2)

Crossed random effects are simply: not nested. individual observations are nested separately within the two or more factors.
lmer (1|group1) + (1|group2)

Linear mixed model covariance decomposition with random slopes- lme4

Daniel

2022-06-17

Linear mixed model covariance decomposition with random slopes- lme4

Load data

Using lmm

Get x,y,z

Get fixed effect and random effect coefficients

Get random effect covariance structure (covariance matrix expandation)

Get residual variance and its identity matrix

Get y covariance matrix

Get fixed effect coefficients covariance matrix

Get random effect coefficients covariance matrix (standard deviation)

Computer fixed effect coefficients

Computer covariance of fixed efffect coefficients

Computer random effect coefficients

Computer predicted values

Using lmm with two grouping factors

Get z2

Using lmm with nested or crossed random effects

Get z3