###################################################
### chunk number 1: nls1
###################################################
#line 159 "Lect12.Rnw"
library(stats)
time<-c(0,1,2,3,4,7,9,11,14,16,18,21,24,29,32,35,38,46)
count<-c(5149,5435,5358,5462,4755,4457,4538,4378,4189,4195,3743,3855,3569,3045,3091,3064,2729,2477)
dataf<-data.frame(cbind(count,time))
names(dataf)<-c("count","time")
# create the data set


###################################################
### chunk number 2: nls2
###################################################
#line 168 "Lect12.Rnw"
par(mfrow=c(1,1))
plot(dataf$time,dataf$count,xlab="time",ylab="count")
# the data, slight indication of curvature
# only slight because we only have one observation after the half-life
# of 40 days


###################################################
### chunk number 3: nls3
###################################################
#line 185 "Lect12.Rnw"
m1<-nls(count~a+b*exp(-cc*time),data=dataf,start=list(a=2000,b=3000,cc=.02))
m1s<-summary(m1)
print(m1s)
# estimation, model summary


###################################################
### chunk number 4: nls3b
###################################################
#line 195 "Lect12.Rnw"
par(mfrow=c(2,2))
ma<-m1s$p[1,1]+m1s$p[2,1]*exp(-m1s$p[3,1]*dataf$time)
plot(dataf$time,dataf$count,xlab="time",ylab="count",main="data")
lines(dataf$time,ma)
qqnorm(m1s$res)
qqline(m1s$res)
plot(ma,m1s$res,main="residuals on fitted")
plot(dataf$time,m1s$res,main="residuals on time")
# residual diagnostics


###################################################
### chunk number 5: nls4
###################################################
#line 218 "Lect12.Rnw"
m2<-nls(count~b*exp(-cc*time),data=dataf,start=list(b=3000,cc=.05))
m2s<-summary(m2)
print(m2s)
# hypothesis testing


###################################################
### chunk number 6: nls5
###################################################
#line 225 "Lect12.Rnw"
print(anova(m1,m2))
# F-test


###################################################
### chunk number 7: nsl6
###################################################
#line 232 "Lect12.Rnw"
m3<-nls(count~exp(a)+exp(b)*exp(-cc*time),data=dataf,start=list(a=log(2000),b=log(3000),cc=.05))
m3s<-summary(m3)
print(m3s)
# constrained NLS


###################################################
### chunk number 8: nls7
###################################################
#line 242 "Lect12.Rnw"
pf<-profile(m1,alphamax=.05)
par(mfrow=c(1,3))
plot(pf)
# profile functions


###################################################
### chunk number 9: nls7d
###################################################
#line 249 "Lect12.Rnw"
print(confint(m1))
# profile based confidence intervals


###################################################
### chunk number 10: glm1
###################################################
#line 314 "Lect12.Rnw"
library(stats)
library(MASS)
#
time<-c(0,1,2,3,4,7,9,11,14,16,18,21,24,29,32,35,38,46)
count<-c(5149,5435,5358,5462,4755,4457,4538,4378,4189,4195,3743,3855,3569,3045,3091,3064,2729,2477)
count2<-c(5121,5263,4930,5266,4521,4434,4590,4290,4005,4123,3659,3785,3538,3159,2946,2956,2643,2500)
dataf<-data.frame(cbind(count,count2,time))
names(dataf)<-c("count","count2","time")


###################################################
### chunk number 11: glm2
###################################################
#line 325 "Lect12.Rnw"
gg<-glm(count~time,data=dataf,family=poisson)
gs<-summary(gg)
print(gs)


###################################################
### chunk number 12: glm2b
###################################################
#line 330 "Lect12.Rnw"
plot(dataf$time,dataf$count)
lines(dataf$time,gg$fit,col="green")


###################################################
### chunk number 13: glm4
###################################################
#line 349 "Lect12.Rnw"
c("1-P(X2>Dev)=",1-pchisq(gg$dev,df=gs$df[2]))


