# Edps 587
# Spring  2019
# c.j.anderson
#
# Riesby data (hedeker) with time varying and invarient predictors
# --- just like any other HLM model
#


library(lme4)
library(lmerTest)
library(texreg)
library(optimx)
library(lattice)
library(ggplot2)

hd <- read.table("C:\\Users\\cja\\Dropbox\\edps587\\lectures\\9 longitudinal\\RIESBY_data_time_varying.txt", header=TRUE)

#
#null
#
summary(model.null <- lmer(hamD ~ 1 + (1 |id),data=hd, REML=FALSE,
                           control = lmerControl(optimizer ="Nelder_Mead")))

#
# with time
summary(model.t1 <- lmer(hamD ~ 1 + week + (1 |id),data=hd, REML=FALSE,
                               control = lmerControl(optimizer ="Nelder_Mead")))


#
# With time, time varying:  lnimi, lndmi
#            not varying:  sex, endog
#
summary(model.t2 <- lmer(hamD ~ 1 + week + lnimi + lndmi + sex + endog + (1 |id),data=hd, REML=FALSE,
                     control = lmerControl(optimizer ="Nelder_Mead")))

#
# Simpler model
#
summary(model.t3 <- lmer(hamD ~ 1 + week  + lndmi + (1 |id),data=hd, REML=FALSE,
                         control = lmerControl(optimizer ="Nelder_Mead")))

screenreg(list(model.null,model.t1,model.t2,model.t3))

#
# A couple of graphs
#
groups <- unique(hd$id)[sample(1:249,25)]
subset <- hd[hd$id %in% groups, ]

xyplot(hamD ~ lndmi | id , data=subset, col.line='black', type=c('p','r'),main='Riesby Hamlton Depression vs Ln(dmi)')

ggplot(hd, aes(x=lndmi, y=hamD, col=id, type='l'),main='Regressions all in one plot') 
               + geom_smooth(method="lm", se=FALSE ) + theme(legend.position="none")
			   
# Plot all regressions in one plot
person.list <- unique(hd$id)
plot(hd$week, hd$hamD,type = 'n', 
     cex.main = 1.5,
         xlab = 'Time', 
     ylab = "Hamliton Depression Index",
     main = "Separate Linear Regression for Person"
     )
# --- put in proper order
	j <- order(hd$week)
for (j in person.list) {
  sub <- hd[ which(hd$id==j),]
  fitted <- fitted(lm(hamD~week,sub))
  lines(sub$week,fitted,col=j)
}

 
# Plot but join
person.list <- unique(hd$id)
hamD.week <- aggregate(hamD ~ week +  id , data=hd, FUN="mean")

plot(hamD.week$week,hamD.week$hamD,type = 'n', 
     cex.main = 1.5,
         xlab = 'Time', 
     ylab = "Hamliton Depression Index",
     main = "Join Mean per Person Person"
     )
# --- put in proper order
	j <- order(hamD.week$week )
for (j in person.list) {
  sub <- hamD.week[ which(hamD.week$id==j),]
  lines(sub$week,sub$hamD,col=j)
 }