################# ESIMERKKI 7.1 ###################################

 a.kuolleet <-
 read.table("http://users.jyu.fi/~junyblom/a-kuolleet.dat",
 header=TRUE)

 attach(a.kuolleet)

 x <- vuosi-1969

 out <- glm(kuolleet~x, family=poisson, offset=log(1e-5*vakiluku))
 coef(out)

 plot(vuosi, 1e5*kuolleet/vakiluku, xlab="Vuosi",
 ylab="Kuolleet/100000 as.")
 curve(exp(coef(out)[1] + coef(out)[2]*(x-1969)), from=1969, to=2005, 
 add = TRUE)

################# ESIMERKKI 7.2 ###################################

doctors <- read.table("http://users.jyu.fi/~junyblom/doctors.dat",
header=T)

attach(doctors)

doctors.glm <- glm(deaths~age*smoker, offset=log(pyears*1e-5), family=poisson)
summary(doctors.glm)
b <- coef(doctors.glm)
b
doctors2.glm <- glm(deaths~age*smoker+ I(age^2), offset=log(pyears*1e-5), 
family=poisson)
summary(doctors2.glm)
b2 <- coef(doctors2.glm)
b2

oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,3,0)) 

 plot(0:4,(deaths/pyears)[smoker=="no"]*1e5,
 xlab=" ", ylab=" ", main =" ")
 points(0:4,(deaths/pyears)[smoker=="yes"]*1e5, ylim=c(1,3000),
 pch=20)
 curve(exp(b[1] + b[2]*x), from=0, to=4, add=TRUE)
 curve(exp(b[1] + b[3]+ (b[2]+b[4])*x), from=0, to=4, lty=2, add=TRUE)

 mtext("Ikäryhmä", font=1, side=1, line=2.5) 
 mtext("Kuolemat 100000 henkilövuotta kohti", font=1, side=2, line=2.5)
 mtext("Ikä lineaarinen", font=1, side=3, line=1.5)

 par(mar=c(3,4,3,0))

 plot(0:4,(deaths/pyears)[smoker=="no"]*1e5,
 xlab=" ", ylab=" ", main =" ")
 points(0:4,(deaths/pyears)[smoker=="yes"]*1e5, ylim=c(1,3000),
 pch=20)
 curve(exp(b2[1]+ b2[2]*x + b2[4]*x^2), 
 from=0, to=4, lty=2, add=TRUE)
 curve(exp(b2[1] + b2[3]+ (b2[2]+b2[5])*x + b2[4]*x^2), 
 from=0, to=4, lty=2, add=TRUE)


 mtext("Ikäryhmä", font=1, side=1, line=2.5) 
 mtext("Kuolemat 100000 henkilövuotta kohti", font=1, side=2, line=2.5)
 mtext("Ikä kvadraattinen", font=1, side=3, line=1.5)

 par(oldpar)

 

############## ESIMERKKI 7.3 #############################
 
 summary(doctors2.glm)

 exp(Smo.d <- b2[3] + b2[5]*(0:4))
 X <- cbind(1,0:4)
 covm <- summary(doctors2.glm)$cov.s                       ## cov(b)
 se <- sqrt((diag(X%*%covm[c(3,5),c(3,5)]%*%t(X))))        ## keskivirheet
 round(exp(cbind(Smo.d-qnorm(.975)*se,Smo.d+qnorm(.975)*se)),2)
 
############### ESIMERKKI 7.5 ############################

 out <- glm(kuolleet~x, family=poisson, offset=log(1e-5*vakiluku))

 summary(out)

 # Pearsonin residuaalit

 res <- (kuolleet - fitted(out))/sqrt(fitted(out))

 # Lag-plot

 n <- length(kuolleet)

 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,3,0)) 

 plot(fitted(out), res, xlab=" ", ylab=" ")

 mtext("Sovite", font=1, side=1, line=2.5) 
 mtext("Jäännös", font=1, side=2, line=2.5)

 par(mar=c(3,4,3,0))

  plot(res[-n], res[-1],, xlab=" ", ylab=" ")

 mtext("Edellinen", font=1, side=1, line=2.5) 
 mtext("Jälkimmäinen", font=1, side=2, line=2.5)

 par(oldpar)

 a <- acf(res)
 a
############## ESIMERKKI 7.7 #################################

aspirin <- read.table("http://users.jyu.fi/~junyblom/aspirin.dat",
header=T)

attach(aspirin)

 aspirin 

 # 2x2x2 taulu
  A <- tapply(freq,list(type,user,ulcer),sum)
  A

 # Kiinteät reunasummat
 A[,1,]+A[,2,]

 type <- relevel(type,ref="control")
 aspirin.glm <- glm(freq~user*ulcer*type, family=poisson)
 summary(aspirin.glm)

# Varsinaiset mahahaavapotilaat
 
 gastric <- glm(freq~user * type, subset=ulcer=="gastric",
 family=poisson)
 summary(gastric)

# Pohjukaissuolipotilaat

 duodenal <- glm(freq~user * type, subset=ulcer=="duodenal",
 family=poisson)
 summary(duodenal)