setwd("D:/Courses/Longitudinal Data Analysis/Data Sets/")
setwd("E:/Courses/Longitudinal Data Analysis/Data Sets/")

sirota.data <- read.table("sirota.txt",header=T,na.strings = "-9")
attach(sirota.data)

sirota.data

panss.mat <- matrix(c(panss,panss.1,panss.2,panss.3,panss.4,panss.5),ncol=6)
sanss.mat <- matrix(c(sanss,sanss.1,sanss.2,sanss.3,sanss.4,sanss.5),ncol=6)
weight.mat <- matrix(c(weight,weight.1,weight.2,weight.3,weight.4,weight.5),ncol=6)

summary(panss.mat)
summary(sanss.mat)
summary(weight.mat)

plot(c(0,1,2,4,8,12),panss.mat[1,],xlab="Time",ylab="PANSS",ylim=c(30,170),
col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(0,1,2,4,8,12),panss.mat[i,],col=Treatment[i])
}

plot(c(0,1,2,4,8,12),sanss.mat[1,],xlab="Time",ylab="SANSS",ylim=c(15,120),
col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(0,1,2,4,8,12),sanss.mat[i,],col=Treatment[i])
}

plot(c(0,1,2,4,8,12),weight.mat[1,],xlab="Time",ylab="Weight",ylim=c(45,110),
col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(0,1,2,4,8,12),weight.mat[i,],col=Treatment[i])
}

table(Gender,Treatment)
plot(Treatment,Age)


# ====================

#  Analysis of PANSS

# ====================


panss.del <- panss.mat[,1]-panss.mat[,2:6]
panss.del.pct <- panss.del/panss.mat[,1]
plot(panss.del[,1],panss.del.pct[,1])
plot(panss.del[,2],panss.del.pct[,2])
plot(panss.del[,3],panss.del.pct[,3])
plot(panss.del[,4],panss.del.pct[,4])
plot(panss.del[,5],panss.del.pct[,5])

summary(panss.del)

plot(c(1,2,4,8,12),panss.del[1,],xlab="Time",ylab="PANSS Change",
ylim=c(-60,80),col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(1,2,4,8,12),panss.del[i,],col=Treatment[i])
}

cov(panss.mat, use = "pairwise.complete.obs")
cov(panss.del, use = "pairwise.complete.obs")
cor(panss.del, use = "pairwise.complete.obs")


# Analysis by Complete Cases

# Look at the response patterns.

is.na(panss.mat)[order(c(apply(is.na(panss.mat),1,sum))),]


cc <- (apply(is.na(panss.del),1,sum)==0)
cc


# Simple analysis comparing last time to first for complete cases.


fit.cc1 <- lm(panss.del[,5] ~ Gender*Treatment, subset=cc)
summary(fit.cc1)

table(Treatment,Gender,cc)

# Arrange all deltas in a vector for repeated measures analysis.

panss.del.1 <- c(panss.del)
panss.del.1

Sex <- rep(Gender,5)-1
Sex

Trt <- rep(Treatment,5)-1
Trt

Time <- c(rep(11,40),rep(10,40),rep(8,40),rep(4,40),rep(0,40))
Time

Subject.1 <- rep(Subject,5)
Subject.1

sel.cc <- rep(cc,5)

library('nlme')

fit.cc2 <- lme(panss.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit, subset=sel.cc )

summary(fit.cc2)


fit.cc3 <- aov(panss.del.1 ~ Sex*Trt*Time+Error(factor(Subject.1)), 
subset=sel.cc)
summary(fit.cc3)



# Analysis by Last Observation Carried Forward

panss.locf <- panss.mat

