# R program for additional exercise 7.6 (VHM 802)
# demonstration of: analysis without outlier
#                   least squares means

sph <- read.csv("r:/hs07_6.csv")
summary (sph.full <- lm(volume ~ as.factor(lot)+as.factor(trough), data=sph))

# analysis without observation equal to 7.5
sph.lm <- lm(volume ~ as.factor(lot)+as.factor(trough), data=sph[sph$volume!=7.5,])
summary(sph.lm)

# least squares means for lots, with SE and 95% CI
# manual construction: note 6 troughs
coef <- coef(sph.lm)
rvcov <- vcov(sph.lm)
w<- matrix(c(6,0,0,0,0,0, 0,0,0,0,0,0, 1,1,1,1,1,
             6,6,0,0,0,0, 0,0,0,0,0,0, 1,1,1,1,1,
             6,0,6,0,0,0, 0,0,0,0,0,0, 1,1,1,1,1,
             6,0,0,6,0,0, 0,0,0,0,0,0, 1,1,1,1,1,
             6,0,0,0,6,0, 0,0,0,0,0,0, 1,1,1,1,1,
             6,0,0,0,0,6, 0,0,0,0,0,0, 1,1,1,1,1,
             6,0,0,0,0,0, 6,0,0,0,0,0, 1,1,1,1,1,
             6,0,0,0,0,0, 0,6,0,0,0,0, 1,1,1,1,1,
             6,0,0,0,0,0, 0,0,6,0,0,0, 1,1,1,1,1,
             6,0,0,0,0,0, 0,0,0,6,0,0, 1,1,1,1,1,
             6,0,0,0,0,0, 0,0,0,0,6,0, 1,1,1,1,1,
             6,0,0,0,0,0, 0,0,0,0,0,6, 1,1,1,1,1), ncol=12)/6; w
est<- t(w)%*%coef
var <- diag(t(w)%*%vcov%*%w)
se <- matrix(sqrt(var), ncol=1)
tstar <- qt(0.05/2, lower.tail=FALSE, df.residual(sph.lm))
low <- est-tstar*se
upp <- est+tstar*se
data.frame(est,se,low,upp)
# using add-on libraries
library(estimability)
library(lsmeans)
lsmeans(sph.lm, "lot")