diff --git a/afq_step3_stats.R b/afq_step3_stats.R
index e80089e..122f699 100644
--- a/afq_step3_stats.R
+++ b/afq_step3_stats.R
@@ -15,6 +15,11 @@ library("lme4")
 dataDir <- "/Users/nmuncy/Projects/emu_AFQ/analyses/"
 privateDir <- "/Users/nmuncy/Projects/emu_private/"
 
+do_lunc = 0
+do_runc = 0
+do_lcing = 1
+do_rcing = 1
+do_latr = 1
 
 ### --- Step 1: Make dataset
 #
@@ -165,7 +170,6 @@ df_afq <- func_makeDF()
 
 
 
-
 ### --- Step 2: Check Memory behavior
 #
 # Run ANOVA for LGI, LDI
@@ -268,181 +272,883 @@ func_memStats()
 
 
 
-### --- Step 3: Model tract, no covariates
-#
-# Plot data, determine distribution,
-#   compare model families.
-#
-# Plot best model
-
-# Get data
-df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
-df_afq$Group <- factor(df_afq$Group)
-df_afq$Sex <- factor(df_afq$Sex)
-
-# Tract
-tract <- "UNC_L"
-df_tract <- df_afq[which(df_afq$tractID == tract), ]
-df_tract$dti_fa <- round(df_tract$dti_fa, 3)
-
-# plot mean data
-ggplot(data = df_tract) +
-  geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
-
-ggplot(data = df_tract) +
-  geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
-  geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
-
-# determine distribution
-descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
-
-fit.beta <- fitdist(df_tract$dti_fa, "beta")
-plot(fit.beta)
-fit.beta$aic
-
-fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
-plot(fit.gamma)
-fit.gamma$aic
-
-
-#  determine k, compare families
-fit_gamma <- bam(dti_fa ~ Group +
-   Sex +
-   s(nodeID, by=Group, k=40) +
-   s(subjectID, bs="re"),
-   data = df_tract,
-   family = Gamma(link = "logit"),
-   method = "REML")
-
-gam.check(fit_gamma, rep = 500)
-
-fit_beta <- bam(dti_fa ~ Group +
-   Sex +
-   s(nodeID, by=Group, k=40) +
-   s(subjectID, bs="re"),
-   data = df_tract,
-   family = betar(link = "logit"),
-   method = "REML")
-
-gam.check(fit_beta, rep = 500)
-
-
-infoMessages('on')
-compareML(fit_gamma, fit_beta)  # fit_gamma recommended
-
-# get stats
-summary(fit_gamma)  # R-sq = 0.819
-
-
-
-### --- Step 4: Model tract, covariates
-
-fit_cov_pds <- bam(dti_fa ~ Group +
-    Sex +
-    s(nodeID, by=Group, k=40) +
-    s(PDS, by=Sex) +
-    s(subjectID, bs="re"),
-    data = df_tract,
-    family = Gamma(link = "logit"),
-    method = "REML")
-
-gam.check(fit_cov_pds, rep = 500)
-
-
-# Test cov model against gamma
-# infoMessages('on')
-compareML(fit_gamma, fit_cov_pds) # PDS wins
-summary(fit_cov_pds) # R-sq = 0.85
-
-
-# plot
-df_pred <- predict.bam(
-  fit_cov_pds,
-  exclude_terms = c("PDS", "Sex","subjectID"),
-  values=list(PDS = NULL, Sex = NULL),
-  se.fit=T,
-  type="response")
-
-df_pred <- data.frame(Group=df_tract$Group,
-    Sex=df_tract$Sex,
-    subjectID=df_tract$subjectID,
-    PDS=df_tract$PDS,
-    nodeID=df_tract$nodeID,
-    fit=df_pred$fit,
-    se.fit=df_pred$se.fit)
-
-ggplot(data = df_pred) +
-  geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
-  ggtitle("GAM Fit of L. Uncinate FA Values") +
-  ylab("Fit FA")
-
-
-
-### --- Step 5: Test for differences
-#
-# Check for group differences in spline.
-
-par(mfrow=c(1,3))
-
-p01 <- plot_diff(fit_cov_pds,
-    view="nodeID",
-    comp=list(Group=c("0", "1")),
-    rm.ranef = T,
-    main = "Difference Scores, Low-Med",
-    ylab = "Est. FA difference",
-    xlab = "Tract Node")
-
-p02 <- plot_diff(fit_cov_pds,
-    view="nodeID",
-    comp=list(Group=c("0", "2")),
-    rm.ranef = T,
-    main = "Difference Scores, Low-High",
-    ylab = "Est. FA difference",
-    xlab = "Tract Node")
-
-p12 <- plot_diff(fit_cov_pds,
-    view="nodeID",
-    comp=list(Group=c("1", "2")),
-    rm.ranef = T,
-    main = "Difference Scores, Med-High",
-    ylab = "Est. FA difference",
-    xlab = "Tract Node")
-
-par(mfrow=c(1,1))
-
-
-
-### --- Step 6: Regress
-#
-# regress node FA value on behavior
+# L. Unc
+if(do_lunc == 1){
+  
+  ### --- Step 3: Model tract, no covariates
+  #
+  # Plot data, determine distribution,
+  #   compare model families.
+  #
+  # Plot best model
+  
+  # Get data
+  df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
+  df_afq$Group <- factor(df_afq$Group)
+  df_afq$Sex <- factor(df_afq$Sex)
+  
+  # Tract
+  tract <- "UNC_L"
+  df_tract <- df_afq[which(df_afq$tractID == tract), ]
+  df_tract$dti_fa <- round(df_tract$dti_fa, 3)
+  
+  # plot mean data
+  ggplot(data = df_tract) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  ggplot(data = df_tract) +
+    geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  # determine distribution
+  descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
+  
+  fit.beta <- fitdist(df_tract$dti_fa, "beta")
+  plot(fit.beta)
+  fit.beta$aic
+  
+  fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
+  plot(fit.gamma)
+  fit.gamma$aic
+  
+  
+  #  determine k, compare families
+  fit_gamma <- bam(dti_fa ~ Group +
+                     Sex +
+                     s(nodeID, by=Group, k=40) +
+                     s(subjectID, bs="re"),
+                   data = df_tract,
+                   family = Gamma(link = "logit"),
+                   method = "REML")
+  
+  gam.check(fit_gamma, rep = 500)
+  
+  fit_beta <- bam(dti_fa ~ Group +
+                    Sex +
+                    s(nodeID, by=Group, k=40) +
+                    s(subjectID, bs="re"),
+                  data = df_tract,
+                  family = betar(link = "logit"),
+                  method = "REML")
+  
+  gam.check(fit_beta, rep = 500)
+  
+  
+  infoMessages('on')
+  compareML(fit_gamma, fit_beta)  # fit_gamma recommended
+  
+  # get stats
+  summary(fit_gamma)  # R-sq = 0.819
+  
+  
+  
+  ### --- Step 4: Model tract, covariates
+  
+  fit_cov_pds <- bam(dti_fa ~ Group +
+                       Sex +
+                       s(nodeID, by=Group, k=40) +
+                       s(PDS, by=Sex) +
+                       s(subjectID, bs="re"),
+                     data = df_tract,
+                     family = Gamma(link = "logit"),
+                     method = "REML")
+  
+  gam.check(fit_cov_pds, rep = 500)
+  
+  
+  # Test cov model against gamma
+  # infoMessages('on')
+  compareML(fit_gamma, fit_cov_pds) # PDS wins
+  summary(fit_cov_pds) # R-sq = 0.85
+  
+  
+  # plot
+  df_pred <- predict.bam(
+    fit_cov_pds,
+    exclude_terms = c("PDS", "Sex","subjectID"),
+    values=list(PDS = NULL, Sex = NULL),
+    se.fit=T,
+    type="response")
+  
+  df_pred <- data.frame(Group=df_tract$Group,
+                        Sex=df_tract$Sex,
+                        subjectID=df_tract$subjectID,
+                        PDS=df_tract$PDS,
+                        nodeID=df_tract$nodeID,
+                        fit=df_pred$fit,
+                        se.fit=df_pred$se.fit)
+  
+  ggplot(data = df_pred) +
+    geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
+    ggtitle("GAM Fit of L. Uncinate FA Values") +
+    ylab("Fit FA")
+  
+  
+  
+  ### --- Step 5: Test for differences
+  #
+  # Check for group differences in spline.
+  
+  par(mfrow=c(1,3))
+  
+  p01 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "1")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-Med",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p02 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p12 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("1", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Med-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  par(mfrow=c(1,1))
+  
+  
+  
+  ### --- Step 6: Regress
+  #
+  # regress node FA value on behavior
+  
+  # find biggest difference
+  df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
+  colnames(df_est) <- colnames(p01)
+  df_est[,1:5] <- rbind(p01, p02, p12)
+  
+  ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
+  node_max <- df_est[ind_max,]$nodeID
+  h_groups <- df_est[ind_max,]$comp
+  groups <- stringr::str_extract_all(h_groups, "\\d+")
+  gA <- as.numeric(groups[[1]][1])
+  gB <- as.numeric(groups[[1]][2])
+  
+  # make df
+  df_max <- as.data.frame(df_tract[which(
+    df_tract$nodeID == node_max &
+      (df_tract$Group == gA | df_tract$Group == gB)
+  ),])
+  
+  # negLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+  # neuLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NeuLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NeuLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+}
 
