## Teht3
miners <- read.table("http://users.jyu.fi/~junyblom/miners.dat",
header=T)
attach(miners)
age<-factor(age)
malli30<-glm(counts~age+wheeze+breath,family=poisson("log"))
summary(malli30)
malli31<-glm(counts~age*wheeze*breath,family=poisson("log"))
summary(malli31)
anova(malli30,malli31,test="Chisq")
1-pchisq(6939.1,df=25)

malli32<-glm(counts~age*breath+age*wheeze+wheeze*breath,family=poisson("log"))
summary(malli32)
anova(malli32,malli31,test="Chisq")
1-pchisq(26.69,df=8)
# Interaktiot merkitseviä

summary(malli31)
kertoimet<-malli31$coef

# wheeze-breath yhdysvaikutukset eri ikäisillä
interak<-kertoimet[28]+c("age1:wheeze:breath"=0,kertoimet[29:36])
interak

OR<-exp(interak)
OR
# Esim, jos wheeze saa arvon 1, niin vedonlyöntisuhde hengästyneisyydelle on
# 25-kertainen verrattuna riskiin, jos wheeze on 0 20-24 ikäisillä

plot(1:9,OR,xlab="ikä",type="p")
# Riippuvuus on pienempää, kun ikä on korkeampi

## teht 4

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

ale<-c(rep(0,times=6),10/12,rep(1,times=3))
aika<-c(rep(0:9))
tulo<-ale*aika

malli<-glm(kuolleet~aika+ale+tulo,
	family=poisson,offset=log(vaki/100000))
summary(malli)
exp(coef(malli))
#Vuonna 98 malli ennustaa 28 kuollutta/100 000 kohden
#Mallin mukaan vuotuinen kuolleisuuden kasvu on 0.6% vuotta kohden
#ennen vuotta 03.
#Malli ennustaa 7.7%kasvun kuolleisuudelle hinnan laskun seurauksena

exp(0.005727  +    0.025877)
#kasvu kuolleisuudessa mallin mukaan on 3.2% vuotta kohden 
#vuodesta 04 alkaen
#a)
1-pchisq(6.3154,6) 
malli riittävä

#b)
#Veron vähennys vaikuttaa sekä tason että kulmakertoimen muutoksena, 
#koska puhutaan suhteellisesta muutoksesta 

#c)
b <- coef(malli)
exp(b[2])              # kasvuvauhti v. 2003 on 0.57 %
exp(b[2] + b[4]*10/12) # kasvuvauhti v. 2004 on 2.8 %
exp(b[2]+b[4])         # kasvuvauhti v. 2007 on 3.2 %

c(exp(confint(malli)[2,1])^aika[6],exp(confint(malli)[2,2])^aika[6])
c(exp(confint(malli)[2,1])^aika[7],exp(confint(malli)[2,2])^aika[7])
c(exp(confint(malli)[2,1])^aika[10],exp(confint(malli)[2,2])^aika[10])

d)
plot(vuosi, 1e5*kuolleet/vaki, ylab="a-kuolleet/100 000 as. kohti",
xlab="Vuosi", main="")
lines(vuosi, 1e5*malli$fit/vaki, type="l")