Introduction to Linear Mixed-Effects Models: nlme Vs lme4 with Rstudio
Introduction to Linear Mixed-Effects Models: nlme Vs lme4 with Rstudio
1 Linear Mixed-Effect Model: Package nlme
2 Restricted Maximum Likelihood (REML)
The way how to find \(\boldsymbol{T}\) and how to estimate \(\boldsymbol{\beta,\Sigma_{U}, \Sigma_{R} }\) can refer to the book of Shayle R. Searle (1992) and the paper of Harville (1977), Laird and Ware (1982), Pinero and Bates (2000), JJ (2006), Crawley (2013) and Satoh (2018)
The advantage of REML approach over ANOVA
- can produce unbiased estimates of variance and covariance parameters;
- can analyze unbalanced designs; and
- has a powerful prediction algorithm that extends the ideas in regression prediction algorithm to cover random as well as fixed effects.
3 Descriptive statistics: ratpupcsv
'data.frame': 322 obs. of 6 variables:
$ pup.id : int 1 2 3 4 5 6 7 8 9 10 ...
$ weight : num 6.6 7.4 7.15 7.24 7.1 6.04 6.98 7.05 6.95 6.29 ...
$ sex : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 1 1 ...
$ litter : int 1 1 1 1 1 1 1 1 1 1 ...
$ litsize : int 12 12 12 12 12 12 12 12 12 12 ...
$ treatment: Factor w/ 3 levels "Control","High",..: 1 1 1 1 1 1 1 1 1 1 ... pup.id weight sex litter litsize
Min. : 1.00 Min. :3.680 Female:151 Min. : 1.00 Min. : 2.00
1st Qu.: 81.25 1st Qu.:5.650 Male :171 1st Qu.: 7.00 1st Qu.:12.00
Median :161.50 Median :6.055 Median :13.50 Median :14.00
Mean :161.50 Mean :6.081 Mean :13.38 Mean :13.33
3rd Qu.:241.75 3rd Qu.:6.397 3rd Qu.:19.75 3rd Qu.:16.00
Max. :322.00 Max. :8.330 Max. :27.00 Max. :18.00
treatment
Control:131
High : 65
Low :126
mysummary <- function(x) c(N = length(x), MEAN = mean(x, na.rm = TRUE), SD = sd(x,
na.rm = TRUE), MIN = min(x, na.rm = TRUE), MAX = max(x, na.rm = TRUE))
pivot_table_01 = summarize(ratpupcsv$weight, by = llist(ratpupcsv$treatment, ratpupcsv$sex),
mysummary)
colnames(pivot_table_01)[1:3] = c("treatment", "sex", "weight")bwplot(weight ~ sex | treatment, data = ratpupcsv, aspect = 2, ylab = "Birth Weights",
xlab = "SEX", main = "Boxplots of Birth Weights by Treatment Levels by Sex")
ggplot(ratpupcsv) + aes(x = sex, y = weight, colour = treatment) + geom_boxplot(fill = "#bdbdbd") +
scale_color_brewer(palette = "Set1") + theme_classic() + facet_wrap(vars(treatment)) +
ggtitle("Boxplots of Birth Weights by Treatment Levels by Sex")
[1] 2.12 2.12 3.23 3.23 3.23 4.03 4.03 4.03 4.03 8.25 8.25 8.25
[13] 8.25 8.25 8.25 8.25 8.25 9.06 9.06 9.06 9.06 9.06 9.06 9.06
[25] 9.06 9.06 9.26 9.26 9.26 9.26 9.26 9.26 9.26 9.26 9.26 9.27
[37] 9.27 9.27 9.27 9.27 9.27 9.27 9.27 9.27 10.19 10.19 10.19 10.19
[49] 10.19 10.19 10.19 10.19 10.19 10.19 10.22 10.22 10.22 10.22 10.22 10.22
[61] 10.22 10.22 10.22 10.22 12.01 12.01 12.01 12.01 12.01 12.01 12.01 12.01
[73] 12.01 12.01 12.01 12.01 12.13 12.13 12.13 12.13 12.13 12.13 12.13 12.13
[85] 12.13 12.13 12.13 12.13 12.24 12.24 12.24 12.24 12.24 12.24 12.24 12.24
[97] 12.24 12.24 12.24 12.24 13.05 13.05 13.05 13.05 13.05 13.05 13.05 13.05
[109] 13.05 13.05 13.05 13.05 13.05 13.10 13.10 13.10 13.10 13.10 13.10 13.10
[ reached getOption("max.print") -- omitted 202 entries ]ranklit <- factor(ranklit)
levels(ranklit) <- c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
"13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25",
"26", "27")bwplot(weight ~ ranklit | treatment * sex, data = ratpupcsv, aspect = 2, ylab = "Birth Weights",
xlab = "Litter ", scales = list(cex = 0.4), main = "Boxplots of Birth Weights by Treatment levels by Sex by Litter")
4 Modeling Specification Strategies:
4.1 Hypthesis Test 01 :
to test the Random Effects (Litter)
Approach: LR Test for Nested linear models
\(H_{0}: \mu_{Litter} = 0\)
\(H_{A}: \mu_{Litter} \neq 0\)
4.1.1 lme: Model.01 with random effects(Litter)
attach(ratpupcsv)
Model.01 <- lme(weight ~ treatment + sex + litsize + treatment * sex, random = ~1 |
litter, method = "REML")
summary(Model.01)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
419.1043 452.8775 -200.5522
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3106722 0.404337
Fixed effects: weight ~ treatment + sex + litsize + treatment * sex
Value Std.Error DF t-value p-value
(Intercept) 7.911652 0.27496390 292 28.773422 0.0000
treatmentHigh -0.799034 0.19429415 23 -4.112495 0.0004
treatmentLow -0.383174 0.15967364 23 -2.399731 0.0249
sexMale 0.411688 0.07315410 292 5.627679 0.0000
litsize -0.128382 0.01875336 23 -6.845819 0.0000
treatmentHigh:sexMale -0.107023 0.13176318 292 -0.812239 0.4173
treatmentLow:sexMale -0.083866 0.10568189 292 -0.793568 0.4281
Correlation:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmentHigh -0.581
treatmentLow -0.336 0.434
sexMale -0.155 0.221 0.269
litsize -0.910 0.372 0.045 -0.001
treatmentHigh:sexMale 0.105 -0.360 -0.150 -0.555 -0.020
treatmentLow:sexMale 0.141 -0.166 -0.345 -0.692 -0.036 0.385
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-7.47250744 -0.50014749 0.02911668 0.57348178 3.00962055
Number of Observations: 322
Number of Groups: 27 numDF denDF F-value p-value
(Intercept) 1 292 9093.772 <.0001
treatment 2 23 5.082 0.0149
sex 1 292 52.602 <.0001
litsize 1 23 47.374 <.0001
treatment:sex 2 292 0.466 0.6282 numDF denDF F-value p-value
(Intercept) 1 292 827.9098 <.0001
treatment 2 23 8.6905 0.0015
sex 1 292 31.6708 <.0001
litsize 1 23 46.8652 <.0001
treatment:sex 2 292 0.4656 0.62824.1.2 gls: Model.01A without random effects(Litter)
library(Hmisc)
attach(ratpupcsv)
Model.01A <- gls(weight ~ treatment + sex + litsize + treatment * sex, data = ratpupcsv)
summary(Model.01A)Generalized least squares fit by REML
Model: weight ~ treatment + sex + litsize + treatment * sex
Data: ratpupcsv
AIC BIC logLik
506.5099 536.5305 -245.255
Coefficients:
Value Std.Error t-value p-value
(Intercept) 7.900732 0.16229542 48.68118 0.0000
treatmentHigh -0.807294 0.12051060 -6.69895 0.0000
treatmentLow -0.381885 0.09273467 -4.11804 0.0000
sexMale 0.339576 0.08902602 3.81435 0.0002
litsize -0.124188 0.01024669 -12.11978 0.0000
treatmentHigh:sexMale -0.178572 0.15325437 -1.16520 0.2448
treatmentLow:sexMale -0.074357 0.12639429 -0.58829 0.5568
Correlation:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmentHigh -0.575
treatmentLow -0.393 0.451
sexMale -0.335 0.439 0.565
litsize -0.907 0.372 0.092 0.014
treatmentHigh:sexMale 0.241 -0.700 -0.333 -0.582 -0.059
treatmentLow:sexMale 0.282 -0.328 -0.733 -0.705 -0.061 0.413
Standardized residuals:
Min Q1 Med Q3 Max
-6.18740247 -0.50909449 -0.04113244 0.61273725 2.93938254
Residual standard error: 0.50151
Degrees of freedom: 322 total; 315 residual4.1.3 Hypthesis Test 01 : results
- results indicated that reject the \(H_{0}: \mu_{Litter} = 0\)
- therefore, \(\mu_{Litter}\) is a significant random factor
Model df AIC BIC logLik Test L.Ratio p-value
Model.01A 1 8 506.5099 536.5305 -245.2550
Model.01 2 9 419.1043 452.8775 -200.5522 1 vs 2 89.40562 <.00014.2 Hypthesis Test 02 :
- to test whether the residual variance differs between treatment groups
4.2.1 lme: Model.02
- with random effects(Litter)
- with a heterogeneous residual variance structure within treatment groups
attach(ratpupcsv)
Model.02 <- lme(weight ~ treatment + sex + litsize + treatment * sex, random = ~1 |
litter, method = "REML", weights = varIdent(form = ~1 | treatment))
summary(Model.02)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
381.8847 423.163 -179.9423
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3134846 0.5147948
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | treatment
Parameter estimates:
Control Low High
1.0000000 0.5649830 0.6394383
Fixed effects: weight ~ treatment + sex + litsize + treatment * sex
Value Std.Error DF t-value p-value
(Intercept) 7.937162 0.27736252 292 28.616565 0.0000
treatmentHigh -0.808610 0.19628137 23 -4.119649 0.0004
treatmentLow -0.390279 0.16301558 23 -2.394121 0.0252
sexMale 0.408131 0.09303486 292 4.386865 0.0000
litsize -0.130007 0.01848708 23 -7.032332 0.0000
treatmentHigh:sexMale -0.094666 0.12919527 292 -0.732737 0.4643
treatmentLow:sexMale -0.076013 0.10811858 292 -0.703053 0.4826
Correlation:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmentHigh -0.616
treatmentLow -0.396 0.498
sexMale -0.197 0.278 0.335
litsize -0.896 0.378 0.070 0.000
treatmentHigh:sexMale 0.155 -0.361 -0.243 -0.720 -0.015
treatmentLow:sexMale 0.189 -0.248 -0.367 -0.860 -0.022 0.620
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.88670114 -0.52493419 0.02123518 0.57307286 2.56409983
Number of Observations: 322
Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-value
Model.01 1 9 419.1043 452.8775 -200.5522
Model.02 2 11 381.8847 423.1630 -179.9423 1 vs 2 41.21964 <.00014.2.2 Hypthesis Test 02 : results
- result indicated that the likelihood ratio testis significant (p<0.0001), so Model.02 with heterogeneous residual variances in the treatment groups has an optimal model fit, implying that the effects of heterogeneous residual variance structure exited within treatment groups.
4.3 Hypthesis Test 03 :
- to test whether the residual variances for the high and low treatment groups are equal
attach(ratpupcsv)
ratpupcsv$trtgrp[treatment == "Control"] <- 1
ratpupcsv$trtgrp[treatment == "Low" | treatment == "High"] <- 24.3.1 lme: Model.03
attach(ratpupcsv)
Model.03 <- lme(weight ~ treatment + sex + litsize + treatment * sex, random = ~1 |
litter, method = "REML", weights = varIdent(form = ~1 | trtgrp))
summary(Model.03)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
381.0807 418.6065 -180.5404
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3145679 0.5147878
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5905488
Fixed effects: weight ~ treatment + sex + litsize + treatment * sex
Value Std.Error DF t-value p-value
(Intercept) 7.942155 0.27838196 292 28.529705 0.0000
treatmentHigh -0.809818 0.19550851 23 -4.142109 0.0004
treatmentLow -0.390200 0.16383300 23 -2.381691 0.0259
sexMale 0.408195 0.09303539 292 4.387529 0.0000
litsize -0.130383 0.01856367 23 -7.023574 0.0000
treatmentHigh:sexMale -0.092026 0.12461723 292 -0.738473 0.4608
treatmentLow:sexMale -0.076397 0.10939797 292 -0.698337 0.4855
Correlation:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmentHigh -0.622
treatmentLow -0.393 0.500
sexMale -0.196 0.279 0.334
litsize -0.897 0.382 0.068 0.000
treatmentHigh:sexMale 0.159 -0.351 -0.250 -0.747 -0.013
treatmentLow:sexMale 0.188 -0.247 -0.369 -0.850 -0.023 0.635
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.88725915 -0.52398492 0.01731335 0.56412339 2.55014733
Number of Observations: 322
Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-value
Model.03 1 10 381.0807 418.6065 -180.5404
Model.02 2 11 381.8847 423.1630 -179.9423 1 vs 2 1.196053 0.27414.3.2 Hypthesis Test 03 : results
- the non-significant result (p=0.2741) indicated that null hypothesis (the residual variances for the high and low dose treatments are equal) is statistically accepted, thereby accepting Model.03 better than Model.02
4.4 Hypthesis Test 04 :
Is the residual variance for the combined high/low treatment group equal to the residual variance for the control group?
with LR Test
Model df AIC BIC logLik Test L.Ratio p-value
Model.01 1 9 419.1043 452.8775 -200.5522
Model.03 2 10 381.0807 418.6065 -180.5404 1 vs 2 40.02358 <.00014.4.1 Hypthesis Test 04 : results
- the very significant result (p<.0001) indicated that Model.03(with the pooled heterogeneous residual variances) is better than Model.01
4.5 Hypthesis Test 05 :
- Can any non-significant fixed effects be moved like the interaction treatment:sex?
numDF denDF F-value p-value
(Intercept) 1 292 813.9441 <.0001
treatment 2 23 8.6436 0.0016
sex 1 292 19.2504 <.0001
litsize 1 23 49.3306 <.0001
treatment:sex 2 292 0.3167 0.72884.5.1 Hypthesis Test 05 : results
- The non-significant interaction (p=0.7288) result indicated that the fixed effects of the \(treatment*sex\) should be removed from Model.03
attach(ratpupcsv)
Model.03M <- lme(weight ~ treatment + sex + litsize, random = ~1 | litter, method = "REML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.03M)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
372.2784 402.3497 -178.1392
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3146374 0.5144324
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5889108
Fixed effects: weight ~ treatment + sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.984202 0.27296903 294 29.249478 0.0000
treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001
treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091
sexMale 0.343431 0.04204323 294 8.168531 0.0000
litsize -0.130681 0.01855194 23 -7.044036 0.0000
Correlation:
(Intr) trtmnH trtmnL sexMal
treatmentHigh -0.613
treatmentLow -0.355 0.464
sexMale -0.051 0.006 0.035
litsize -0.910 0.403 0.064 -0.043
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992
Number of Observations: 322
Number of Groups: 27 4.6 Hypthesis Test 06 :
- Is there an effect of treatment on birth weight?
- Remember to use maximum likelihood (ML) for mixed models with different fixed effects structures if anova() is used for nested models comparsion
4.6.1 REML
attach(ratpupcsv)
Model.06_WTrm <- lme(weight ~ treatment + sex + litsize, random = ~1 | litter, method = "REML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.06_WTrm)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
372.2784 402.3497 -178.1392
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3146374 0.5144324
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5889108
Fixed effects: weight ~ treatment + sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.984202 0.27296903 294 29.249478 0.0000
treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001
treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091
sexMale 0.343431 0.04204323 294 8.168531 0.0000
litsize -0.130681 0.01855194 23 -7.044036 0.0000
Correlation:
(Intr) trtmnH trtmnL sexMal
treatmentHigh -0.613
treatmentLow -0.355 0.464
sexMale -0.051 0.006 0.035
litsize -0.910 0.403 0.064 -0.043
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992
Number of Observations: 322
Number of Groups: 27 Model.06_WNOTrm <- lme(weight ~ sex + litsize, random = ~1 | litter, method = "REML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.06_WNOTrm)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
381.7586 404.3497 -184.8793
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.4381779 0.5150813
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5881025
Fixed effects: weight ~ sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.211043 0.28139850 294 25.625733 0e+00
sexMale 0.348955 0.04213828 294 8.281194 0e+00
litsize -0.099343 0.02212882 25 -4.489285 1e-04
Correlation:
(Intr) sexMal
sexMale -0.043
litsize -0.947 -0.035
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.81957536 -0.51656194 0.02117273 0.55878848 2.53064896
Number of Observations: 322
Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-value
Model.06_WNOTrm 1 6 381.7586 404.3497 -184.8793
Model.06_WTrm 2 8 372.2784 402.3497 -178.1392 1 vs 2 13.48015 0.00124.6.1.1 ML
attach(ratpupcsv)
Model.06_WTrm <- lme(weight ~ treatment + sex + litsize, random = ~1 | litter, method = "ML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.06_WTrm)Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
353.7734 383.9698 -168.8867
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.2884174 0.5137273
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5881703
Fixed effects: weight ~ treatment + sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.986093 0.25810690 294 30.941030 0.0000
treatmentHigh -0.864775 0.17148477 23 -5.042866 0.0000
treatmentLow -0.434407 0.14242656 23 -3.050046 0.0057
sexMale 0.341896 0.04223590 294 8.094909 0.0000
litsize -0.130701 0.01750951 23 -7.464595 0.0000
Correlation:
(Intr) trtmnH trtmnL sexMal
treatmentHigh -0.616
treatmentLow -0.355 0.467
sexMale -0.054 0.006 0.037
litsize -0.911 0.407 0.065 -0.046
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.98845212 -0.52758888 0.02174839 0.55158792 2.59709135
Number of Observations: 322
Number of Groups: 27 Model.06_WNOTrm <- lme(weight ~ sex + litsize, random = ~1 | litter, method = "ML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.06_WNOTrm)Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
368.3706 391.0179 -178.1853
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.4209366 0.5146485
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5871372
Fixed effects: weight ~ sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.208088 0.27310606 294 26.392998 0e+00
sexMale 0.348624 0.04223474 294 8.254428 0e+00
litsize -0.099158 0.02146517 25 -4.619464 1e-04
Correlation:
(Intr) sexMal
sexMale -0.044
litsize -0.947 -0.036
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.81693672 -0.52012175 0.01677768 0.56092789 2.54095453
Number of Observations: 322
Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-value
Model.06_WNOTrm 1 6 368.3706 391.0179 -178.1853
Model.06_WTrm 2 8 353.7734 383.9698 -168.8867 1 vs 2 18.59723 1e-044.6.2 Hypthesis Test 06 : results
The results with “ML” and “REML” indicated that the significant treatment effects is very significant in Model.06_WTrm (p< 1e-04)
It should be aware that since Model.06_WTrm and Model.06_WNOTrm are under nested models specfication with one fixed parameter difference(treatment), therefore, the result of testing significance for both two approach are same. However, if their models ar not in nested specification, ML approach is more valid than the REML for two unnested model specification.
5 The first best model should be : Model.FB
attach(ratpupcsv)
Model.FB <- lme(weight ~ treatment + sex + litsize, random = ~1 | litter, method = "REML",
weights = varIdent(form = ~1 | trtgrp))
summary(Model.FB)Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
372.2784 402.3497 -178.1392
Random effects:
Formula: ~1 | litter
(Intercept) Residual
StdDev: 0.3146374 0.5144324
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | trtgrp
Parameter estimates:
1 2
1.0000000 0.5889108
Fixed effects: weight ~ treatment + sex + litsize
Value Std.Error DF t-value p-value
(Intercept) 7.984202 0.27296903 294 29.249478 0.0000
treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001
treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091
sexMale 0.343431 0.04204323 294 8.168531 0.0000
litsize -0.130681 0.01855194 23 -7.044036 0.0000
Correlation:
(Intr) trtmnH trtmnL sexMal
treatmentHigh -0.613
treatmentLow -0.355 0.464
sexMale -0.051 0.006 0.035
litsize -0.910 0.403 0.064 -0.043
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992
Number of Observations: 322
Number of Groups: 27 numDF denDF F-value p-value
(Intercept) 1 294 855.5320 <.0001
treatment 2 23 11.3870 4e-04
sex 1 294 66.7249 <.0001
litsize 1 23 49.6184 <.00016 Model diagnostics
6.1 Histogram
attach(ratpupcsv)
res <- resid(Model.FB)
ratred <- data.frame(ratpupcsv, res)
histogram(~res | factor(trtgrp), data = ratred, layout = c(2, 1), aspect = 2, xlab = "Residual",
strip = strip.custom(factor.levels = c("Control", "High/Low")))
6.2 Quantile-quantile plots
attach(ratpupcsv)
qqnorm(Model.FB, ~resid(.) | factor(trtgrp), layout = c(2, 1), aspect = 2, id = 0.05,
strip = strip.custom(factor.levels = c("Control", "High/Low")))
6.3 Homogeneity of variance plot
- fitted values vs. residuals
- the plot shows that residual variance is similar for small and large fitted values
attach(ratpupcsv)
plot(Model.FB, resid(.) ~ fitted(.) | factor(trtgrp), layout = c(2, 1), aspect = 2,
abline = 0, strip = strip.custom(factor.levels = c("Control", "High/Low")))
6.4 Post-hoc test between all treatment levels.
library(multcomp)
Model.FB <- lme(weight ~ treatment + sex + litsize, random = ~1 | litter, method = "REML",
ratpupcsv, weights = varIdent(form = ~1 | trtgrp))
summary(glht(Model.FB, linfct = mcp(treatment = "Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = weight ~ treatment + sex + litsize, data = ratpupcsv,
random = ~1 | litter, weights = varIdent(form = ~1 | trtgrp),
method = "REML")
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
High - Control == 0 -0.8623 0.1829 -4.714 <0.001 ***
Low - Control == 0 -0.4337 0.1523 -2.848 0.0120 *
Low - High == 0 0.4286 0.1755 2.442 0.0384 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
7 Linear Mixed-Effect Model: Package lme4
lmer.model <- lmer(weight ~ treatment + sex + litsize + treatment * sex + (1 | litter),
REML = FALSE, ratpupcsv)
summary(lmer.model)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: weight ~ treatment + sex + litsize + treatment * sex + (1 | litter)
Data: ratpupcsv
AIC BIC logLik deviance df.resid
395.8 429.8 -188.9 377.8 313
Scaled residuals:
Min 1Q Median 3Q Max
-7.5188 -0.5045 0.0254 0.5880 3.0233
Random effects:
Groups Name Variance Std.Dev.
litter (Intercept) 0.0807 0.2841
Residual 0.1617 0.4021
Number of obs: 322, groups: litter, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 7.91056 0.25784 41.07091 30.680 < 2e-16 ***
treatmentHigh -0.79972 0.18207 40.29948 -4.392 7.92e-05 ***
treatmentLow -0.38343 0.14911 35.83419 -2.572 0.0144 *
sexMale 0.41055 0.07273 299.07067 5.645 3.84e-08 ***
litsize -0.12821 0.01755 38.34969 -7.305 9.05e-09 ***
treatmentHigh:sexMale -0.11001 0.13083 308.82323 -0.841 0.4011
treatmentLow:sexMale -0.08414 0.10502 302.81404 -0.801 0.4236
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmntHgh -0.580
treatmentLw -0.335 0.433
sexMale -0.165 0.234 0.286
litsize -0.911 0.372 0.045 -0.001
trtmntHgh:M 0.112 -0.382 -0.160 -0.556 -0.022
trtmntLw:sM 0.149 -0.176 -0.368 -0.693 -0.037 0.386model1.fit.lmer <- lmer(weight ~ treatment + sex + litsize + treatment * sex + (1 |
litter), ratpupcsv, REML = T)
summary(model1.fit.lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: weight ~ treatment + sex + litsize + treatment * sex + (1 | litter)
Data: ratpupcsv
REML criterion at convergence: 401.1
Scaled residuals:
Min 1Q Median 3Q Max
-7.4725 -0.5001 0.0291 0.5735 3.0096
Random effects:
Groups Name Variance Std.Dev.
litter (Intercept) 0.09652 0.3107
Residual 0.16349 0.4043
Number of obs: 322, groups: litter, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 7.91165 0.27496 33.70992 28.773 < 2e-16 ***
treatmentHigh -0.79903 0.19429 32.88474 -4.112 0.000245 ***
treatmentLow -0.38317 0.15967 29.56842 -2.400 0.022919 *
sexMale 0.41169 0.07315 295.30143 5.628 4.24e-08 ***
litsize -0.12838 0.01875 31.79751 -6.846 9.94e-08 ***
treatmentHigh:sexMale -0.10702 0.13176 303.74417 -0.812 0.417291
treatmentLow:sexMale -0.08387 0.10568 298.58691 -0.794 0.428077
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmnH trtmnL sexMal litsiz trtH:M
treatmntHgh -0.581
treatmentLw -0.336 0.434
sexMale -0.155 0.221 0.269
litsize -0.910 0.372 0.045 -0.001
trtmntHgh:M 0.105 -0.360 -0.150 -0.555 -0.020
trtmntLw:sM 0.141 -0.166 -0.345 -0.692 -0.036 0.385Marginal Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
treatment 2.8416 1.4208 2 31.278 8.6905 0.0009965 ***
sex 5.1778 5.1778 1 295.301 31.6708 4.244e-08 ***
litsize 7.6619 7.6619 1 31.798 46.8652 9.937e-08 ***
treatment:sex 0.1522 0.0761 2 302.303 0.4656 0.6282199
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.1 Hypthesis Test 01 : with lem4
ANOVA-like table for random-effects: Single term deletions
Model:
weight ~ treatment + sex + litsize + (1 | litter) + treatment:sex
npar logLik AIC LRT Df Pr(>Chisq)
<none> 9 -200.55 419.10
(1 | litter) 8 -245.25 506.51 89.406 1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of Deviance Table (Type II Wald chisquare tests)
Response: weight
Chisq Df Pr(>Chisq)
treatment 23.3610 2 8.457e-06 ***
sex 56.9064 1 4.571e-14 ***
litsize 46.8652 1 7.604e-12 ***
treatment:sex 0.9312 2 0.6278
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model2.fit.lmer <- lmer(weight ~ treatment + sex + litsize + (1 | litter), ratpupcsv,
REML = F)
model3.fit.lmer <- lmer(weight ~ treatment + sex + (1 | litter), ratpupcsv, REML = F)
model4.fit.lmer <- lmer(weight ~ treatment + litsize + (1 | litter), ratpupcsv, REML = F)
model5.fit.lmer <- lmer(weight ~ sex + litsize + (1 | litter), ratpupcsv, REML = F)
anova(model1.fit.lmer, model2.fit.lmer) #testing interaction effectData: ratpupcsv
Models:
model2.fit.lmer: weight ~ treatment + sex + litsize + (1 | litter)
model1.fit.lmer: weight ~ treatment + sex + litsize + treatment * sex + (1 | litter)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model2.fit.lmer 7 392.79 419.21 -189.39 378.79
model1.fit.lmer 9 395.81 429.78 -188.91 377.81 0.9721 2 0.615Data: ratpupcsv
Models:
model3.fit.lmer: weight ~ treatment + sex + (1 | litter)
model2.fit.lmer: weight ~ treatment + sex + litsize + (1 | litter)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model3.fit.lmer 6 422.98 445.62 -205.49 410.98
model2.fit.lmer 7 392.79 419.21 -189.39 378.79 32.19 1 1.398e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Data: ratpupcsv
Models:
model4.fit.lmer: weight ~ treatment + litsize + (1 | litter)
model2.fit.lmer: weight ~ treatment + sex + litsize + (1 | litter)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model4.fit.lmer 6 442.05 464.69 -215.02 430.05
model2.fit.lmer 7 392.79 419.21 -189.39 378.79 51.262 1 8.084e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Data: ratpupcsv
Models:
model5.fit.lmer: weight ~ sex + litsize + (1 | litter)
model2.fit.lmer: weight ~ treatment + sex + litsize + (1 | litter)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model5.fit.lmer 5 407.48 426.35 -198.74 397.48
model2.fit.lmer 7 392.79 419.21 -189.39 378.79 18.695 2 8.716e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
7.3 check residuals for normality
model6.fit.lmer <- lmer(weight ~ treatment + sex + litsize + (1 | litter), ratpupcsv,
REML = T)
summary(glht(model6.fit.lmer, linfct = mcp(treatment = "Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = weight ~ treatment + sex + litsize + (1 | litter),
data = ratpupcsv, REML = T)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
High - Control == 0 -0.8587 0.1818 -4.723 <0.001 ***
Low - Control == 0 -0.4285 0.1504 -2.849 0.0119 *
Low - High == 0 0.4302 0.1804 2.385 0.0446 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)
Reference
Crawley, Michael J. 2013. “The R Book Second Edition.” John Wiley & Sons.
Harville, David A. 1977. “Maximum Likelihood Approaches to Variance Component Estimation and to Related Problems.” Journal of the American Statistical Association 72 (358): 320–38. www.jstor.org/stable/2286796.
JJ, Faraway. 2006. “Binomial Data. Extending the Linear Model with R.” Chapman & Hall/CRC.
Laird, Nan M., and James H. Ware. 1982. “Random-Effects Models for Longitudinal Data.” Biometrics 38 (4): 963–74. www.jstor.org/stable/2529876.
Pinero, Jose, and Douglas Bates. 2000. “Mixed-Effects Models in S and S-Plus (Statistics and Computing).” Springer, New York.
Satoh, Masahiro. 2018. “An Alternative Derivation Method of Mixed Model Equations from Best Linear Unbiased Prediction (Blup) and Restricted Blup of Breeding Values Not Using Maximum Likelihood.” Animal Science Journal 89 (6): 876–79.
Shayle R. Searle, Charles E. McCulloch, George Casella. 1992. Variance Components. NewJersey: Wiley-Interscience.
2020-02-28