###########################################################
## Beispiel vorlesung zu SOEP-Daten 
## HK 25.6.2015
###########################################################

# Daten laden
library(foreign)
soep <- read.dta(file="http://www.stat.uni-muenchen.de/~abender/daten/soep_lebensz.dta")
head(soep)
boxplot(soep$bildung)
table(soep$lebensz_org)
table(soep$jahr)

# Benötigte Pakete laden
library(statmod)
library(nlme)


### Variable Zufriedenheit metrisch 0-10 
### die ersten 3 Kategorein werden gestrichen 
soep$zuf<- as.numeric(soep$lebensz_org) -4

##  Modell mit geschlecht als Einflussgroesse
table(soep$zuf)

### Bildung zentrieren ###
soep$bildungz <- soep$bildung- mean(soep$bildung, na.rm=TRUE)
summary(soep$bildungz)
soepc <- subset(soep,bildungz!="NA" & zuf!="NA")
summary(soepc$zuf)

### Modell mit zufaelligem Personeneffekt
zmod1 <- lme(zuf ~ sex + as.factor(jahr), random = ~1|persnr, data=soepc)
summary(zmod1)

### Jetzt mit AR(1)
### Nicht ganz korrekt, da Jahre im Korrelationssturktur noch fehlen
zmod2 <- lme(zuf ~ sex + bildungz + jahr, random=~1|persnr, data=soepc,
             corr=corAR1(form=~1|persnr) )
summary(zmod2)

### Jetzt mit stetigem AR1
zmod3<-lme(zuf~sex+bildungz+jahr,random=~1|persnr, data=soepc,  
           corr=corCAR1(form=~jahr|persnr) )
summary(zmod3)