# Leser inn dataene:

ctsib=read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/data/ctsib.txt",header=TRUE)

# Punkt a
# Tilpasser en vanlig logistisk regresjonsmodell

fit.a=glm(stable~Sex+Age+Height+Weight+Surface+Vision,family=binomial,data=ctsib)
summary(fit.a)


# Punkt b
# Tilpasser en logistisk regresjonsmodell med "Subject" som kategorisk kovariat (faktor)

fit.b=glm(stable~Sex+Age+Height+Weight+Surface+Vision+factor(Subject),family=binomial,data=ctsib)

# Merk warnings!

summary(fit.b)

# Merk menigsløse estimater

anova(fit.a,fit.b,test="Chisq")

# Merk at vanlig likelihood ratio testing er problematisk her (siden antall parametere øker med antall observasjoner)


# Punkt c
# Tilpasser så modell med tilfeldig effekt for "Subject".
# Det er flere mulige kommandoer her: glmmPQL, glmmML og lmer.
# glmmPQL maksimerer  en "penalized quasi-likelihood" 
# glmmML og lmer maksimerer (en tilnærmelse til) den marginale likelihooden (dvs likelihooden hvor tilfeldige effekter er integrert ut).


# glmmPQL
library(MASS)
fit.c.PQL=glmmPQL(stable~Sex+Age+Height+Weight+Surface+Vision, random=~1|Subject,family=binomial,data=ctsib)
summary(fit.c.PQL)


# Laplace tilnærmelse med glmmML
library(glmmML)
fit.c.ML.Laplace=glmmML(stable~Sex+Age+Height+Weight+Surface+Vision, cluster=Subject,family=binomial,data=ctsib)
summary(fit.c.ML.Laplace)


# Gauss-Hermite tilnærmelse med glmmML
fit.c.ML.GH=glmmML(stable~Sex+Age+Height+Weight+Surface+Vision, cluster=Subject,family=binomial,data=ctsib,method="ghq")
summary(fit.c.ML.GH)


# Laplace tilnærmelse med lmer
library(lme4)
fit.c.lme.Laplace=lmer(stable~Sex+Age+Height+Weight+Surface+Vision+(1|Subject),family=binomial,data=ctsib)
summary(fit.c.lme.Laplace)

# Gauss-Hermite tilnærmelse med lmer
fit.c.lme.GH=lmer(stable~Sex+Age+Height+Weight+Surface+Vision+(1|Subject),family=binomial,data=ctsib, nAGQ=10)
summary(fit.c.lme.GH)


# Punkt d
# Vi undersøker så hvilke variable som er signifikanet.
# Vi nøyer oss med å bruke lmer med Gauss-Hermite tilnærmelse

# Tar utgangspunkt i siste modell i forrige punkt
fit1= fit.c.lme.GH

# Prøver alle modeller der en kovariat er kuttet ut:
drop1(fit1,test="Chisq")

# Kutter ut alder. Estimerer modell uten alder
fit2=lmer(stable~Sex+Height+Weight+Surface+Vision+(1|Subject),family=binomial,data=ctsib, nAGQ=10)

# Prøver modeller der en kovariat til kuttes ut:
drop1(fit2,test="Chisq")


# Kutter ut vekt. Estimerer modell uten vekt
fit3=lmer(stable~Sex+Height+Surface+Vision+(1|Subject),family=binomial,data=ctsib, nAGQ=10)


# Prøver modeller der en kovariat til kuttes ut:
drop1(fit3,test="Chisq")


# Kutter ut høyde. Estimerer modell uten høyde
fit4=lmer(stable~Sex+Surface+Vision+(1|Subject),family=binomial,data=ctsib, nAGQ=10)

# Prøver modeller der en kovariat til kuttes ut:
drop1(fit4,test="Chisq")


# Kutter ut kjønn. Estimerer modell uten kjønn
fit5=lmer(stable~Surface+Vision+(1|Subject),family=binomial,data=ctsib, nAGQ=10)

drop1(fit5,test="Chisq")

fit6=lmer(stable~Surface+(1|Subject),family=binomial,data=ctsib, nAGQ=10)
fit7=lmer(stable~1+(1|Subject),family=binomial,data=ctsib, nAGQ=10)


# Sammenligner modellene:

anova(fit7,fit6,fit5,fit4,fit3,fit2,fit1,test="Chisq")


# Vi ender opp modellen med Surface og Vision:
summary(fit5)