alcohol <- read.delim("alcohol.txt")
attach(alcohol)

## Slide 3
summary(fit <- lm(Metabol~Sex))
mean(Metabol)
by(Metabol,Sex,mean)
VarB <- vcov(fit)
SE <- sqrt(diag(VarB))
2*(1-pt(coef(fit)[2]/SE[2],fit$df.residual))

## Equivalent to t-test
t.test(Metabol~Sex,var.equal=TRUE)

## Set up indicator variables
Male <- 1*(Sex=="Male")
Female <- 1*(Sex=="Female")
Alcoholic <- 1*(Alcohol=="Alcoholic")
Nonalc <- 1*(Alcohol=="Non-alcoholic")

## Slide 8
summary(fit2 <- lm(Metabol~0+Female+Male))

## Slide 10
summary(lm(Metabol~Gastric))
CGastric <- Gastric - mean(Gastric)
summary(lm(Metabol~CGastric))
summary(lm(Metabol~Alcohol))

## Slide 12-13
summary(lm(Metabol~Sex+Alcohol))
summary(lm(Metabol~Sex+Gastric))

## Slide 15
summ <- summary(fit <- lm(Metabol~Male*Alcoholic))
lambda <- c(0,0,1,1)

## Writing a function to do what an ESTIMATE statement in SAS does
estimate <- function(lambda,fit)
  {
    Estimate <- crossprod(lambda,fit$coef)
    summ <- summary(fit)
    SE <- sqrt(summ$sigma^2*t(lambda)%*%summ$cov.unscaled%*%lambda)
    t <- Estimate/SE
    p <- 2*pt(abs(t),fit$df.resid,low=FALSE)
    val <- c(Estimate,SE,t,p)
    names(val) <- c("Estimate","SE","t","p")
    return(val)
  }
estimate(lambda,fit)

## Slide 22
summary(lm(Metabol~Female*Nonalc))
summary(lm(Metabol~Male*Gastric))

## Trellis plot
xyplot(Metabol~Gastric|Sex,type=c("p","r"))

## Slides 30-31
summary(lm(Metabol~Alcoholic+Male+Gastric))
summary(lm(Metabol~Alcoholic+Male*Gastric))

## Slide 34
trellis.par.set(plot.symbol=list(pch=19))
xyplot(Metabol~Gastric|Sex+Alcohol,type=c("p","r"))

## Slide 35
summary(lm(Metabol~Alcoholic*Male*Gastric))