##########################################################################################################
######    Cleaned-up version of code
##########################################################################################################
source("http://www.stat.uchicago.edu/~pmcc/projects/regress.R")
homedir <- "http://www.stat.uchicago.edu/~pmcc/projects/ch01/"
#
Y <- scan(paste0(homedir, "rats.dat"))	# unstructured numeric read
nobs <- 120;  nsites <- 5;  nrats <- 24
site <- gl(nsites, 1, nobs)	# site factor
rat <- gl(nrats, nsites, nobs)	# rat indicator
treat <- gl(2, 70, 120)		# treatment level
levels(treat) <- c("HypO2", "Cntrl")
treat <- relevel(treat, ref="Cntrl")	# make "Cntrl" the reference level
w <- (Y > 0)			# non-missing
logY <- log(Y + 1-w)		# use log(Y) as response
ratdata <- data.frame(logY=logY, site=site, rat=rat, treat=treat)[w,]

mainfit <- regress(logY~site+treat, ~rat, data=ratdata)		# typical syntax for regress() with block factor
summary(mainfit)

###  REML LR for treat effect	(H0: no treatment effect)
Ks <- model.matrix(~site, data=ratdata)
fit0 <- regress(logY~site, ~rat, data=ratdata, kernel=Ks)	# default kernel for REML is Ks
fit1 <- regress(logY~site+treat, ~rat, data=ratdata, kernel=Ks)	# overrides the default
2*(fit1$llik - fit0$llik)	# on one d.f.   Essentially the same as (betahat/se)^2

###  REML LR statistic for site effects	(H0: no site effects)
Kt <- model.matrix(~treat, data=ratdata)
fit0 <- regress(logY~treat, ~rat, data=ratdata, kernel=Kt)	# default kernel is Kt
fit1 <- regress(logY~site+treat, ~rat, data=ratdata, kernel=Kt)	# overrides the default
2*(fit1$llik - fit0$llik)	# on nsites-1 d.f.

### ... or if you insist on ML rather than REML use the zero kernel
K <- 0	#
fit0 <- regress(logY~treat, ~rat, data=ratdata, kernel=K)
fit1 <- regress(logY~site+treat, ~rat, data=ratdata, kernel=K)
2*(fit1$llik - fit0$llik)       # on nsites-1 d.f.


##########################################################################################################
######    Earlier working version of code (archival use only)
##########################################################################################################
homedir <- "/Users/pmcc/projects/ch01/"
source("http://www.stat.uchicago.edu/~pmcc/courses/regress.R")
y <-scan(paste(homedir, "rats.dat", sep=""))
n <- 120
w <- (y != 0)	# missing values coded as zero
site <- gl(5, 1, 120)[w > 0]	# as a five-level factor
rat <- gl(24, 5, 120)[w > 0]	# 24 rats
treat <- gl(2, 70, 120)[w > 0]	# treatment level
levels(treat) <- c("hypO2", "cntrl")
ly <- log(y[w > 0])		# use log(y) as response

summary(y[w>0])
lybar <- tapply(ly, rat, mean)
trt <- c(rep("hyp02", 14), rep("cntrl", 10))
tapply(lybar, trt, mean)
tapply(lybar, trt, var)
round(tapply(y[w], rat, mean), 3)
round(tapply(y[w], site, mean), 3)


K <- model.matrix(~site)
K1 <- model.matrix(~site+treat)
fit0 <- regress(ly~1, V=~rat)
fit0a <- regress(ly~site, V=~rat, kernel=1)
fit1 <- regress(ly~site, V=~rat)
fit2 <- regress(ly~site, V=~rat, kernel=K)
fit2a <- regress(ly~site+treat, V=~rat)
fit3 <- regress(ly~site+treat, V=~rat, kernel=K1)
fit3a <- regress(ly~site*treat, V=~rat, kernel=K1)

vsite <- as.numeric(site); vsite <- abs(outer(vsite, vsite, "-"))
vsite <- exp(-vsite/100)
fit <- regress(ly~treat, V=~rat+site+vsite);  fit$llik

ratC <- outer(rat, rat, "==")*outer(treat=="cntrl", treat=="cntrl")
ratT <- outer(rat, rat, "==")*outer(treat=="hypO2", treat=="hypO2")
fit0 <- regress(ly~site+treat, ~rat)
fit1 <- regress(ly~site+treat, ~ratC+ratT)
c(fit0$llik, fit1$llik)

fit0 <- lm(y[w]~rat+site)
fv2 <- fit0$fitted^2
fit1 <- lm(y[w]~rat+site+fv2)
summary(fit1)
1 - 2*fit1$beta[29]*mean(y[w])

fit <- regress(ly~site+treat, ~rat)
install.packages("lme4",
   repos=c("http://lme4.r-forge.r-project.org/repos",
      getOption("repos")[["CRAN"]]))
library(lme4)
fitlme <- lmer(ly~site+treat+(1|rat))

site1 <- -2*(site==1) - 1*(site==2) + 0*(site==3) + 1*(site==4) + 2*(site==5)
site2 <-  2*(site==1) - 1*(site==2) - 2*(site==3) - 1*(site==4) + 2*(site==5)
site3 <- -1*(site==1) + 2*(site==2) - 0*(site==3) - 2*(site==4) + 1*(site==5)
site4 <-  1*(site==1) - 4*(site==2) + 6*(site==3) - 4*(site==4) + 1*(site==5)
fit2 <- regress(ly~site1+site2+site3+site4+treat, ~rat)
fitlme2 <- lmer(ly~site1+site2+site3+site4+treat+(1|rat))

###  For Exercise 1.14
D0 <- diag(treat=="cntrl");  D1 <- diag(treat=="hypO2")
fit0 <- regress(ly~treat+site, ~rat)
fit1 <- regress(ly~treat+site, ~rat+D1)
fit2 <- regress(ly~treat+site, ~rat+D0+D1, identity=FALSE, start=c(3,2,2))