Extra assignment for session 6:

For whiteside data, fit simple linear regressions and check the 
normality of standardized residuals using Shapiro-Wilks test 
(shapiro.test). By simulation, compute a P-value that takes into 
account the dependence between residuals. Use n=1000
simulated datasets with normal errors.

library(MASS)
summary(whiteside)
gasB <- lm(Gas ~ Temp, whiteside, subset = Insul=="Before")
gasA <- update(gasB, subset = Insul=="After")
qqnorm(studres(gasB))
shapiro.test(studres(gasB))
qqnorm(studres(gasA))
shapiro.test(studres(gasA))
hist(studres(gasA))
# The residuals for the analysis of Insul==After data are left-skewed.
# We restrict our analysis to this subdataset

Gas0 <- whiteside$Gas[whiteside$Insul=="After"]
Temp0 <- whiteside$Temp[whiteside$Insul=="After"]
summary(lm(Gas0 ~ Temp0))
set.seed(210)
Gas <- rnorm(30,4.7238-0.2779*Temp0,0.355)
w<-c(shapiro.test(studres(lm(Gas ~ Temp0)))$statistic)
for (i in 2:1000)
{ 
  Gas <- rnorm(30,4.7238-0.2779*Temp0,0.355)
  w<-c(w,shapiro.test(studres(lm(Gas ~ Temp0)))$statistic)
}
names(w)<-NULL
pvec<-seq(0.005,0.10,0.005)
quantile(w,pvec)
z<-(w<=0.9205)
mean(z)
# The simulation P-value is 0.071, contrasting the given value of 0.02764.