Example 8.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup (p-250)
Source:R/Exam8.1.R
Exam8.1.RdExam8.1 Nested factorial structure
References
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
Author
Muhammad Yaseen (myaseen208@gmail.com)
Adeela Munawar (adeela.uaf@gmail.com)
Examples
data(DataSet8.1)
DataSet8.1$block <- factor(x = DataSet8.1$block)
DataSet8.1$set <- factor(x = DataSet8.1$set)
DataSet8.1$trt <- factor(x = DataSet8.1$trt)
library(lmerTest)
Exam8.1Lmer <- lmer(y ~ set + trt %in% set + (1|set/block), DataSet8.1)
#> fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients
#> Warning: Model is nearly unidentifiable: large eigenvalue ratio
#> - Rescale variables?
#> Warning: Model may not have converged with 1 eigenvalue close to zero: 8.9e-16
summary(Exam8.1Lmer)
#> Linear mixed model fit by REML. t-tests use Satterthwaite's method [
#> lmerModLmerTest]
#> Formula: y ~ set + trt %in% set + (1 | set/block)
#> Data: DataSet8.1
#>
#> REML criterion at convergence: 170.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -1.5721 -0.2644 -0.1207 0.3136 2.0303
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> block:set (Intercept) 60.5499 7.7814
#> set (Intercept) 0.1599 0.3999
#> Residual 22.7506 4.7698
#> Number of obs: 30, groups: block:set, 10; set, 2
#>
#> Fixed effects:
#> Estimate Std. Error df t value Pr(>|t|)
#> (Intercept) 8.200 4.101 11.863 1.999 0.069000 .
#> set2 8.960 5.800 11.863 1.545 0.148635
#> set1:trt2 0.700 3.017 16.000 0.232 0.819445
#> set1:trt3 0.540 3.017 16.000 0.179 0.860180
#> set2:trt4 -12.720 3.017 16.000 -4.217 0.000655 ***
#> set2:trt5 -9.880 3.017 16.000 -3.275 0.004762 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Correlation of Fixed Effects:
#> (Intr) set2 st1:t2 st1:t3 st2:t4
#> set2 -0.707
#> set1:trt2 -0.368 0.260
#> set1:trt3 -0.368 0.260 0.500
#> set2:trt4 0.000 -0.260 0.000 0.000
#> set2:trt5 0.000 -0.260 0.000 0.000 0.500
#> fit warnings:
#> fixed-effect model matrix is rank deficient so dropping 6 columns / coefficients
#> optimizer (nloptwrap) convergence code: 0 (OK)
#> Model is nearly unidentifiable: large eigenvalue ratio
#> - Rescale variables?
#>
anova(Exam8.1Lmer)
#> Missing cells for: set2:trt1, set2:trt2, set2:trt3, set1:trt4, set1:trt5, set1:trt6.
#> Interpret type III hypotheses with care.
#> Type III Analysis of Variance Table with Satterthwaite's method
#> Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
#> set 54.29 54.294 1 11.863 2.3865 0.148635
#> set:trt 447.14 111.786 4 16.000 4.9135 0.008898 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
library(emmeans)
emmeans(object = Exam8.1Lmer, specs = ~trt|set)
#> NOTE: A nesting structure was detected in the fitted model:
#> trt %in% set
#> set = 1:
#> trt emmean SE df lower.CL upper.CL
#> 1 8.20 4.1 11.9 -0.7446 17.1
#> 2 8.90 4.1 11.9 -0.0446 17.8
#> 3 8.74 4.1 11.9 -0.2046 17.7
#>
#> set = 2:
#> trt emmean SE df lower.CL upper.CL
#> 4 4.44 4.1 11.9 -4.5046 13.4
#> 5 7.28 4.1 11.9 -1.6646 16.2
#> 6 17.16 4.1 11.9 8.2154 26.1
#>
#> Degrees-of-freedom method: kenward-roger
#> Confidence level used: 0.95
contrast(emmeans(object = Exam8.1Lmer, specs = ~trt|set), method = "pairwise", by = "set")
#> NOTE: A nesting structure was detected in the fitted model:
#> trt %in% set
#> set = 1:
#> contrast estimate SE df t.ratio p.value
#> trt1 - trt2 -0.70 3.02 16 -0.232 0.9999
#> trt1 - trt3 -0.54 3.02 16 -0.179 1.0000
#> trt2 - trt3 0.16 3.02 16 0.053 1.0000
#>
#> set = 2:
#> contrast estimate SE df t.ratio p.value
#> trt4 - trt5 -2.84 3.02 16 -0.941 0.9295
#> trt4 - trt6 -12.72 3.02 16 -4.217 0.0071
#> trt5 - trt6 -9.88 3.02 16 -3.275 0.0452
#>
#> Degrees-of-freedom method: kenward-roger
#> P value adjustment: tukey method for comparing a family of 6 estimates