## initialization
library(lme4)
beer <- read.table('beer.txt', header=T)
beer$tmt <- factor(beer$tmt)
beer$temp <- factor(beer$temp)
beer$press <- factor(beer$press)
beer$block <- factor(beer$block)

## 1 - b
beer.lme <- lmer( density ~ temp*press  + (1|block), data=beer)
anova(beer.lme)
#
#Analysis of Variance Table
#           Df  Sum Sq Mean Sq F value
#temp        1  33.676  33.676  9.5095
#press       1 181.693 181.693 51.3066
#temp:press  1   1.315   1.315  0.3713

pf(9.509, df1=1, df2=53, lower.tail=F)
# 0.0032
pf(51.307, df1=1, df2=53, lower.tail=F)
# 2.466e-09
pf(0.371, df1=1, df2=53, lower.tail=F)
#0.5448

## 1 - c
aov.beer <- aov( density ~ temp*press  + Error(block), data=beer)
summary(aov.beer)
#Error: block
#          Df  Sum Sq Mean Sq F value Pr(>F)
#Residuals  7 299.896  42.842               
#
#Error: Within
#           Df  Sum Sq Mean Sq F value    Pr(>F)    
#temp        1  33.676  33.676  9.5095  0.003243 ** 
#press       1 181.693 181.693 51.3066 2.466e-09 ***
#temp:press  1   1.315   1.315  0.3713  0.544881    
#Residuals  53 187.690   3.541                      
#---
# same as 1 - b.

#2
#initialization
milk <- read.table("milk.txt", header=T)
milk$day <- factor(milk$day)
milk$cow <- factor(milk$cow)

#2-b
mle.milk <- lmer(level ~ treatment + (1|day)+ (1 | day:cow), data=milk)
anova(mle.milk)#
#Linear mixed model fit by REML 
#Formula: level ~ treatment + (1 | day) + (1 | day:cow) 
#   Data: milk 
#   AIC   BIC logLik deviance REMLdev
# 267.7 280.8 -126.8    263.8   253.7
#Random effects:
# Groups   Name        Variance   Std.Dev.  
# day:cow  (Intercept) 6.9759e-14 2.6412e-07 <-Notice that this is almost zero!
# day      (Intercept) 4.7687e+01 6.9056e+00
# Residual             9.8081e+00 3.1318e+00
#Number of obs: 48, groups: day:cow, 24; day, 6
#
#Fixed effects:
#            Estimate Std. Error t value
#(Intercept)  61.9905     2.9604  20.940
#treatmentB   -0.1588     1.2785  -0.124
#treatmentC   -6.2191     1.2785  -4.864
#treatmentD   11.0150     1.2785   8.615
pf( 63.073, df1=3, df2=15, lower.tail=F)
#[1] 9.799215e-09

#2-c
aov.milk <- aov(level ~ treatment + Error(day/cow), data=milk)
summary(aov.milk)
#Error: day
#          Df Sum Sq Mean Sq F value Pr(>F)
#Residuals  5 1956.5   391.3               
#
#Error: day:cow
#          Df  Sum Sq Mean Sq F value    Pr(>F)    
#treatment  3 1855.88  618.63  78.686 2.096e-09 ***
#Residuals 15  117.93    7.86                      
#
#Error: Within
#          Df  Sum Sq Mean Sq F value Pr(>F)
#Residuals 24 264.588  11.024               
# not exactly same as 2-b. This happens when
# the estimated variance of (day:cow) is numerically zero.

#2-d
t.025 <- qt(0.025, df=15, lower.tail=F)

std.e <- 1.2785

m <- c(61.9905, 61.9905-0.1588 , 61.9905-6.2191 , 61.9905+11.0150)

names(m) <-c("A","B","C","D")

M <- outer(m, m, "-")

abs(M)