panss.locf[is.na(panss.locf[,2]),2] <- panss.locf[is.na(panss.locf[,2]),1]
panss.locf[is.na(panss.locf[,3]),3] <- panss.locf[is.na(panss.locf[,3]),2]
panss.locf[is.na(panss.locf[,4]),4] <- panss.locf[is.na(panss.locf[,4]),3]
panss.locf[is.na(panss.locf[,5]),5] <- panss.locf[is.na(panss.locf[,5]),4]
panss.locf[is.na(panss.locf[,6]),6] <- panss.locf[is.na(panss.locf[,6]),5]
panss.locf

fit.locf <- lm(panss.locf[,1]-panss.locf[,6] ~ factor(Gender)*factor(Treatment))
summary(fit.locf)


#  Analysis by the linear mixed model on all subjects

fit <- lme(panss.del.1 ~ Sex+Trt+Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit)

fit.1 <- lme(panss.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit.1)

plot(Subject.1[is.na(panss.del.1)==F],fit.1$resid[,1])
plot(Time[is.na(panss.del.1)==F],fit.1$resid[,2])
plot(Sex[is.na(panss.del.1)==F],fit.1$resid[,2])

qqnorm(fit.1$resid[,2])


fit.2 <- lme(panss.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.2)

anova(fit.1,fit.2)


Age.1 <- rep(Age,5)

fit.3 <- lme(panss.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.3)


# To get a likelihood ratio test for Age, fit models 2 and 3 by ML.

fit.2ML <- lme(panss.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

fit.3ML <- lme(panss.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.3ML)

anova(fit.2ML,fit.3ML)


# Consider dropping the sex by treatment interaction to get a "clean"
# assessment of the treatment effect.

fit.4ML <- lme(panss.del.1 ~ Age.1+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.4ML)

anova(fit.4ML,fit.3ML)


# Is it possible that the treatment effect is related to age?

fit.5ML <- lme(panss.del.1 ~ Age.1*Trt+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.5ML)

anova(fit.4ML,fit.5ML)


# Might we need a qudratic trend in time?

Time.sq <- Time^2

fit.6ML <- lme(panss.del.1 ~ Age.1+Trt+Sex*(Time+Time.sq)+Trt*(Time+Time.sq), 
random = ~1+Time|Subject.1, na.action=na.omit, method="ML" )

summary(fit.6ML)

anova(fit.4ML,fit.6ML)



# As a final model, include the age adjustment and fit by REML.

fit.4 <- lme(panss.del.1 ~ Age.1+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="REML" )

summary(fit.4)

plot(Subject.1[is.na(panss.del.1)==F],fit.4$resid[,1])
plot(Time[is.na(panss.del.1)==F],fit.4$resid[,2])
plot(Sex[is.na(panss.del.1)==F],fit.4$resid[,2])

qqnorm(fit.4$resid[,2])


# ====================

#  Analysis of SANSS

# ====================


sanss.del <- sanss.mat[,1]-sanss.mat[,2:6]
summary(sanss.del)

plot(c(1,2,4,8,12),sanss.del[1,],xlab="Time",ylab="SANSS Change",
ylim=c(-30,80),col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(1,2,4,8,12),sanss.del[i,],col=Treatment[i])
}

cov(sanss.mat, use = "pairwise.complete.obs")
cov(sanss.del, use = "pairwise.complete.obs")
cor(sanss.del, use = "pairwise.complete.obs")


# Analysis by Complete Cases

# Look at the response patterns.

is.na(sanss.mat)[order(c(apply(is.na(sanss.mat),1,sum))),]


cc <- (apply(is.na(sanss.del),1,sum)==0)
cc


# Simple analysis comparing last time to first for complete cases.


fit.cc1 <- lm(sanss.del[,5] ~ Gender*Treatment, subset=cc)
summary(fit.cc1)


# Arrange all deltas in a vector for repeated measures analysis.

sanss.del.1 <- c(sanss.del)
sanss.del.1


fit.cc2 <- lme(sanss.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit, subset=sel.cc )

summary(fit.cc2)


