# ------------------------------------------------------------------------------- # calculating adjusted means and testing for pairwise differences in ANCOVA # with unequal slopes in R using the glht function in the multcomp package # ------------------------------------------------------------------------------- # for the pie data, we have: > anova(dessert.lm1) Analysis of Variance Table Response: taste Df Sum Sq Mean Sq F value Pr(>F) hours 1 12.5939 12.5939 738.798 < 2.2e-16 *** pie 2 2.8430 1.4215 83.389 5.587e-13 *** hours:pie 2 0.8372 0.4186 24.557 4.820e-07 *** Residuals 30 0.5114 0.0170 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # indicating a lack of parallelism # a common approach in this situation is to compare the grouping factor (pies) # separately at a few different x values, similar to using a slice approach # for understanding interaction in factorial models # ------------------------------------------------------------------------------- # In the following approach, we have to understand the ANCOVA model that R is using: # ------------------------------------------------------------------------------- > summary(dessert.lm1) Call: lm(formula = taste ~ hours + pie + hours:pie, data = dessert) Residuals: Min 1Q Median 3Q Max -0.22114 -0.08820 0.01494 0.09436 0.21474 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.50397 0.17191 20.383 < 2e-16 *** hours 0.47164 0.05201 9.069 4.24e-10 *** piemincemeat -1.11895 0.26608 -4.205 0.000217 *** piepumpkin 0.01687 0.28999 0.058 0.953994 hours:piemincemeat 0.42528 0.07157 5.942 1.64e-06 *** hours:piepumpkin 0.00846 0.06934 0.122 0.903707 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1306 on 30 degrees of freedom Multiple R-squared: 0.9695, Adjusted R-squared: 0.9645 F-statistic: 190.9 on 5 and 30 DF, p-value: < 2.2e-16 # The fitted model is (predicted taste, rounding to three digits): yhat = 3.50 + .472*hours -1.12*I(pie = mincemeat) +.0169*I(pie = pumpkin) +.425*hours*I(pie = mincemeat) +.00846*hours*I(pie = pumpkin), so, for example, if a pie is in the custard group, yhat = 3.50 + .472*hours, but if a pie is in the pumpkin group, yhat = 3.50 + .472*hours +.0169 +.00846*hours = 3.52 + .480*hours # ------------------------------------------------------------------------------- # We can calculate adjusted means and perform pairwise multiple comparison tests # in ANCOVA with unequal slopes using the glht function in the multcomp package, # but we have to carefully specify the comparisons of interest using contrasts # ------------------------------------------------------------------------------- # To estimate the mean taste for custard pies at x = 3.2 (hours), we use the # following contrast and the glht function: c.cust.at.3.2 <- rbind("custard at 3.2" <- c(1,3.2,0,0,0,0)) summary(glht(dessert.lm1, linfct= c.cust.at.3.2)) # We could estimate all three pie means at 3.2 hours with: contr3.2 <- rbind("cust at 3.2" <- c(1,3.2,0,0,0,0), "mince at 3.2" <- c(1,3.2,1,0,3.2,0), "pump at 3.2" <- c(1,3.2,0,1,0,3.2)) summary(glht(dessert.lm1, linfct= contr3.2)) # Notice how all three contrasts start with a '1' and '3.2', because all three # lines use the first two terms above: yhat = 3.50 + .472*hours + .... (4 other terms) # and 3.2 is the requested value for hours # Note that the third and fourth coefficients are only nonzero if referring to the # corresponding group (mince at 3.2 has a '1' for the third term, which is 'piemincemeat') # Similarly, the fifth and sixth coefficients have the requested hour value if referring # to the corresponding group # Now that we have the coefficients needed to estimate means, to estimate differences of means, # we take differences of the corresponding coefficients: "cust at 3.2" <- c(1,3.2,0,0, 0, 0), "mince at 3.2" <- c(1,3.2,1,0,3.2,0), so, cust - mince at 3.2 is c(1,3.2,0,0, 0, 0) - c(1,3.2,1,0,3.2,0) = c(0, 0,-1,0,-3.2,0) # We now calculate these differences for each pairwise comparison and again use the # glht function to get: > contrdiff3.2 <- rbind("cust -mince at 3.2" <- c(0,0,-1,0,-3.2,0), "mince - pump at 3.2" <- c(0,0,1,-1,3.2,-3.2), "cust - pump at 3.2" <- c(0,0,0,-1,0,-3.2)) > > summary(glht(dessert.lm1, linfct= contrdiff3.2)) Simultaneous Tests for General Linear Hypotheses Fit: lm(formula = taste ~ hours + pie + hours:pie, data = dessert) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) 1 == 0 -0.24194 0.06800 -3.558 0.00368 ** 2 == 0 0.19800 0.10784 1.836 0.17004 3 == 0 -0.04394 0.09925 -0.443 0.89540 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method) # These three P values are nearly the same as those shown by SAS. The glht function # uses a different multiple comparison adjustment than Tukey, so they will not have # identical P values # ------------------------------------------------------------------------------- # ------------------------------------------------------------------------------- # For the other two time points, use this for x = 4.2: contrdiff4.2 <- rbind("cust -mince at 4.2" <- c(0,0,-1,0,-4.2,0), "mince - pump at 4.2" <- c(0,0,1,-1,4.2,-4.2), "cust - pump at 4.2" <- c(0,0,0,-1,0,-4.2)) summary(glht(dessert.lm1, linfct= contrdiff4.2)) # And this for x = 4.8: contrdiff4.8 <- rbind("cust -mince at 4.8" <- c(0,0,-1,0,-4.8,0), "mince - pump at 4.8" <- c(0,0,1,-1,4.8,-4.8), "cust - pump at 4.8" <- c(0,0,0,-1,0,-4.8)) summary(glht(dessert.lm1, linfct= contrdiff4.8))