###################################################
### chunk number 1: 
###################################################
library(lattice)
options(width=65,show.signif.stars=FALSE)
options(SweaveHooks=list(fig=function() {
ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
ltheme$strip.background$col <- "transparent" ## change strip bg
ltheme$plot.symbol$col <- "black"
ltheme$axis.text$cex <- 1.2
ltheme$par.xlab.text$cex <- 1.2
ltheme$par.ylab.text$cex <- 1.2
lattice.options(default.theme = ltheme)      ## set as default
}))
source("my.boxcox.R")
library(xtable)
library(car)
#library(lattice)
#ltheme <- canonical.theme(color = FALSE)     ## in-built B&W theme
#ltheme$strip.background$col <- "transparent" ## change strip bg
#lattice.options(default.theme = ltheme)      ## set as default


###################################################
### chunk number 2: 
###################################################
solv <- read.csv('Data/solve.csv')
solv$group <- relevel(solv$group,ref="row")


###################################################
### chunk number 3:  eval=FALSE
###################################################
## data(solv,package="SMIR")


###################################################
### chunk number 4:  eval=FALSE
###################################################
## library(lattice)
## print(xyplot(time~eft|group,data=solv,type=c("p","r")))


###################################################
### chunk number 5: solve.plot eval=FALSE
###################################################
## print(xyplot(time ~ eft | group, data=solv,
##        panel=function(x, y, ...) {
##            panel.xyplot(x, y, ...)
##            mydata.lm <- lm(y~x)
##            llines(x,fitted(mydata.lm),...)
##        }))


###################################################
### chunk number 6: 
###################################################
print(xyplot(time ~ eft | group, data=solv,
       panel=function(x, y, ...) {
           panel.xyplot(x, y, ...)
           mydata.lm <- lm(y~x)
           llines(x,fitted(mydata.lm),...)
       }))


###################################################
### chunk number 7:  eval=FALSE
###################################################
## solv$group <- gl(2,12,labels=c('row','corner'))


###################################################
### chunk number 8: solve.models
###################################################
solv$group <- relevel(solv$group,'row') 
solv.glm5 <- glm(time~1,data=solv) 
solv.glm4 <- update(solv.glm5,.~group)
solv.glm3 <- update(solv.glm5,.~eft)
solv.glm2 <- update(solv.glm4,.~.+eft) 
solv.glm1 <- update(solv.glm2,.~.+eft:group) 


###################################################
### chunk number 9: solv.model.figures
###################################################
# Make one large data frame with data set repeated (once per model)
# and fitted values added. Thanks to Paul Murrell for this one
all <- rbind(cbind(solv,model="Model 5", fitted=fitted(solv.glm5)),
             cbind(solv, model="Model 4", fitted=fitted(solv.glm4)),
             cbind(solv, model="Model 3", fitted=fitted(solv.glm3)),
             cbind(solv, model="Model 2", fitted=fitted(solv.glm2)),
             cbind(solv, model="Model 1", fitted=fitted(solv.glm1)))
print(xyplot(time ~ eft | group * model, data=all, subscripts=TRUE,
#             scales=list(y=list(labels=NULL)),
       panel=function(x, y, subscripts, ...) {
           panel.xyplot(x, y, ...)
           llines(x, all$fitted[subscripts], ...)
       }))


###################################################
### chunk number 10: 
###################################################
# Make one large data frame with data set repeated (once per model)
# and fitted values added. Thanks to Paul Murrell for this one
all <- rbind(cbind(solv,model="Model 5", fitted=fitted(solv.glm5)),
             cbind(solv, model="Model 4", fitted=fitted(solv.glm4)),
             cbind(solv, model="Model 3", fitted=fitted(solv.glm3)),
             cbind(solv, model="Model 2", fitted=fitted(solv.glm2)),
             cbind(solv, model="Model 1", fitted=fitted(solv.glm1)))
print(xyplot(time ~ eft | group * model, data=all, subscripts=TRUE,
#             scales=list(y=list(labels=NULL)),
       panel=function(x, y, subscripts, ...) {
           panel.xyplot(x, y, ...)
           llines(x, all$fitted[subscripts], ...)
       }))


###################################################
### chunk number 11:  eval=FALSE
###################################################
## glm(time~eft+group,data=solv)


###################################################
### chunk number 12: solve.models eval=FALSE
###################################################
## solv.glm5 <- glm(time~1,data=solv)
## solv.glm4 <- update(solv.glm5,.~.+group)
## solv.glm2 <- update(solv.glm4,.~.+eft)
## solv.glm1 <- update(solv.glm2,.~.+eft:group)


###################################################
### chunk number 13: second.ANOVA.table eval=FALSE
###################################################
## (anova2 <- anova(solv.glm5,solv.glm4,solv.glm2,solv.glm1,test="F"))


