pupilDS_AX-CPT.rmd

---
title: "Pupillometry and AX-CPT"
author: "Jeremy Elman"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: readable
  pdf_document:
    latex_engine: xelatex
  word_document: default
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

Load libraries
```{r results='hide', message=FALSE, warning=FALSE}
library(ggplot2)
library(gridExtra)
library(gtools)
library(lme4)
library(arm)
library(lmerTest)
library(dplyr)
library(magrittr)
library(sjPlot)
library(sjmisc)
library(multcomp)
source("~/netshare/K/code/misc/summarySE.R")
source("~/netshare/K/code/misc/normDataWithin.R")
source("~/netshare/K/code/misc/summarySEwithin.R")
source("~/netshare/K/code/misc/get_legend.R")

sjp.setTheme(base=theme_bw())
```

Load data.
```{r}
pupilAXCPT = read.csv("~/netshare/K/Projects/Pupil_AX-CPT/data/pupil_AX-CPT.csv", 
                         stringsAsFactors = FALSE)
pupilAXCPT = pupilAXCPT[!is.na(pupilAXCPT$PCA) & 
                        !is.na(pupilAXCPT$CPTDPRIME_V2_nasp),]
pupilAXCPT = subset(pupilAXCPT, (VETSAGRP=='V1V2') | (VETSAGRP=='V2AR'))
str(pupilAXCPT)
```

Recode two-level factors to 0/1. Center continuous data.
*TRANSFORM*
```{r}
## Create factors
# Set 6 as reference level to contrast adjacent levels. Flip sign of 
# coefficient contrast 3 vs 6 to interpret slope.
pupilAXCPT$facLoad = factor(pupilAXCPT$Load)
pupilAXCPT$facLoad = relevel(pupilAXCPT$facLoad, ref="6")
loads = unique(pupilAXCPT$Load)
pupilAXCPT$MZ = ifelse(pupilAXCPT$zygos==1, 1, 0)
pupilAXCPT$site_v2rev = ifelse(pupilAXCPT$site_v2rev==1, 1, 0)
pupilAXCPT$Device = factor(pupilAXCPT$LR_Device)
pupilAXCPT$acamedbin = ifelse(pupilAXCPT$acamedtot==0, 0, 1)

#Bin MCI data into 0 = No impairment, 1 = Single domain MCI, 2 = Multi-domain MCI.
pupilAXCPT$rMCI_cons_3grp = ifelse(pupilAXCPT$rMCI_cons_v2pe==0, 0, 
                                        ifelse(pupilAXCPT$rMCI_cons_v2pe==1 | 
                                               pupilAXCPT$rMCI_cons_v2pe==2, 1, 2))
pupilAXCPT$rMCI_cons_3grp = factor(pupilAXCPT$rMCI_cons_3grp, 
                                   labels=c("CN","SMCI","MMCI"))

#Bin MCI data into 0 = No impairment, 1 = Non-amnestic S-MCI, 2 = Amnestic S-MCI, 3 = Multi-domain MCI.
pupilAXCPT$rMCI_cons_4grp = pupilAXCPT$rMCI_cons_v2pe
pupilAXCPT$rMCI_cons_4grp = gsub(4, 3, pupilAXCPT$rMCI_cons_4grp)
pupilAXCPT$rMCI_cons_4grp = factor(pupilAXCPT$rMCI_cons_4grp,
                                   labels=c("CN","naSMCI","aSMCI","MMCI"))

# Original MCI categories: 0 = No impairment, 1 = Non-amnestic S-MCI, 2 = Amnestic S-MCI, 3 = Non-amnestic M-MCI, 4 = Amnestic M-MCI
pupilAXCPT$rMCI_cons_v2pe = factor(pupilAXCPT$rMCI_cons_v2pe, 
                                   labels=c("CN","naSMCI","aSMCI","naMMCI","aMMCI"))

# Center continuous variables to interpret intercept at mean value (rather than 0).
# Subject specific variables centered at subject level, trial specific variables centered at trial level.
contTrialVars = c("PCA","pctPCA","adjPCA")
contSubjVars = c("age_v2","nas201tran","DSFMAX_V2_nasp","CPTDPRIME_V2_nasp")

# Center trial specific variables
for (x in contTrialVars) {
  newcol = paste0('c.',x)
  pupilAXCPT[[newcol]] = as.numeric(scale(pupilAXCPT[[x]], center=TRUE, scale=FALSE))
}

# Center subject specific variables
subjDF = pupilAXCPT[c("vetsaid",contSubjVars)]
subjDF %<>% group_by(vetsaid) %>% dplyr::summarise_each(funs(mean))
nums <- sapply(subjDF, is.numeric)
c.subjDF = as.data.frame(apply(subjDF[,nums],2, function(y) y - mean(y, na.rm=TRUE)))
names(c.subjDF) = paste0("c.",names(c.subjDF))
c.subjDF$vetsaid = subjDF$vetsaid
pupilAXCPT = left_join(pupilAXCPT, c.subjDF, by="vetsaid")

# Create quantile groups data
qVars = c("DSFMAX_V2_nasp","CPTDPRIME_V2_nasp")
new.qVars = qVars
names(new.qVars) = paste0("q.",qVars)
pupilAXCPT %<>%
  group_by(vetsaid) %>%
  dplyr::summarise_each_(funs(mean),qVars) %>%
  mutate_each_(funs(quantcut(., labels=seq(1,4))),qVars) %>%
  dplyr::select_("vetsaid"="vetsaid",.dots=new.qVars) %>%
  inner_join(pupilAXCPT, by="vetsaid")

# Create residPCA by regressing BaselineDiameter from PCA and taking the residuals.
lme.PCA.resid = lmer(PCA ~ BaselineDiameter + (1 | case/vetsaid),data=pupilAXCPT)
pupilAXCPT$residPCA = resid(lme.PCA.resid)
```