-# find biggest difference
-df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
-colnames(df_est) <- colnames(p01)
-df_est[,1:5] <- rbind(p01, p02, p12)
 
-ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
-node_max <- df_est[ind_max,]$nodeID
-h_groups <- df_est[ind_max,]$comp
-groups <- stringr::str_extract_all(h_groups, "\\d+")
-gA <- as.numeric(groups[[1]][1])
-gB <- as.numeric(groups[[1]][2])
+# R. UNC
+if(do_runc == 1){
+  
+  # Get data
+  df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
+  df_afq$Group <- factor(df_afq$Group)
+  df_afq$Sex <- factor(df_afq$Sex)
+  
+  # Tract
+  tract <- "UNC_R"
+  df_tract <- df_afq[which(df_afq$tractID == tract), ]
+  df_tract$dti_fa <- round(df_tract$dti_fa, 3)
+  
+  # plot mean data
+  ggplot(data = df_tract) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  ggplot(data = df_tract) +
+    geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  # determine distribution
+  descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
+  
+  fit.beta <- fitdist(df_tract$dti_fa, "beta")
+  plot(fit.beta)
+  fit.beta$aic
+  
+  fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
+  plot(fit.gamma)
+  fit.gamma$aic
+  
+  
+  #  determine k, compare families
+  fit_gamma <- bam(dti_fa ~ Group +
+                     Sex +
+                     s(nodeID, by=Group, k=40) +
+                     s(subjectID, bs="re"),
+                   data = df_tract,
+                   family = Gamma(link = "logit"),
+                   method = "REML")
+  
+  gam.check(fit_gamma, rep = 500)
+  
+  fit_beta <- bam(dti_fa ~ Group +
+                    Sex +
+                    s(nodeID, by=Group, k=40) +
+                    s(subjectID, bs="re"),
+                  data = df_tract,
+                  family = betar(link = "logit"),
+                  method = "REML")
+  
+  gam.check(fit_beta, rep = 500)
+  
+  
+  infoMessages('on')
+  compareML(fit_gamma, fit_beta)  # fit_gamma recommended
+  
+  # get stats
+  summary(fit_gamma)  # R-sq = 0.89
+  
+  # covariates
+  fit_cov_pds <- bam(dti_fa ~ Group +
+                       Sex +
+                       s(nodeID, by=Group, k=40) +
+                       s(PDS, by=Sex) +
+                       s(subjectID, bs="re"),
+                     data = df_tract,
+                     family = Gamma(link = "logit"),
+                     method = "REML")
+  
+  gam.check(fit_cov_pds, rep = 500)
+  
+  
+  # Test cov model against gamma
+  # infoMessages('on')
+  compareML(fit_gamma, fit_cov_pds) # PDS wins
+  summary(fit_cov_pds) # R-sq = 0.91
+  
+  
+  # plot
+  df_pred <- predict.bam(
+    fit_cov_pds,
+    exclude_terms = c("PDS", "Sex","subjectID"),
+    values=list(PDS = NULL, Sex = NULL),
+    se.fit=T,
+    type="response")
+  
+  df_pred <- data.frame(Group=df_tract$Group,
+                        Sex=df_tract$Sex,
+                        subjectID=df_tract$subjectID,
+                        PDS=df_tract$PDS,
+                        nodeID=df_tract$nodeID,
+                        fit=df_pred$fit,
+                        se.fit=df_pred$se.fit)
+  
+  ggplot(data = df_pred) +
+    geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
+    ggtitle("GAM Fit of R. Uncinate FA Values") +
+    ylab("Fit FA")
+  
+  
+  # Check for group differences in spline.
+  
+  par(mfrow=c(1,3))
+  
+  p01 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "1")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-Med",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p02 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p12 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("1", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Med-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  par(mfrow=c(1,1))
+  
+  
+  
+  # find biggest difference
+  df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
+  colnames(df_est) <- colnames(p01)
+  df_est[,1:5] <- rbind(p01, p02, p12)
+  
+  ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
+  node_max <- df_est[ind_max,]$nodeID
+  h_groups <- df_est[ind_max,]$comp
+  groups <- stringr::str_extract_all(h_groups, "\\d+")
+  gA <- as.numeric(groups[[1]][1])
+  gB <- as.numeric(groups[[1]][2])
+  
+  # make df
+  df_max <- as.data.frame(df_tract[which(
+    df_tract$nodeID == node_max &
+      (df_tract$Group == gA | df_tract$Group == gB)
+  ),])
+  
+  # negLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+  # neuLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NeuLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NeuLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+}
 
