# R script for lecture 10 in STK4900/9900

# Slide 7: Pairwise comparison - Fecal fat -------------------------------------

fecfat <- read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/data/fecfat.txt", header=T)
x <- fecfat$fecfat[fecfat$pilltype==1]
y <- fecfat$fecfat[fecfat$pilltype==2]
t.test(x,y,paired=T)



# Slide 8: Pairwise comparison - Birth weight ----------------------------------

babies <- read.table("http://www.uio.no/studier/emner/matnat/math/STK4900/data/gababies.txt",header=T)
first <- babies$bweight[babies$birthord==1]
fifth <- babies$bweight[babies$birthord==5]
diff <- fifth-first
t.test(diff)

t.test(fifth,first,paired=T)


# Slide 9: Fixed effects models ------------------------------------------------


fit <- lm(fecfat~factor(pilltype)+factor(subject),data=fecfat)
anova(fit)

fit <- lm(bweight~factor(birthord)+factor(momid),data=babies)
anova(fit)



# Slide 11: Random effects model - Fecal fat -----------------------------------

library(nlme)

fit.fecfat <- lme(fecfat~factor(pilltype), random=~1|subject, data=fecfat)
summary(fit.fecfat)



# Slide 14: Random effects model - Birth weight --------------------------------

fit.babies <- lme(bweight~birthord+initage, random=~1|momid, data=babies)
summary(fit.babies)



# Slide 20: GEE - logistic model -----------------------------------------------

geefit <- gee(lowbrth~birthord+initage,id=momid,family=binomial,data=babies,corstr="unstructured")
summary(geefit)

or.ci(geefit)



# Slide 28: Time series -----------------------------------------------

par(mfrow=c(2,2))
time <- seq(1,100,1)
a <- c(0.1, 0.5, 0.9,-0.9)
for(i in 1:4){
  y_t <- arima.sim(n=100,model=list(ar=a[i]),sd=1)
  plot(time, y_t, type="l", main=a[i])
}


y_t <- arima.sim(n=100,model=list(ar=a[2]),sd=1)

arima(y_t, order=c(1,0,0))