# Basic sample descriptives
```{r, include=FALSE}
#Create dataset with subject level variables to assess relationships with BOLD variance.

subjDatAXCPT = pupilAXCPT %>% 
            dplyr::select(-contains("PCA"),-contains("Load")) %>%
            group_by(vetsaid) %>%
            summarise_each(funs(first))
#write.csv(subjDatAXCPT, "~/netshare/K/Projects/Pupil_AX-CPT/data/AX-CPT_subjDat.csv")
```

How many subjects?
```{r}
n_distinct(subjDatAXCPT$vetsaid)
```

How many attrition replacements?
```{r}
dplyr::count(subjDatAXCPT, VETSAGRP)
```

How many twin pairs vs unpaired twins? 
```{r}
subjDatAXCPT %>%
  group_by(case) %>%
  dplyr::summarise(n_twins = n_distinct(vetsaid)) %>%
  dplyr::count(n_twins)
```


How many MZ and DZ pairs? (excludes unpaired subjects)
```{r}
subjDatAXCPT %>% 
  group_by(case) %>%
  mutate(n_twins = n_distinct(vetsaid)) %>%
  filter(n_twins > 1) %>%
  dplyr::summarise(zyg = mean(zygos)) %>%
  group_by(zyg) %>%
  dplyr::count(zyg)
```

How many levels of load did subjects complete? Shows number of subjects who completed 1, 2, or 3 levels (corresponding to digit spans of 3, 6, and 9).
```{r}
pupilAXCPT %>%
  group_by(vetsaid) %>%
  dplyr::summarise(load = n()) %>%
  dplyr::count(load)
```

How many of each MCI group?
```{r}
dplyr::count(subjDatAXCPT, rMCI_cons_v2pe)
```

How many of each MCI group binned by single or multi-doman MCI?
```{r}
dplyr::count(subjDatAXCPT, rMCI_cons_3grp)
```

How many of each MCI group binned by single amnestic/non-amnestic or multi-doman MCI?
```{r}
dplyr::count(subjDatAXCPT, rMCI_cons_4grp)
```

Plot frequencies of MCI measures
```{r}
sjp.frq(subjDatAXCPT$rMCI_cons_v2pe, 
          title="MCI groups",
          axis.labels=c("Normal","Single Domain\nNon-amnestic",
                         "Single Domain\nAmnestic","Multi Domain\nNon-Amnestic",
                         "Multi Domain Amnestic"))
```

```{r, include=FALSE}
# png("~/netshare/K/Projects/Pupil_AX-CPT/results/MCI_groups.png", width=8, height=6)
# grid.arrange(mci.p1$plot,mci.p2$plot)
# dev.off()
```


Plot pctPCA and residPCA data to check for ceiling effects on pupil dilation
```{r}
pctPcaPlot = ggplot(pupilAXCPT, aes(x=pctPCA)) + 
  geom_histogram(fill="steelblue3",color="black",size=.4) + 
  facet_wrap(~ Load, nrow=3) + 
  ggtitle("Histogram of pctPCA by Load") + 
  xlab("% PCA") + theme_bw() 

residPcaPlot = ggplot(pupilAXCPT, aes(x=residPCA)) + 
  geom_histogram(fill="steelblue3",color="black",size=.4) + 
  facet_wrap(~ Load, nrow=3) + 
  ggtitle("Histogram of residPCA by Load") + 
  xlab("PCA (adjusted for baseline)") + theme_bw() 

grid.arrange(pctPcaPlot, residPcaPlot, ncol=2)
```

