Example 7.1 from Generalized Linear Mixed Models: Modern Concepts, Methods and Applications by Walter W. Stroup (p-213)
Exam7.1.RdExam7.1 explains multifactor models with all factors qualitative
References
Stroup, W. W. (2012). Generalized Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press.
See also
DataSet7.1
Examples
library(lsmeans) library(car) data(DataSet7.1) DataSet7.1$a <- factor(x = DataSet7.1$a) DataSet7.1$b <- factor(x = DataSet7.1$b) Exam7.1.lm1 <- lm( formula = y ~ a + b + a*b , data = DataSet7.1 # , subset # , weights # , na.action , method = "qr" , model = TRUE # , x = FALSE # , y = FALSE , qr = TRUE , singular.ok = TRUE , contrasts = NULL # , offset # , ... ) summary( Exam7.1.lm1 )#> #> Call: #> lm(formula = y ~ a + b + a * b, data = DataSet7.1, method = "qr", #> model = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL) #> #> Residuals: #> Min 1Q Median 3Q Max #> -10.50 -2.50 0.00 3.25 7.75 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 50.500 2.561 19.718 1.23e-13 *** #> a1 17.250 3.622 4.763 0.000156 *** #> b1 9.250 3.622 2.554 0.019936 * #> b2 21.000 3.622 5.798 1.71e-05 *** #> a1:b1 -9.750 5.122 -1.904 0.073089 . #> a1:b2 -17.750 5.122 -3.465 0.002761 ** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #> #> Residual standard error: 5.122 on 18 degrees of freedom #> Multiple R-squared: 0.7352, Adjusted R-squared: 0.6617 #> F-statistic: 9.997 on 5 and 18 DF, p-value: 0.0001049 #>anova(Exam7.1.lm1)#> Analysis of Variance Table #> #> Response: y #> Df Sum Sq Mean Sq F value Pr(>F) #> a 1 392.04 392.04 14.9428 0.0011330 ** #> b 2 603.25 301.62 11.4966 0.0006068 *** #> a:b 2 316.08 158.04 6.0238 0.0099348 ** #> Residuals 18 472.25 26.24 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1##---Result obtained as in SLICE statement in SAS for a0 & a1 library(phia) a0 <- list(a=c("0"=1)) phia::testInteractions(Exam7.1.lm1, custom=a0, across="b")#> F Test: #> P-value adjustment method: holm #> b1 b2 Df Sum of Sq F Pr(>F) #> a1 -21 -11.75 2 886.17 16.888 7.417e-05 *** #> Residuals 18 472.25 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1#> F Test: #> P-value adjustment method: holm #> b1 b2 Df Sum of Sq F Pr(>F) #> a1 -3.25 -3.75 2 33.17 0.6321 0.5429 #> Residuals 18 472.25##---Interaction plot lsmip( object = Exam7.1.lm1 , formula = a~b , ylab = "y Lsmeans" , main = "Lsmeans for a*b" )#------------------------------------------------------------- ## Individula least squares treatment means #------------------------------------------------------------- Lsm7.1 <- lsmeans::lsmeans( object = Exam7.1.lm1 , specs = ~a*b # , ... ) Lsm7.1#> a b lsmean SE df lower.CL upper.CL #> 0 0 50.50 2.56106 18 45.11941 55.88059 #> 1 0 67.75 2.56106 18 62.36941 73.13059 #> 0 1 59.75 2.56106 18 54.36941 65.13059 #> 1 1 67.25 2.56106 18 61.86941 72.63059 #> 0 2 71.50 2.56106 18 66.11941 76.88059 #> 1 2 71.00 2.56106 18 65.61941 76.38059 #> #> Confidence level used: 0.95##---Simpe effects comparison of interaction by a ## (but it doesn't give the same p-value as in article 7.4.2 page#215) SimpleEff7.1 <- lsmeans::lsmeans( object = Exam7.1.lm1 , specs = pairwise~b|a # , ... )$contrasts SimpleEff7.1#> a = 0: #> contrast estimate SE df t.ratio p.value #> 0 - 1 -9.25 3.621886 18 -2.554 0.0498 #> 0 - 2 -21.00 3.621886 18 -5.798 <.0001 #> 1 - 2 -11.75 3.621886 18 -3.244 0.0119 #> #> a = 1: #> contrast estimate SE df t.ratio p.value #> 0 - 1 0.50 3.621886 18 0.138 0.9896 #> 0 - 2 -3.25 3.621886 18 -0.897 0.6488 #> 1 - 2 -3.75 3.621886 18 -1.035 0.5649 #> #> P value adjustment: tukey method for comparing a family of 3 estimates##---Alternative method of pairwise comparisons by applying contrast ## coefficient (gives the same p-value as in 7.4.2) ContrastLsm7.1 <- lsmeans::contrast( Lsm7.1 , list ( c1 = c(1,0,-1,0,0,0) , c2 = c(1,0,0,0,-1,0) , c3 = c(0,0,1,0,-1,0) , c4 = c(0,1,0,-1,0,0) , c5 = c(0,1,0,0,0,-1) , c6 = c(0,1,0,0,-1,0) ) ) ContrastLsm7.1#> contrast estimate SE df t.ratio p.value #> c1 -9.25 3.621886 18 -2.554 0.0199 #> c2 -21.00 3.621886 18 -5.798 <.0001 #> c3 -11.75 3.621886 18 -3.244 0.0045 #> c4 0.50 3.621886 18 0.138 0.8917 #> c5 -3.25 3.621886 18 -0.897 0.3814 #> c6 -3.75 3.621886 18 -1.035 0.3142 #>##---Nested Model (page 216)---- Exam7.1.lm2 <- lm( formula = y ~ a + a %in% b , data = DataSet7.1 # , subset # , weights # , na.action , method = "qr" , model = TRUE # , x = FALSE # , y = FALSE , qr = TRUE , singular.ok = TRUE , contrasts = NULL # , offset # , ... ) summary( Exam7.1.lm2 )#> #> Call: #> lm(formula = y ~ a + a %in% b, data = DataSet7.1, method = "qr", #> model = TRUE, qr = TRUE, singular.ok = TRUE, contrasts = NULL) #> #> Residuals: #> Min 1Q Median 3Q Max #> -10.50 -2.50 0.00 3.25 7.75 #> #> Coefficients: #> Estimate Std. Error t value Pr(>|t|) #> (Intercept) 50.500 2.561 19.718 1.23e-13 *** #> a1 17.250 3.622 4.763 0.000156 *** #> a0:b1 9.250 3.622 2.554 0.019936 * #> a1:b1 -0.500 3.622 -0.138 0.891734 #> a0:b2 21.000 3.622 5.798 1.71e-05 *** #> a1:b2 3.250 3.622 0.897 0.381392 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #> #> Residual standard error: 5.122 on 18 degrees of freedom #> Multiple R-squared: 0.7352, Adjusted R-squared: 0.6617 #> F-statistic: 9.997 on 5 and 18 DF, p-value: 0.0001049 #>anova(Exam7.1.lm2)#> Analysis of Variance Table #> #> Response: y #> Df Sum Sq Mean Sq F value Pr(>F) #> a 1 392.04 392.04 14.9428 0.0011330 ** #> a:b 4 919.33 229.83 8.7602 0.0004147 *** #> Residuals 18 472.25 26.24 #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1#> Linear hypothesis test #> #> Hypothesis: #> a0:b1 - a0:b2 = 0 #> #> Model 1: restricted model #> Model 2: y ~ a + a %in% b #> #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 19 748.38 #> 2 18 472.25 1 276.12 10.525 0.004503 ** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1#> Linear hypothesis test #> #> Hypothesis: #> a1:b1 - a1:b2 = 0 #> #> Model 1: restricted model #> Model 2: y ~ a + a %in% b #> #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 19 500.38 #> 2 18 472.25 1 28.125 1.072 0.3142##---Bonferroni's adjusted p-values SimpleEff7.1B <- lsmeans::lsmeans( object = Exam7.1.lm2 , specs = pairwise~b|a , adjust = "bonferroni" )$contrasts SimpleEff7.1B#> a = 0: #> contrast estimate SE df t.ratio p.value #> 0 - 1 -9.25 3.621886 18 -2.554 0.0598 #> 0 - 2 -21.00 3.621886 18 -5.798 0.0001 #> 1 - 2 -11.75 3.621886 18 -3.244 0.0135 #> #> a = 1: #> contrast estimate SE df t.ratio p.value #> 0 - 1 0.50 3.621886 18 0.138 1.0000 #> 0 - 2 -3.25 3.621886 18 -0.897 1.0000 #> 1 - 2 -3.75 3.621886 18 -1.035 0.9426 #> #> P value adjustment: bonferroni method for 3 tests##---Alternative method of pairwise comparisons by applying contrast coefficient with Bonferroni adjustment Bonferroni7.1 <- lsmeans::contrast( Lsm7.1 , list( c1 = c(1,0,-1,0,0,0) , c2 = c(1,0,0,0,-1,0) , c3 = c(0,0,1,0,-1,0) , c4 = c(0,1,0,-1,0,0) , c5 = c(0,1,0,0,0,-1) , c6 = c(0,1,0,0,-1,0) ) , adjust="bonferroni" ) Bonferroni7.1#> contrast estimate SE df t.ratio p.value #> c1 -9.25 3.621886 18 -2.554 0.1196 #> c2 -21.00 3.621886 18 -5.798 0.0001 #> c3 -11.75 3.621886 18 -3.244 0.0270 #> c4 0.50 3.621886 18 0.138 1.0000 #> c5 -3.25 3.621886 18 -0.897 1.0000 #> c6 -3.75 3.621886 18 -1.035 1.0000 #> #> P value adjustment: bonferroni method for 6 tests