lsd <- t.025 * std.e
abs(M) > lsd
2*pt(abs(M/std.e), df=15, lower.tail=F)#
#             A            B            C            D
#A 1.000000e+00 9.027997e-01 2.061971e-04 3.416025e-07
#B 9.027997e-01 1.000000e+00 2.630609e-04 2.850173e-07
#C 2.061971e-04 2.630609e-04 1.000000e+00 8.685723e-10
#D 3.416025e-07 2.850173e-07 8.685723e-10 1.000000e+00

#3
news <-read.table('news.txt', header=T)
news$city <- factor(news$city)
news$street <- factor(news$street)
news$station <- factor(news$station)
lm.news = lm(sale ~ cover + city + city:street + city:station, data=news)

anova(lm.news)
#Analysis of Variance Table

#Response: sale
#             Df Sum Sq Mean Sq F value  Pr(>F)  
#cover         2 6.6023  3.3011  7.4321 0.02378 *
#city          1 2.1376  2.1376  4.8125 0.07071 .
#city:street   4 1.8586  0.4647  1.0461 0.45651  
#city:station  4 2.0491  0.5123  1.1534 0.41627  
#Residuals     6 2.6650  0.4442                  


#4
chem = read.table("chem.txt", header=T)
chem$machine <- factor(chem$machine)
str(chem)


#4-a numerical summary and interaction plot
with(chem, table(catalyst,machine,tmt)) 

with(chem, interaction.plot(tmt, catalyst, response))


#4-b
lm.chem = lm(response ~ machine + tmt*catalyst, data=chem)
anova(lm.chem)

#Analysis of Variance Table
#Response: response
#             Df  Sum Sq Mean Sq F value   Pr(>F)   
#machine       3 10.5682  3.5227  4.0940 0.026158 * 
#tmt           2  9.6145  4.8073  5.5869 0.015371 * 
#catalyst      1 11.5278 11.5278 13.3973 0.002321 **
#tmt:catalyst  2 12.6832  6.3416  7.3700 0.005897 **
#Residuals    15 12.9069  0.8605     

summary(lm.chem)

#Coefficients:
#               Estimate Std. Error t value Pr(>|t|)    
#(Intercept)      3.2849     0.5680   5.783 3.61e-05 ***
#machine2         1.2841     0.5356   2.398  0.02996 *  
#machine3        -0.4742     0.5356  -0.885  0.38987    
#machine4         0.6463     0.5356   1.207  0.24623    
#tmtB            -0.0907     0.6559  -0.138  0.89186    
#tmtC             0.2054     0.6559   0.313  0.75851    
#catalystY        0.1873     0.6559   0.286  0.77909    
#tmtB:catalystY  -1.2139     0.9276  -1.309  0.21037    
#tmtC:catalystY  -3.5065     0.9276  -3.780  0.00182 ** 


lmer.chem = lmer(response ~ tmt*catalyst + (1|machine), data=chem)
anova(lmer.chem)

#Analysis of Variance Table
#             Df  Sum Sq Mean Sq F value
#tmt           2  9.6145  4.8073  5.5869
#catalyst      1 11.5278 11.5278 13.3973
#tmt:catalyst  2 12.6832  6.3416  7.3700
summary(lmer.chem)
#Random effects:
# Groups   Name        Variance Std.Dev.
# machine  (Intercept) 0.44371  0.66612 
# Residual             0.86046  0.92761 
#Number of obs: 24, groups: machine, 4
#
#Fixed effects:
#               Estimate Std. Error t value
#(Intercept)      3.6489     0.5710   6.390
#tmtB            -0.0907     0.6559  -0.138
#tmtC             0.2054     0.6559   0.313
#catalystY        0.1873     0.6559   0.286
#tmtB:catalystY  -1.2139     0.9276  -1.309
#tmtC:catalystY  -3.5065     0.9276  -3.780

#4-c model assumption

#for fixed effect model
par(mfrow=c(1,2))
plot(lm.chem, which=1)
plot(lm.chem, which=2)

#for mixed effect model
par(mfrow=c(1,3))
plot(resid(lmer.chem)~ fitted(lmer.chem))
abline(h=0)
qqnorm(resid(lmer.chem))
qqnorm(ranef(lmer.chem)$machine$"(Intercept)")