```{r, include=FALSE}
# png("~/netshare/K/Projects/Pupil_AX-CPT/results/PCAxLoad.png")
# grid.arrange(pctPcaPlot, residPcaPlot, ncol=2)
# dev.off()
```

--------------------------------------------------

# AX-CPT performance 

## Error Rates

d' by MCI group
```{r}
dprime.p1 = ggplot(subjDatAXCPT, aes(x=CPTDPRIME_V2_nasp)) + 
  geom_histogram(fill="steelblue3",color="black",size=.4) + 
  theme(plot.title=element_text(face="bold")) + theme_bw(14)

dprime.p2 = ggplot(subjDatAXCPT, aes(x=rMCI_cons_4grp,y=CPTDPRIME_V2_nasp)) + 
              geom_boxplot(fill="firebrick",alpha=.75) +
              geom_jitter(alpha=0.3,position=position_jitter(width = 0.1)) +
              theme(plot.title=element_text(face="bold")) +
              theme_bw(14)
g.dprime = grid.arrange(dprime.p1,dprime.p2, ncol=2, top="d' by MCI status")
summary(lm(CPTDPRIME_V2_nasp ~ rMCI_cons_4grp, data=subjDatAXCPT))
```

```{r,include=FALSE}
# png("~/netshare/K/Projects/Pupil_AX-CPT/results/AX-CPTperformance_X_MCI.png", width=16,height=12)
# grid.arrange(g.dprime,g.ax,g.bx,g.ay, ncol=2,top="AX-CPT Performance by MCI status")
# dev.off()
```

# Mixed Effects Models

Run model testing interaction of Load with d prime, digit span and MCI
```{r}
lme.facLoad.dprime = lmer(residPCA ~ facLoad*(CPTDPRIME_V2_nasp + 
                                                  c.DSFMAX_V2_nasp) + 
                                    apoe4 + c.age_v2 + site_v2rev + Device + acamedbin +  
                                    (1 | case/vetsaid), data=pupilAXCPT)
summary(lme.facLoad.dprime)
anova(lme.facLoad.dprime)
```

Test for differences based on d prime within each level of Load
```{r}
lme.dprime.loads <- lapply(loads, function(x) {
    lmer(residPCA ~ CPTDPRIME_V2_nasp + apoe4 + c.age_v2 + site_v2rev + 
                    Device + acamedbin + c.DSFMAX_V2_nasp  + 
                    (1 | case), 
         data=filter(pupilAXCPT, Load==x))
})
lapply(lme.dprime.loads, summary)
```

Plot pupil dilation by d' quantile
```{r, echo=FALSE}

# Plot raw values
summaryResidPCA.dprime = summarySEwithin(pupilAXCPT, 
                measurevar="residPCA", 
                idvar="vetsaid", 
                withinvars="Load",
                betweenvars="q.CPTDPRIME_V2_nasp",
                na.rm=TRUE)
p.Draw = ggplot(summaryResidPCA.dprime, aes(x=Load,y=residPCA,
                                        color=q.CPTDPRIME_V2_nasp,
                                        group=q.CPTDPRIME_V2_nasp)) +
  geom_line(size=1) +
  geom_errorbar(width=.2,size=1, aes(ymin=residPCA-ci,ymax=residPCA+ci)) +
  theme_bw(20) + ylab("Change in Pupil Diameter") +
  scale_color_brewer(name=" d' quartile", palette='RdBu') + 
          theme(plot.title = element_text(size=20),
                axis.title = element_text(size=18),
                axis.text.x = element_text(size=16),
                legend.background = element_rect(fill = "transparent"),
                legend.justification=c(1,0), legend.position=c(1,0), 
                legend.text.align=1,
                legend.title.align=0,
                legend.title = element_text(size=16),
                legend.text = element_text(size=16))

# Plot predicted values
dfDpred = sjp.int(lme.facLoad.dprime, type="eff", int.term="facLoad:CPTDPRIME_V2_nasp", 
               swap.pred=T, mdrt.values = "quart", prnt.plot=F)$data.list[[1]]
levels(dfDpred$grp) = seq(length(levels(dfDpred$grp)))
p.Dpred = ggplot(data=subset(dfDpred, !grp==3), aes(x=x, y=y, color=grp, group=grp)) + 
          geom_line(size=1) +
          geom_errorbar(width=.2,size=1, aes(ymin=conf.low,ymax=conf.high)) +
          ylab(expression(paste("Predicted ", Delta, " Pupil Diameter"))) + xlab("Load") +
          scale_color_brewer(name="d' quartile",palette='RdBu',direction=-1,labels=seq(4)) + 
          scale_x_continuous(breaks=c(3,6,9)) + theme_bw(16) +
          theme(plot.title = element_text(size=20),
                axis.title = element_text(size=18),
                axis.text.x = element_text(size=16),
                legend.background = element_rect(fill = "transparent"),
                legend.justification=c(1,0), legend.position=c(1,0), 
                legend.text.align=1,
                legend.title.align=0,
                legend.title = element_text(size=16),
                legend.text = element_text(size=16))

gDprime = arrangeGrob(p.Draw, p.Dpred,ncol=2)

ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_dprime_All.png", gDprime, height=6, width=12, dpi=300)
ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_dprime_Raw.png", p.Draw, height=6, width=8, dpi=300)
ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_dprime_Predicted.png", p.Dpred, height=6, width=8, dpi=300)
```

