# brief R program for additional exercise 10.7 (VHM 802)
# focus is on new modeling features

salso <- read.csv("h:/vhm/vhm802/data_csv/hs10_7.csv")
salso$id <- 100*salso$group + salso$person
salso$group <- as.factor(salso$group)

# profile plots using nlme library
library(nlme)
salso.grp <- groupedData( salsolinol ~ day|id, outer= ~group, data=salso)
plot(salso.grp, outer = ~ group)
lnsalso.grp <- groupedData( log(salsolinol) ~ day|id, outer= ~group, data=salso)
plot(lnsalso.grp, outer = ~ group)

# Box-Cox analysis for fixed effects version of hierarchical model
library(MASS)
bc <- boxcox(salsolinol ~ group*as.factor(day)+as.factor(id), data=salso, lambda=seq(-1, 1, 0.01))
bc

# single time point (day=1) analysis by restriction
summary( lm( log(salsolinol) ~ group, data=salso[salso$day==1,]))
# same thing after reshape to wide data format
salso.wide <- reshape(salso, direction="wide", timevar="day", idvar="id", v.names="salsolinol")
summary( lm( log(salsolinol.1) ~ group, data=salso.wide))

# no repeated measures ANOVA in standard libraries

# hierarchical models, using nlme library
salso.hier <- lme( log(salsolinol) ~ group*as.factor(day), random=~1|id, data=salso.grp)
summary(salso.hier)
anova(salso.hier)
summary (lme( log(salsolinol) ~ group*day, random=~1|id, data=salso.grp))

# hierarchical models with additional correlation structure, using nlme library
salso.hier.ar1 <- update( salso.hier, correlation=corAR1())

# models with correlation structure only, using nlme library
# compound symmetry structure
salso.cs <- gls( log(salsolinol) ~ group*as.factor(day), correlation=corCompSymm (form = ~1|id), data=salso.grp)
summary(salso.cs)
anova(salso.cs)
# ar(1) structure
salso.ar1 <- gls( log(salsolinol) ~ group*as.factor(day), correlation=corAR1 (form = ~1|id), data=salso.grp)
summary(salso.ar1)
anova(salso.ar1)