---
title: "MLM Chapter 6"
author: "Doug Luke"
date: "12/29/2019"
output: 
  html_document:
    toc: true
    toc_float: true
    number_sections: false
    theme: lumen
    highlight: textmate
---

# Longitudinal Models
{#chap:longitudinal} 

```{r setuplong,message=FALSE}
library(LukeMLM)
library(lme4)
library(xtable)
library(ggplot2)
```


## Longitudinal Data as Hierarchical: Time Nested Within Person
{#sec:longtimeperson}


## Intra-individual Change
{#sec:longintra}


```{r longintramodels,warning=FALSE,message=FALSE}
library(lmerTest)

# Fixing dataset name
data(lsoa)
lsoadat <- lsoa

lm_null <- lmer(numadl ~ 1 + (1|caseid), data = lsoadat, REML = FALSE)
lm_lin <- lmer(numadl ~ wave + (1|caseid), data = lsoadat, REML = FALSE)
lm_quad <- lmer(numadl ~ wave + wave2 + (1|caseid), data = lsoadat, REML = FALSE)

# Same models with age as the time variable
lsoadat$newage <- lsoadat$age84 + 2*lsoadat$wave

lma_lin <- lmer(numadl ~ newage + (1|caseid), data = lsoadat, REML = FALSE)
lma_quad <- lmer(numadl ~ newage + I(newage^2) + (1|caseid), data = lsoadat, REML = FALSE)

summary(lm_null)
summary(lm_lin)
summary(lm_quad)

```



```{r plotquad}

# Fixes for later: put into ggplot 
x <- seq(0,6,.01)

# Note that we divide x by 2 to get back to the scaling of wave and wave2
ypred <- (fixef(lm_quad)[1] + fixef(lm_quad)[2]*(x/2) + fixef(lm_quad)[3]*((x/2)^2))


plot(x,ypred,type="l",col="black",lty=1,
	xlim=c(0,6.5),ylim=c(0,2),
	ylab="Predicted number of ADLs",xlab="Years after study enrollment")

```

```{r plotquadage}

# Fixes for later: put into ggplot 
x <- seq(70,103,.01)

ypred <- (fixef(lma_quad)[1] + fixef(lma_quad)[2]*(x) + fixef(lma_quad)[3]*((x)^2))


plot(x,ypred,type="l",col="black",lty=1,
	xlim=c(65,105),ylim=c(0,7),
	ylab="Predicted number of ADLs",xlab="Age")

```

## Inter-individual Change
{#sec:longinter}


```{r longintermods}
# Use these in earlier sections
lma_0 <- lmer(numadl ~ poly(newage,2) + (1|caseid), data = lsoadat, REML = FALSE)
summary(lma_0)

lma_1 <- lmer(numadl ~ poly(newage,2) + sex + (1|caseid), data = lsoadat, REML = FALSE)
summary(lma_1)

lma_2 <- lmer(numadl ~ poly(newage,2) + sex + marstat + (1|caseid), data = lsoadat, REML = FALSE)
summary(lma_2)

lma_3 <- lmer(numadl ~ poly(newage,2)*sex + marstat + (1|caseid), data = lsoadat, REML = FALSE)
summary(lma_3)

```

```{r plotlonginterresults}

# Fixes for later: put into ggplot and don't use hand-coded values
# First, get predicted values for men and women
# Assume not married

adlmen <- predict(lma_3,newdata = data.frame('newage' = seq(70,105,0.5),
                                   'sex' = rep("Male",71),
                                   'marstat' = rep("No",71)),
                  re.form=NA)
adlwomen <- predict(lma_3,newdata = data.frame('newage' = seq(70,105,0.5),
                                   'sex' = rep("Female",71),
                                   'marstat' = rep("No",71)),
                  re.form=NA)


x <- seq(70,105,.5)


plot(x,adlmen,type="l",col="black",lty=1,
	xlim=c(65,110),ylim=c(0,7.1),
	ylab="Predicted number of ADLs",xlab="Age")

points(x,adlwomen,type="l",lty=2,col="gray20")

legend(c(70,80),c(6.5,4), 
       c("Unmarried Women",
         "Unmarried Men"),
	col=c("gray20","black"),lty=c(2,1),bty="n",cex=1.3)
```

```{r plotlonginterresultsggplot}

# Keeping the code here for now, but not using it in the book
# First, get predicted values for men and women
# Assume not married

adlmen <- predict(lma_3,newdata = data.frame('newage' = seq(70,105,0.5),
                                   'sex' = rep("Male",71),
                                   'marstat' = rep("No",71)),
                  re.form=NA)
adlwomen <- predict(lma_3,newdata = data.frame('newage' = seq(70,105,0.5),
                                   'sex' = rep("Female",71),
                                   'marstat' = rep("No",71)),
                  re.form=NA)

df.men <- data.frame(
  'adlpred' = adlmen,
  'newage' = seq(70,105,0.5),
  'Gender' = rep("Male",71)
)

df.women <- data.frame(
  'adlpred' = adlwomen,
  'newage' = seq(70,105,0.5),
  'Gender' = rep("Female",71)
)

df.adl <- rbind(df.men,df.women)

ggplot(data = df.adl, aes(x = newage, y = adlpred, group = Gender)) +
  geom_line(aes(linetype=Gender, color=Gender),size=1) +
  labs(x = "Age", y = "Predicted Number of ADLs") +
  xlim(65,112) +
  
  scale_y_continuous(breaks=0:7,labels=c("0","1","2","3","4","5","6","7")) +
  scale_linetype_manual(values=c("longdash", "solid")) +
  scale_color_manual(values=c('darkblue','red'))+
  
  theme_minimal() +
  theme(legend.position = c(.8,.3),
        legend.justification = c("left","top"),
        legend.title=element_text(size=14),
        legend.text=element_text(size=12))
```


## Alternative Covariance Structures
{#sec:longcovar}


```{r longaltcov,warning=FALSE,message=FALSE}
# Now some alternative correlations

library(nlme)


# No additional covariances fitted
lma_1_cor1 <- lme(numadl ~ poly(newage,2) + sex, random = ~1|caseid, data = lsoadat, correlation=NULL)
summary(lma_1_cor1)

# Compound symmetry
lma_1_cor2 <- lme(numadl ~ poly(newage,2) + sex, random = ~1|caseid, data = lsoadat, correlation=corCompSymm(form=~1|caseid))
summary(lma_1_cor2)

# Autoregressive
lma_1_cor3 <- lme(numadl ~ poly(newage,2) + sex, random = ~1|caseid, data = lsoadat, correlation=corAR1(form=~1|caseid))
summary(lma_1_cor3)

# Unstructured
lma_1_cor4 <- lme(numadl ~ poly(newage,2) + sex, random = ~1|caseid, data = lsoadat, 
                  correlation=corSymm(form=~ 1|caseid))
summary(lma_1_cor4)


```


