#Example 1: Phlebitis

library(ggplot2)
library(nlme) 

setwd("C:/Users/Mark.Williamson.2/Desktop/Williamson Data/Example Datasets")
phlebitisdata <- read.table("phlebitis.csv", header=T, sep=",")
attach(phlebitisdata)
phlebitisdata #This isn’t necessary, but you might want to see the data structure

g1 <- ggplot(data = phlebitisdata, aes(x = Time, y = Y, group = Animal))
g1 + geom_point()
g1 + geom_line() #spaghetti plot

means<-aggregate(x=Y, by=list(Treatment,Time), FUN=mean)

names(means)[names(means)=="Group.1"] <- "Treatment"
names(means)[names(means)=="Group.2"] <- "Time"
names(means)[names(means)=="x"] <- "Y_mean"
means$Treatment<-factor(means$Treatment)
print(means)

g2 <-ggplot(data=means, aes(x=Time, y=Y_mean, group=Treatment, color=Treatment))
g2 + geom_line(size=2) 

interaction.plot(Time, factor(Treatment), Y, lty=c(1:3),lwd=2,ylab="mean of Y", xlab="time", trace.label="Treatment")
 
#Model:   Y = Treatment + Animal(Treatment) + Time + Treatment*Time + Error
aov.p = aov(Y~(factor(Treatment)*factor(Time))+Error(factor(Animal)),phlebitisdata )
summary(aov.p)


nestinginfo <- groupedData(Y ~ Treatment | Animal, data= phlebitisdata)

fit.compsym <- gls(Y ~ factor(Treatment)*factor(Time), data=nestinginfo, corr=corCompSymm(, form= ~ 1 | Animal))

fit.nostruct <- gls(Y ~ factor(Treatment)*factor(Time), data=nestinginfo, corr=corSymm(, form= ~ 1 | Animal), 
                    weights = varIdent(form = ~ 1 | Time))

fit.ar1 <- gls(Y ~ factor(Treatment)*factor(Time), data=nestinginfo, corr=corAR1(, form= ~ 1 | Animal))

fit.ar1het <- gls(Y ~ factor(Treatment)*factor(Time), data=nestinginfo, corr=corAR1(, form= ~ 1 | Animal), 
                  weights=varIdent(form = ~ 1 | Time))

anova(fit.compsym, fit.nostruct, fit.ar1, fit.ar1het) #compares the models

anova(fit.compsym)
anova(fit.ar1)

detach(phlebitisdata)
################################################################################################################################


#Example 2: Beat the Blues (Linear Mixed Effects Model)

library("lme4")
library('HSAUR')
head(BtheB)

data("BtheB", package="HSAUR")

BtheB$subject <- factor(rownames(BtheB))
nobs <-nrow(BtheB)

BtheB_long <- reshape(BtheB, idvar='subject',
                      varying=c('bdi.2m','bdi.4m','bdi.6m','bdi.8m'),
                      direction='long')
BtheB_long$time<-rep(c(2,4,6,8, rep(nobs, 4)))

subset(BtheB_long, subject %in% c("1","2","3"))

layout(matrix(1:2, nrow=1))

ylim <- range(BtheB[, grep("bdi", names(BtheB))], na.rm = TRUE)

boxplot(subset(BtheB, treatment == "TAU")[, grep("bdi",
           names(BtheB))], main = "Treated as usual", ylab = "BDI",
           xlab = "Time (in months)", names = c(0, 2, 4, 6, 8), ylim = ylim)

boxplot(subset(BtheB, treatment == "BtheB")[, grep("bdi",
           names(BtheB))], main = "Beat the Blues", ylab = "BDI",
           xlab = "Time (in months)", names = c(0, 2, 4, 6, 8), ylim = ylim)

BtheB_lmer1 <- lmer(bdi ~ bdi.pre + time + treatment + drug + length + (1 | subject),
                    data = BtheB_long, REML=FALSE, na.action = na.omit) #Random intercept
          
BtheB_lmer2 <- lmer(bdi ~ bdi.pre + time + treatment + drug + length + (time | subject),
                    data = BtheB_long, REML=FALSE, na.action = na.omit) #Random intercept and slope


anova(BtheB_lmer1, BtheB_lmer2)

summary(BtheB_lmer1)