-# make df
-df_max <- as.data.frame(df_tract[which(
-  df_tract$nodeID == node_max &
-  (df_tract$Group == gA | df_tract$Group == gB)
-),])
 
-# negLGI x group
-ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
-  geom_point() +
-  geom_smooth(method = "lm") +
-  facet_wrap(~ Group)
+# L. Cing
+if(do_lcing == 1){
+  
+  # Get data
+  df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
+  df_afq$Group <- factor(df_afq$Group)
+  df_afq$Sex <- factor(df_afq$Sex)
+  
+  # Tract
+  tract <- "CGC_L"
+  df_tract <- df_afq[which(df_afq$tractID == tract), ]
+  df_tract$dti_fa <- round(df_tract$dti_fa, 3)
+  
+  # plot mean data
+  ggplot(data = df_tract) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  ggplot(data = df_tract) +
+    geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  # determine distribution
+  descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
+  
+  fit.beta <- fitdist(df_tract$dti_fa, "beta")
+  plot(fit.beta)
+  fit.beta$aic
+  
+  fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
+  plot(fit.gamma)
+  fit.gamma$aic
+  
+  
+  #  determine k, compare families
+  fit_gamma <- bam(dti_fa ~ Group +
+                     Sex +
+                     s(nodeID, by=Group, k=40) +
+                     s(subjectID, bs="re"),
+                   data = df_tract,
+                   family = Gamma(link = "logit"),
+                   method = "REML")
+  
+  gam.check(fit_gamma, rep = 500)
+  
+  fit_beta <- bam(dti_fa ~ Group +
+                    Sex +
+                    s(nodeID, by=Group, k=40) +
+                    s(subjectID, bs="re"),
+                  data = df_tract,
+                  family = betar(link = "logit"),
+                  method = "REML")
+  
+  gam.check(fit_beta, rep = 500)
+  
+  
+  infoMessages('on')
+  compareML(fit_gamma, fit_beta)  # fit_gamma recommended
+  
+  # get stats
+  summary(fit_gamma)  # R-sq = 0.67
+  
+  # covariates
+  fit_cov_pds <- bam(dti_fa ~ Group +
+                       Sex +
+                       s(nodeID, by=Group, k=40) +
+                       s(PDS, by=Sex) +
+                       s(subjectID, bs="re"),
+                     data = df_tract,
+                     family = Gamma(link = "logit"),
+                     method = "REML")
+  
+  gam.check(fit_cov_pds, rep = 500)
+  
+  
+  # Test cov model against gamma
+  # infoMessages('on')
+  compareML(fit_gamma, fit_cov_pds) # PDS wins
+  summary(fit_cov_pds) # R-sq = 0.69
+  
+  
+  # plot
+  df_pred <- predict.bam(
+    fit_cov_pds,
+    exclude_terms = c("PDS", "Sex","subjectID"),
+    values=list(PDS = NULL, Sex = NULL),
+    se.fit=T,
+    type="response")
+  
+  df_pred <- data.frame(Group=df_tract$Group,
+                        Sex=df_tract$Sex,
+                        subjectID=df_tract$subjectID,
+                        PDS=df_tract$PDS,
+                        nodeID=df_tract$nodeID,
+                        fit=df_pred$fit,
+                        se.fit=df_pred$se.fit)
+  
+  ggplot(data = df_pred) +
+    geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
+    ggtitle("GAM Fit of L. Cingulate FA Values") +
+    ylab("Fit FA")
+  
+  
+  # Check for group differences in spline.
+  
+  par(mfrow=c(1,3))
+  
+  p01 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "1")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-Med",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p02 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p12 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("1", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Med-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  par(mfrow=c(1,1))
+  
+  
+  
+  # find biggest difference
+  df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
+  colnames(df_est) <- colnames(p01)
+  df_est[,1:5] <- rbind(p01, p02, p12)
+  
+  ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
+  node_max <- df_est[ind_max,]$nodeID
+  h_groups <- df_est[ind_max,]$comp
+  groups <- stringr::str_extract_all(h_groups, "\\d+")
+  gA <- as.numeric(groups[[1]][1])
+  gB <- as.numeric(groups[[1]][2])
+  
+  # make df
+  df_max <- as.data.frame(df_tract[which(
+    df_tract$nodeID == node_max &
+      (df_tract$Group == gA | df_tract$Group == gB)
+  ),])
+  
+  # negLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+  # neuLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NeuLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NeuLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+}
 
-fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
-summary(fit)
 
+# R. Cing
+if(do_lcing == 1){
+  
+  # Get data
+  df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
+  df_afq$Group <- factor(df_afq$Group)
+  df_afq$Sex <- factor(df_afq$Sex)
+  
+  # Tract
+  tract <- "CGC_R"
+  df_tract <- df_afq[which(df_afq$tractID == tract), ]
+  df_tract$dti_fa <- round(df_tract$dti_fa, 3)
+  
+  # plot mean data
+  ggplot(data = df_tract) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  ggplot(data = df_tract) +
+    geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  # determine distribution
+  descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
+  
+  fit.beta <- fitdist(df_tract$dti_fa, "beta")
+  plot(fit.beta)
+  fit.beta$aic
+  
+  fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
+  plot(fit.gamma)
+  fit.gamma$aic
+  
+  
+  #  determine k, compare families
+  fit_gamma <- bam(dti_fa ~ Group +
+                     Sex +
+                     s(nodeID, by=Group, k=40) +
+                     s(subjectID, bs="re"),
+                   data = df_tract,
+                   family = Gamma(link = "logit"),
+                   method = "REML")
+  
+  gam.check(fit_gamma, rep = 500)
+  
+  fit_beta <- bam(dti_fa ~ Group +
+                    Sex +
+                    s(nodeID, by=Group, k=40) +
+                    s(subjectID, bs="re"),
+                  data = df_tract,
+                  family = betar(link = "logit"),
+                  method = "REML")
+  
+  gam.check(fit_beta, rep = 500)
+  
+  
+  infoMessages('on')
+  compareML(fit_gamma, fit_beta)  # fit_gamma recommended
+  
+  # get stats
+  summary(fit_gamma)  # R-sq = 0.48
+  
+  # covariates
+  fit_cov_pds <- bam(dti_fa ~ Group +
+                       Sex +
+                       s(nodeID, by=Group, k=40) +
+                       s(PDS, by=Sex) +
+                       s(subjectID, bs="re"),
+                     data = df_tract,
+                     family = Gamma(link = "logit"),
+                     method = "REML")
+  
+  gam.check(fit_cov_pds, rep = 500)
+  
+  
+  # Test cov model against gamma
+  # infoMessages('on')
+  compareML(fit_gamma, fit_cov_pds) # PDS wins
+  summary(fit_cov_pds) # R-sq = 0.51
+  
+  
+  # plot
+  df_pred <- predict.bam(
+    fit_cov_pds,
+    exclude_terms = c("PDS", "Sex","subjectID"),
+    values=list(PDS = NULL, Sex = NULL),
+    se.fit=T,
+    type="response")
+  
+  df_pred <- data.frame(Group=df_tract$Group,
+                        Sex=df_tract$Sex,
+                        subjectID=df_tract$subjectID,
+                        PDS=df_tract$PDS,
+                        nodeID=df_tract$nodeID,
+                        fit=df_pred$fit,
+                        se.fit=df_pred$se.fit)
+  
+  ggplot(data = df_pred) +
+    geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
+    ggtitle("GAM Fit of R. Cingulate FA Values") +
+    ylab("Fit FA")
+  
+  
+  # Check for group differences in spline.
+  
+  par(mfrow=c(1,3))
+  
+  p01 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "1")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-Med",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p02 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p12 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("1", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Med-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  par(mfrow=c(1,1))
+  
+  
+  
+  # find biggest difference
+  df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
+  colnames(df_est) <- colnames(p01)
+  df_est[,1:5] <- rbind(p01, p02, p12)
+  
+  ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
+  node_max <- df_est[ind_max,]$nodeID
+  h_groups <- df_est[ind_max,]$comp
+  groups <- stringr::str_extract_all(h_groups, "\\d+")
+  gA <- as.numeric(groups[[1]][1])
+  gB <- as.numeric(groups[[1]][2])
+  
+  # make df
+  df_max <- as.data.frame(df_tract[which(
+    df_tract$nodeID == node_max &
+      (df_tract$Group == gA | df_tract$Group == gB)
+  ),])
+  
+  # negLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+  # neuLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NeuLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NeuLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+}
 
 
+# L. ATR
+if(do_latr == 1){
+  
+  # Get data
+  df_afq <- read.csv(paste0(dataDir, "Master_dataframe.csv"))
+  df_afq$Group <- factor(df_afq$Group)
+  df_afq$Sex <- factor(df_afq$Sex)
+  
+  # Tract
+  tract <- "ATR_L"
+  df_tract <- df_afq[which(df_afq$tractID == tract), ]
+  df_tract$dti_fa <- round(df_tract$dti_fa, 3)
+  
+  # plot mean data
+  ggplot(data = df_tract) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  ggplot(data = df_tract) +
+    geom_point(mapping = aes(x=nodeID, y=dti_fa, color=Group),size=0.3) +
+    geom_smooth(mapping = aes(x=nodeID, y=dti_fa, color=Group))
+  
+  # determine distribution
+  descdist(df_tract$dti_fa, discrete=F) # Could be beta or gamma
+  
+  fit.beta <- fitdist(df_tract$dti_fa, "beta")
+  plot(fit.beta)
+  fit.beta$aic
+  
+  fit.gamma <- fitdist(df_tract$dti_fa, "gamma")
+  plot(fit.gamma)
+  fit.gamma$aic
+  
+  
+  #  determine k, compare families
+  fit_gamma <- bam(dti_fa ~ Group +
+                     Sex +
+                     s(nodeID, by=Group, k=40) +
+                     s(subjectID, bs="re"),
+                   data = df_tract,
+                   family = Gamma(link = "logit"),
+                   method = "REML")
+  
+  gam.check(fit_gamma, rep = 500)
+  
+  fit_beta <- bam(dti_fa ~ Group +
+                    Sex +
+                    s(nodeID, by=Group, k=45) +
+                    s(subjectID, bs="re"),
+                  data = df_tract,
+                  family = betar(link = "logit"),
+                  method = "REML")
+  
+  gam.check(fit_beta, rep = 500)
+  
+  
+  infoMessages('on')
+  compareML(fit_gamma, fit_beta)  # fit_gamma recommended
+  
+  # get stats
+  summary(fit_gamma)  # R-sq = 0.864
+  
+  # covariates
+  fit_cov_pds <- bam(dti_fa ~ Group +
+                       Sex +
+                       s(nodeID, by=Group, k=40) +
+                       s(PDS, by=Sex) +
+                       s(subjectID, bs="re"),
+                     data = df_tract,
+                     family = Gamma(link = "logit"),
+                     method = "REML")
+  
+  gam.check(fit_cov_pds, rep = 500)
+  
+  
+  # Test cov model against gamma
+  # infoMessages('on')
+  compareML(fit_gamma, fit_cov_pds) # PDS wins
+  summary(fit_cov_pds) # R-sq = 0.864
+  
+  
+  # plot
+  df_pred <- predict.bam(
+    fit_cov_pds,
+    exclude_terms = c("PDS", "Sex","subjectID"),
+    values=list(PDS = NULL, Sex = NULL),
+    se.fit=T,
+    type="response")
+  
+  df_pred <- data.frame(Group=df_tract$Group,
+                        Sex=df_tract$Sex,
+                        subjectID=df_tract$subjectID,
+                        PDS=df_tract$PDS,
+                        nodeID=df_tract$nodeID,
+                        fit=df_pred$fit,
+                        se.fit=df_pred$se.fit)
+  
+  ggplot(data = df_pred) +
+    geom_smooth(mapping = aes(x=nodeID, y=fit, color=Group)) +
+    ggtitle("GAM Fit of R. Cingulate FA Values") +
+    ylab("Fit FA")
+  
+  
+  # Check for group differences in spline.
+  
+  par(mfrow=c(1,3))
+  
+  p01 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "1")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-Med",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p02 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("0", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Low-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  p12 <- plot_diff(fit_cov_pds,
+                   view="nodeID",
+                   comp=list(Group=c("1", "2")),
+                   rm.ranef = T,
+                   main = "Difference Scores, Med-High",
+                   ylab = "Est. FA difference",
+                   xlab = "Tract Node")
+  
+  par(mfrow=c(1,1))
+  
+  
+  
+  # find biggest difference
+  df_est <- as.data.frame(matrix(NA, nrow=3*dim(p01)[1], ncol=dim(p01)[2]))
+  colnames(df_est) <- colnames(p01)
+  df_est[,1:5] <- rbind(p01, p02, p12)
+  
+  ind_max <- which(abs(df_est$est) == max(abs(df_est$est)))
+  node_max <- df_est[ind_max,]$nodeID
+  h_groups <- df_est[ind_max,]$comp
+  groups <- stringr::str_extract_all(h_groups, "\\d+")
+  gA <- as.numeric(groups[[1]][1])
+  gB <- as.numeric(groups[[1]][2])
+  
+  # make df
+  df_max <- as.data.frame(df_tract[which(
+    df_tract$nodeID == node_max &
+      (df_tract$Group == gA | df_tract$Group == gB)
+  ),])
+  
+  # negLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NegLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NegLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+  # neuLGI x group
+  ggplot(df_max, aes(x=dti_fa, y=NeuLGI)) +
+    geom_point() +
+    geom_smooth(method = "lm") +
+    facet_wrap(~ Group)
+  
+  fit <- lmList(NeuLGI ~ dti_fa | Group, data = df_max)
+  summary(fit)
+  
+}