###################################################
### chunk number 14: 
###################################################
(anova2 <- anova(solv.glm5,solv.glm4,solv.glm2,solv.glm1,test="F"))
#library(xtable)
#xtable(anova2,
#caption="ANOVA table for the \texttt{solv} dataset.",label="tab:solv.anova",
#digits=c(0,0,0,0,0,3,3))


###################################################
### chunk number 15: reorder.anova.table eval=FALSE
###################################################
## solv.glm3 <- update(solv.glm5,.~.+eft)
## (anova3 <- anova(solv.glm5,solv.glm3,solv.glm2,solv.glm1,test="F"))


###################################################
### chunk number 16: 
###################################################
solv.glm3 <- update(solv.glm5,.~.+eft)
(anova3 <- anova(solv.glm5,solv.glm3,solv.glm2,solv.glm1,test="F"))


###################################################
### chunk number 17: solv.fitted eval=FALSE
###################################################
## plot(time~eft,data=solv)
## lines(solv$eft,fitted(solv.glm3))


###################################################
### chunk number 18: 
###################################################
solv.glm2$call
coef(summary(solv.glm2))


###################################################
### chunk number 19: confint.solve.glm2
###################################################
confint(solv.glm2)


###################################################
### chunk number 20: 
###################################################
round(coef(summary(solv.glm1)),4)


###################################################
### chunk number 21: solve.glm6
###################################################
solv$eft.row <- solv$eft*(solv$group=='row')
solv.glm6 <- glm(time~eft.row+group,data=solv)
coef(summary(solv.glm6))
anova(solv.glm6,solv.glm1,test="F")


###################################################
### chunk number 22:  eval=FALSE
###################################################
## summary(solv.glm6)$sigma^2*summary(solv.glm6)$cov.unscaled


###################################################
### chunk number 23: 
###################################################
Qx <- function(x) 3.3315*x^2 -2*39.7825*x-17064
uniroot(Qx,lower=0,upper=200)$root


###################################################
### chunk number 24: solve.sex
###################################################
solv$sex <- gl(2,6,labels=c('girl','boy'))


###################################################
### chunk number 25: solve.sex.group.eft
###################################################
round(coef(summary(glm(time~sex*group*eft,data=solv))),4)


###################################################
### chunk number 26: group.sex.coplot eval=FALSE
###################################################
## #coplot(time~eft|group*sex,data=solv,panel=panel.smooth)
## print(xyplot(time~eft|group*sex,data=solv,type=c("p","r")))
## #panel=function(...){
## #panel.xyplot(...)
## #panel.lmline(...)}))


###################################################
### chunk number 27: group.sex.coplot eval=FALSE
###################################################
## print(xyplot(time ~ eft | group * sex, data=solv,
## #             scales=list(x=list(relation="free"),y=list(relation="free")),
##        panel=function(x, y, ...) {
##            panel.xyplot(x, y, ...)
##            mydata.lm <- lm(y~x)
##            llines(x,fitted(mydata.lm),...)
## }))


###################################################
### chunk number 28: 
###################################################
print(xyplot(time ~ eft | group * sex, data=solv,
#             scales=list(x=list(relation="free"),y=list(relation="free")),
       panel=function(x, y, ...) {
           panel.xyplot(x, y, ...)
           mydata.lm <- lm(y~x)
           llines(x,fitted(mydata.lm),...)
}))


###################################################
### chunk number 29: 
###################################################
source("NPL.bands.R")


###################################################
### chunk number 30:  eval=FALSE
###################################################
## library(SMIR)


###################################################
### chunk number 31: solve.bounds eval=FALSE
###################################################
## #source('NPL.bands.R')
## #env <- NPL.bands(resid(solv.glm1))
## print(xyplot(lower+upper~x,
##              data=NPL.bands(resid(solv.glm1)),
##              pch=20,
##              xlab="residual",ylab="cumulative proportion",
##              panel=function(x,y,...){
##                panel.xyplot(x,y,...)
##                panel.curve(pnorm(x,0,sqrt(summary(solv.glm1)$dispersion)))
## }))


###################################################
### chunk number 32: 
###################################################
#source('NPL.bands.R')
#env <- NPL.bands(resid(solv.glm1))
print(xyplot(lower+upper~x,
             data=NPL.bands(resid(solv.glm1)),
             pch=20,
             xlab="residual",ylab="cumulative proportion",
             panel=function(x,y,...){
               panel.xyplot(x,y,...)
               panel.curve(pnorm(x,0,sqrt(summary(solv.glm1)$dispersion)))
}))


###################################################
### chunk number 33: 
###################################################
library(MASS)
round(stdres(solv.glm1),3)


###################################################
### chunk number 34: tab.solve.influence eval=FALSE
###################################################
## sort(round(hatvalues(solv.glm1),2))


###################################################
### chunk number 35: 
###################################################
sort(round(hatvalues(solv.glm1),2))


###################################################
### chunk number 36: 
###################################################
sort(round(hatvalues(glm(time~sex*group*eft,data=solv)),2))


