Example 5.2 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup(p-164)
Source:R/Exam5.2.R
Exam5.2.RdExam5.2 three factor main effects only design
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
DataSet5.2$a <- factor( x = DataSet5.2$a)
DataSet5.2$b <- factor( x = DataSet5.2$b)
DataSet5.2$c <- factor(x = DataSet5.2$c)
##---first adding factor a in model
Exam5.2.lm1 <- lm(formula = y~ a, data = DataSet5.2)
summary(Exam5.2.lm1)
#>
#> Call:
#> lm(formula = y ~ a, data = DataSet5.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.5833 -1.2833 0.4333 0.9333 1.1167
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 9.7667 0.6913 14.13 2.11e-06 ***
#> a1 1.6167 0.8466 1.91 0.0978 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.197 on 7 degrees of freedom
#> Multiple R-squared: 0.3425, Adjusted R-squared: 0.2486
#> F-statistic: 3.646 on 1 and 7 DF, p-value: 0.09783
#>
library(parameters)
model_parameters(Exam5.2.lm1)
#> Parameter | Coefficient | SE | 95% CI | t(7) | p
#> ------------------------------------------------------------------
#> (Intercept) | 9.77 | 0.69 | [ 8.13, 11.40] | 14.13 | < .001
#> a [1] | 1.62 | 0.85 | [-0.39, 3.62] | 1.91 | 0.098
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
library(emmeans)
##---A first
emmeans(object = Exam5.2.lm1, specs = ~a)
#> a emmean SE df lower.CL upper.CL
#> 0 9.77 0.691 7 8.13 11.4
#> 1 11.38 0.489 7 10.23 12.5
#>
#> Confidence level used: 0.95
contrast(emmeans(object = Exam5.2.lm1, specs = ~a), method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> a0 - a1 -1.62 0.847 7 -1.910 0.0978
#>
anova(object = Exam5.2.lm1)
#> Analysis of Variance Table
#>
#> Response: y
#> Df Sum Sq Mean Sq F value Pr(>F)
#> a 1 5.2272 5.2272 3.6463 0.09783 .
#> Residuals 7 10.0350 1.4336
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##---then adding factor b in model
Exam5.2.lm2 <- lm(formula = y~ a + b, data = DataSet5.2)
summary(Exam5.2.lm2)
#>
#> Call:
#> lm(formula = y ~ a + b, data = DataSet5.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.8667 -0.7833 0.5667 0.7667 1.1833
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 10.2667 0.9083 11.303 2.87e-05 ***
#> a1 1.3667 0.9083 1.505 0.183
#> b1 -0.7500 0.8617 -0.870 0.418
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.219 on 6 degrees of freedom
#> Multiple R-squared: 0.4162, Adjusted R-squared: 0.2216
#> F-statistic: 2.139 on 2 and 6 DF, p-value: 0.199
#>
model_parameters(Exam5.2.lm2)
#> Parameter | Coefficient | SE | 95% CI | t(6) | p
#> ------------------------------------------------------------------
#> (Intercept) | 10.27 | 0.91 | [ 8.04, 12.49] | 11.30 | < .001
#> a [1] | 1.37 | 0.91 | [-0.86, 3.59] | 1.50 | 0.183
#> b [1] | -0.75 | 0.86 | [-2.86, 1.36] | -0.87 | 0.418
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
emmeans(object = Exam5.2.lm2, specs = ~b)
#> b emmean SE df lower.CL upper.CL
#> 0 10.9 0.609 6 9.46 12.4
#> 1 10.2 0.609 6 8.71 11.7
#>
#> Results are averaged over the levels of: a
#> Confidence level used: 0.95
contrast(emmeans(object = Exam5.2.lm2, specs = ~b), method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> b0 - b1 0.75 0.862 6 0.870 0.4175
#>
#> Results are averaged over the levels of: a
anova(object = Exam5.2.lm2)
#> Analysis of Variance Table
#>
#> Response: y
#> Df Sum Sq Mean Sq F value Pr(>F)
#> a 1 5.2272 5.2272 3.5200 0.1097
#> b 1 1.1250 1.1250 0.7576 0.4175
#> Residuals 6 8.9100 1.4850
##---then adding factor c in model
Exam5.2.lm3 <- lm(formula = y~ a + b + c, data = DataSet5.2)
summary(Exam5.2.lm3)
#>
#> Call:
#> lm(formula = y ~ a + b + c, data = DataSet5.2)
#>
#> Residuals:
#> 1 2 3 4 5 6 7
#> -1.500e-01 1.500e-01 3.000e-01 1.000e-01 -8.000e-01 4.000e-01 -2.500e-01
#> 8 9
#> 2.500e-01 4.857e-17
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 10.5000 0.3469 30.264 7.39e-07 ***
#> a1 1.6000 0.3469 4.612 0.00578 **
#> b1 -0.0500 0.3469 -0.144 0.89104
#> c1 -2.1000 0.3469 -6.053 0.00178 **
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 0.4626 on 5 degrees of freedom
#> Multiple R-squared: 0.9299, Adjusted R-squared: 0.8878
#> F-statistic: 22.11 on 3 and 5 DF, p-value: 0.002584
#>
model_parameters(Exam5.2.lm3)
#> Parameter | Coefficient | SE | 95% CI | t(5) | p
#> ------------------------------------------------------------------
#> (Intercept) | 10.50 | 0.35 | [ 9.61, 11.39] | 30.26 | < .001
#> a [1] | 1.60 | 0.35 | [ 0.71, 2.49] | 4.61 | 0.006
#> b [1] | -0.05 | 0.35 | [-0.94, 0.84] | -0.14 | 0.891
#> c [1] | -2.10 | 0.35 | [-2.99, -1.21] | -6.05 | 0.002
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
emmeans(object = Exam5.2.lm3, specs = ~c)
#> c emmean SE df lower.CL upper.CL
#> 0 11.28 0.200 5 10.76 11.8
#> 1 9.18 0.283 5 8.45 9.9
#>
#> Results are averaged over the levels of: a, b
#> Confidence level used: 0.95
contrast(emmeans(object = Exam5.2.lm3, specs = ~c), method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> c0 - c1 2.1 0.347 5 6.053 0.0018
#>
#> Results are averaged over the levels of: a, b
anova(object = Exam5.2.lm3)
#> Analysis of Variance Table
#>
#> Response: y
#> Df Sum Sq Mean Sq F value Pr(>F)
#> a 1 5.2272 5.2272 24.426 0.004314 **
#> b 1 1.1250 1.1250 5.257 0.070401 .
#> c 1 7.8400 7.8400 36.636 0.001775 **
#> Residuals 5 1.0700 0.2140
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##---Now Change the order and add b first in model
Exam5.2.lm4 <- lm(formula = y ~ b, data = DataSet5.2)
summary(Exam5.2.lm4)
#>
#> Call:
#> lm(formula = y ~ b, data = DataSet5.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.96 -0.10 0.00 0.84 1.14
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 11.3600 0.5921 19.185 2.6e-07 ***
#> b1 -1.1600 0.8882 -1.306 0.233
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.324 on 7 degrees of freedom
#> Multiple R-squared: 0.1959, Adjusted R-squared: 0.08105
#> F-statistic: 1.706 on 1 and 7 DF, p-value: 0.2328
#>
model_parameters(Exam5.2.lm4)
#> Parameter | Coefficient | SE | 95% CI | t(7) | p
#> ------------------------------------------------------------------
#> (Intercept) | 11.36 | 0.59 | [ 9.96, 12.76] | 19.18 | < .001
#> b [1] | -1.16 | 0.89 | [-3.26, 0.94] | -1.31 | 0.233
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
emmeans(object = Exam5.2.lm4, specs = ~b)
#> b emmean SE df lower.CL upper.CL
#> 0 11.4 0.592 7 9.96 12.8
#> 1 10.2 0.662 7 8.63 11.8
#>
#> Confidence level used: 0.95
contrast(emmeans(object = Exam5.2.lm4, specs = ~b), method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> b0 - b1 1.16 0.888 7 1.306 0.2328
#>
anova(object = Exam5.2.lm4)
#> Analysis of Variance Table
#>
#> Response: y
#> Df Sum Sq Mean Sq F value Pr(>F)
#> b 1 2.9902 2.9902 1.7056 0.2328
#> Residuals 7 12.2720 1.7531
##---then adding factor a in model
Exam5.2.lm5 <- lm(formula = y ~ b + a, data = DataSet5.2)
summary(Exam5.2.lm5)
#>
#> Call:
#> lm(formula = y ~ b + a, data = DataSet5.2)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.8667 -0.7833 0.5667 0.7667 1.1833
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 10.2667 0.9083 11.303 2.87e-05 ***
#> b1 -0.7500 0.8617 -0.870 0.418
#> a1 1.3667 0.9083 1.505 0.183
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 1.219 on 6 degrees of freedom
#> Multiple R-squared: 0.4162, Adjusted R-squared: 0.2216
#> F-statistic: 2.139 on 2 and 6 DF, p-value: 0.199
#>
model_parameters(Exam5.2.lm5)
#> Parameter | Coefficient | SE | 95% CI | t(6) | p
#> ------------------------------------------------------------------
#> (Intercept) | 10.27 | 0.91 | [ 8.04, 12.49] | 11.30 | < .001
#> b [1] | -0.75 | 0.86 | [-2.86, 1.36] | -0.87 | 0.418
#> a [1] | 1.37 | 0.91 | [-0.86, 3.59] | 1.50 | 0.183
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
emmeans(object = Exam5.2.lm5, specs = ~a)
#> a emmean SE df lower.CL upper.CL
#> 0 9.89 0.718 6 8.13 11.6
#> 1 11.26 0.518 6 9.99 12.5
#>
#> Results are averaged over the levels of: b
#> Confidence level used: 0.95
contrast(emmeans(object = Exam5.2.lm5, specs = ~a), method = "pairwise")
#> contrast estimate SE df t.ratio p.value
#> a0 - a1 -1.37 0.908 6 -1.505 0.1831
#>
#> Results are averaged over the levels of: b
anova(object = Exam5.2.lm5)
#> Analysis of Variance Table
#>
#> Response: y
#> Df Sum Sq Mean Sq F value Pr(>F)
#> b 1 2.9902 2.9902 2.0136 0.2057
#> a 1 3.3620 3.3620 2.2640 0.1831
#> Residuals 6 8.9100 1.4850