fit.cc3 <- aov(sanss.del.1 ~ Sex*Trt*Time+Error(factor(Subject.1)), 
subset=sel.cc)
summary(fit.cc3)



# Analysis by Last Observation Carried Forward

sanss.locf <- sanss.mat

sanss.locf[is.na(sanss.locf[,2]),2] <- sanss.locf[is.na(sanss.locf[,2]),1]
sanss.locf[is.na(sanss.locf[,3]),3] <- sanss.locf[is.na(sanss.locf[,3]),2]
sanss.locf[is.na(sanss.locf[,4]),4] <- sanss.locf[is.na(sanss.locf[,4]),3]
sanss.locf[is.na(sanss.locf[,5]),5] <- sanss.locf[is.na(sanss.locf[,5]),4]
sanss.locf[is.na(sanss.locf[,6]),6] <- sanss.locf[is.na(sanss.locf[,6]),5]
sanss.locf

fit.locf <- lm(sanss.locf[,1]-sanss.locf[,6] ~ factor(Gender)*factor(Treatment))
summary(fit.locf)

fit.locf.1 <- lm(sanss.locf[,1]-sanss.locf[,6] ~ factor(Treatment))
summary(fit.locf.1)



#  Analysis by the linear mixed model on all subjects

fit <- lme(sanss.del.1 ~ Sex+Trt+Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit)

fit.1 <- lme(sanss.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit.1)

plot(Subject.1[is.na(sanss.del.1)==F],fit.1$resid[,1])
plot(Time[is.na(sanss.del.1)==F],fit.1$resid[,2])
plot(Sex[is.na(sanss.del.1)==F],fit.1$resid[,2])

qqnorm(fit.1$resid[,2])


fit.2 <- lme(sanss.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.2)

anova(fit.1,fit.2)


Age.1 <- rep(Age,5)

plot(Age.1[is.na(sanss.del.1)==F],fit.2$resid[,2])

fit.3 <- lme(sanss.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.3)


# To get a likelihood ratio test for Age, fit models 2 and 3 by ML.