###################################################
### chunk number 37: 
###################################################
model <- glm(time~sex*group*eft,data=solv)
summary(update(model,subset=(hatvalues(model)<(2*model$rank/length(solv$time)))))


###################################################
### chunk number 38: 
###################################################
trees <- read.csv('../Old.Book/smir/Rversion/Data/trees.csv')


###################################################
### chunk number 39:  eval=FALSE
###################################################
## data(trees,package="SMIR")


###################################################
### chunk number 40:  eval=FALSE
###################################################
## plot(~.,data=trees)


###################################################
### chunk number 41: trees.plot eval=FALSE
###################################################
## library(car)
## scatterplot.matrix(~d+h+v,data=trees,smooth=FALSE,col=c(1,1,1))


###################################################
### chunk number 42: 
###################################################
library(car)
scatterplot.matrix(~d+h+v,data=trees,smooth=FALSE,col=c(1,1,1))


###################################################
### chunk number 43: 
###################################################
trees.glm1 <- glm(v~h+d,data=trees)
summary.lm(trees.glm1)


###################################################
### chunk number 44: 
###################################################
round(hatvalues(trees.glm1),3)


###################################################
### chunk number 45: tree.glm1.residuals.plot eval=FALSE
###################################################
## scatter.smooth(trees$d,residuals(trees.glm1,type="working"),xlab='d',
## ylab=expression(paste('residual from model  ',v,' ~ ',h+d))) 


###################################################
### chunk number 46: 
###################################################
scatter.smooth(trees$d,residuals(trees.glm1,type="working"),xlab='d',
ylab=expression(paste('residual from model  ',v,' ~ ',h+d))) 


###################################################
### chunk number 47:  eval=FALSE
###################################################
## require(MASS)


###################################################
### chunk number 48: tree.boxcox
###################################################
tree.bc <- boxcox(trees.glm1,plotit=FALSE,lambda=seq(-2,2,by=0.5))
names(tree.bc$y) <- tree.bc$x
print(as.vector(tree.bc$y),digits=3)


###################################################
### chunk number 49: 
###################################################
tree.boxcox <- boxcox(trees.glm1,plotit=FALSE,lambda=seq(from=0,to=1,by=0.1))
tree.boxcox$x[which.max(tree.boxcox$y)]
lambda.max <- tree.boxcox$x
logl <- tree.boxcox$y
names(logl) <- tree.boxcox$x
print(logl,digits=3)
rbind(lambda.max,logl)
logl <- logl-max(logl)
a <- seq(from=0,to=1,by=0.1)[logl>log(0.15)]


###################################################
### chunk number 50: trees.cuberoot.plots eval=FALSE
###################################################
## trees$vth <- trees$v^(1/3)
## par(mfrow=c(1,2))
## plot(vth~d+h,data=trees,ask=FALSE)
## #plot(vth~d,data=trees)
## #plot(vth~h,data=trees)


###################################################
### chunk number 51: 
###################################################
trees$vth <- trees$v^(1/3)
par(mfrow=c(1,2))
plot(vth~d+h,data=trees,ask=FALSE)
#plot(vth~d,data=trees)
#plot(vth~h,data=trees)


###################################################
### chunk number 52: 
###################################################
trees.vth.glm1 <- glm(vth~h + d,data=trees)
summary.lm(trees.vth.glm1)


###################################################
### chunk number 53: 
###################################################
round(rbind(fitted=fitted(trees.vth.glm1)^3,observed=trees$v),2)


###################################################
### chunk number 54: trees.log.model.plots eval=FALSE
###################################################
## trees$lv <- log(trees$v)
## trees$ld <- log(trees$d)
## trees$lh <- log(trees$h)
## par(mfrow=c(1,2))
## plot(v~d+h,data=trees,ask=FALSE,log="xy")
## #
## #plot(lv~ld,data=trees,xlab='log diameter',ylab='log volume')
## #plot(lv~lh,data=trees,xlab='log height',ylab='log volume')


###################################################
### chunk number 55: 
###################################################
trees$lv <- log(trees$v)
trees$ld <- log(trees$d)
trees$lh <- log(trees$h)
par(mfrow=c(1,2))
plot(v~d+h,data=trees,ask=FALSE,log="xy")
#
#plot(lv~ld,data=trees,xlab='log diameter',ylab='log volume')
#plot(lv~lh,data=trees,xlab='log height',ylab='log volume')


###################################################
### chunk number 56: 
###################################################
trees.lv.glm1 <- glm(lv~lh+ld,data=trees)
summary.lm(trees.lv.glm1)


###################################################
### chunk number 57: 
###################################################
library(MASS)
boxcox(v~lh+ld,data=trees,plotit=FALSE,lambda=seq(from=-2,to=2,by=0.5))


###################################################
### chunk number 58: 
###################################################
trees <- transform(trees,df=d/12)
trees <- transform(trees,ldf=log(df))


###################################################
### chunk number 59: 
###################################################
(trees.log.glm1 <- glm(lv~ldf+lh,data=trees))


