Linear mixed model covariance decomposition (nlme)
Daniel
2022-06-17
Linear mixed model covariance decomposition (nlme)
Load data
data(Orthodont, package="nlme")
Data <- Orthodontlibrary(ggplot2)
ggplot(Data, aes(y = distance, x = age, color = Sex)) +
geom_smooth(method = "lm",se=F) + geom_point() + theme_classic()## `geom_smooth()` using formula 'y ~ x'
library(ggplot2)
ggplot(Data, aes(y = distance, x = age, fill = Subject, color = Sex)) +
geom_smooth(method = "lm",se=F) + geom_point() + theme_classic()## `geom_smooth()` using formula 'y ~ x'
Using lme (nlme)
library(nlme)
fit.lmer2 <- lme(distance ~ age + Sex , random = ~1 | Subject, data = Data)
summary(fit.lmer2)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 447.5125 460.7823 -218.7563
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.807425 1.431592
##
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.706713 0.8339225 80 21.233044 0.0000
## age 0.660185 0.0616059 80 10.716263 0.0000
## SexFemale -2.321023 0.7614168 25 -3.048294 0.0054
## Correlation:
## (Intr) age
## age -0.813
## SexFemale -0.372 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.74889609 -0.55034466 -0.02516628 0.45341781 3.65746539
##
## Number of Observations: 108
## Number of Groups: 27library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmList# lme4
fit.lmer <- lmer(distance ~ age + Sex +(1 | Subject), data = Data)
# compound symmetry- t test
tstat <- fixef(fit.lmer2)/sqrt(diag(vcov(fit.lmer2)))
pval <- 2*pt(-abs(tstat), df = fit.lmer2$fixDF$X)
Tresult <- data.frame(t = tstat, p = round(pval, 4))
print(Tresult)## t p
## (Intercept) 21.233044 0.0000
## age 10.716263 0.0000
## SexFemale -3.048294 0.0054- F test
anova(fit.lmer2 ) ## numDF denDF F-value p-value
## (Intercept) 1 80 4123.156 <.0001
## age 1 80 114.838 <.0001
## Sex 1 25 9.292 0.0054- restricted model test
anova(fit.lmer2, L = c(0, 1, 0)) ## F-test for linear combination(s)
## [1] 1
## numDF denDF F-value p-value
## 1 1 80 114.8383 <.0001- likelihood test
fit.lmer2.2 <- lme(distance ~ Sex , random = ~1 | Subject, data = Data)
anova(fit.lmer2, fit.lmer2.2)## Warning in anova.lme(fit.lmer2, fit.lmer2.2): fitted objects with different
## fixed effects. REML comparisons are not meaningful.## Model df AIC BIC logLik Test L.Ratio p-value
## fit.lmer2 1 5 447.5125 460.7823 -218.7562
## fit.lmer2.2 2 4 513.8718 524.5255 -252.9359 1 vs 2 68.35926 <.0001And now by hand
logLik <- logLik(fit.lmer2)
logLik0 <- logLik(fit.lmer2.2)
LR <- 2 * (logLik - logLik0)
pval <- pchisq(LR, df = 2, lower.tail = FALSE)
LRresult <- data.frame(LR = LR, p = round(pval, 4), row.names = "age")
print(LRresult)## LR p
## age 68.35926 0Model diagnosing
library(nlme)
plot.lme(fit.lmer2)
Get x, y, z matrixs (using lme4)
X=(getME(fit.lmer, "X"))
y=getME(fit.lmer, "y")
Z <- getME(fit.lmer, "Z")
Z2 <- model.matrix(~Subject-1,data=Data)check z matrix
dummyz<- as.matrix (Z)==Z2
table(dummyz)## dummyz
## TRUE
## 2916dim(Z)## [1] 108 27check y and X matrixs
dummyx<- Data[,c(2)]==X[,2]
table(dummyx)## dummyx
## TRUE
## 108Fixed effect coefficient
bhat <- fit.lmer2$coef$fixed
bhat## (Intercept) age SexFemale
## 17.7067130 0.6601852 -2.3210227# fixef()- fixed effect coefficient confidence intervals
intervals(fit.lmer2)## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 16.0471544 17.7067130 19.3662716
## age 0.5375855 0.6601852 0.7827849
## SexFemale -3.8891901 -2.3210227 -0.7528554
##
## Random Effects:
## Level: Subject
## lower est. upper
## sd((Intercept)) 1.31035 1.807425 2.493061
##
## Within-group standard error:
## lower est. upper
## 1.226102 1.431592 1.671522Plot fixed effect coefficients and theri ci
tab <- cbind(Est = intervals(fit.lmer2)[["fixed"]][,2], LL = intervals(fit.lmer2)[["fixed"]][,1], UL=
intervals(fit.lmer2)[["fixed"]][,3])
#Extracting the fixed effect estimates and manually calculating the 95% confidence limits
#qnorm extracts the standard normal critical value for 1-alpha/2= 1-0.05/2=0.975
Odds= (tab)
round(Odds,4)## Est LL UL
## (Intercept) 17.7067 16.0472 19.3663
## age 0.6602 0.5376 0.7828
## SexFemale -2.3210 -3.8892 -0.7529Odds=as.data.frame(Odds)
Odds$Label=rownames(Odds)
ggplot(Odds[-1,],aes(x=Est,y=Label))+geom_point()+
geom_errorbarh(aes(xmax = UL, xmin = LL))+
theme_bw()+geom_vline(xintercept=0,col="darkgray",size=1.2,linetype=2)+
theme(axis.title.y = element_blank(),axis.text = element_text(size=10),
axis.title.x = element_text(size=12))+labs(x="Beta coefficients")
And now by hand
print(tquants <- qt(0.975, df = fit.lmer2$fixDF$X))## (Intercept) age SexFemale
## 1.990063 1.990063 2.059539Low <- fixef(fit.lmer2) - tquants * sqrt(diag(vcov(fit.lmer2)))
Upp <- fixef(fit.lmer2) + tquants * sqrt(diag(vcov(fit.lmer2)))
EstInt <- data.frame(Lower = Low, Estimate = fixef(fit.lmer2), Upper = Upp)
print(EstInt)## Lower Estimate Upper
## (Intercept) 16.0471544 17.7067130 19.3662716
## age 0.5375855 0.6601852 0.7827849
## SexFemale -3.8891901 -2.3210227 -0.7528554Random effect coefficient
fit.lmer2$coef$random## $Subject
## (Intercept)
## M16 -1.70183357
## M05 -1.70183357
## M02 -1.37767479
## M11 -1.16156894
## M07 -1.05351602
## M08 -0.94546309
## M03 -0.62130432
## M12 -0.62130432
## M13 -0.62130432
## M14 -0.08103969
## M09 0.13506616
## M15 0.78338371
## M06 1.21559540
## M04 1.43170125
## M01 2.40417758
## M10 3.91691853
## F10 -3.58539251
## F09 -1.31628108
## F06 -1.31628108
## F01 -1.10017523
## F05 -0.01964599
## F07 0.30451279
## F02 0.30451279
## F08 0.62867156
## F03 0.95283034
## F04 1.92530666
## F11 3.22194176Fixed effect coefficients covariance
fit.lmer2$varFix## (Intercept) age SexFemale
## (Intercept) 0.69542669 -4.174818e-02 -2.361967e-01
## age -0.04174818 3.795289e-03 -2.492877e-17
## SexFemale -0.23619673 -2.492877e-17 5.797556e-01# vcov(fit.lmer2)
sqrt(diag(vcov(fit.lmer2)))## (Intercept) age SexFemale
## 0.83392247 0.06160592 0.76141685fixed effect coefficients correlation
-0.04174818/prod(sqrt( diag(fit.lmer2$varFix)[-3] ))## [1] -0.8126236Random effect covariance and correlation
pdMatrix(fit.lmer2$modelStruct$reStruct[[1]])*fit.lmer2$sigma**2## (Intercept)
## (Intercept) 3.266784getVarCov(fit.lmer2) #using default ## Random effects variance covariance matrix
## (Intercept)
## (Intercept) 3.2668
## Standard Deviations: 1.8074getVarCov(fit.lmer2, individual = "F01", type = "marginal")## Subject F01
## Marginal variance covariance matrix
## 1 2 3 4
## 1 5.3162 3.2668 3.2668 3.2668
## 2 3.2668 5.3162 3.2668 3.2668
## 3 3.2668 3.2668 5.3162 3.2668
## 4 3.2668 3.2668 3.2668 5.3162
## Standard Deviations: 2.3057 2.3057 2.3057 2.3057getVarCov(fit.lmer2, type = "conditional")## Subject M01
## Conditional variance covariance matrix
## 1 2 3 4
## 1 2.0495 0.0000 0.0000 0.0000
## 2 0.0000 2.0495 0.0000 0.0000
## 3 0.0000 0.0000 2.0495 0.0000
## 4 0.0000 0.0000 0.0000 2.0495
## Standard Deviations: 1.4316 1.4316 1.4316 1.4316vc <- as.numeric(as.matrix(VarCorr(fit.lmer2))[,1])
vc/sum(vc)## [1] 0.6144914 0.3855086VarCorr(fit.lmer2)## Subject = pdLogChol(1)
## Variance StdDev
## (Intercept) 3.266784 1.807425
## Residual 2.049456 1.431592Get y covariance directly
= Z REcov Zt + sigma**2
require(mgcv)## Loading required package: mgcv## This is mgcv 1.8-40. For overview type 'help("mgcv-package")'.library(nlme)
cov <- extract.lme.cov(fit.lmer2,Data)
head(cov)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.316240 3.266784 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [2,] 3.266784 5.316240 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [3,] 3.266784 3.266784 5.316240 3.266784 0.000000 0.000000 0.000000 0.000000
## [4,] 3.266784 3.266784 3.266784 5.316240 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.316240 3.266784 3.266784 3.266784
## [6,] 0.000000 0.000000 0.000000 0.000000 3.266784 5.316240 3.266784 3.266784
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov)## [1] 108 108Residual variance
fit.lmer2$sigma**2## [1] 2.049456Compute fixed effect coefficients
library(matlib)
inv(t(as.matrix(X))%*%inv(as.matrix(cov))%*%as.matrix(X))%*%t(as.matrix(X))%*%inv(as.matrix(cov))%*%y## [,1]
## [1,] 17.7067076
## [2,] 0.6601871
## [3,] -2.3210230Compute covariance of fixed efffect coefficients
inv(t(as.matrix(X))%*%inv(as.matrix(cov))%*%as.matrix(X))##
## [1,] 0.69542669 -0.04174818 -0.2361967
## [2,] -0.04174818 0.00379529 0.0000000
## [3,] -0.23619674 0.00000000 0.5797556Compute random effect coefficients
Lambda_new <-as.numeric(VarCorr(fit.lmer2)[1])*diag(length(levels(Data$Subject)))
# head(Lambda_new)uhat <- fit.lmer2$coef$random
comput_uhat <- (as.matrix(Lambda_new))%*%t(Z)%*%inv(as.matrix(cov))%*%(y-as.matrix(X)%*%(bhat))
cbind((comput_uhat@x),(uhat[["Subject"]]))## (Intercept)
## M16 -1.70183367 -1.70183357
## M05 -1.70183367 -1.70183357
## M02 -1.37767488 -1.37767479
## M11 -1.16156901 -1.16156894
## M07 -1.05351608 -1.05351602
## M08 -0.94546315 -0.94546309
## M03 -0.62130436 -0.62130432
## M12 -0.62130436 -0.62130432
## M13 -0.62130436 -0.62130432
## M14 -0.08103970 -0.08103969
## M09 0.13506616 0.13506616
## M15 0.78338375 0.78338371
## M06 1.21559548 1.21559540
## M04 1.43170134 1.43170125
## M01 2.40417773 2.40417758
## M10 3.91691877 3.91691853
## F10 -3.58539273 -3.58539251
## F09 -1.31628117 -1.31628108
## F06 -1.31628117 -1.31628108
## F01 -1.10017530 -1.10017523
## F05 -0.01964599 -0.01964599
## F07 0.30451281 0.30451279
## F02 0.30451281 0.30451279
## F08 0.62867160 0.62867156
## F03 0.95283040 0.95283034
## F04 1.92530678 1.92530666
## F11 3.22194196 3.22194176Compute predicted values
usually computing mean and prediction ci by using bootstrapping and simulated methods
approximate at below:
This webpage calculates the sd like above, but I think the following is also correct.
the conditional distribution
\[ \mathbf{y}_{} \mid \mathbf{b}_{} \sim \mathbf{N}\left(\mathbf{X}_{} \boldsymbol{\beta}+\mathbf{Z}_{} \mathbf{b}_{}, \boldsymbol{\Sigma}_{}\right) \] therefore
conditional mean prediction\[ \begin{align} E(\hat{Y}_0)&=E\mathbf{(X_0\hat{\beta} + Z_0\hat{u})}=\mathbf{X_0\beta + Z_0{u}}=E\mathbf{(Y_0)}\\ var(\hat{Y}_0)&=E\mathbf{((X_0\hat{\beta}-X_0\beta)+ (Z_0\hat{u} - Z_0{u}))}^2\\&=E\mathbf{((X_0\hat{\beta}-X_0\beta)^2+ (Z_0\hat{u} - Z_0{u})^2+ 2 ((X_0\hat{\beta}-X_0\beta)* (Z_0\hat{u} - Z_0{u})) )} \\ &=E\mathbf{\left( X_0(\hat{\beta}-\beta)(\hat{\beta}-\beta)'X_0' \right)}+ E\mathbf{\left( Z_0(\hat{u}-u)(\hat{u}-u)'Z_0' \right)}+2E ((X_0(\hat{\beta}-\beta))* (Z_0(\hat{u} - {u}))) ? \\ &=\sigma^2\mathbf{X_0\left( X'X \right)^{-1}X_0'}+ \sigma^2\mathbf{Z_0\left( Z'Z \right)^{-1}Z_0'} + 2 *0*0 \quad \text {if independent}\\ &= \mathbf{X_0 cov (X) X_0'}+ \mathbf{Z_0 cov (Z) Z_0'} \\ \end{align} \]generally\[ \begin{align} Var(X+Y) &= Cov(X+Y,X+Y) \\ &= E((X+Y)^2)-E(X+Y)E(X+Y) \\ &\text{by expanding,} \\ &= E(X^2) - (E(X))^2 + E(Y^2) - (E(Y))^2 + 2(E(XY) - E(X)E(Y)) \\ &= Var(X) + Var(Y) + 2(E(XY)) - E(X)E(Y)) \\ \end{align} \]conditional individual prediction\[ \begin{align} E(\hat{Y}_0)&= E\mathbf{(Y_0)}\\ var(\hat{Y}_0) &= \mathbf{X_0 cov (X) X_0'}+ \mathbf{Z_0 cov (Z) Z_0'} + \sigma^2 \\ \end{align} \]
\[ \begin{align} \hat{Y}_0& \sim N(\mu_{\hat{Y}_0},\sigma^2_{\hat{Y}_0})\\ \end{align} \]
lower level (conditional prediction) in R
# how to calculate predicted values
yhat <- X%*%(fit.lmer2$coef$fixed)+Z%*% as.numeric ( uhat[["Subject"]])
head(cbind (yhat,predict(fit.lmer2),y)) #create individual trajectory curve ## 6 x 3 Matrix of class "dgeMatrix"
## y
## 1 25.39237 25.39237 26.0
## 2 26.71274 26.71274 25.0
## 3 28.03311 28.03311 29.0
## 4 29.35348 29.35348 31.0
## 5 21.61052 21.61052 21.5
## 6 22.93089 22.93089 22.5#compute standard error for marginal predictions
predvar_rand <- diag(X %*% fit.lmer2$varFix %*% t(X)) + diag(Z %*% diag(getVarCov(fit.lmer2)[1] ,27) %*% t(Z))
SE_rand <- sqrt (predvar_rand) #mean prediction
SE_rand2 <- sqrt(predvar_rand+fit.lmer2$sigma^2) #individual prediction
head(SE_rand,20)## 1 2 3 4 5 6 7 8
## 1.880728 1.872639 1.872639 1.880728 1.880728 1.872639 1.872639 1.880728
## 9 10 11 12 13 14 15 16
## 1.880728 1.872639 1.872639 1.880728 1.880728 1.872639 1.872639 1.880728
## 17 18 19 20
## 1.880728 1.872639 1.872639 1.880728head(SE_rand2,20)## 1 2 3 4 5 6 7 8
## 2.363598 2.357166 2.357166 2.363598 2.363598 2.357166 2.357166 2.363598
## 9 10 11 12 13 14 15 16
## 2.363598 2.357166 2.357166 2.363598 2.363598 2.357166 2.357166 2.363598
## 17 18 19 20
## 2.363598 2.357166 2.357166 2.363598see below
The marginal distribution
\[ \mathbf{y}_{} \sim \mathrm{N}\left(\mathbf{X}_{} \boldsymbol{\beta}, \mathbf{V}_{}\right) \]
\[ \operatorname{Cov}[\mathbf{y}]=\mathbf{V}=\mathbf{Z D Z}^{\prime}+\Sigma \]
\(\mathbf{Z}_{} \mathbf{D Z}_{}^{\prime}\) represents the random-effects structure
therefore
marginal mean prediction
\[ \begin{align} E(\hat{Y}_0)&=E\mathbf{(X_0\hat{\beta})}=\mathbf{X_0\beta}=E\mathbf{(Y_0)}\\ var(\hat{Y}_0)&=E\mathbf{(X_0\hat{\beta}-X_0\beta)}^2\\ &=E\mathbf{\left( X_0(\hat{\beta}-\beta)(\hat{\beta}-\beta)'X_0' \right)}\\ &=E\mathbf{X_0\left( (\hat{\beta}-\beta)(\hat{\beta}-\beta)' \right)X_0'}\\ &=\mathbf{X_0Cov(X)X_0'}\\ \end{align} \]
marginal individual prediction\[ \begin{align} E(\hat{Y}_0)&= E\mathbf{(Y_0)}\\ var(\hat{Y}_0) &= \mathbf{X_0 cov (X) X_0'}+ \sigma^2 \\ \end{align} \]higher level (marginal)in R
#compute standard error for marginal predictions
predvar <- diag(X %*% fit.lmer2$varFix %*% t(X))
SE <- sqrt (predvar) #mean prediction
SE2 <- sqrt(predvar+fit.lmer2$sigma^2) #individual prediction
head(SE,10)## 1 2 3 4 5 6 7 8
## 0.5199561 0.4898898 0.4898898 0.5199561 0.5199561 0.4898898 0.4898898 0.5199561
## 9 10
## 0.5199561 0.4898898head(SE2,10)## 1 2 3 4 5 6 7 8
## 1.523092 1.513092 1.513092 1.523092 1.523092 1.513092 1.513092 1.523092
## 9 10
## 1.523092 1.513092up=predict(fit.lmer2, newdata=Data, level=0) +1.96 *SE #mean prediction
up2=predict(fit.lmer2, newdata=Data, level=0) +1.96 *SE2 #individual prediction
head(up)## 1 2 3 4 5 6
## 24.00731 25.26875 26.58912 27.96842 24.00731 25.26875head(up2)## 1 2 3 4 5 6
## 25.97346 27.27423 28.59460 29.93457 25.97346 27.27423library(tidyverse)
library(ggeffects)
ggpredict(fit.lmer2,terms=c("age" ))## # Predicted values of distance
##
## age | Predicted | 95% CI
## --------------------------------
## 8 | 22.99 | [21.97, 24.01]
## 10 | 24.31 | [23.35, 25.27]
## 12 | 25.63 | [24.67, 26.59]
## 14 | 26.95 | [25.93, 27.97]
##
## Adjusted for:
## * Sex = Male
## * Subject = 0 (population-level)ggpredict(fit.lmer2,terms=c("age" )) %>% plot(rawdata = T, dot.alpha = 0.2)
# ggpredict(fit.lmer2,"age", type = "re" ) %>% plot(rawdata = T, dot.alpha = 0.2)Compute SE by hand and compare with above SE:
0.520408163 0.489795918 0.489795918 0.520408163
head(SE,10)## 1 2 3 4 5 6 7 8
## 0.5199561 0.4898898 0.4898898 0.5199561 0.5199561 0.4898898 0.4898898 0.5199561
## 9 10
## 0.5199561 0.4898898# head(predict(fit.lmer2))
# head (predict(fit.lmer2, newdata=Data, level=1)) #predicted value
# head (predict(fit.lmer2, newdata=Data, level=0)) #mean and ci Using lme with Gaussian
it means the correlation in a group but not correlation in variables (different covariances)
library(nlme)
fit.lmer.gaus <- lme(distance ~ age + Sex , random = ~1 | Subject, correlation=corGaus(form= ~ age|Subject, nugget=TRUE), data = Data)
summary(fit.lmer.gaus)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 451.4311 470.0089 -218.7156
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.795504 1.446192
##
## Correlation Structure: Gaussian spatial correlation
## Formula: ~age | Subject
## Parameter estimate(s):
## range nugget
## 1.1324270 0.1039537
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.717033 0.8448813 80 20.969848 0.0000
## age 0.659620 0.0628344 80 10.497756 0.0000
## SexFemale -2.325355 0.7612563 25 -3.054629 0.0053
## Correlation:
## (Intr) age
## age -0.818
## SexFemale -0.367 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.70618679 -0.54458678 -0.01435195 0.45606804 3.62915585
##
## Number of Observations: 108
## Number of Groups: 27- Plotting of the empirical covariation (of residuals) versus “age”
when describing the spatial pattern of a measured variable, a common way of visualizing the spatial autocorrelation of a variable is a variogram plot.
checking the semi-variance change by the distance of spatial autocorrelation.
print(plot(Variogram(fit.lmer.gaus, form =~as.numeric(age)|Subject,
data = Data)))
require(mgcv)
library(nlme)
cov.gaus <- extract.lme.cov(fit.lmer.gaus,Data)
head(cov.gaus)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.315307 3.306657 3.223843 3.223836 0.000000 0.000000 0.000000 0.000000
## [2,] 3.306657 5.315307 3.306657 3.223843 0.000000 0.000000 0.000000 0.000000
## [3,] 3.223843 3.306657 5.315307 3.306657 0.000000 0.000000 0.000000 0.000000
## [4,] 3.223836 3.223843 3.306657 5.315307 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.315307 3.306657 3.223843 3.223836
## [6,] 0.000000 0.000000 0.000000 0.000000 3.306657 5.315307 3.306657 3.223843
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov.gaus)## [1] 108 108- Compute fixed effect coefficients using Gaussian
library(matlib)
inv(t(as.matrix(X))%*%inv(as.matrix(cov.gaus))%*%as.matrix(X))%*%t(as.matrix(X))%*%inv(as.matrix(cov.gaus))%*%y## [,1]
## [1,] 17.7170272
## [2,] 0.6596201
## [3,] -2.3253552fit.lmer.gaus$coef$fixed## (Intercept) age SexFemale
## 17.717033 0.659620 -2.325355the same below
Using lme with autoregressive
fit.lmer.AR1 <- lme(distance ~ age + Sex , random = ~1 | Subject, correlation=corAR1(form= ~ as.numeric(age)|Subject ), data = Data)
summary(fit.lmer.AR1)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 449.5125 465.4363 -218.7563
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.807425 1.431592
##
## Correlation Structure: ARMA(1,0)
## Formula: ~as.numeric(age) | Subject
## Parameter estimate(s):
## Phi1
## 0
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.706713 0.8339225 80 21.233044 0.0000
## age 0.660185 0.0616059 80 10.716263 0.0000
## SexFemale -2.321023 0.7614168 25 -3.048294 0.0054
## Correlation:
## (Intr) age
## age -0.813
## SexFemale -0.372 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.74889609 -0.55034466 -0.02516628 0.45341781 3.65746539
##
## Number of Observations: 108
## Number of Groups: 27- Plotting of the empirical covariation (of residuals) versus “age”
print(plot(Variogram(fit.lmer.AR1, form =~as.numeric(age)|Subject,
data = Data)))
require(mgcv)
library(nlme)
cov.AR1 <- extract.lme.cov(fit.lmer.AR1,Data)
head(cov.AR1)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.316240 3.266784 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [2,] 3.266784 5.316240 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [3,] 3.266784 3.266784 5.316240 3.266784 0.000000 0.000000 0.000000 0.000000
## [4,] 3.266784 3.266784 3.266784 5.316240 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.316240 3.266784 3.266784 3.266784
## [6,] 0.000000 0.000000 0.000000 0.000000 3.266784 5.316240 3.266784 3.266784
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov.AR1)## [1] 108 108Using lme with exponential
fit.lmer.Exp <- lme(distance ~ age + Sex , random = ~1 | Subject, correlation=corExp(form= ~ as.numeric(age)|Subject ), data = Data)
summary(fit.lmer.Exp)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 449.3968 465.3206 -218.6984
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.788899 1.454494
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~as.numeric(age) | Subject
## Parameter estimate(s):
## range
## 0.7045117
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.721416 0.8500193 80 20.848250 0.0000
## age 0.659405 0.0634074 80 10.399499 0.0000
## SexFemale -2.327485 0.7611852 25 -3.057711 0.0053
## Correlation:
## (Intr) age
## age -0.821
## SexFemale -0.365 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.683026667 -0.540915318 -0.008097445 0.461167542 3.612579065
##
## Number of Observations: 108
## Number of Groups: 27- Plotting of the empirical covariation (of residuals) versus “age”
print(plot(Variogram(fit.lmer.Exp, form =~as.numeric(age)|Subject,
data = Data)))
require(mgcv)
library(nlme)
cov.Exp <- extract.lme.cov(fit.lmer.Exp,Data)
head(cov.Exp)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.315710 3.323904 3.207397 3.200582 0.000000 0.000000 0.000000 0.000000
## [2,] 3.323904 5.315710 3.323904 3.207397 0.000000 0.000000 0.000000 0.000000
## [3,] 3.207397 3.323904 5.315710 3.323904 0.000000 0.000000 0.000000 0.000000
## [4,] 3.200582 3.207397 3.323904 5.315710 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.315710 3.323904 3.207397 3.200582
## [6,] 0.000000 0.000000 0.000000 0.000000 3.323904 5.315710 3.323904 3.207397
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov.Exp)## [1] 108 108Using lme with unstructured
fit.lmer.corSymm <- lme(distance ~ age + Sex , random = ~1 | (Subject), correlation=corSymm(form= ~ 1| (Subject) ),
weights = varIdent(form=~1|age),
data = Data)
summary(fit.lmer.corSymm)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 456.6945 493.85 -214.3473
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.783141 1.481553
##
## Correlation Structure: General
## Formula: ~1 | (Subject)
## Parameter estimate(s):
## Correlation:
## 1 2 3
## 2 -0.260
## 3 0.238 -0.149
## 4 -0.251 -0.007 0.426
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | age
## Parameter estimates:
## 8 10 12 14
## 1.0000000 0.6868412 1.1990772 1.0004135
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.417592 0.8657149 80 20.119316 0.0000
## age 0.674651 0.0702289 80 9.606462 0.0000
## SexFemale -2.045167 0.7361374 25 -2.778241 0.0102
## Correlation:
## (Intr) age
## age -0.840
## SexFemale -0.346 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.33886286 -0.56222407 0.01374192 0.52235930 3.99060002
##
## Number of Observations: 108
## Number of Groups: 27- Plotting of the empirical covariation (of residuals) versus “age”
print(plot(Variogram(fit.lmer.corSymm, form =~as.numeric(age)|Subject,
data = Data)))
require(mgcv)
library(nlme)
cov.corSymm <- extract.lme.cov(fit.lmer.corSymm,Data)
head(cov.corSymm)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.374591 2.786933 3.807033 2.628359 0.000000 0.000000 0.000000 0.000000
## [2,] 2.786933 4.215084 2.909655 3.168383 0.000000 0.000000 0.000000 0.000000
## [3,] 3.807033 2.909655 6.335531 4.301450 0.000000 0.000000 0.000000 0.000000
## [4,] 2.628359 3.168383 4.301450 5.376406 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.374591 2.786933 3.807033 2.628359
## [6,] 0.000000 0.000000 0.000000 0.000000 2.786933 4.215084 2.909655 3.168383
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov.corSymm)## [1] 108 108Using lme with compound symm
fit.lmer.Symm <- lme(distance ~ age + Sex , random = ~1 | (Subject), correlation=corCompSymm(form= ~ as.numeric(age)| (Subject) ), data = Data)
summary(fit.lmer.Symm)## Linear mixed-effects model fit by REML
## Data: Data
## AIC BIC logLik
## 449.5125 465.4363 -218.7563
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 1.807425 1.431592
##
## Correlation Structure: Compound symmetry
## Formula: ~as.numeric(age) | (Subject)
## Parameter estimate(s):
## Rho
## 0
## Fixed effects: distance ~ age + Sex
## Value Std.Error DF t-value p-value
## (Intercept) 17.706713 0.8339225 80 21.233044 0.0000
## age 0.660185 0.0616059 80 10.716263 0.0000
## SexFemale -2.321023 0.7614168 25 -3.048294 0.0054
## Correlation:
## (Intr) age
## age -0.813
## SexFemale -0.372 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.74889609 -0.55034466 -0.02516628 0.45341781 3.65746539
##
## Number of Observations: 108
## Number of Groups: 27- Plotting of the empirical covariation (of residuals) versus “age”
print(plot(Variogram(fit.lmer.Symm, form =~as.numeric(age)|Subject,
data = Data)))
require(mgcv)
library(nlme)
cov.Symm <- extract.lme.cov(fit.lmer.Symm,Data)
head(cov.Symm)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 5.316240 3.266784 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [2,] 3.266784 5.316240 3.266784 3.266784 0.000000 0.000000 0.000000 0.000000
## [3,] 3.266784 3.266784 5.316240 3.266784 0.000000 0.000000 0.000000 0.000000
## [4,] 3.266784 3.266784 3.266784 5.316240 0.000000 0.000000 0.000000 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 5.316240 3.266784 3.266784 3.266784
## [6,] 0.000000 0.000000 0.000000 0.000000 3.266784 5.316240 3.266784 3.266784
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [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
## [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [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
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [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
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [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
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## [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
## [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [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
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92]
## [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
## [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101] [,102] [,103]
## [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
## [,104] [,105] [,106] [,107] [,108]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 0 0 0 0dim(cov.Symm)## [1] 108 108Using gls with unstructured, corSymm
gls.corsymm <- gls(distance ~ Sex +age, Orthodont,
correlation = corSymm(form = ~ 1 | Subject),
weights = varIdent(form=~1|age))
summary(gls.corsymm)## Generalized least squares fit by REML
## Model: distance ~ Sex + age
## Data: Orthodont
## AIC BIC logLik
## 454.6945 489.196 -214.3473
##
## Correlation Structure: General
## Formula: ~1 | Subject
## Parameter estimate(s):
## Correlation:
## 1 2 3
## 2 0.586
## 3 0.652 0.563
## 4 0.489 0.666 0.737
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | age
## Parameter estimates:
## 8 10 12 14
## 1.0000000 0.8855835 1.0857261 1.0001634
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 17.417594 0.8657135 20.119350 0.0000
## SexFemale -2.045174 0.7361405 -2.778238 0.0065
## age 0.674651 0.0702285 9.606514 0.0000
##
## Correlation:
## (Intr) SexFml
## SexFemale -0.346
## age -0.840 0.000
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.508189030 -0.584409277 0.003646472 0.531072913 2.179754376
##
## Residual standard error: 2.318327
## Degrees of freedom: 108 total; 105 residualUsing gls with compound symm, which = lmer model default = lme default = lme compound symm =lme auto reg in this example
gls.comsys <- gls(distance ~ Sex +age, Orthodont,
correlation = corCompSymm(form = ~ 1 | Subject))
summary(gls.comsys)## Generalized least squares fit by REML
## Model: distance ~ Sex + age
## Data: Orthodont
## AIC BIC logLik
## 447.5125 460.7823 -218.7563
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | Subject
## Parameter estimate(s):
## Rho
## 0.6144914
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 17.706713 0.8339225 21.233044 0.0000
## SexFemale -2.321023 0.7614169 -3.048294 0.0029
## age 0.660185 0.0616059 10.716263 0.0000
##
## Correlation:
## (Intr) SexFml
## SexFemale -0.372
## age -0.813 0.000
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.59712955 -0.64544226 -0.02540005 0.51680604 2.32947531
##
## Residual standard error: 2.305697
## Degrees of freedom: 108 total; 105 residuallme4 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.