##Load the nlme package ## library(nlme) install.packages('lsmeans') ## I manually inputted the data ## content <- c(40.89,37.99,37.18,34.98,34.89,42.07,41.22,49.42,45.85,50.15, 41.99,46.69,44.57,52.68,37.61,36.94,46.65,40.23,41.90,39.20,43.29, 40.45,42.91,39.97) block <- as.factor(sort(rep(1:4,6))) nitrogen <- as.factor(c(2,5,4,1,6,3,1,3,4,6,5,2,6,3,5,1,2,4,2,4,6,5,3,1)) ## Made a data frame called wheat wheat <- data.frame(content,block,nitrogen) ## Always a good idea to make a boxplot or visualize the data before performing a formal analysis boxplot(content ~ nitrogen,data = wheat) ## Our model ## model <- lme(content ~ nitrogen , data = wheat, random = ~ 1|block, method = 'REML') anova.lme(model) ## So we see that not all treatments are the same. We now ask : What treatments do differ? ## ## This is where pairwise contrasts become useful ## ## Compare each treatment to one another to see what differs ## library(lsmeans) lsmeans(model,pairwise ~ nitrogen , adjust = 'Tukey') ## Cotton Data ## residue <- c(120,110,120,100,140,130,71,71,70,76,63,68) trt <- as.factor(sort( rep(c("A","B"),6) )) batch <- as.factor(sort(rep(1:6,2))) sample <- as.factor(1:12) cotton <- data.frame(residue,trt,batch,sample) View(cotton) boxplot(residue ~ trt) out <- lme(residue ~ trt , data = cotton , random = ~ 1|sample(batch) , method = "REML" ) anova.lme(out) lsmeans(out,pairwise~trt) ## Rabbit Data ## temp <- c(-.3,-.2,1.2,3.1,-.5,2.2,3.3,3.7,-1.1,2.4,2.2,2.7,1.0,1.7,2.1,2.5,-.3,.8,.6,.9,-1.1,-2.2,.2,.3,-1.4,-.2,-.5,-.1,-.1,-.1,-.5,-.3,-.2, .1,-.2,.4,-.1,-.2,.7,-.3,-1.8,.2,.1,.6,-.5,0.0,1.0,.5,-1.0,-.3,-2.1,.6,.4,.4,-.7,-.3,-.5,.9,-.4,-.3) treat <- as.factor(sort(rep(c('A','B','C'),20))) time <- as.factor(rep(c(0,30,60,90),15)) rabbit <- as.factor(sort(rep(1:15,4))) data <- data.frame(rabbit,temp,time,treat) ## These are called Profile Plots -- Very useful for longitudinal data analysis## interaction.plot(time,treat,temp,mean) ##model## o <- lme(temp ~ treat + time + treat*time, random = ~time|rabbit , correlation = corCompSymm(form = ~time|rabbit)) anova.lme(o)