###################################################
### chunk number 60: trees.offset.model
###################################################
trees$z <- trees$lh + 2*trees$ldf
trees.log.glm2 <- glm(lv~offset(z),data=trees)
summary.lm(trees.log.glm2)


###################################################
### chunk number 61: 
###################################################
1-(sum(resid(trees.log.glm2)^2))/(sum((trees$lv-mean(trees$lv))^2))


###################################################
### chunk number 62: 
###################################################
log(c(pi/4,pi/12))


###################################################
### chunk number 63: 
###################################################
fv <- exp(fitted(trees.log.glm2))
round(rbind("f median"=fv,"observed"=trees$v),1)
#plot(trees$v,fv)


###################################################
### chunk number 64: 
###################################################
(trees.glm1 <- glm(v~h+d,data=trees,
                  family=quasi(power(lambda=1/3),variance="constant")))
#summary.lm(trees.glm1)


###################################################
### chunk number 65: 
###################################################
qqnorm(resid(trees.glm1))
qqline(resid(trees.glm1))


###################################################
### chunk number 66: 
###################################################
trees.glm2 <- glm(v~lh+ld,data=trees,quasi(link='log',variance='constant'))
summary.glm(trees.glm2)


###################################################
### chunk number 67: 
###################################################
trees.log.glm3 <- glm(lv~lh+ld,data=trees)
new <- data.frame(lh=log(80),ld=log(15))
predict.lm(trees.log.glm3,new,interval='prediction')


###################################################
### chunk number 68: PRESS
###################################################
PRESS <- function(x)sum((resid(x)^2/(1-hatvalues(x))^2))


###################################################
### chunk number 69: CVE
###################################################
CVE <- function(x)sum((resid(x)^2/ 
                       (1-hatvalues(x))^2))/length(fitted(x))
RSS <- function(x)sum(resid(x)^2)


###################################################
### chunk number 70:  eval=FALSE
###################################################
## data(trees,package="SMIR")


###################################################
### chunk number 71: 
###################################################
trees <- read.csv('Data/trees.csv')
options(digits=4)


###################################################
### chunk number 72: trees.glm
###################################################
trees <- transform(trees, 
                   lv=log(v),
                   lh=log(h),
                   ld=log(d))
#trees$lv <- log(trees$v)
#trees$lh <- log(trees$h)
#trees$ld <- log(trees$d)
trees.glm <- glm(lv~ld+lh,data=trees)
#summary(trees.glm)
PRESS(trees.glm)
CVE(trees.glm)


###################################################
### chunk number 73: 
###################################################
R2 <- function(model) {
  yhat <- fitted(model)
  ehat <- residuals(model)
  y <- yhat + ehat
  cor(yhat,y)^2
}  
R2CV <- function(model){
  yhat <- fitted(model)
  ehat <- residuals(model)
  y <- yhat + ehat
  hi <- influence(model)$hat
  yhatcve <- (yhat-hi*y)/(1-hi)
  cor(y,yhatcve)^2
}


###################################################
### chunk number 74: 
###################################################
R2 <- function(x) {f <- fitted(x); 1 - (sum(resid(x)^2))/(sum((x$y-mean(f))^2))}
R2CV <- function(x){yj <- (fitted(x)-hatvalues(x)*x$y)/(1-hatvalues(x));cor(x$y,yj)^2}


###################################################
### chunk number 75: trees.glm1
###################################################
trees$z <- 2*trees$ld+trees$lh
trees.glm1 <- glm(lv~offset(z),data=trees)
#summary(trees.glm1)
RSS(trees.glm1)
summary(trees.glm1)$dispersion
PRESS(trees.glm1)# <- sum((resid(trees.glm1)^2/(1-hatvalues
R2(trees.glm1)
CVE(trees.glm1)# <- sum((resid(trees.glm1)^2/(1-hatvalues(trees.glm1))^2))/
#  trees.glm1$df.residual)
R2CV(trees.glm1)


###################################################
### chunk number 76: 
###################################################
car <- read.csv('Data/car.csv')


###################################################
### chunk number 77:  eval=FALSE
###################################################
## data(car, package="SMIR")


###################################################
### chunk number 78: car
###################################################
row.names(car) <- abbreviate(car$model,minlength=14)
car$model <- NULL
format(car)


###################################################
### chunk number 79: car.plots eval=FALSE
###################################################
## par(mfrow=c(3,3),mar=c(4,4,2,1.5))
## plot(mpg~s+c+t+g+disp+hp+cb+drat+wt,data=car,ask=FALSE)


###################################################
### chunk number 80: 
###################################################
par(mfrow=c(3,3),mar=c(4,4,2,1.5))
plot(mpg~s+c+t+g+disp+hp+cb+drat+wt,data=car,ask=FALSE)


###################################################
### chunk number 81: car.log.plots eval=FALSE
###################################################
## car <- transform(car,lmpg=log(mpg),
##                  ldis=log(disp))
## plot(lmpg~ldis,data=car,xlab='log displacement',ylab='log mpg')


