# Leser inn dataene:
claims=read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/h11/data/claims.txt",header=T,colClasses=c("factor","factor","factor","numeric","numeric"))


# Punkt b)

# Modeller med alle faktorene som hovedeffekter:
fit.avd=glm(antskader~offset(log(antforsikret))+alder+motorvolum+distrikt,family=poisson,data=claims)


# Sammenligning av modellene (flere modeller har vært prøvd ut og at de som gis her svarer til baklengs eliminasjon)
fit.av=glm(antskader~offset(log(antforsikret))+alder+motorvolum,family=poisson,data=claims)
fit.v=glm(antskader~offset(log(antforsikret))+motorvolum,family=poisson,data=claims)
fit.0=glm(antskader~offset(log(antforsikret)),family=poisson,data=claims)
anova(fit.0,fit.v,fit.av,fit.avd,test="Chisq")

# Prøver også modeller med interaksjon, for eksempel


fit.avd.av=glm(antskader~offset(log(antforsikret))+alder+motorvolum+distrikt+alder:motorvolum,family=poisson,data=claims)
anova(fit.avd,fit.avd.av,test="Chisq")

# Finner ingen signifikante interaksjoner


# Punkt c)

# Definerer alder og motorvolum som numeriske kovariater. (Vi bruker verdiene 1,2,3,4 som numeriske kovariater. Alternativt kunne vi brukt midtpunktet i hvert intervall, med passende modifikasjoner for første og siste intervall.)

claims$numalder=as.numeric(claims$alder)
claims$numvolum=as.numeric(claims$motorvolum)


# Tilpasser alder som numerisk variabel og sammenligner med modellen der alder er en faktor:

fit.navd=glm(antskader~offset(log(antforsikret))+numalder+motorvolum+distrikt,family=poisson,data=claims)
anova(fit.navd,fit.avd,test="Chisq")

# Finner at alder kan modelleres numerisk.

# Tilpasser motorvolum som numerisk variabel og sammenligner med modellen der motorvolum er en faktor:

fit.anvd=glm(antskader~offset(log(antforsikret))+alder+numvolum+distrikt,family=poisson,data=claims)
anova(fit.anvd,fit.avd,test="Chisq")

# Finner at motorvolum kan modelleres numerisk.


# For distrikt er det ikke meningsfylt å modellere numerisk, men vi kan slå sammen distriktene 1-3 (fortsatt definert som faktor)

claims$nydistrikt=1*(as.numeric(claims$distrikt)<=3)+4*(claims$distrikt==4)
claims$nydistrikt=factor(claims$nydistrikt)

fit.avnd=glm(antskader~offset(log(antforsikret))+alder+motorvolum+nydistrikt,family=poisson,data=claims)
anova(fit.avnd,fit.avd,test="Chisq")

# Finner at vi kan slå sammen distriktene 1-3.


# Tilpasser modell med alle de nye variablene:
fit.nanvnd=glm(antskader~offset(log(antforsikret))+numalder+numvolum+nydistrikt,family=poisson,data=claims)


# d) 

# Plotter devians- og Pearsonresidualene mot alder (forskyer Pearsonresidualene litt mot høyre for å få bedre oversikt)

plot(claims$numalder,residuals(fit.nanvnd))
points(claims$numalder+0.03,residuals(fit.nanvnd,type="pearson"),pch=3)

# Det er her liten forskjell mellom de to typene av residualer (så vi nøyer oss med deviansresidualene heretter)

# Plotter residualene mot motorvolum og distrikt
plot(claims$numvolum,residuals(fit.nanvnd))
plot(as.numeric(claims$distrikt),residuals(fit.nanvnd)) 

# Lager normalfordelingplott av residualene (rimelig her siden vi har grupperte data)

qqnorm(residuals(fit.nanvnd))


# d) 

# Leser inn funksjon som beregner RR med konfidensintervall (jf forelesningene):
RRCItab<-function(glmfit){
sumglm<-summary(glmfit)$coef
RR<-exp(sumglm[,1])
RRL<-exp(sumglm[,1]-1.96*sumglm[,2])
RRU<-exp(sumglm[,1]+1.96*sumglm[,2])
cbind(RR,RRL,RRU)
}

# Bestemmer rate ratioene med konfidensintervall
RRCItab(fit.nanvnd)


# f)

# Rate med konfidensintervall for forsikret i alder 25-29 år, motorvolum 1.5-2 liter og bosatt i London (jf avsnitt 5.6 i læreboka)
x=matrix(c(1,2,2,1),4,1)
b=matrix(fit.nanvnd$coef,4,1)
kov=summary(fit.nanvnd)$cov.unscaled
se= sqrt(t(x)%*%kov%*%x)
rate=exp(t(x)%*%b)
lower=exp(t(x)%*%b-1.96*se)
upper= exp(t(x)%*%b+1.96*se)
c(rate, lower, upper)