hvac <- read.csv("http://www.stat.ufl.edu/~winner/data/HVAC_perform.csv") attach(hvac); names(hvac) colMeans(hvac[,-1]) plot(hvac[,2:6]) cor(hvac[,2:6]) n <- length(powerp) X0 <- rep(1,n) X <- cbind(X0,airflux,wheelspd,regtemp,humid,drybulb) Y <- powerp (XPXI <- solve(t(X) %*% X)) (XPY <- t(X) %*% Y) (beta_hat <- XPXI %*% XPY) hvac.mod1 <- lm(powerp ~ airflux+wheelspd+regtemp+humid+drybulb) summary(hvac.mod1) anova(hvac.mod1) drop1(hvac.mod1, test="F") hvac.mod2 <- lm(powerp ~ airflux+humid+drybulb) summary(hvac.mod2) anova(hvac.mod2) drop1(hvac.mod2, test="F")