###################################################
### chunk number 82: 
###################################################
car <- transform(car,lmpg=log(mpg),
                 ldis=log(disp))
plot(lmpg~ldis,data=car,xlab='log displacement',ylab='log mpg')


###################################################
### chunk number 83: car.log.plots2 eval=FALSE
###################################################
## car <- transform(car,lhp=log(hp),
##                  lwt=log(wt*1000))
## par(mfrow=c(1,2),mar=c(4,4,2,1.5))
## plot(lmpg~lhp+lwt,data=car,ask=FALSE)


###################################################
### chunk number 84: 
###################################################
car <- transform(car,lhp=log(hp),
                 lwt=log(wt*1000))
par(mfrow=c(1,2),mar=c(4,4,2,1.5))
plot(lmpg~lhp+lwt,data=car,ask=FALSE)


###################################################
### chunk number 85: 
###################################################
car.bc <- boxcox(mpg~s+c+t+g+ldis+lhp+cb+drat+lwt,data=car,
plotit=FALSE,lambda=seq(-2,2,by=0.5))


###################################################
### chunk number 86: 
###################################################
names(car.bc$y) <- car.bc$x
round(car.bc$y,1)


###################################################
### chunk number 87: 
###################################################
car.bc2 <- boxcox(mpg~s+c+t+g+disp+hp+cb+drat+wt,data=car,
plotit=FALSE,lambda=seq(-2,2,by=0.5))


###################################################
### chunk number 88: 
###################################################
names(car.bc2$y) <- car.bc2$x
round(car.bc2$y,1)


###################################################
### chunk number 89: 
###################################################
#car$d <- ifelse(rownames(car)=="MERCEDES 240D",1,0)
car$d <- c(rep(0,7),1,rep(0,24))
summary(glm(lmpg~s+c+t+g+ldis+lhp+cb+drat+lwt+d,data=car))


###################################################
### chunk number 90: car.glm
###################################################
summary(car.glm <- glm(lmpg~s+c+t+g+ldis+lhp+cb+drat+lwt,data=car))


###################################################
### chunk number 91: 
###################################################
round(sort(hatvalues(car.glm)),3)


###################################################
### chunk number 92: 
###################################################
round(PRESS(car.glm),4)# <- sum((resid(car.glm)^2/(1-hatvalues(car.glm))^2)))
round(CVE(car.glm),5)# <- sum((resid(car.glm)^2/(1-hatvalues(car.glm))^2))/
#  car.glm$df.residual)


###################################################
### chunk number 93: 
###################################################
boxcox(qmt~s+c+t+g+disp+hp+cb+drat+wt,data=car,plotit=TRUE)


###################################################
### chunk number 94: 
###################################################
boxcox(qmt~s+c+t+g+disp+hp+cb+drat+wt,data=car[rownames(car)!="MERCEDES 230",])


###################################################
### chunk number 95: 
###################################################
car <- transform(car, sp=900/qmt)
summary.lm(glm(sp~s+c+t+g+disp+hp+cb+drat+wt,data=car,subset=rownames(car)!="MERCEDES 230"))


###################################################
### chunk number 96: 
###################################################
car <- transform(car, g2 = g^2)
summary.lm(glm(sp~s+c+t+g+disp+hp+cb+drat+wt+g2,data=car,subset=rownames(car)!="MERCEDES 230"))


###################################################
### chunk number 97: 
###################################################
summary.lm(glm(sp~s+hp+wt+g,data=car,subset=rownames(car)!="MERCEDES 230"))


###################################################
### chunk number 98: 
###################################################
#car$lhp <- log(car$hp)
#car$lwt <- log(car$wt)
#car$ldis <- log(car$disp)
car <- transform(car, lhp=log(hp),lwt=log(wt),ldis=log(disp))
boxcox(sp~s+c+t+g+ldis+lhp+cb+drat+lwt,data=car[-9,],plotit=FALSE)


###################################################
### chunk number 99: 
###################################################
summary.lm(glm(qmt~s+g+lhp+lwt,data=car,subset=rownames(car)!="MERCEDES 230"))


###################################################
### chunk number 100: 
###################################################
#solv <- read.csv('Data/solve.csv')
options(digits=8)
solv.glm <- glm(time~eft,data=solv)
coef(summary(solv.glm))
summary(solv.glm)$dispersion
summary(solv.glm)$dispersion*summary(solv.glm)$cov.unscaled
qt(0.975,solv.glm$df.residual)


