setwd("D:/Courses/Longitudinal Data Analysis/Data Sets/") setwd("E:/Courses/Longitudinal Data Analysis/Data Sets/") rats.data <- read.table("rats.txt",header=T,na.strings = "-9") attach(rats.data) plot(Time,Response) sel <- (is.na(Response)==F) G1mean <- tapply(Response[(Group==1)&(sel==T)],Time[(Group==1)&(sel==T)],mean) G2mean <- tapply(Response[(Group==2)&(sel==T)],Time[(Group==2)&(sel==T)],mean) G3mean <- tapply(Response[(Group==3)&(sel==T)],Time[(Group==3)&(sel==T)],mean) print(cbind(G1mean,G2mean,G3mean)) plot(c(50,60,70,80,90,100,110),G1mean,col=1,ylim=c(70,85),pch=16) points(c(50,60,70,80,90,100,110),G2mean,col=2,pch=16) points(c(50,60,70,80,90,100,110),G3mean,col=3,pch=16) G1N <- tapply(Group[(Group==1)&(sel==T)],Time[(Group==1)&(sel==T)],sum) G2N <- tapply(Group[(Group==2)&(sel==T)]/2,Time[(Group==2)&(sel==T)],sum) G3N <- tapply(Group[(Group==3)&(sel==T)]/3,Time[(Group==3)&(sel==T)],sum) print(cbind(G1N,G2N,G3N)) X <- log(1+(Time-45)/10) X[1:7] plot(X[1:7],G1mean,col=1,ylim=c(70,85),pch=16) points(X[1:7],G2mean,col=2,pch=16) points(X[1:7],G3mean,col=3,pch=16) # TWO-STAGE ANALYSIS subj <- tapply(Subject,Subject,mean) grp <- tapply(Group,Subject,mean) b <- NULL for (i in 1:length(subj)){ fit <- lm(Response[(Subject==subj[i])&(sel==T)] ~ X[(Subject==subj[i])&(sel==T)]) b <- rbind(b,fit$coef) } print(cbind(grp,b)) plot(grp,b[,1]) plot(grp,b[,2]) sel1 <- (is.na(b[,2])==F) fit <- aov(b[,1] ~ factor(grp), na.action=na.omit, subset=sel1) summary(fit) fit <- aov(b[,2] ~ factor(grp), na.action=na.omit, subset=sel1) summary(fit) plot(fit) library('nlme') # The following 3 variables will let us fit separate slopes for # each group of mice, with no differences in intercept. X1 <- X*(Group==1) X2 <- X*(Group==2) X3 <- X*(Group==3) # First model ignores the nesting of data within mice. fit <- lm(Response ~ X1+X2+X3, subset=(sel==T)) summary(fit) # Second model has a random effect for the intercept of each mouse. fit <- lme(Response ~ X1+X2+X3, random = ~1|Subject, na.action=na.omit ) summary(fit) # Third model adds a random effect for slopes of each mouse. fit <- lme(Response ~ X1+X2+X3, random = ~1+X|Subject, na.action=na.omit ) summary(fit) # This model has a single common slope for each group. Comparing its # log likelihood to the last model lets us test whether there is # convincing evidence that the slope differs across groups. fit <- lme(Response ~ X, random = ~1+X|Subject, na.action=na.omit ) summary(fit) # An alternative approach is to define contrast variables for the groups. G1 <- ((Group==1)-(Group==3)) G2 <- ((Group==2)-(Group==3)) # The next model includes main effects for groups (i.e. different intercepts # in each group) and group by time (X) interaction (i.e. slopes vary by group). fit <- lme(Response ~ X*(G1+G2), random = ~1+X|Subject, na.action=na.omit ) summary(fit) # The following way to write the model adds the time main effect and the # time by group interactions but WITHOUT the group main effects. Thanks # to Amir for pointing this out. fit <- lme(Response ~ X+X:(G1+G2), random = ~1+X|Subject, na.action=na.omit ) summary(fit)