###################################################
### chunk number 14: glm4b
###################################################
#line 354 "Lect12.Rnw"
par(mfrow=c(2,2))
plot(gg)


###################################################
### chunk number 15: glm5
###################################################
#line 369 "Lect12.Rnw"
c("Z-interval")
c("(Intercept)",round(gs$coef[1,1]-1.96*gs$coef[1,2],8), round(gs$coef[1,1]+1.96*gs$coef[1,2],8))
c("(time)",round(gs$coef[2,1]-1.96*gs$coef[2,2],8), round(gs$coef[2,1]+1.96*gs$coef[2,2],8))


###################################################
### chunk number 16: glm6
###################################################
#line 375 "Lect12.Rnw"
plot(profile(gg))


###################################################
### chunk number 17: glm7
###################################################
#line 388 "Lect12.Rnw"
print(confint(gg))


###################################################
### chunk number 18: glm8
###################################################
#line 393 "Lect12.Rnw"
c("phi-hat=",round(sum(residuals(gg,'pearson')^2)/gg$df.r,5))


###################################################
### chunk number 19: glmc1
###################################################
#line 399 "Lect12.Rnw"
A<-as.factor(c(rep(0,length(time)),rep(1,length(time))))
Time<-rep(time,2)
Count<-c(count,count2)
Dataf<-data.frame(cbind(Count,Time,A))
names(Dataf)<-c("Count","Time","A")
gg2<-glm(Count~A*Time,data=Dataf,family=poisson)
gs2<-summary(gg2)
par(mfrow=c(1,1))
plot(Time,Count)
points(time,count,col="red")
points(time,count2,col="green")
lines(Time[A==0],gg2$fit[A==0],col="red")
lines(Time[A==1],gg2$fit[A==1],col="green")


###################################################
### chunk number 20: glmc1b
###################################################
#line 414 "Lect12.Rnw"
print(gs2)


###################################################
### chunk number 21: glmc2
###################################################
#line 431 "Lect12.Rnw"
c("Dev/df=",round(gg2$dev/gs2$df[2],5))


###################################################
### chunk number 22: glmc4
###################################################
#line 437 "Lect12.Rnw"
print(confint(gg2))


###################################################
### chunk number 23: glmc5
###################################################
#line 442 "Lect12.Rnw"
gg2b<-update(gg2,.~.-A:Time)
anova(gg2b,gg2)


###################################################
### chunk number 24: glmc5b
###################################################
#line 447 "Lect12.Rnw"
c("Pvalue", round((1-pchisq(2.131,1)),6))


###################################################
### chunk number 25: glmc6
###################################################
#line 453 "Lect12.Rnw"
gg2c<-update(gg2,.~.-A-A:Time)
anova(gg2c,gg2)
c("Pvalue", round((1-pchisq(23.462,1)),6))


###################################################
### chunk number 26: glmc7
###################################################
#line 461 "Lect12.Rnw"
FF<-(gg2c$dev-gg2$dev)/(2*(gg2$dev/gs2$df[2]))
c("F-statistic=",round(FF,6))
c("Pvalue", round((1-pf(FF,2,gs2$df[2])),6))


###################################################
### chunk number 27: glmc8
###################################################
#line 469 "Lect12.Rnw"
gg3<-glm.nb(Count~A*Time,data=Dataf)
gs3<-summary(gg3)
print(gs3)


###################################################
### chunk number 28: glmc9
###################################################
#line 476 "Lect12.Rnw"
print(confint(gg3))


###################################################
### chunk number 29: glmc10
###################################################
#line 480 "Lect12.Rnw"
gg3c<-update(gg3,.~.-A-A:Time)
anova(gg3,gg3c)


###################################################
### chunk number 30: glm11
###################################################
#line 487 "Lect12.Rnw"
step(gg2)


###################################################
### chunk number 31: glm12
###################################################
#line 491 "Lect12.Rnw"
step(gg3)