# Covariate Plot 

Plot pupil dilation by digit span max quantile
```{r, echo=FALSE}
summaryResidPCA.DSFMAX_v2 = summarySEwithin(pupilAXCPT, 
                measurevar="residPCA", 
                idvar="vetsaid", 
                withinvars="Load",
                betweenvars="q.DSFMAX_V2_nasp",
                na.rm=TRUE)
  
pDSFMAX_v2raw = summaryResidPCA.DSFMAX_v2 %>%
      ggplot(., aes(x=Load,y=residPCA,
                    color=q.DSFMAX_V2_nasp,
                    group=q.DSFMAX_V2_nasp)) +
        theme_bw(16) +
        ylab("Change in Pupil Diameter\n(Adjusted for Baseline)") +
        geom_line(size=1) +
        geom_errorbar(width=.2,size=1, aes(ymin=residPCA-ci,ymax=residPCA+ci)) +
        theme_bw(20) + ylab("Change in Pupil Diameter") +
        scale_color_brewer(name="d' quartile",palette='RdBu') + 
        theme(plot.title = element_text(size=20),
              axis.title = element_text(size=18),
              axis.text.x = element_text(size=16),
              legend.background = element_rect(fill = "transparent"),
              legend.justification=c(1,0), legend.position=c(1,0), 
              legend.text.align=1,
              legend.title.align=0,
              legend.title = element_text(size=16),
              legend.text = element_text(size=16))

# Plot predicted values
dfDSFRAWpred = sjp.int(lme.facLoad.dprime, type="eff", int.term="facLoad:c.DSFMAX_V2_nasp", 
               swap.pred=T, mdrt.values = "quart", prnt.plot=F)$data.list[[1]]
levels(dfDSFRAWpred$grp) = seq(length(levels(dfDSFRAWpred$grp)))
pDSFMAX_v2pred = ggplot(data=subset(dfDSFRAWpred, !grp==3), aes(x=x, y=y, color=grp, group=grp)) + 
          geom_line(size=1) +
          geom_errorbar(width=.2,size=1, aes(ymin=conf.low,ymax=conf.high)) +
          ylab(expression(paste("Predicted ", Delta, " Pupil Diameter"))) + xlab("Load") +
          scale_color_brewer(name="Max Digit Span quartile",palette='RdBu',direction=-1,labels=seq(4)) + 
          scale_x_continuous(breaks=c(3,6,9)) + theme_bw(16) +
          theme(plot.title = element_text(size=20),
                axis.title = element_text(size=18),
                axis.text.x = element_text(size=16),
                legend.background = element_rect(fill = "transparent"),
                legend.justification=c(1,0), legend.position=c(1,0), 
                legend.text.align=1,
                legend.title.align=0,
                legend.title = element_text(size=16),
                legend.text = element_text(size=16))

gDSFMAX = arrangeGrob(pDSFMAX_v2raw, pDSFMAX_v2pred,ncol=2)

ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_DSFMAX_All.png", gDSFMAX, height=6, width=12, dpi=300)
ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_DSFMAX_All_Raw.png", pDSFMAX_v2raw, height=6, width=8, dpi=300)
ggsave("~/netshare/K/Projects/Pupil_AX-CPT/results/pupil_DSFMAX_All_Predicted.png", pDSFMAX_v2pred, height=6, width=8, dpi=300)

plot(gDSFMAX)
```


----------------------------
```{r, include=FALSE}
#Save out analysis dataset
write.csv(pupilAXCPT, "~/netshare/K/Projects/Pupil_AX-CPT/data/pupilDS_AX-CPT_AnalysisDataset.csv")
```

```{r}
print(sessionInfo(), locale = FALSE)
```