layout(matrix(1:2, ncol=2))
qint<-ranef(BtheB_lmer1)$subject[["(Intercept)"]]
qres<-residuals(BtheB_lmer1)

qqnorm(qint, ylab="Estimated random intercepts",
       xlim=c(-3,3), ylim=c(-20,20), main="Random intercepts")
qqline(qint)

qqnorm(qres, ylab="Estimated residuals",
       xlim=c(-3,3), ylim=c(-20,20), main="Residuals")
qqline(qres)

#DCAR, DAR, Non-Ignorable
################################################################################################################################


#Example 3: Respitory (Logistic Generalized Estimating Equation)

library("gee")

head(respiratory) #binary response; would use Logistic regression if not for the repeated measures

resp <- subset(respiratory, month > "0")
resp$baseline <- rep(subset(respiratory, month =="0")$status, rep(4, 111))

resp$nstat <- as.numeric(resp$status == "good") #nstat is a dummy coding for a poor respiratory status

resp_glm <- glm(status ~ centre + treatment + sex + baseline + age, data = resp, family = "binomial")

resp_gee1 <- gee(nstat ~ centre + treatment + sex + baseline + age, data = resp, family = "binomial",
                 id = subject, corstr = "independence", scale.fix = TRUE, scale.value = 1)

resp_gee2 <- gee(nstat ~ centre + treatment + sex + baseline + age, data = resp, family = "binomial",
                 id = subject, corstr = "exchangeable", scale.fix = TRUE, scale.value = 1)

summary(resp_glm)
summary(resp_gee1)
summary(resp_gee2)

#log-odds of treatment 1.299

se <-summary(resp_gee2)$coefficients["treatmenttreatment", "Robust S.E."]
coef(resp_gee2)["treatmenttreatment"] + c(-1,1)*se*qnorm(0.975)

exp(1.299)
exp(coef(resp_gee2)["treatmenttreatment"] + c(-1,1)*se*qnorm(0.975))
################################################################################################################################


#Example 4: Epilepsy (Poisson Generalized Estimating Equation)

head(epilepsy) #count response; would use Poisson regression if not for the repeated measures

itp <- interaction(epilepsy$treatment, epilepsy$period)
tapply(epilepsy$seizure.rate, itp, mean)
tapply(epilepsy$seizure.rate, itp, var)
per <- rep(log(2), nrow(epilepsy))
epilepsy$period <- as.numeric(epilepsy$period)

layout(matrix(1:2, nrow = 1))
ylim <- range(epilepsy$seizure.rate)
placebo <- subset(epilepsy, treatment == "placebo")
progabide <- subset(epilepsy, treatment == "Progabide")

boxplot(seizure.rate ~ period, data = placebo, ylab = "Number of seizures", xlab = "Period", ylim = ylim, main = "Placebo")
boxplot(seizure.rate ~ period, data = progabide, main = "Progabide", ylab = "Number of seizures", xlab = "Period", ylim = ylim)


layout(matrix(1:2, nrow = 1))
ylim <- range(log(epilepsy$seizure.rate + 1))

boxplot(log(seizure.rate + 1) ~ period, data = placebo, main = "Placebo", ylab = "Log number of seizures", xlab = "Period", ylim = ylim)
boxplot(log(seizure.rate + 1) ~ period, data = progabide, main = "Progabide", ylab = "Log number of seizures", xlab = "Period", ylim = ylim)


epilepsy_glm <- glm(seizure.rate ~ base + age + treatment + offset(per), 
                    data = epilepsy, family = "poisson")

epilepsy_gee1 <- gee(seizure.rate ~ base + age + treatment + offset(per), 
                     data = epilepsy, family = "poisson", id = subject, corstr = "independence", 
                     scale.fix = TRUE, scale.value = 1)

epilepsy_gee2 <- gee(seizure.rate ~ base + age + treatment + offset(per), 
                     data = epilepsy, family = "poisson", id = subject, corstr = "exchangeable", 
                     scale.fix = TRUE, scale.value = 1)

epilepsy_gee3 <- gee(seizure.rate ~ base + age + treatment + offset(per),
                     data = epilepsy, family = "poisson", id = subject, corstr = "exchangeable", 
                     scale.fix = FALSE, scale.value = 1)

summary(epilepsy_glm)
summary(epilepsy_gee1)
summary(epilepsy_gee2)
summary(epilepsy_gee3)