###################################################
### chunk number 32: chd1
###################################################
#line 499 "Lect12.Rnw"
SA<-data.frame(read.table('SA.dat',sep='\t',header=T))
SA$ldl<-log(SA$ldl)
SA$age<-log(SA$age)
SA$famhist<-SA$famhist-1
gg<-glm(chd~ldl+age+sbp+alcohol+alcind+tobacco+tobind+famhist+typea+adiposity+obesity,'binomial',data=SA)
print(summary(gg))


###################################################
### chunk number 33: chd2
###################################################
#line 508 "Lect12.Rnw"
sum(residuals(gg,'pearson')^2/gg$df.r)


###################################################
### chunk number 34: chd3
###################################################
#line 512 "Lect12.Rnw"
plot(profile(gg))


###################################################
### chunk number 35: chd4
###################################################
#line 525 "Lect12.Rnw"
confint(gg)


###################################################
### chunk number 36: chd5
###################################################
#line 531 "Lect12.Rnw"
frac<-2/3
ii<-sample(seq(1,dim(SA)[1]),round(dim(SA)[1]*frac))
yy<-SA[ii,10]
xx<-as.matrix(SA[ii,-10])
yyt<-SA[-ii,10]
xxt<-as.matrix(SA[-ii,-10])
#
library(leaps)
rleaps<-regsubsets(xx,yy,int=T,nbest=100,nvmax=dim(SA)[2],really.big=T,method=c("ex"))
cleaps<-summary(rleaps,matrix=T)
tt<-apply(cleaps$which,1,sum)
BICvec<-rep(0,dim(cleaps$which)[1])
AICvec<-BICvec
y<-SA[,10]
x<-as.matrix(SA[,-10])
data1<-as.data.frame(cbind(y,x))
for (zz in (1:dim(cleaps$which)[1])) {
  gg<-glm(y~x[,cleaps$which[zz,2:dim(cleaps$which)[2]]==T],'binomial',subset=ii,data=data1)
  AICvec[zz]<-gg$dev+2*tt[zz]
  BICvec[zz]<-gg$dev+tt[zz]*log(length(yy)) }


###################################################
### chunk number 37: chd5b
###################################################
#line 553 "Lect12.Rnw"
plot(tt,AICvec,xlab="modelsize",ylab="AIC")


###################################################
### chunk number 38: chd5c
###################################################
#line 556 "Lect12.Rnw"
plot(tt,BICvec,xlab="modelsize",ylab="BIC")


###################################################
### chunk number 39: chd6
###################################################
#line 569 "Lect12.Rnw"
print(aicmod<-cleaps$which[AICvec==min(AICvec),])
print(bicmod<-cleaps$which[BICvec==min(BICvec),])


###################################################
### chunk number 40: chd7
###################################################
#line 574 "Lect12.Rnw"
mmaic<-glm(y~x[,aicmod[2:length(aicmod)]==T],'binomial',subset=ii)
mmbic<-glm(y~x[,bicmod[2:length(bicmod)]==T],'binomial',subset=ii)
ppaic<-predict(mmaic,as.data.frame(x[-ii,]),'response')
ppbic<-predict(mmbic,as.data.frame(x[-ii,]),'response')
print(PEaicK<-sum(yyt!=round(ppaic[-ii]))/length(yyt))
print(PEbicK<-sum(yyt!=round(ppbic[-ii]))/length(yyt))


###################################################
### chunk number 41: chd8
###################################################
#line 582 "Lect12.Rnw"
source('randomsplitsBin.R')
rr<-randomsplitsBin(SA,10,2/3,100)
print(rr$modseltab)


###################################################
### chunk number 42: chd9
###################################################
#line 587 "Lect12.Rnw"
boxplot(cbind(rr$PEa,rr$PEb),names=c("AIC","BIC"),main="Prediction errors")


###################################################
### chunk number 43: chd10
###################################################
#line 590 "Lect12.Rnw"
boxplot(cbind(rr$modsizea,rr$modsizeb),names=c("AIC","BIC"),main="Modelsizes")