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.
'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")# Visualisations: ratpupcsv
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.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.5522Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3106722 0.404337Fixed effects: weight ~ treatment + sex + litsize + treatment * sex Value Std.Error DF t-value p-value(Intercept) 7.911652 0.27496390 292 28.773422 0.0000treatmentHigh -0.799034 0.19429415 23 -4.112495 0.0004treatmentLow -0.383174 0.15967364 23 -2.399731 0.0249sexMale 0.411688 0.07315410 292 5.627679 0.0000litsize -0.128382 0.01875336 23 -6.845819 0.0000treatmentHigh:sexMale -0.107023 0.13176318 292 -0.812239 0.4173treatmentLow:sexMale -0.083866 0.10568189 292 -0.793568 0.4281 Correlation: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmentHigh -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.385Standardized Within-Group Residuals: Min Q1 Med Q3 Max -7.47250744 -0.50014749 0.02911668 0.57348178 3.00962055 Number of Observations: 322Number of Groups: 27 numDF denDF F-value p-value(Intercept) 1 292 9093.772 <.0001treatment 2 23 5.082 0.0149sex 1 292 52.602 <.0001litsize 1 23 47.374 <.0001treatment:sex 2 292 0.466 0.6282 numDF denDF F-value p-value(Intercept) 1 292 827.9098 <.0001treatment 2 23 8.6905 0.0015sex 1 292 31.6708 <.0001litsize 1 23 46.8652 <.0001treatment: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.255Coefficients: Value Std.Error t-value p-value(Intercept) 7.900732 0.16229542 48.68118 0.0000treatmentHigh -0.807294 0.12051060 -6.69895 0.0000treatmentLow -0.381885 0.09273467 -4.11804 0.0000sexMale 0.339576 0.08902602 3.81435 0.0002litsize -0.124188 0.01024669 -12.11978 0.0000treatmentHigh:sexMale -0.178572 0.15325437 -1.16520 0.2448treatmentLow:sexMale -0.074357 0.12639429 -0.58829 0.5568 Correlation: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmentHigh -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.413Standardized 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-valueModel.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.9423Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3134846 0.5147948Variance 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.0000treatmentHigh -0.808610 0.19628137 23 -4.119649 0.0004treatmentLow -0.390279 0.16301558 23 -2.394121 0.0252sexMale 0.408131 0.09303486 292 4.386865 0.0000litsize -0.130007 0.01848708 23 -7.032332 0.0000treatmentHigh:sexMale -0.094666 0.12919527 292 -0.732737 0.4643treatmentLow:sexMale -0.076013 0.10811858 292 -0.703053 0.4826 Correlation: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmentHigh -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.620Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.88670114 -0.52493419 0.02123518 0.57307286 2.56409983 Number of Observations: 322Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-valueModel.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"] <- 1ratpupcsv$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.5404Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3145679 0.5147878Variance 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.0000treatmentHigh -0.809818 0.19550851 23 -4.142109 0.0004treatmentLow -0.390200 0.16383300 23 -2.381691 0.0259sexMale 0.408195 0.09303539 292 4.387529 0.0000litsize -0.130383 0.01856367 23 -7.023574 0.0000treatmentHigh:sexMale -0.092026 0.12461723 292 -0.738473 0.4608treatmentLow:sexMale -0.076397 0.10939797 292 -0.698337 0.4855 Correlation: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmentHigh -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.635Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.88725915 -0.52398492 0.01731335 0.56412339 2.55014733 Number of Observations: 322Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-valueModel.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-valueModel.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 <.0001treatment 2 23 8.6436 0.0016sex 1 292 19.2504 <.0001litsize 1 23 49.3306 <.0001treatment: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.1392Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3146374 0.5144324Variance 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.0000treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091sexMale 0.343431 0.04204323 294 8.168531 0.0000litsize -0.130681 0.01855194 23 -7.044036 0.0000 Correlation: (Intr) trtmnH trtmnL sexMaltreatmentHigh -0.613 treatmentLow -0.355 0.464 sexMale -0.051 0.006 0.035 litsize -0.910 0.403 0.064 -0.043Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992 Number of Observations: 322Number 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.1392Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3146374 0.5144324Variance 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.0000treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091sexMale 0.343431 0.04204323 294 8.168531 0.0000litsize -0.130681 0.01855194 23 -7.044036 0.0000 Correlation: (Intr) trtmnH trtmnL sexMaltreatmentHigh -0.613 treatmentLow -0.355 0.464 sexMale -0.051 0.006 0.035 litsize -0.910 0.403 0.064 -0.043Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992 Number of Observations: 322Number 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.8793Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.4381779 0.5150813Variance 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+00sexMale 0.348955 0.04213828 294 8.281194 0e+00litsize -0.099343 0.02212882 25 -4.489285 1e-04 Correlation: (Intr) sexMalsexMale -0.043 litsize -0.947 -0.035Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.81957536 -0.51656194 0.02117273 0.55878848 2.53064896 Number of Observations: 322Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-valueModel.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.8867Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.2884174 0.5137273Variance 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.0000treatmentHigh -0.864775 0.17148477 23 -5.042866 0.0000treatmentLow -0.434407 0.14242656 23 -3.050046 0.0057sexMale 0.341896 0.04223590 294 8.094909 0.0000litsize -0.130701 0.01750951 23 -7.464595 0.0000 Correlation: (Intr) trtmnH trtmnL sexMaltreatmentHigh -0.616 treatmentLow -0.355 0.467 sexMale -0.054 0.006 0.037 litsize -0.911 0.407 0.065 -0.046Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.98845212 -0.52758888 0.02174839 0.55158792 2.59709135 Number of Observations: 322Number 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.1853Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.4209366 0.5146485Variance 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+00sexMale 0.348624 0.04223474 294 8.254428 0e+00litsize -0.099158 0.02146517 25 -4.619464 1e-04 Correlation: (Intr) sexMalsexMale -0.044 litsize -0.947 -0.036Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.81693672 -0.52012175 0.01677768 0.56092789 2.54095453 Number of Observations: 322Number of Groups: 27 Model df AIC BIC logLik Test L.Ratio p-valueModel.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.
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.1392Random effects: Formula: ~1 | litter (Intercept) ResidualStdDev: 0.3146374 0.5144324Variance 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.0000treatmentHigh -0.862268 0.18293359 23 -4.713556 0.0001treatmentLow -0.433663 0.15226167 23 -2.848140 0.0091sexMale 0.343431 0.04204323 294 8.168531 0.0000litsize -0.130681 0.01855194 23 -7.044036 0.0000 Correlation: (Intr) trtmnH trtmnL sexMaltreatmentHigh -0.613 treatmentLow -0.355 0.464 sexMale -0.051 0.006 0.035 litsize -0.910 0.403 0.064 -0.043Standardized Within-Group Residuals: Min Q1 Med Q3 Max -5.97351495 -0.53695365 0.01508652 0.54234475 2.58286992 Number of Observations: 322Number of Groups: 27 numDF denDF F-value p-value(Intercept) 1 294 855.5320 <.0001treatment 2 23 11.3870 4e-04sex 1 294 66.7249 <.0001litsize 1 23 49.6184 <.00016.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 hom*ogeneity 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 HypothesesMultiple Comparisons of Means: Tukey ContrastsFit: 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)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, 27Fixed 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 ' ' 1Correlation of Fixed Effects: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmntHgh -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: ratpupcsvREML criterion at convergence: 401.1Scaled 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, 27Fixed 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 ' ' 1Correlation of Fixed Effects: (Intr) trtmnH trtmnL sexMal litsiz trtH:MtreatmntHgh -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 deletionsModel: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: ratpupcsvModels: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: ratpupcsvModels: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: ratpupcsvModels: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: ratpupcsvModels: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 ' ' 17.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 HypothesesMultiple Comparisons of Means: Tukey ContrastsFit: 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)
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.