###################################################
### chunk number 101: 
###################################################
predictx <-  function(model,y,alpha=0.05){
    a <- coefficients(model)[1]
    b <- coefficients(model)[2]
s2 <- if (class(model)=="lm"){ 
summary(model)$sigma^2} else {
if ((class(model)[1]=="glm")&(family(model)$family=="gaussian")){
summary(model)$dispersion} else {stop("only available for normal models")}}
    V <- summary(model)$cov.unscaled*s2
    v11 <- V[1,1]
    v12 <- V[1,2]
    v22 <- V[2,2]
    df <- model$df.residual
    c <- qt(1-alpha/2,df)^2
    z3 <- (b^2-c*v22)
    z2 <- (b*(y-a)+c*v12)
    z1 <- ((y-a)^2-c*(s2+v11))
    z2 <- z2/z3*2
    z1 <- z1/z3
    z <- -as.real(polyroot(c(z1,z2,1)))
    z <- c(z[1],((y-a)/b),z[2])
    names(z) <- c("ll","x","ul")
    z
  }


###################################################
### chunk number 102: 
###################################################
predictx(lm(time~eft,data=solv),y=400)


###################################################
### chunk number 103: 
###################################################
poison <- read.csv('Data/poison.csv')
poison$type <- factor(poison$type,label=c("I","II","III"))
poison$treat <- factor(poison$treat,label=LETTERS[1:4])


###################################################
### chunk number 104:  eval=FALSE
###################################################
## data(poison,package="SMIR")


###################################################
### chunk number 105: poison.glm
###################################################
poison.glm1 <- glm(time~1,data=poison)
poison.glm2 <- update(poison.glm1,.~.+type)
poison.glm3 <- update(poison.glm2,.~.+treat)
poison.glm4 <- update(poison.glm3,.~.+type:treat)
round(coef(summary(poison.glm4)),4)


###################################################
### chunk number 106: poison.ANOVA
###################################################
library(xtable)
xtable(anova(poison.glm4,test="F"),
caption="ANOVA table, time.", label="tab:poison.anova",caption.placement="top")


###################################################
### chunk number 107: 
###################################################
round(with(poison,tapply(time,list(type=type,treatment=treat),mean)),4)


###################################################
### chunk number 108: 
###################################################
round(with(poison,tapply(time,treat,mean)),4)


###################################################
### chunk number 109: 
###################################################
round(tvar <- with(poison,tapply(time,list(type=type,treatment=treat),var)),7)


###################################################
### chunk number 110: log.quant.log.var eval=FALSE
###################################################
## tvar <- as.vector(tvar)
## q <- qchisq(rank(tvar)/13,3)
## print(xyplot(q~tvar,ylab='quantile (log scale)',
## xlab="variance (log scale)",ylim=c(0.02,20),
## scales=list(y=list(log="e"),x=list(log="e")),
## panel=function(...){
## panel.xyplot(...)
## panel.curve(x+log(3)-log(0.0222))}))


###################################################
### chunk number 111: 
###################################################
tvar <- as.vector(tvar)
q <- qchisq(rank(tvar)/13,3)
print(xyplot(q~tvar,ylab='quantile (log scale)',
xlab="variance (log scale)",ylim=c(0.02,20),
scales=list(y=list(log="e"),x=list(log="e")),
panel=function(...){
panel.xyplot(...)
panel.curve(x+log(3)-log(0.0222))}))


###################################################
### chunk number 112: 
###################################################
poison.bc <- boxcox(poison.glm4,plotit=FALSE,lambda=seq(-1.5,-0.2,by=0.1))
names(poison.bc$y) <- poison.bc$x
round(poison.bc$y,2)


###################################################
### chunk number 113: 
###################################################
poison.rate.glm1 <- glm(1/time~1,data=poison)
poison.rate.glm2 <- update(poison.rate.glm1,.~.+type)
poison.rate.glm3 <- update(poison.rate.glm2,.~.+treat)
poison.rate.glm4 <- update(poison.rate.glm3,.~.+type:treat)


###################################################
### chunk number 114: 
###################################################
round(anova(poison.rate.glm4),4)


###################################################
### chunk number 115: 
###################################################
round(with(poison,tapply(1/time,list(type=type,treatment=treat),mean)),3)
round(tvar <- with(poison,tapply(1/time,list(type=type,treatment=treat),var)),5)


###################################################
### chunk number 116: poisson.recip.quantile eval=FALSE
###################################################
## tvar <- as.vector(tvar)
## q <- qchisq(rank(tvar)/13,3)
## print(xyplot(log(q)~log(tvar),ylab='quantile (log scale)',
## xlab="variance (log scale)",
## panel=function(...){
## panel.xyplot(...)
## panel.curve(x+log(3)-log(0.240))}))


###################################################
### chunk number 117: 
###################################################
tvar <- as.vector(tvar)
q <- qchisq(rank(tvar)/13,3)
print(xyplot(log(q)~log(tvar),ylab='quantile (log scale)',
xlab="variance (log scale)",
panel=function(...){
panel.xyplot(...)
panel.curve(x+log(3)-log(0.240))}))


###################################################
### chunk number 118:  eval=FALSE
###################################################
## data(hostility,package="SMIR")


