###########################################
# QUESTION 1
###########################################
data(Aids2)
attach(Aids2)
time<-death-diag+1
c<-codes(status)-1
rate<-c/time
summary(glm(rate~state+sex+diag+T.categ+age, family=poisson, weight=time))
glm(rate~state+sex+diag+age, family=poisson, weight=time)
glm(rate~sex+diag+T.categ+age, family=poisson, weight=time)
glm(rate~state+diag+T.categ+age, family=poisson, weight=time)
glm(rate~diag+T.categ+age, family=poisson, weight=time)
glm(rate~sex+diag+age, family=poisson, weight=time)
glm(rate~state+diag+age, family=poisson, weight=time)
summary(glm(rate~diag+age, family=poisson, weight=time))

###########################################
# QUESTION 2
###########################################
data(biopsy)
attach(biopsy)
Y<-codes(class)-1
fV1<-factor(V1)

par(mfrow=c(2,2))
glm0<-glm(Y~V1)
summary(glm0)
glm1<-glm(Y~fV1)
summary(glm1)
plot(V1,fitted(glm1))
points(V1,fitted(glm0),pch=2)
title('Figure 2.1: Normal family')

glm0<-glm(Y~V1,family=binomial)
summary(glm0)
glm1<-glm(Y~fV1, family=binomial)
summary(glm1)
plot(V1,fitted(glm1))
points(V1,fitted(glm0),pch=2)
title('Figure 2.2: Binomial family')

summary(glm(Y~V1+V2+V3+V4+V5+V6+V7+V8+V9, family=binomial))
summary(glm(Y~V1+V4+V6+V7, family=binomial))

###########################################
# QUESTION 3
###########################################

freq<-c(scan("C:/keith/teaching/datasets/happy"))
(t(matrix(freq,5,12)))
happ<-gl(3,20,60)
years<-gl(4,5,60)
sibs<-gl(5,1,60)
glm(freq~happ+years+sibs, family=poisson)
glm(freq~happ+years+sibs+happ:years, family=poisson)
glm(freq~happ+years+sibs+happ:sibs, family=poisson)
glm(freq~happ+years+sibs+years:sibs, family=poisson)
glm(freq~happ+years+sibs+happ:years+happ:sibs, family=poisson)
glm(freq~happ+years+sibs+happ:years+years:sibs, family=poisson)
glm(freq~happ+years+sibs+happ:sibs+years:sibs, family=poisson)
glm(freq~happ+years+sibs+happ:years+happ:sibs+years:sibs, family=poisson)
xyears<-codes(years)
xsibs<-codes(sibs)
summary(glm(freq~happ+years+sibs+happ:xyears+happ:xsibs+years:sibs, family=poisson))
glm(freq~happ+years+sibs+happ:xsibs+years:sibs, family=poisson)
glm(freq~happ+years+sibs+happ:xyears+years:sibs, family=poisson)
glm(freq~happ+years+sibs+years:sibs, family=poisson)

###########################################
# QUESTION 4
###########################################

m<-t(matrix(scan("C:/keith/teaching/datasets/steel"),2,40))
stress<-m[,1]
time<-m[,2]
ltime<-log(time)
plot(stress,ltime)
title('Figure 4.1: Time vs. stress')
summary(glm(time~stress, family=Gamma(link="log")))
fstress<-factor(stress)
summary(glm(time~stress+fstress, family=Gamma(link="log")))

###########################################
# QUESTION 5
###########################################

m<-t(matrix(scan("C:/keith/teaching/datasets/turbines"),4,45))
loc<-factor(m[,1])
temp<-m[,2]
stren<-m[,3]
grow<-m[,4]
plot(temp,grow)
title('Figure 5.1: grow vs. temp')
lgrow<-log(grow)
plot(temp,lgrow)
title('Figure 5.2: log(grow) vs. temp')
summary(glm(lgrow~loc+temp+stren))
summary(glm(lgrow~loc+temp+stren+loc:temp))
summary(glm(lgrow~loc+temp+stren+loc:temp+loc:stren))
plot(codes(loc),temp)
title('Figure 5.3: temp vs. loc')
plot(codes(loc),stren)
title('Figure 5.4: stren vs. loc')
glm0<-glm(lgrow~loc+temp+stren+loc:temp)
plot(fitted(glm0),resid(glm0))
title('Figure 5.5: fitted vs. resid')

vl<-predict.glm(glm0,se.fit=T)$se.fit^2
sc<-summary(glm0)$dispersion
r<-resid(glm0)
z<-r/sqrt(sc-vl)
df<-glm0$df.residual
tstat<-z/sqrt((df-z^2)/(df-1))
cbind(r,z,tstat)
z[7]=0
u<-pnorm(z)
i<-1:45
uhat<-(i-0.5)/45
u<-sort(u)
u-uhat

###########################################
# QUESTION 6
###########################################

y<-(-1000:1000)/1000*10
mu<-5
s<-0
for(alpha in 1:5) {
phi=1/alpha
j<-0:1000
p<-rowSums(log(1+outer(y^2,1/(1+2*phi*j)^2)))
f0<-2^(alpha-2)*exp(-p-lgamma(alpha))*alpha
# plot(y,f0,type="l")
s=c(s,sum(exp(-p))/100)
}
s<-s[-1]
alpha<-1:5
s*2^(alpha-2)
sqrt(pi)








#