# Leser inn dataene (jf exercise 15 i oppgaveheftet, men bruk extension "txt")
car=read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/h11/data/car.txt",header=T,sep=",")

# En oversikt over dataene er gitt i Tabell 1.6, og en nøyaktig beskrivelse av kodingen av dataene er gitt på bokas web-side: http://www.afas.mq.edu.au/acst_docs/glms_for_insurance_data/data/vehicle.txt

# Vi vil lage en modell for hvordan antall skader (numclaims) avhenger av bilens verdi (veh_value), type bil (veh_body), bilens alder (veh_age), kjønn til føreren (gender), område der føreren bor (area) og alderen til føreren (agecat)

# Vi bruker "SEDAN" som referanse for biltype
car$veh_body=relevel(car$veh_body,ref="SEDAN")

# Lager en trunkert variabel for bilens verdi (som er veldig skjevt fordelt -- bare vel 2% har verdi over 50 000)
car$veh_value_trunc=car$veh_value*(car$veh_value<=5)+5*(car$veh_value>5)

# Tilpasser modell med alle variablene:

fit.alle=glm(numclaims~offset(log(exposure))+veh_value_trunc+veh_body+factor(veh_age)+gender+area+factor(agecat),family=poisson, data=car)


# Kutter ut variable ved baklengs eliminasjon. 

# Kjønn er minst signifikant og kuttes ut først:
fit.uten.gender=glm(numclaims~offset(log(exposure))+veh_value_trunc+veh_body+factor(veh_age)+area+factor(agecat),family=poisson, data=car)
anova(fit.uten.gender,fit.alle,test="Chisq")

# Kutter så ut bilens alder:
fit.uten.gender.vehage=glm(numclaims~offset(log(exposure))+veh_value_trunc+veh_body+area+factor(agecat),family=poisson, data=car)
anova(fit.uten.gender.vehage,fit.uten.gender,test="Chisq")

# Kutte endelig ut område:
fit.uten.gender.vehage.area=glm(numclaims~offset(log(exposure))+veh_value_trunc+veh_body+factor(agecat),family=poisson, data=car)
anova(fit.uten.gender.vehage.area,fit.uten.gender.vehage,test="Chisq")

# Nå er alle variable signifikante, så vi kutter ikke ut flere.
# Vi burde også ha sjekket for interaksjoner, men vi stopper denne illustrasjonen her.


# Endelig modell:
summary(fit.uten.gender.vehage.area)


# 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)
}

round(RRCItab(fit.uten.gender.vehage.area),3)