###################################################
### chunk number 119: 
###################################################
hostility <- read.csv('Data/hostility.csv')
hostility$g <- factor(hostility$group,
                          labels=c('overdoses','F controls','T controls'))
hostility$nation <- factor(hostility$nationality,
                                labels=c('Australian','British'))
hostility$po <- factor(hostility$po,labels=c('none','previous'))


###################################################
### chunk number 120: hostility.table
###################################################
ftable(xtabs(~g+po+nation,data=hostility),
       row.vars=c(1),col.vars=c(3,2))


###################################################
### chunk number 121: 
###################################################
host.glm1 <- glm(positivity~g*nation*po,data=hostility)
host.glm2 <- glm(positivity~po*g*nation,data=hostility)
host.glm3 <- glm(positivity~nation*po*g,data=hostility)


###################################################
### chunk number 122: 
###################################################
#options(digits=4)
host.glm4 <- glm(positivity~(g+nation)^2+po,data=hostility)
print(summary(host.glm4),digits=3)


###################################################
### chunk number 123: host.glm5
###################################################
#options(digits=4)
hostility$g23 <- factor(hostility$group)
levels(hostility$g23) <- c('overdoses','control','control')
host.glm5 <- update(host.glm4,.~.-g:nation+g23+g23:nation)
print(summary(host.glm5),digits=4)


###################################################
### chunk number 124: host.glm6
###################################################
#options(digits=4)
host.glm6 <- update(host.glm5,.~.-g)
round(coef(summary(host.glm6)),4)


###################################################
### chunk number 125: host.glm7
###################################################
#options(digits=4)
hostility$nation.for.controls <- with(hostility,
(g23=="control")&(nation=="British"))
host.glm7 <- update(host.glm6,.~po+g23+nation.for.controls)
round(coef(summary(host.glm7)),4)


###################################################
### chunk number 126: 
###################################################
anova(host.glm7,host.glm1,test="F")


###################################################
### chunk number 127: host.glm8
###################################################
#options(digits=4)
hostility$comp <- with(hostility,
(g23=="control")+((nation=="British")*(g23=="control"))-(po=='previous'))
host.glm8 <- update(host.glm7,.~comp,data=hostility)
summary(host.glm8)


###################################################
### chunk number 128:  eval=FALSE
###################################################
## round(ftable(with(hostility,tapply(fitted(host.glm8),
##                              list(g,po,nation),mean)),
##        row.vars=c(1),col.vars=c(3,2)),2)
## round(ftable(with(hostility,tapply(positivity,
##                              list(g,po,nation),mean)),
##        row.vars=c(1),col.vars=c(3,2)),2)


###################################################
### chunk number 129: 
###################################################
round(ftable(with(hostility,tapply(fitted(host.glm8),
                             list(g,po,nation),mean)),
       row.vars=c(1),col.vars=c(3,2)),2)
round(ftable(with(hostility,tapply(positivity,
                             list(g,po,nation),mean)),
       row.vars=c(1),col.vars=c(3,2)),2)


###################################################
### chunk number 130: 
###################################################
round(coef(summary(host.glm1))[,1:2],4)


###################################################
### chunk number 131: 
###################################################
host.glm2 <- update(host.glm1,.~.-nation:po:g)
anova(host.glm2,host.glm1,test="F")
round(coef(summary(host.glm2))[,1:3],4)


###################################################
### chunk number 132: 
###################################################
host.glm3 <- update(host.glm2,.~.-g:po-nation:po)
anova(host.glm3,host.glm2,test="F")
round(coef(summary(host.glm3))[,1:3],4)


###################################################
### chunk number 133: 
###################################################
round(with(hostility,tapply(positivity,list(g=g,po=po),mean,na.rm=TRUE)),2)


###################################################
### chunk number 134:  eval=FALSE
###################################################
## data(poison,package="SMIR")


###################################################
### chunk number 135: 
###################################################
library(dglm)
poison$rate <- 1/poison$time
poison.dglm <- dglm(rate~type*treat,dformula=~type*treat,data=poison)


###################################################
### chunk number 136: 
###################################################
print(summary(poison.dglm,dispersion=1),digits=4)


###################################################
### chunk number 137: 
###################################################
poison.dglm2 <- dglm(rate~type*treat,dformula=~type+treat,data=poison)


###################################################
### chunk number 138: poison.dglm2
###################################################
print(summary(poison.dglm2,dispersion=1),digits=4)


###################################################
### chunk number 139: 
###################################################
poison.dglm3 <- dglm(rate~type*treat,dformula=~type,data=poison)


###################################################
### chunk number 140: poison.dglm3
###################################################
summary(poison.dglm3$dispersion)


###################################################
### chunk number 141: poison.dmod4
###################################################
poison.dglm4 <- dglm(rate~type*treat,dformula=~1,data=poison)
summary(poison.dglm4$dispersion)


###################################################
### chunk number 142: 
###################################################
poison$type2 <- ifelse(poison$type=="II",1,0)
poison.dglm5 <- dglm(rate~type*treat,dformula=~type2,data=poison)


