PJ
PSY710
Lab Bennett, #7
Statistics Lab #7
PSY710 November 12, 2008
Initialize R by entering the following command at the prompt. You must type the command exactly as shown. > source(url("http://psycserv.mcmaster.ca/bennett/psy710/Rlab7/Rlab7_init.R") ) > options(contrasts=c("contr.sum","contr.poly"))
1
ANCOVA on heating data
The following passage is copied directly from Rs description of the whiteside data in the MASS package: Mr Derek Whiteside of the UK Building Research Station recorded the weekly gas consumption and average external temperature at his own house in south-east England for two heating seasons, one of 26 weeks before, and one of 30 weeks after cavity-wall insulation was installed. The object of the exercise was to assess the eect of the insulation on gas consumption. I have created a data frame, heating, that contains fake data that matches the design of Whitesides experiment. The variables are Insul, a factor with the levels Before and After; Temp, which contains the weekly average external temperature (C); and Gas, which is the weekly gas consumption in units of 1000 cubic feet. We are going to perform an analysis of covariance (ANCOVA) on these data. The dependent variable is Gas, the independent variable is Insul, and the covariate is Temp.
1.1 1.2
tting linear models to the data graphing the data
The analysis of covariance examines the eect of a factor (Insul) after controlling for the linear association between a covariate (Temp) and dependent variable (Gas). So, before we start, I want to see if a linear model is appropriate for these data. I begin by extracting the data from the Before and After conditions. > heating.before <- subset(heating, Insul == "Before") Next, I use aov to compute the best-tting line that relates Gas and Temp in that condition: > lm.before <- lm(Gas ~ 1 + Temp, data = heating.before) Finally, I graph the data points and best-tting line. The linear t to the data appears adequate: there is no obvious evidence of nonlinearities over the range of the covariate. > with(heating.before, plot(Temp, Gas, ylim = c(0, 10), xlim = c(-5, + 25), pch = 21)) > abline(lm.before) > title(main = "Heating Data (Before)") 1
Bennett, PJ
PSY710
Lab #7
Heating Data (Before)
10
q
8
q q
q
6 Gas
q q q qq
q q q q
4
q
q q q
q q q q
2
q
q
q q q q q q q
0 5
0
5
10 Temp
15
20
25
We can examine the parameters of the linear model with the summary command. summary.lm prints the values of the intercept (5.59) and slope (-0.31) of the tted line, t tests that determine if each parameter diers from zero, and R2 (0.62) indicates the proportion of variance in the dependent variable that is accounted for by the model. > summary(lm.before) Call: lm(formula = Gas ~ 1 + Temp, data = heating.before) Residuals: Min 1Q Median 3Q Max -2.25641 -1.19450 -0.05102 0.76520 3.61773 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.59125 0.40635 13.760 5.53e-14 *** Temp -0.30784 0.04559 -6.752 2.48e-07 *** --Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Residual standard error: 1.441 on 28 degrees of freedom Multiple R-squared: 0.6195, Adjusted R-squared: 0.606 F-statistic: 45.6 on 1 and 28 DF, p-value: 2.479e-07 Task: Extract the data from the After condition, compute the best-tting line, and then graph the data and line. Does the linear model appear to t the data adequately? How do the intercepts and slopes dier between the Before and After conditions?
2
Bennett, PJ
PSY710
Lab #7
1.3
calculating the ANCOVA
Next, we compute the sum of squares associated with Insul after correcting for the linear association between Temp and Gas: > heating.ancova.1 <- lm(Gas ~ 1 + Temp + Insul, data = heating) > anova(heating.ancova.1) Analysis of Variance Table Response: Gas Df Sum Sq Mean Sq F value Pr(>F) Temp 1 148.691 148.691 96.223 7.578e-14 *** Insul 1 18.218 18.218 11.790 0.001116 ** Residuals 57 88.081 1.545 --Signif. codes: 0 *** 0.001 0.01 ** * 0.05 . 0.1 1 The eect of insulation after controlling for the linear eect of temperature is signicant (F (1, 57) = 11.779, p = .001).
1.4
homogeneity of slopes assumption
The previous section used the analysis of covariance to show that the main eect of insulation was signicant after controlling for the linear eect of temperature. However, the linear model used in that analysis did not include an interaction term: it assumed that the interaction was zero. This assumption is equivalent to assuming that the slope of the line relating Gas to temperature had the same slope in the two conditions. We can test the homogeneity of slopes assumption by testing the signicance of the interaction term: > heating.ancova.2 <- lm(Gas ~ 1 + Temp + Insul + Temp:Insul, data = heating) > anova(heating.ancova.2) Analysis of Variance Table Response: Gas Df Sum Sq Mean Sq Temp 1 148.691 148.691 Insul 1 18.218 18.218 Temp:Insul 1 0.003 0.003 Residuals 56 88.078 1.573 --Signif. codes: 0 *** 0.001
F value Pr(>F) 94.5377 1.253e-13 *** 11.5833 0.001236 ** 0.0018 0.966513
** 0.01 * 0.05 . 0.1 1
The interaction is not signicant (F (1, 56) = 0.0018, p = 0.96), so the homogeneity of slopes assumption appears to be valid for these data. Task: Perform the ANCOVA using the aov command.
1.5
inspecting the residuals
ANCOVA assumes that the errors are distributed normally, and that error variance is constant across the range of the covariate. We can test the constant variance assumption by plotting the residuals against the temperature. If the constant variance assumption is true, then the residuals should be scattered about the zero line and the 3
Bennett, PJ
PSY710
Lab #7
spread of the scatter should be constant across the covariate. The graph indicates that the magnitude of the residuals is roughly constant across the range of the covariate, so the constant variance assumption appears reasonable. > plot(heating$Temp, residuals(heating.ancova.1)) > abline(h = 0) > title(main = "ancova: residuals vs. covariate")
ancova: residuals vs. covariate
q
3
q
residuals(heating.ancova.1)
2
q q
q q
q
1
q
qq q q q q qq q q q
q q q qq q q q q
q
q q q q q q qq qq q q q q q
0
q
q
q
1
q q q q q q
q q q q q
q
2
q
5
0
5
10 heating$Temp
15
20
We can check the normality assumption by using a qq plot to inspect the residuals > qqnorm(residuals(heating.ancova...