> ratspeed = read.table(file="D:\\stat141\\ratspeed.dat", header = T) # fit regression line to amphetamine data ex 12.1 > plot(dose, food, pch = 20) > tapply( food,dose, mean) 0 2.5 5 99.9875 75.5000 54.9500 > plot(c(0,2.5,5),tapply( food,dose, mean), pch = 15) #fig 12.10 plot of ave(Y|X) for X = 0, 2.5,5 > ratreg = lm(food ~ dose) # object for the food on dose straight-line fit > summary(ratreg) Call:lm(formula = food ~ dose) Residuals: Min 1Q Median 3Q Max -21.513 -7.031 1.528 7.448 27.006 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 99.331 3.680 26.99 < 2e-16 *** dose -9.008 1.140 -7.90 7.27e-08 *** --- Residual standard error: 11.4 on 22 degrees of freedom Multiple R-Squared: 0.7394, Adjusted R-squared: 0.7275 F-statistic: 62.41 on 1 and 22 DF, p-value: 7.265e-08 > plot(dose, food, pch = 20) > abline(ratreg) > anova(ratreg) Analysis of Variance Table Response: food Df Sum Sq Mean Sq F value Pr(>F) dose 1 8113.5 8113.5 62.415 7.265e-08 *** Residuals 22 2859.9 130.0 > sqrt(130.0) [1] 11.40175 > coef(ratreg) (Intercept) dose 99.33125 -9.00750 #confidence interval for slope > coef(ratreg)[2] - qt(.975,22)*1.140; coef(ratreg)[2] + qt(.975,22)*1.140 dose -11.37172 dose -6.643285 > qqnorm(resid(ratreg)) # not perfect but not awful normality graphic > plot(fitted(ratreg), residuals(ratreg), pch = 20) #resid vs fit > predict.lm(ratreg, se.fit=T, interval = "confidence", level=.95) $fit fit lwr upr 1 99.33125 91.69979 106.96271 9 76.81250 71.98594 81.63906 17 54.29375 46.66229 61.92521 $se.fit [1] 3.679811 3.679811 3.679811 3.679811 3.679811 3.679811 [7] 3.679811 3.679811 2.327317 2.327317 2.327317 2.327317 [13] 2.327317 2.327317 2.327317 2.327317 3.679811 3.679811 [19] 3.679811 3.679811 3.679811 3.679811 3.679811 3.679811 > predict.lm(ratreg, se.fit=T, interval = "prediction", level=.95) $fit fit lwr upr 1 99.33125 74.48502 124.17748 9 76.81250 52.67970 100.94530 17 54.29375 29.44752 79.13998 Revisit Regression ex # fit regression line to amphetamine data ex 12.1 > ratspeed = read.table(file="D:\\stat141\\ratspeed.dat", header = T) > ratreg = lm(food ~ dose) # object for the food on dose straight-line fit > summary(ratreg) Call: lm(formula = food ~ dose) Residuals: Min 1Q Median 3Q Max -21.513 -7.031 1.528 7.448 27.006 Coefficients: Estimate Std. Error t value Pr(>|t|) # fit is 99.3 - 9.0*dose (Intercept) 99.331 3.680 26.99 < 2e-16 *** # slope is change in food dose -9.008 1.140 -7.90 7.27e-08 *** # consumption per unit --- increase in dose Residual standard error: 11.4 on 22 degrees of freedom Multiple R-Squared: 0.7394, Adjusted R-squared: 0.7275 F-statistic: 62.41 on 1 and 22 DF, p-value: 7.265e-08 R-square equivalents > cor(food,dose) > cor(food,dose)^2 #squared mult corr Rsq [1] -0.859873 [1] 0.7393816 > cor(food,fitted(ratreg))^2 #squared mult corr Rsq [1] 0.7393816 > anova(ratreg) # show overhead, decomposition STT = SSReg + SSE Analysis of Variance Table Response: food Df Sum Sq Mean Sq F value Pr(>F) dose 1 8113.5 8113.5 62.415 7.265e-08 *** Residuals 22 2859.9 130.0 > rsq = 8113.5/(8113.5 + 2859.9) [1] 0.7393789 > adjrsq = 1 - (23/22)*2859.8/(8113.5 + 2859.9) [1] 0.727542 > ssqydotx = (23/22)*var(food)*(1 - cor(food,dose)^2) > ssqydotx [1] 129.9937 > sqrt(ssqydotx) [1] 11.40148 > var(residuals(ratreg))*(23/22) [1] 129.9937 > sum(residuals(ratreg)) > quantile(residuals(ratreg)) [1] 5.32907e-15 0% 25% 50% 75% 100% -21.512500 -7.031250 1.528125 7.448437 27.006250 > predict.lm(ratreg, se.fit=T, interval = "confidence", level=.95) fit lwr upr # carryover from 11/29 1 99.33125 91.69979 106.96271 9 76.81250 71.98594 81.63906 17 54.29375 46.66229 61.92521 $se.fit [1] 3.679811 [9] 2.327317 [17] 3.679811 > predict.lm(ratreg, se.fit=T, interval = "prediction", level=.95) fit lwr upr 1 99.33125 74.48502 124.17748 9 76.81250 52.67970 100.94530 17 54.29375 29.44752 79.13998 > predict.lm(ratreg, newdata = data.frame(dose=4),se.fit=T, interval = "confidence", level=.95) fit lwr upr # CI for for at interpolated value dose = 4 [1,] 63.30125 57.31165 69.29085 $se.fit [1] 2.888124 > predict.lm(ratreg, newdata = data.frame(dose=4),se.fit=T, interval = "prediction", level=.95) fit lwr upr [1,] 63.30125 38.90921 87.69329