# clean the workspace

rm(list=ls())

# load R packages

library(nlme)

library(lattice)

library(multcomp)

# choose your .csv data file

dat <- read.csv(file.choose(), header=T)

# data formulation: create variable dog, eye

dat$dog <- as.factor(rep(rep(c('Morse','Agnes','Anjou','Alice','Arlo','Sherlock','Clem','Clyde'),c(6,3,rep(6,6))), 3))

dat$eye <- as.factor(rep(rep(1:15, each=3), 3))

dat$region <- factor(dat$region, levels=c('AC','VS','SPN'))

# remove missing data(rows)

dat <- data.frame(na.omit(dat))

# data formulation: create variable injection

dat$x <- rep('injected', nrow(dat))

dat$x[which(dat$injected==F)] <- 'uninjected'

dat$x <- factor(dat$x, levels=c('uninjected', 'injected'))

# mixed-effects model with nested random effects and just main effects

summary(fit0 <- lme(scone ~ region + x, random=~1|dog/eye/region, data=dat, method = 'ML'))

# F-tests for main effects of region and x

anova(fit0, type = 'marginal')

dat$fit <- as.vector(fitted(fit0))

#1)Calculate the actual means for injected vs. uninjected areas (these don?°•t need a formal post-hoc comparison because the xinjected coefficient and associated p-value directly tell you whether these 2 means differ).

# the descriptive statistics: minimum, lower quartile (25%th), median (50%th),mean, upper quartile(75%th), maximum

summary(dat$fit[dat$injected==TRUE]) #descriptive statistics for injected area

summary(dat$fit[dat$injected==FALSE]) #descriptive statistics for uninjected area

#2)Calculate the actual means for each level of region. Because there are 3 levels in this variable, you need to do step 3 below to fully describe the pattern in the main effect.

summary(dat$fit[dat$injected==TRUE & dat$region=='AC']) #descriptive statistics for injected area, region AC

summary(dat$fit[dat$injected==TRUE & dat$region=='VS']) #descriptive statistics for injected area, region VS

summary(dat$fit[dat$injected==TRUE & dat$region=='SPN']) #descriptive statistics for injected area, region SPN

summary(dat$fit[dat$injected==FALSE & dat$region=='AC']) #descriptive statistics for uninjected area, region AC

summary(dat$fit[dat$injected==FALSE & dat$region=='VS']) #descriptive statistics for uninjected area, region VS

summary(dat$fit[dat$injected==FALSE & dat$region=='SPN']) #descriptive statistics for uninjected area, region SPN

# Pairwise comparisons illustrating main effect of region (differences between regions).

summary(glht(fit0, linfct = mcp(region = "Tukey")))

# mixed-effects model with nested random effects, plus both main effects and an interaction

summary(fit1 <- lme(scone ~ region*x, random=~1|dog/eye/region, data=dat, method = 'ML'))

# Likelihood ratio test for overall significance of the interaction term region*x

# if this is significant, then you want to do the multiple comparison post-hoc tests that Zhen illustrated,

# but now using fit1 instead of fit0

anova(fit0, fit1)

# multiple comparison post-hoc analysis to LMM

# if you want to know whether there are differences between injection vs. uninjection

# within each region, and whether this effect differs across

# regions (interactions).

c1 <- rbind('AC: injected vs uninjected effect:' = c(0,0,0,1,0,0),

'SPN: injected vs uninjected effect:' = c(0,0,0,1,1,0),

'VS: injected vs uninjected effect:' = c(0,0,0,1,0,1),

'injected vs uninjected effect, AC vs SPN' = c(0,0,0,0,1, 0),

'injected vs uninjected effect, AC vs VS' = c(0,0,0,0,0, 1),

'injected vs uninjected effect, SPN vs VS' = c(0,0,0,0,-1,1)

)

summary(glht(fit1, c1))

# if you only want to know whether there are differences between injection vs. uninjection

# within each region

c2 <- rbind('AC: injected vs uninjected effect:' = c(0,0,0,1,0,0),

'SPN: injected vs uninjected effect:' = c(0,0,0,1,1,0),

'VS: injected vs uninjected effect:' = c(0,0,0,1,0,1))

summary(glht(fit1, c2))