fit.2ML <- lme(sanss.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

fit.3ML <- lme(sanss.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.3ML)

anova(fit.2ML,fit.3ML)


# Consider dropping the sex by treatment interaction to get a "clean"
# assessment of the treatment effect.

fit.4ML <- lme(sanss.del.1 ~ Age.1+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.4ML)

anova(fit.4ML,fit.3ML)


# Is it possible that the treatment effect is related to age?

fit.5ML <- lme(sanss.del.1 ~ Age.1*Trt+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

fit.5aML <- lme(sanss.del.1 ~ Age.1*Trt+Sex*Time+Trt*Time, random = ~1|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.5aML)

anova(fit.4ML,fit.5aML)


# Might we need a qudratic trend in time?

Time.sq <- Time^2

fit.6ML <- lme(sanss.del.1 ~ Age.1+Trt+Sex*(Time+Time.sq)+Trt*(Time+Time.sq), 
random = ~1+Time|Subject.1, na.action=na.omit, method="ML" )

summary(fit.6ML)

anova(fit.4ML,fit.6ML)



# As a final model, include the age adjustment and fit by REML.

fit.4 <- lme(sanss.del.1 ~ Age.1+Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="REML" )

summary(fit.4)

plot(Subject.1[is.na(sanss.del.1)==F],fit.4$resid[,1])
plot(Time[is.na(sanss.del.1)==F],fit.4$resid[,2])
plot(Sex[is.na(sanss.del.1)==F],fit.4$resid[,2])

qqnorm(fit.4$resid[,2])




# ====================

#  Analysis of WEIGHT

# ====================

weight.del <- weight.mat[,1]-weight.mat[,2:6]
summary(weight.del)

plot(c(1,2,4,8,12),weight.del[1,],xlab="Time",ylab="Weight Change",
ylim=c(-20,15),col=Treatment[1],type="l")
for (i in 2:40){
  lines(c(1,2,4,8,12),weight.del[i,],col=Treatment[i])
}

cov(weight.mat, use = "pairwise.complete.obs")
cov(weight.del, use = "pairwise.complete.obs")
cor(weight.del, use = "pairwise.complete.obs")


# Analysis by Complete Cases

# Look at the response patterns.

# Simple analysis comparing last time to first for complete cases.

fit.cc1 <- lm(weight.del[,5] ~ Gender*Treatment, subset=cc)
summary(fit.cc1)


# Arrange all deltas in a vector for repeated measures analysis.

weight.del.1 <- c(weight.del)
weight.del.1


library('nlme')

fit.cc2 <- lme(weight.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit, subset=sel.cc )

summary(fit.cc2)


fit.cc3 <- aov(weight.del.1 ~ Sex*Trt*Time+Error(factor(Subject.1)), 
subset=sel.cc)
summary(fit.cc3)



# Analysis by Last Observation Carried Forward

weight.locf <- weight.mat

weight.locf[is.na(weight.locf[,2]),2] <- weight.locf[is.na(weight.locf[,2]),1]
weight.locf[is.na(weight.locf[,3]),3] <- weight.locf[is.na(weight.locf[,3]),2]
weight.locf[is.na(weight.locf[,4]),4] <- weight.locf[is.na(weight.locf[,4]),3]
weight.locf[is.na(weight.locf[,5]),5] <- weight.locf[is.na(weight.locf[,5]),4]
weight.locf[is.na(weight.locf[,6]),6] <- weight.locf[is.na(weight.locf[,6]),5]
weight.locf

fit.locf <- lm(weight.locf[,1]-weight.locf[,6] ~ factor(Gender)*factor(Treatment))
summary(fit.locf)


#  Analysis by the linear mixed model on all subjects

fit <- lme(weight.del.1 ~ Sex+Trt+Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit)

fit.1 <- lme(weight.del.1 ~ Sex*Trt*Time, random = ~1|Subject.1, 
na.action=na.omit )

summary(fit.1)

plot(Subject.1[is.na(weight.del.1)==F],fit.1$resid[,1])
plot(Time[is.na(weight.del.1)==F],fit.1$resid[,2])
plot(Sex[is.na(weight.del.1)==F],fit.1$resid[,2])

qqnorm(fit.1$resid[,2])


fit.2 <- lme(weight.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.2)

anova(fit.1,fit.2)


Age.1 <- rep(Age,5)

fit.3 <- lme(weight.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit )

summary(fit.3)


# To get a likelihood ratio test for Age, fit models 2 and 3 by ML.

fit.2ML <- lme(weight.del.1 ~ Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

fit.3ML <- lme(weight.del.1 ~ Age.1+Sex*Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.3ML)

anova(fit.2ML,fit.3ML)


# Consider dropping the sex by treatment interaction to get a "clean"
# assessment of the treatment effect.

fit.4ML <- lme(weight.del.1 ~ Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="ML" )

summary(fit.4ML)

anova(fit.4ML,fit.2ML)


# Might we need a qudratic trend in time?

Time.sq <- Time^2

fit.6ML <- lme(weight.del.1 ~ Trt+Sex*(Time+Time.sq)+Trt*(Time+Time.sq), 
random = ~1+Time|Subject.1, na.action=na.omit, method="ML" )

summary(fit.6ML)

anova(fit.4ML,fit.6ML)



# As a final model, drop the age adjustment and fit by REML.

fit.4 <- lme(weight.del.1 ~ Sex*Time+Trt*Time, random = ~1+Time|Subject.1, 
na.action=na.omit, method="REML" )

summary(fit.4)

plot(Subject.1[is.na(weight.del.1)==F],fit.4$resid[,1])
plot(Time[is.na(weight.del.1)==F],fit.4$resid[,2])
plot(Sex[is.na(weight.del.1)==F],fit.4$resid[,2])

qqnorm(fit.4$resid[,2])