# edps 587
# spring 2020
# c.j.anderson
#
# Riesby depression data from Hedeker:  SERIAL CORRELATION
#
# Note the R code uses REML so you should not compare models 
# that have different fixed effects structures using anything
# based on the max of the likelikehood (e.g., AIC, BIC). 
#
#
#
library(nlme)       # new one
library(lme4)
library(lmerTest)
library(texreg)

source("D:\\Dropbox\\edps587\\lectures\\9 longitudinal\\All_functions.txt")

hd <- read.table("D:\\Dropbox\\edps587\\lectures\\9 longitudinal\\riesby_data_time_invarient_for_r.txt", header=TRUE)

# Remove missing --- not ideal but for lecture....
hd <- na.omit(hd)

hd$hamd <- as.numeric(hd$hamD)
hd$wksq <- hd$week*hd$week
		


#		   
# empty model, no serial correlation  
# (use REML to get lmer & nlme.0 to be exactly equal)
#
lmer.0 <- lmer(hamd ~ 1 + (1|id), data=hd)

#-- note "lme" is command in nlme
nlme.0 <- lme(hamd ~ 1,
              random= ~1 | id, 
			  data=hd)

# more fixed effects, no serial correlation 			  
lmer.1 <- lmer(hamd ~ 1 + week + wksq + endog + (1 |id), data=hd)

summary( nlme.1 <- lme(hamd ~ 1+ week + wksq + endog,
              random= ~ 1 | id, 
			  data=hd) )

			  
#
# AR(1): Add serial to model with random intercept & slope
#
# ~ page 			  
summary( nlme.2 <- lme(hamd ~ 1+ week + wksq + endog,
              random= ~ 1 + week + wksq | id, 
			  correlation=corAR1 (form= ~ 1 + week |id),
			  data=hd))
			  
#
# I used SAS to check that I programed this correctly),
#
# The results from R code below == SAS results (except some small differences
#                              variances & covariances between random effects)
#
#  SAS: I changed method=REML
#  R:  I re-coded endog so that used same dummy as SES (make comparison easier)
#
# ~ page 58
hd$nodog <- ifelse(hd$endog==1,0,1)

summary( nlme.3 <- lme(hamd ~ 1 + week + wksq + nodog,
             random= ~ 1 + week + wksq | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 

#
# Movging averge using ARMA(0,1) ??  I don't think so.  Below gives same result as nlme3, corAR1)
#
summary( nlme.4 <- lme(hamd ~ 1 + week + wksq + nodog,
             random= ~ 1 + week + wksq | id,
			 correlation=corARMA(form= ~ 1 + week|id, p=0,q=1),
			 data=hd) )

#
#ARMA(1,1) -- doesn't converge (in SAS "estimated G matrix not positive definite")
#			
#R:  Error in lme.formula(hamd ~ 1 + week + nodog, random = ~1 + week | id,  : 
#    nlminb problem, convergence error code = 1
#    message = iteration limit reached without convergence (10)
# 
nlme.5 <- lme(hamd ~ 1 + week  + nodog,
             random= ~ 1 + week  | id,
			 correlation=corARMA(form= ~ 1+ week|id, p=1,q=1),
			 data=hd)
			 
#
# Drop random wksq 
#
summary( nlme.6 <- lme(hamd ~ 1 + week + wksq + nodog,
             random= ~ 1 + week  | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 
			 
# Changed random effects (need to compute mixture of p-values)
anova(nlme.6,nlme.3)
p1 <- 1-pchisq(3.40522,3)
p2 <- 1-pchisq(3.40522,2)
(p.value <- .5*(p1+p2))

intervals(nlme.6,method=profile)
#
# Drop fixed wksq 
#
summary( nlme.6b <- lme(hamd ~ 1 + week + nodog,
             random= ~ 1 + week  | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 

#
# Drop random wksq & random week, so just random intercept with AR1
#
summary( nlme.7 <- lme(hamd ~ 1 + week + wksq + nodog,
             random= ~ 1   | id, 
			 correlation=corAR1 (form= ~ 1 |id),
			 data=hd) ) 
anova(nlme.7,nlme.6)
p1 <- 1-pchisq(24.26085,2)
p2 <- 1-pchisq(24.26085,1)
(p.value <- .5*(p1+p2))

			 
screenreg(list(nlme.0,nlme.1,nlme.2,nlme.6,nlme.7))

# Now getting rid of fixed effect for "endog" from nlme.6
#  because it is not significant
summary( nlme.8 <- lme(hamd ~ 1 + week + wksq ,
             random= ~ 1 + week  | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 
# NO LR test between nlme.6 and nlme.8 because we're using REML

# Now getting rid of fixed effect for week**2 from model nlme.8
#  because it is not significant
summary( nlme.8 <- lme(hamd ~ 1 + week ,
             random= ~ 1 + week  | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 
# NO LR test between nlme.6 and nlme.8 because we're using REML
	

screenreg(list(nlme.0,nlme.1,nlme.2,nlme.6,nlme.7,nlme.8))

###############################################################
# Final Model: No U0j but random week, AR(1)                  #
#                                                             #
# ~ page 68-70                                                #
###############################################################
summary( nlme.9 <- lme(hamd ~ 1 + week + nodog,
             random= ~ -1 + week  | id, 
			 correlation=corAR1 (form= ~ 1 + week|id),
			 data=hd) ) 
#
#  The variance covariance matrix for the fixed effects
#
Cov.beta9 <- (nlme.9$varFix)
# and SE of the paramter estimates
se.beta9 <- sqrt(diag(nlme.9$varFix))
#
# and the Hessian:   Cov.beta9. = to (−H)−1.
# "solve" is base R function to get inverase of a matrix

H <- -solve(Cov.beta9)

#############################################################
#  Get model based covariance & correlation matrix for Y    #
#  This does require some linear alegbra (matrix opertations#
#############################################################
# Part to between indivuals (random slope but not random intercept)
tau2 <-1.214405**2
week <- matrix(c(0,1,2,3,4,5),nrow=6,ncol=1)
timeXtime <- week %*% t(week)
(between <- tau2*timeXtime)

# Part due to serial correlation (i.e., AR(1))
rho <- .473371
sigma2 <- nlme.9$sigma**2
omega <- diag(6)
for (i in 2:6) { 
  for (j in 1:5) { 
    power <- i-1
	omega[i,j] <- rho**power
	omega[j,i] <- omega[i,j]
  }
 }
 
(within <- sigma2*omega)
(covY <- between + within)

# Correlation matrix for Y
isd <- solve(diag(sqrt(diag(covY))))
Ry <- isd %*% covY %*% isd