source("http://www.uvm.edu/~rsingle/stat221/data/scripts-221.R")

dat <- bookdata("ch05q02.txt")
head(dat)

mod <- lm(SBP~QUET+AGE, data=dat)

hatvalues(mod)[1:5]      #leverage or "hat" values

leverages <- hatvalues(mod) 
plot(dat$SBP, leverages)
mean(leverages)      #(2+1)/32 = 0.09375
abline(h= 2*mean(leverages) )
identify(dat$SBP, leverages, dat$PERSON)
dat$PERSON[leverages>2*.09375]

high.lev = leverages>2*.09375
plot(dat$QUET,dat$AGE)
points(dat$QUET[high.lev],dat$AGE[high.lev],pch="X",col=2)
#------------------------------------------------
# The hat() function works on the design matrix and gives the
# hatvalues() for the model with that design matrix

# Recode the data in terms of Y & X
Y  <- dat$SBP
X1 <- dat$QUET
X2 <- dat$AGE
X0 <- rep(1, times=length(X1))

# Create the design matrix X
X <- cbind(X0,X1,X2)  # bind column vectors together
rm(X0,X1,X2)
head(X)

hatvalues(mod)[1:5] #first 5 leverage values
hat(X)[1:5]         #X is the design matrix
hat(X[,-1], intercept = TRUE)[1:5]  #adds back a column of 1s
hat(X[,-1], intercept = FALSE)[1:5] #wrong

 H = X %*% inverse(t(X) %*% X) %*% t(X) #H = X (X'X)^-1 X'
 diag(H)[1:5] 
 dim(H)

#------------------------------------------------
#Residuals
 rstandard(mod)[1:5]      #internally studentized residuals
 rstudent(mod)[1:5]       #externally studentized residuals (a.k.a. jackknife)
 hatvalues(mod)[1:5]      #leverage or "hat" values
 cooks.distance(mod)[1:5] #Cook's distances

 plot(mod, which=1) #resid vs fit
 plot(mod, which=2) #QQ of std res
 plot(mod, which=3) #scale location plot (sqrt(std.res) vs fit)
 plot(mod, which=4) #cooks plot

 max(rstandard(mod)); abs( qt(.05,29  ))
 max(rstudent(mod));  abs( qt(.05,29-1))

#------------------------------------------------
 dat2 = dat
 dat2$SBP[3] = dat2$SBP[3]+50 #Create an outlier at obs#3
 mod2 <- lm(SBP~QUET+AGE, data=dat2)
 rstandard(mod2)[1:5]      #internally studentized residuals
 rstudent(mod2)[1:5]     

#------------------------------------------------
#Computations
 summary(mod)
 mod.summary = summary(mod)
 s=mod.summary$sigma
 mod$resid[1:5]
 rstandard(mod)[1:5]      
 
 h.ii = hatvalues(mod)
 my.rstandard = (mod$resid - 0)/(s*sqrt(1-h.ii))
 my.rstandard[1:5]
 rstandard(mod)[1:5]