###################################################
### chunk number 143: poison.dmod5
###################################################
summary(poison.dglm5$dispersion)
#deviance(poison.dmod5)


###################################################
### chunk number 144: 
###################################################
poison.dglm6 <- dglm(rate~type+treat,data=poison,dformula=~type2)


###################################################
### chunk number 145: poison.dmod6
###################################################
print(summary(poison.dglm6),digits=4)


###################################################
### chunk number 146: 
###################################################
poison$d34 <- (poison$type=="III")*(poison$treat=="D")
poison.dglm7 <- dglm(rate~type+treat+d34,data=poison,dformula=~type2)
#poison.dmod7 <- lm.disp(rate~type+treat+d34,data=poison,
#                        var.formula~type2)


###################################################
### chunk number 147: poison.dglm7
###################################################
print(summary(poison.dglm7,disp=1),digits=4)


###################################################
### chunk number 148: fc2
###################################################
fc2 <- function(x) formatC(x,digits=2,format="f")


###################################################
### chunk number 149: 
###################################################
#newpoison <- data.frame(type=rep(c(3,2,1),4),treat=rep(c(1,3,4,2),c(3,3,3,3)),
#d34=c(rep(0,6),1,rep(0,5)))
#predict(poison.dmod7$mean,new=newpoison)
poison.fitted <- predict(poison.dglm7)
tab.means <- with(poison,tapply(rate,list(Type=type,Treat=treat),
mean,na.rm=TRUE))
tab.fitted <- tapply(poison.fitted,
list(Type=poison$type,Treat=poison$treat),mean,na.rm=TRUE)


###################################################
### chunk number 150: 
###################################################
fc3 <- function(x)formatC(x,digits=3,format="f")


###################################################
### chunk number 151: 
###################################################
poison.var.fitted <- predict(poison.dglm7$dispersion.fit)
tab.var <- with(poison,tapply(rate,list(Type=type,Treat=treat),
var,na.rm=TRUE))
tab.var.fitted <- tapply(poison.var.fitted,
list(Type=poison$type,Treat=poison$treat),mean,na.rm=TRUE)


###################################################
### chunk number 152:  eval=FALSE
###################################################
## data(trees,package="SMIR")


###################################################
### chunk number 153: tree.dglm
###################################################
library(dglm)
trees <- transform(trees,
                   lh2 = lh^2,
                   ld2 = ld^2,
                   lhld = lh*ld)
trees.dglm <- dglm(lv~lh+ld,data=trees,
                       dformula~lh+ld+lh2+lhld+ld2)


###################################################
### chunk number 154: 
###################################################
summary(trees.dglm,disp=1)


###################################################
### chunk number 155: 
###################################################
trees.dglm2 <- dglm(lv~lh+ld,data=trees,
                    dformula=~ld+ld2)
trees$z <- 2*trees$ld+trees$lh
trees.dglm3 <- dglm(lv-z~1,data=trees,
                    dformula=~ld+ld2)
trees.dglm4 <- dglm(lv-z~1,data=trees,dformula=~1)


###################################################
### chunk number 156: 
###################################################
fc0 <- function(x) formatC(x,digits=0,format="f" )


###################################################
### chunk number 157: tree.var.plot eval=FALSE
###################################################
## plot(log(fitted(trees.dglm2$dispersion.fit))~log(d),data=trees,type="l",las=1,ylim=c(-13,-3),
## xlab="log d",ylab="fitted log variance")
## points(log(trees$d),log(trees.dglm2$dispersion.fit$y))


###################################################
### chunk number 158: 
###################################################
plot(log(fitted(trees.dglm2$dispersion.fit))~log(d),data=trees,type="l",las=1,ylim=c(-13,-3),
xlab="log d",ylab="fitted log variance")
points(log(trees$d),log(trees.dglm2$dispersion.fit$y))


###################################################
### chunk number 159: 
###################################################
plot(fitted(trees.dglm2$dispersion.fit)~log(d),data=trees,type="l",las=1,ylim=c(0,0.030),
xlab="log d",ylab="fitted variance")
points(log(trees$d),trees.dglm2$dispersion.fit$y)
#plot(trees$ld,log(fitted(trees.dmod2$var)),xlab=expression("log",d),ylab="fitted variance",
#type="l")
#points(trees$ld,log(resid(trees.dmod2$mean)^2))


###################################################
### chunk number 160: 
###################################################
trees.dglm2 <- dglm(lv~lh+ld,data=trees,dformula=~ld+ld2)


###################################################
### chunk number 161: tree_quadratic_influence eval=FALSE
###################################################
## plot(trees$ld,hatvalues(trees.dglm2$dispersion.fit),xlab='log d',
##      ylab='variance influence',las=1)


###################################################
### chunk number 162: 
###################################################
plot(trees$ld,hatvalues(trees.dglm2$dispersion.fit),xlab='log d',
     ylab='variance influence',las=1)