################################################################################
## Chapter 12
## April 30, 2009
## Alan T. Arnholt
options(width=65)
library(PASWR)
################################################################################
# Example 12.6
# part a.
attach(Grades)
plot(sat,gpa)
# part bi.
Y <- gpa
x <- sat
# Solving using summation notation as in (12.18) (b1) and (12.14) (b0):
b1 <- sum((x-mean(x))*(Y-mean(Y))) / sum((x-mean(x))^2)
b0 <- mean(Y)-b1*mean(x)
c(b0,b1)
# part bii. Solving using matrix notation as in (12.23):
X <- cbind(rep(1,200),x)
Y <- matrix(Y,nrow=200)
betahat <- solve(t(X)%*%X)%*%t(X)%*%Y
beta0hat <- betahat[1,1]
beta1hat <- betahat[2,1]
c(beta0hat,beta1hat)
# part biii. Solving use lm()
model <- lm(Y~x)
model$coefficients
# part c.
b1*50
detach(Grades)
################################################################################
# Example 12.7
attach(Grades)
Y <- gpa
x <- sat
simple.model <- lm(Y~x)
X <- cbind(rep(1,200),x)
XTX <- t(X)%*%X
solve(XTX)
# The preferred way to compute (X'X)^-1 with S is now shown:
XTXI <- summary(simple.model)$cov.unscaled
XTXI
MSE <- sum(resid(simple.model)^2)/(200-2)
MSE
# A more direct method of obtaining the MSE is
summary(simple.model)$sigma^2
#
var.cov.b <- MSE*XTXI
var.cov.b
# Computing above directly (only R)
vcov(simple.model)
# part b.
summary(simple.model)$coef
# part c.
b0 <- summary(simple.model)$coef[1,1]
b1 <- summary(simple.model)$coef[2,1]
s.b0 <- summary(simple.model)$coef[1,2]
s.b1 <- summary(simple.model)$coef[2,2]
ct <- qt(1-.10/2,198) # alpha = 0.10
ct
CI.B0 <- c(b0 - ct*s.b0, b0 + ct*s.b0)
CI.B0
CI.B1 <- c(b1 - ct*s.b1, b1 + ct*s.b1)
CI.B1
# Or, if working in R only, a method requiring less typing is:
confint(lm(Y~x),level=.90)
detach(Grades)
################################################################################
# Example 12.8
attach(Grades)
model.lm <- lm(gpa~sat)
anova(model.lm)
detach(Grades)
################################################################################
# Example 12.10
attach(HSwrestler)
Y <- HWFAT
x1 <- ABS
x2 <- TRICEPS
x3 <- SUBSCAP
n <- 78
X <- cbind(rep(1,n),x1,x2,x3)
Y <- matrix(HWFAT,nrow=n)
J <- matrix(rep(1,n*n),nrow=n)
H <- X%*%solve(t(X)%*%X)%*%t(X)
b <- solve(t(X)%*%X)%*%t(X)%*%Y
SSR <- t(b)%*%t(X)%*%Y-(1/n)*t(Y)%*%J%*%Y
SSE <- t(Y)%*%Y - t(b)%*%t(X)%*%Y
SSTO <- t(Y)%*%Y -(1/n)*t(Y)%*%J%*%Y
SS <- rbind(SSR,SSE,SSTO)
r.names <- c("SSR","SSE","SSTO")
c.names <- c("Sum of Squares")
dimnames(SS) <- list(r.names,c.names)
SS
# Computing SSR, SSE, and SST with (12.42), (12.40), and (12.41)
ID <- diag(nrow=n) # n*n Identity matrix
SSRa <- t(Y)%*%(H-(1/n)*J)%*%Y
SSEa <- t(Y)%*%(ID-H)%*%Y
SSTOa <- t(Y)%*%(ID-(1/n)*J)%*%Y
SSa <- rbind(SSRa,SSEa,SSTOa)
r.names <- c("SSR","SSE","SSTO")
c.names <- c("Sum of Squares")
dimnames(SSa) <- list(r.names,c.names)
SSa
#
mod123 <- lm(Y~x1+x2+x3)
anova(mod123)
#
mod321 <- lm(Y~x3+x2+x1)
anova(mod321)
detach(HSwrestler)
################################################################################
# Example 12.11
attach(HSwrestler)
Y <- HWFAT
x1 <- ABS
x2 <- TRICEPS
x3 <- SUBSCAP
mod132 <- lm(Y~x1+x3+x2)
anova(mod132)
anova(mod132)[3,4]                  # Fobs value
summary(mod132)
summary(mod132)$coefficients[4,3]^2 # tobs value squared
detach(HSwrestler)
#
drop1(mod132, test="F")
################################################################################
# Example 12.12
attach(HSwrestler)
Y <- HWFAT
x1 <- ABS
x2 <- TRICEPS
x3 <- SUBSCAP
mod312 <- lm(Y~x3+x1+x2)
anova(mod312)
#
pf(19.98, 2, 74,lower.tail=FALSE)  # Only R
# Another approach to test H0 : Beta1 = Beta2 = 0 is to use anova()
# on the reduced and full models
mod3 <- lm(Y~x3)
anova(mod3,mod312)
detach(HSwrestler)
################################################################################
# General Linear Hypothesis
# Example 12.13 - General Linear Model
attach(HSwrestler)
Y <- HWFAT
x1 <- ABS
x2 <- TRICEPS
x3 <- SUBSCAP
mod312 <- lm(Y~x3+x1+x2)
K <- matrix(c(0,0,1,0,0,0,0,1),byrow=TRUE,nrow=2)
Q <- qr(K)$rank
b <- matrix(coef(mod312),ncol=1)
m <- matrix(c(0,0),byrow=TRUE,nrow=2)
XTXI <- summary(mod312)$cov.unscaled
NUM <- t(K%*%b-m)%*%solve(K%*%XTXI%*%t(K))%*%(K%*%b-m)
MSE <- anova(mod312)[4,3]
Fobs <- NUM/(Q*MSE)
c(Fobs,1-pf(Fobs,Q,74))
# Using glht() from the multcomp package.
library(multcomp)
summary(glht(mod312, linfct=K, rhs=c(0,0)), test=Ftest())
# part b.
K <- matrix(c(0,2,1,-1,0,-5,0,1),byrow=TRUE,nrow=2)
summary(glht(mod312, linfct=K, rhs=c(0,0.2)), test=Ftest())
# part c.
K <- matrix(c(0,0,1,-1),byrow=TRUE,nrow=1)
summary(glht(mod312, linfct=K, rhs=0), test=Ftest())
detach(HSwrestler)
################################################################################
# Example 12.14 - Model Selection with HSwrestler
attach(HSwrestler)
# Backward elimination showing all steps
reg.all <- lm(HWFAT ~ AGE + HT + WT + ABS + TRICEPS + SUBSCAP)
summary(reg.all)$coefficients
# Note that HT has the largest p-value of 6.762255e-01, so it is eliminated
# from the model.
reg.m1 <- lm(HWFAT ~ AGE + WT + ABS + TRICEPS + SUBSCAP)
summary(reg.m1)$coefficients
# Note that SUBSCAP has the largest p-value of 4.306601e-01, so it is
# eliminated from the model.
reg.m2 <- lm(HWFAT ~ AGE + WT + ABS + TRICEPS)
summary(reg.m2)$coefficients
# Note that WT has the largest p-value of 3.869075e-01, so it is eliminated
# from the model.
reg.m3 <- lm(HWFAT ~ AGE + ABS + TRICEPS)
summary(reg.m3)
# Backward elimination using drop1() in R
drop1(lm(HWFAT ~ AGE + HT +WT + ABS +TRICEPS +SUBSCAP),test="F")
# Note that HT has the largest p-value of 0.676225, so it is eliminated
# from the model.
drop1(lm(HWFAT ~ AGE + WT + ABS +TRICEPS +SUBSCAP),test="F")
# Note that SUBSCAP has the largest p-value of 0.430660, so it is eliminated
# from the model.
drop1(lm(HWFAT ~ AGE + WT + ABS +TRICEPS),test="F")
# Note that WT has the largest p-value of 0.3869, so it is eliminated from
# the model.
drop1(lm(HWFAT ~ AGE + ABS +TRICEPS),test="F")
# part b. --- Forward Selection
# Forward selection showing all steps using add1() in R
add1(lm(HWFAT~1), scope=(~.+ AGE + HT + WT + ABS + TRICEPS + SUBSCAP),
test="F")
# The variable ABS has the most significant (smallest) p-value = 2.2e - 16
# with the largest F value= 407.9929, so it is added to the model.
add1(lm(HWFAT~ABS),scope=(~.+AGE +HT +WT +TRICEPS +SUBSCAP),test="F")
# The variable TRICEPS has the most significant (smallest)
# p-value = 2.639e - 06 with the largest F value= 25.8443, so it is
# added to the model.
add1(lm(HWFAT~ABS+TRICEPS),scope=(~.+ AGE + HT + WT + SUBSCAP),test="F")
# The variable AGE has the most significant (smallest) p-value = 0.04441
# with the largest F value= 4.1823, so it is added to the model.
add1(lm(HWFAT~ABS+TRICEPS+AGE), scope=(~.+ HT + WT + SUBSCAP), test="F")
########################
summary(lm(HWFAT~ABS+TRICEPS+AGE+HT))$coefficients
#
summary(lm(HWFAT~ABS+TRICEPS+AGE+WT))$coefficients
#
summary(lm(HWFAT~ABS+TRICEPS+AGE+SUBSCAP))$coefficients
# part c.
library(leaps)
a <- regsubsets(as.matrix(HSwrestler[,-c(7,8,9)]),HSwrestler[,7])
summary(a)
#
summary(a)$adjr2
# part d.
summary(a)$cp
par(pty="s")
plot(2:7,summary(a)$cp,ylim=c(2,7),xlab="p",ylab="Cp")
abline(a=0,b=1)
par(pty="m")
# part e.
library(MASS) # For function stepAIC()
reg.all <- lm(HWFAT ~ AGE + HT + WT + ABS + TRICEPS + SUBSCAP)
mod.aic <- stepAIC(reg.all, direction="both",
scope=(~.+SUBSCAP+TRICEPS+ABS+WT+HT+AGE),k=2)
# part f.
mod.bic <- stepAIC(reg.all, direction="both",
scope=(~.+SUBSCAP+TRICEPS+ABS+WT+HT+AGE),k=log(length(HWFAT)))
mod.bic
detach(HSwrestler)
################################################################################
#  Example 12.15
attach(HSwrestler)
mod.3 <- lm(HWFAT~AGE+ABS+TRICEPS)
# Getting the values manually...
X <- model.matrix(mod.3)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
hii <- diag(H)
hii[1:5]
# Extracts hatvalues in R and S-PLUS
influence(mod.3)$hat[1:5]
# Verifying that sum(h_ii)=p
sum(hii)
detach(HSwrestler)
################################################################################
# Example 12.16
attach(HSwrestler)
library(MASS)
mod.2 <- lm(HWFAT~ABS+TRICEPS)
# Figure 12.9
par(mfrow=c(2,2))
plot(fitted(mod.2),resid(mod.2),ylim=c(-10,10),main="")
title(main="Residuals vs Fitted")
abline(h=0,lty=2)
plot(fitted(mod.2),stdres(mod.2),ylim=c(-3.5,3.5),main="")
title(main="Standardized Residuals vs Fitted")
abline(h=0,lty=2)
plot(fitted(mod.2),studres(mod.2),ylim=c(-3.5,3.5),main="")
title(main="Studentized Residuals vs Fitted")
abline(h=0,lty=2)
plot(mod.2,which=1,main="Default Graph 1")
par(mfrow=c(1,1))
sort(abs(resid(mod.2)))[76:78]    # Extract three largest values
sort(abs(stdres(mod.2)))[76:78]   # Extract three largest values
sort(abs(studres(mod.2)))[76:78]  # Extract three largest values
qt(1-.2/(2*78),78-3-1)            # Critical value
detach(HSwrestler)
################################################################################
# Example 12.17
attach(Kinder)
# part a. Figure 12.10
plot(wt,ht)
# part b.
mod <- lm(ht~wt)
hii <- lm.influence(mod)$hat
hii
#
X <- model.matrix(mod)
n <- nrow(X)
p <- ncol(X)
H <- X%*%solve(t(X)%*%X)%*%t(X)
hi <- diag(H)
plot(hi,type="h",ylab="leverage")
abline(h=2*p/n)
#
modk <- lm(ht[-c(19,20)]~wt[-c(19,20)])
# part c.
library(MASS)     # Need for function studres()
library(car)      # Need for function cookd
modk19 <- lm(ht[-19]~wt[-19])
n <- 19
p <- 2
par(mfrow=c(2,3))
hiik19 <- lm.influence(modk19)$hat # extracting hii values
plot(hiik19,ylab="Leverage")
cv <- 2*p/n
abline(h=cv,lty=2)
plot(lm.influence(modk19)$coefficients[,2],
ylab="Difference in Coefficients")
plot(dfbetas(modk19)[,2],ylab="DFBETAS")
cv <- 2/sqrt(n)               # Critical value for DFBETAS
abline(h=c(-cv,cv),lty=2)
plot(studres(modk19),ylab="Studentized Residuals")
cv <- qt(1-.10/(2*n), n-p-1)    # Critical value
abline(h=c(-cv,cv),lty=2)
DFFITS <- studres(modk19)*(hiik19/(1-hiik19))^.5 #See *
plot(DFFITS,ylab="DFFITS")
cv <- 2*sqrt(p/n)             # Critical value for DFITS
abline(h=c(-cv,cv),lty=2)
cd <- cookd(modk19)           # Cook’s distance
plot(cd,ylab="Cook’s Distance")
CF <- qf(.50,p,n-p)           # Critical value for Cook’s Distance
abline(h=CF,lty=2)
par(mfrow=c(1,1))
# part d.
# Figure 12.13
modk20 <- lm(ht[-20]~wt[-20])
n <- 19
p <- 2
par(mfrow=c(2,3))
hiik20 <- lm.influence(modk20)$hat
plot(hiik20,ylab="Leverage")
cv <- 2*p/n
abline(h=cv,lty=2)
plot(lm.influence(modk20)$coefficients[,2],ylab="Difference in Coefficients")
plot(dfbetas(modk20)[,2],ylab="DFBETAS")
cv <- 2/sqrt(n)             # Critical value for DFBETAS
abline(h=c(-cv,cv),lty=2)
plot(studres(modk20),ylab="Studentized Residuals")
cv <- qt(1-.10/(2*n), n-p-1)  # Critical value
abline(h=c(-cv,cv),lty=2)
DFFITS <- studres(modk20)*(hiik20/(1-hiik20))^.5
plot(DFFITS,ylab="DFFITS")
cv <- 2*sqrt(p/n)           # Critical value for DFITS
abline(h=c(-cv,cv),lty=2)
cd <- cookd(modk20)         # Cook’s distance
plot(cd,ylab="Cook’s Distance")
CF <- qf(.50,p,n-p)         # Critical value for Cook’s Distance
abline(h=CF,lty=2)
par(mfrow=c(1,1))
# part e.
# Figure 12.14
# Scatterplot of height versus weight for data from Kinder with four
# superimposed regression lines: The solid line is for model mod
# (all observations); the dot dash line is for model modk20 (case 20 omitted);
# the dashed line is for model modk (case 19 and 20 omitted); and the dotted
# line is for model modk19 (case 19 omitted).
plot(wt[-c(19,20)],ht[-c(19,20)],cex=2,xlim=c(28,62),ylim=c(36,52),
xlab="Weight in Pounds",ylab="Height in Inches")
abline(mod,lty=1,lwd=2)
abline(modk,lty=2,lwd=2)
abline(modk19,col="red",lty=3,lwd=2)
abline(modk20,col="blue",lty=4,lwd=2)
abline(mod,lty=4,lwd=2)
points(wt[19],ht[19],pch=16,cex=2,col="red")
points(wt[20],ht[20],pch=17,cex=2,col="blue")
detach(Kinder)
################################################################################
rm(list=ls(all=TRUE))
# Example 12.18
# Transform x_1
# Figure 12.15
attach(SimDataXT)
par(mfrow=c(2,3))
plot(x1,Y)
lines(x1,x1^.5)               # function Y = x1^.5
plot(lm(Y~x1),which=c(1,2))   # Residual and Q-Q normal plots
plot(x1^.5,Y)
mod1 <- lm(Y~I(x1^.5))        # Works in R for S-PLUS see *
abline(mod1)
plot(mod1,which=c(1,2))       # Q-Q plot in S-PLUS is 4 not 2
par(mfrow=c(1,1))
# Transform x_2
# Figure 12.16
par(mfrow=c(2,3))
plot(x2,Y)
lines(x2,x2^2)                # function Y = x2^2
plot(lm(Y~x2),which=c(1,2))   # Residuals and Q-Q normal plots
plot(x2^2,Y)
mod2 <- lm(Y~I(x2^2))
abline(mod2)
plot(mod2,which=c(1,2))
par(mfrow=c(1,1))
# Transform x_3
# Figure 12.17
par(mfrow=c(2,3))
plot(x3,Y)
lines(x3,x3^(-1))             # function Y = 1/x3
plot(lm(Y~x3),which=c(1,2))   # Residuals and Q-Q normal plots
plot(x3^(-1),Y)
mod3 <- lm(Y~I(x3^(-1)))
abline(mod3)
plot(mod3,which=c(1,2))
par(mfrow=c(1,1))
detach(SimDataXT)
################################################################################
# Example 12.19
attach(SimDataXT)
modC <- lm(Y~I(x1^.5)+I(x2^2)+I(x3^(-1)))
summary(modC)
cbind(x2,x3^(-.5))[1:5,]
# part b.
modB <- lm(Y~I(x1^.5)+I(x2^2))
X <- model.matrix(modB)
eigen(t(X)%*%X, only.values=TRUE)$values # extracting eigenvalues
lambda.max <- max(eigen(t(X)%*%X, only.values=TRUE)$values)
lambda.min <- min(eigen(t(X)%*%X, only.values=TRUE)$values)
condition.number <- (lambda.max/lambda.min)^.5
condition.number
#
kappa(X,exact=TRUE)
# Computing the VIF with (12.69):
1/(1-summary(lm(I(x1^.5)~I(x2^2)))$r.square)
# Computing the variance inflation factors with the function vif()
# from the car package.
library(car) # For function vif()
vif(modB)
# Verifying that the standard error for betahat_1 from a model where Y is
# regressed solely on x_1^0.5 to modB increases by approximately sqrt(VIF_1).
mod1 <- lm(Y~I(x1^.5))
summary(mod1)$coefficients
se.x1.mod1 <- summary(mod1)$coefficients[2,2]
se.x1.modB <- summary(modB)$coefficients[2,2]
ratio <- se.x1.modB/se.x1.mod1
ratio
vif(modB)^.5
detach(SimDataXT)
################################################################################
# Example 12.20
# Box-Cox Transformations
# Figure 12.18
attach(SimDataST)
library(MASS)
par(mfrow=c(1,2)) # 1 row by 2 columns
modx1 <- lm(Y1~x1)
boxcox(modx1)
boxcox(modx1,lambda=seq(-.3,.3,.01))
par(mfrow=c(1,1))
detach(SimDataST)
################################################################################
# Example 12.21
# Transforming Y with ln
# Figure 12.19
attach(SimDataST)
library(MASS)
par(mfrow=c(2,3))
plot(x1,Y1)
modx1 <- lm(Y1~x1)
plot(modx1,which=1)
boxcox(modx1,lambda=seq(-.3,.3,.01))
plot(x1,log(Y1))
plot(lm(log(Y1)~x1),which=c(1,2)) #Q-Q plot in S-PLUS is 4 not 2
par(mfrow=c(1,1))
detach(SimDataST)
################################################################################
# Example 12.22
# Transforming with a reciprocal
# Figure 12.20
attach(SimDataST)
par(mfrow=c(2,3))
plot(x2,Y2)
modx2 <- lm(Y2~x2)
plot(modx2,which=1)
boxcox(modx2,lambda=seq(-1.4,-.6,.01))
plot(x2,Y2^(-1))
plot(lm(Y2^(-1)~x2),which=c(1,2)) #Q-Q plot in S-PLUS is 4 not 2
par(mfrow=c(1,1))
# Figure 12.21
library(MASS)
par(mfrow=c(2,3))
plot(x3,Y1)
plot(x3,log(Y1))
mod <- lm(log(Y1)~x3)
abline(mod)
plot(x3^.5,log(Y1))
mod2 <- lm(log(Y1)~I(x3^.5))
abline(mod2)
boxcox(lm(Y1~x3),lambda=seq(-.3,.3,.01))
plot(mod,which=1)
plot(mod2,which=1)
par(mfrow=c(1,1))
detach(SimDataST)
################################################################################
# Interpreting a Logarithmically Transformed Model
library(MASS)
attach(Animals)
SA <- Animals[order(body),]
NoDINO <- SA[-c(28:26),]
detach(Animals)
attach(NoDINO)
Y <- log(brain)
x <- log(body)
simple.model <- lm(Y~x)
(simple.model)$coef
detach(NoDINO)
################################################################################
# Qualitative Predictor
contr.treatment(4)
################################################################################
# Example 12.23
attach(EPIDURAL)
contrasts(Ease)
levels(Ease) <- c("Easy", "Difficult", "Impossible")
levels(Ease)
contrasts(Ease)
detach(EPIDURAL)
################################################################################
# Example 12.24 - Elevators
# part a.
attach(vit2005)
Elevator <- as.factor(elevator)
contrasts(Elevator)
modSimpl <- lm(totalprice~area)
summary(modSimpl)
# part bi.
modTotal <- lm(totalprice~area+Elevator+area:Elevator)
modSimpl <- lm(totalprice~area)
anova(modSimpl,modTotal)
# part bii.
anova(modTotal)
# part biii.
modTotal <- lm(totalprice~area+Elevator+area:Elevator)
modInter <- lm(totalprice~area+area:Elevator)
anova(modInter,modTotal)
# Since the p-value for testing the null hypothesis is 0.1133, one fails to
# reject H0 and should conclude that the two lines have the same intercept
# but different slopes.
summary(modInter)$coef
detach(vit2005)
################################################################################
# Example 12.25 - Grades
# part a.
attach(Grades)
Y <- gpa
x <- sat
simple.model <- lm(Y~x)
betahat <- simple.model$coef
betahat
Xh <- matrix(c(1,1300),nrow=1)
Yhath <- Xh%*%betahat
Yhath
# A method which requires less typing is:
predict(simple.model,newdata=data.frame(x=1300))
# part b.
MSE <- anova(simple.model)[2,3]
MSE
XTXI <- summary(simple.model)$cov.unscaled
XTXI
var.cov.b <- MSE*XTXI
var.cov.b
s2yhath <- Xh%*%var.cov.b%*%t(Xh)
s2yhath
syhath <- sqrt(s2yhath)
syhath
crit.t <- qt(.95,198)
crit.t
CI.EYh <- c(Yhath - crit.t*syhath, Yhath + crit.t*syhath)
CI.EYh
# Or in R
predict(lm(Y~x),data.frame(x=1300),interval="conf",level=.90)
# For S-SPLUS
# predict(lm(Y~x),data.frame(x=1300),ci.fit=TRUE,conf.level=.90)$ci.fit
# part c.
s2yhathnew <- MSE + s2yhath
syhathnew <- sqrt(s2yhathnew)
syhathnew
PI <- c(Yhath - crit.t*syhathnew, Yhath + crit.t*syhathnew)
PI
# Or only in R
predict(lm(Y~x),data.frame(x=1300),interval="pred",level=.90)
detach(Grades)
#
# Only in S-PLUS
# predict(lm(Y~x),data.frame(x=1300),pi.fit=TRUE,conf.level=.90)$pi.fit
################################################################################
# Example 12.26
# part a.
attach(HSwrestler)
alpha <- 0.10
mult.model <- lm(HWFAT~AGE+ABS+TRICEPS)
summary(mult.model)$coef
b <- summary(mult.model)$coef[2:4,1]
s.b <- summary(mult.model)$coef[2:4,2]
g <- 3
B <- qt((1-alpha/(2*g)),78-4)
B
BonSimCI.b <- matrix(c(b-B*s.b, b+B*s.b),ncol=2)
conf <- c("5%","95%")
bnam <- c("AGE","ABS","TRICEPS")
dimnames(BonSimCI.b) <- list(bnam,conf)
BonSimCI.b
Q <- 3
S <- sqrt(Q*qf(.9,Q,78-4))
S
SchSimCI.b <- matrix(c(b-S*s.b, b+S*s.b),ncol=2)
dimnames(SchSimCI.b) <- list(bnam,conf)
SchSimCI.b
# part b. Bonferroni
confidence.ellipse(lm(HWFAT~AGE+ABS+TRICEPS),level=.90,
which.coef=c(3,4),Scheffe=FALSE,main="")
title(main="Bonferroni Confidence Region")
abline(v=BonSimCI.b[2,])
abline(h=BonSimCI.b[3,])
# Scheffe
confidence.ellipse(lm(HWFAT~AGE+ABS+TRICEPS),level=.90,
which.coef=c(3,4),Scheffe=TRUE,main="")
title(main="Scheffe Confidence Region")
abline(v=SchSimCI.b[2,])
abline(h=SchSimCI.b[3,])
# part c.
g <- 3
alpha <- 0.10
SC <- sqrt(g*qf(1-alpha,3,74))
TC <- qt(1-alpha/(2*g),74)
c(TC,SC)
RES <- predict(mult.model, newdata=data.frame(AGE=c(16,17,18),
ABS=c(10,11,8), TRICEPS=c(9,11,8)), se.fit=TRUE)
Yhath <- RES$fit
Syhath <- RES$se.fit
ll <- Yhath - TC*Syhath
ul <- Yhath + TC*Syhath
BCI <- cbind(Yhath,Syhath,ll,ul)
BCI
round(BCI,3)
# part d.
g <- 3
alpha <- 0.10
SC <- sqrt(g*qf(1-alpha,3,74))
TC <- qt(1-alpha/(2*g),74)
c(SC,TC)
# Use TC with equation 12.69
MSE <- anova(mult.model)[4,3]
MSE
s2yhathnew <- MSE + Syhath^2
Syhathnew <- sqrt(s2yhathnew)
ll <- Yhath - TC*Syhathnew
ul <- Yhath + TC*Syhathnew
SPI <- cbind(Yhath,Syhathnew,ll,ul)
SPI
detach(HSwrestler)
################################################################################