Load Packages

library(mgcv)
library(voxel)
library('knitr')
library("dplyr")
library(cowplot)

Read in r2star & group data

r2star <- read.csv("r2star_finalSample_03152017.csv")
groups <- read.csv("n1601_go1_diagnosis_dxpmr_20161014.csv")
r2star_groups <- merge(r2star, groups, by=c("bblid","scanid"))

# include only TD & PS subjects
r2star_groups2 <- subset(r2star_groups, select=c("bblid","scanid","sex","ageAtScan1","Right.Pallidum","Left.Pallidum","Right.Accumbens.Area","Left.Accumbens.Area","Right.Caudate","Left.Caudate","Right.Putamen","Left.Putamen","goassessDxpmr6"))
r2star_final <- subset(r2star_groups2, goassessDxpmr6 == "TD" | goassessDxpmr6 == "PS")

Set groups

td <- subset(r2star_final, goassessDxpmr6 == "TD")
ps <- subset(r2star_final, goassessDxpmr6 == "PS")

Average brain regions

# TD
wpallidum_td <- (td$Left.Pallidum + td$Right.Pallidum)/2
waccumbens_td <- (td$Left.Accumbens.Area + td$Right.Accumbens.Area)/2
wcaudate_td <- (td$Left.Caudate + td$Right.Caudate)/2
wputamen_td <- (td$Left.Putamen + td$Right.Putamen)/2

# PS
wpallidum_ps <- (ps$Left.Pallidum + ps$Right.Pallidum)/2
waccumbens_ps <- (ps$Left.Accumbens.Area + ps$Right.Accumbens.Area)/2
wcaudate_ps <- (ps$Left.Caudate + ps$Right.Caudate)/2
wputamen_ps <- (ps$Left.Putamen + ps$Right.Putamen)/2

Set age and sex variables

# TD
td$sex <- ordered(td$sex)
sex_td <- td$sex
age_td <- (td$ageAtScan1)/12

# PS
ps$sex <- ordered(ps$sex)
sex_ps <- ps$sex
age_ps <- (ps$ageAtScan1)/12

Group Models

Pallidum

# right pallidum
rpallidum_td <- gam(Right.Pallidum ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(rpallidum_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

rpall_pval_td <- as.data.frame(summary(rpallidum_td)$s.table[,4])
colnames(rpall_pval_td) <- "rpall_pval_td"
print(rpall_pval_td)
##   rpall_pval_td
## 1  3.534849e-29
rpallidum_ps <- gam(Right.Pallidum ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(rpallidum_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

rpall_pval_ps <- as.data.frame(summary(rpallidum_ps)$s.table[,4])
colnames(rpall_pval_ps) <- "rpall_pval_ps"
print(rpall_pval_ps)
##   rpall_pval_ps
## 1  9.104527e-15
# left pallidum
lpallidum_td <- gam(Left.Pallidum ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(lpallidum_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

lpall_pval_td <- as.data.frame(summary(lpallidum_td)$s.table[,4])
colnames(lpall_pval_td) <- "lpall_pval_td"
print(lpall_pval_td)
##   lpall_pval_td
## 1  2.586817e-25
lpallidum_ps <- gam(Left.Pallidum ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(lpallidum_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

lpall_pval_ps <- as.data.frame(summary(lpallidum_ps)$s.table[,4])
colnames(lpall_pval_ps) <- "lpall_pval_ps"
print(lpall_pval_ps)
##   lpall_pval_ps
## 1  5.521619e-15
# whole pallidum
pallidum_td <- gam(wpallidum_td ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(pallidum_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

pall_pval_td <- as.data.frame(summary(pallidum_td)$s.table[,4])
colnames(pall_pval_td) <- "pall_pval_td"
print(pall_pval_td)
##   pall_pval_td
## 1 4.064404e-32
pallidum_ps <- gam(wpallidum_ps ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(pallidum_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

pall_pval_ps <- as.data.frame(summary(pallidum_ps)$s.table[,4])
colnames(pall_pval_ps) <- "pall_pval_ps"
print(pall_pval_ps)
##   pall_pval_ps
## 1 2.597176e-18

Accumbens

# right accumbens
raccumbens_td <- gam(Right.Accumbens.Area ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(raccumbens_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

racc_pval_td <- as.data.frame(summary(raccumbens_td)$s.table[,4])
colnames(racc_pval_td) <- "racc_pval_td"
print(racc_pval_td)
##   racc_pval_td
## 1  8.87645e-08
raccumbens_ps <- gam(Right.Accumbens.Area ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(raccumbens_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

racc_pval_ps <- as.data.frame(summary(raccumbens_ps)$s.table[,4])
colnames(racc_pval_ps) <- "racc_pval_ps"
print(racc_pval_ps)
##   racc_pval_ps
## 1 0.0005362273
#left accumbens
laccumbens_td <- gam(Left.Accumbens.Area ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(laccumbens_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

lacc_pval_td <- as.data.frame(summary(laccumbens_td)$s.table[,4])
colnames(lacc_pval_td) <- "lacc_pval_td"
print(lacc_pval_td)
##   lacc_pval_td
## 1 9.421719e-08
laccumbens_ps <- gam(Left.Accumbens.Area ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(laccumbens_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

lacc_pval_ps <- as.data.frame(summary(laccumbens_ps)$s.table[,4])
colnames(lacc_pval_ps) <- "lacc_pval_ps"
print(lacc_pval_ps)
##   lacc_pval_ps
## 1    0.0157165
# whole accumbens
accumbens_td <- gam(waccumbens_td ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(accumbens_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

acc_pval_td <- as.data.frame(summary(accumbens_td)$s.table[,4])
colnames(acc_pval_td) <- "acc_pval_td"
print(acc_pval_td)
##    acc_pval_td
## 1 8.217211e-10
accumbens_ps <- gam(waccumbens_ps ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(accumbens_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

acc_pval_ps <- as.data.frame(summary(accumbens_ps)$s.table[,4])
colnames(acc_pval_ps) <- "acc_pval_ps"
print(acc_pval_ps)
##   acc_pval_ps
## 1 0.001440166

Caudate

# right caudate
rcaudate_td <- gam(Right.Caudate ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(rcaudate_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

rcau_pval_td <- as.data.frame(summary(rcaudate_td)$s.table[,4])
colnames(rcau_pval_td) <- "rcau_pval_td"
print(rcau_pval_td)
##   rcau_pval_td
## 1 5.204701e-07
rcaudate_ps <- gam(Right.Caudate ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(rcaudate_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

rcau_pval_ps <- as.data.frame(summary(rcaudate_ps)$s.table[,4])
colnames(rcau_pval_ps) <- "rcau_pval_ps"
print(rcau_pval_ps)
##   rcau_pval_ps
## 1 2.426957e-05
# left caudate
lcaudate_td <- gam(Left.Caudate ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(lcaudate_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

lcau_pval_td <- as.data.frame(summary(lcaudate_td)$s.table[,4])
colnames(lcau_pval_td) <- "lcau_pval_td"
print(lcau_pval_td)
##   lcau_pval_td
## 1 0.0002420983
lcaudate_ps <- gam(Left.Caudate ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(lcaudate_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

lcau_pval_ps <- as.data.frame(summary(lcaudate_ps)$s.table[,4])
colnames(lcau_pval_ps) <- "lcau_pval_ps"
print(lcau_pval_ps)
##   lcau_pval_ps
## 1  0.007110724
# whole caudate
caudate_td <- gam(wcaudate_td ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(caudate_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

cau_pval_td <- as.data.frame(summary(caudate_td)$s.table[,4])
colnames(cau_pval_td) <- "cau_pval_td"
print(cau_pval_td)
##    cau_pval_td
## 1 2.079886e-06
caudate_ps <- gam(wcaudate_ps ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(caudate_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

cau_pval_ps <- as.data.frame(summary(caudate_ps)$s.table[,4])
colnames(cau_pval_ps) <- "cau_pval_ps"
print(cau_pval_ps)
##   cau_pval_ps
## 1 5.52171e-05

Putamen

# right putamen
rputamen_td <- gam(Right.Putamen ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(rputamen_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

rput_pval_td <- as.data.frame(summary(rputamen_td)$s.table[,4])
colnames(rput_pval_td) <- "rput_pval_td"
print(rput_pval_td)
##   rput_pval_td
## 1  2.98818e-13
rputamen_ps <- gam(Right.Putamen ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(rputamen_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

rput_pval_ps <- as.data.frame(summary(rputamen_ps)$s.table[,4])
colnames(rput_pval_ps) <- "rput_pval_ps"
print(rput_pval_ps)
##   rput_pval_ps
## 1 4.080813e-07
# left putamen
lputamen_td <- gam(Left.Putamen ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(lputamen_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

lput_pval_td <- as.data.frame(summary(lputamen_td)$s.table[,4])
colnames(lput_pval_td) <- "lput_pval_td"
print(lput_pval_td)
##   lput_pval_td
## 1 5.915704e-16
lputamen_ps <- gam(Left.Putamen ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(lputamen_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

lput_pval_ps <- as.data.frame(summary(lputamen_ps)$s.table[,4])
colnames(lput_pval_ps) <- "lput_pval_ps"
print(lput_pval_ps)
##   lput_pval_ps
## 1 1.461125e-06
# whole putamen
putamen_td <- gam(wputamen_td ~ s(age_td,k=4) + sex_td, data=td, method="REML")
plotGAM(putamen_td,"age_td","sex_td")
## [1] "Working with a GAM Model"

put_pval_td <- as.data.frame(summary(putamen_td)$s.table[,4])
colnames(put_pval_td) <- "put_pval_td"
print(put_pval_td)
##    put_pval_td
## 1 2.043481e-17
putamen_ps <- gam(wputamen_ps ~ s(age_ps,k=4) + sex_ps, data=ps, method="REML")
plotGAM(putamen_ps,"age_ps","sex_ps")
## [1] "Working with a GAM Model"

put_pval_ps <- as.data.frame(summary.gam(putamen_ps)$s.table[,4])
colnames(put_pval_ps) <- "put_pval_ps"
print(put_pval_ps)
##    put_pval_ps
## 1 2.434185e-08

P values

# TD p values
pvalues_td = list(rpall_pval_td, lpall_pval_td, pall_pval_td, racc_pval_td, lacc_pval_td, acc_pval_td, rcau_pval_td, lcau_pval_td, cau_pval_td, rput_pval_td, lput_pval_td, put_pval_td)
pvals_td <- as.data.frame(pvalues_td)
pvals_td_t <- t(pvals_td)

fdr_td <- as.data.frame(p.adjust(pvals_td[1,], "fdr"))
names(fdr_td) <- "fdr_td"
pvals_fdr_td <- as.data.frame(cbind(pvals_td_t, fdr_td))
pvals_fdr_td <- t(pvals_fdr_td)

# PS p values
pvalues_ps = list(rpall_pval_ps, lpall_pval_ps, pall_pval_ps, racc_pval_ps, lacc_pval_ps, acc_pval_ps, rcau_pval_ps, lcau_pval_ps, cau_pval_ps, rput_pval_ps, lput_pval_ps, put_pval_ps)
pvals_ps <- as.data.frame(pvalues_ps)
pvals_ps_t <- t(pvals_ps)

fdr_ps <- as.data.frame(p.adjust(pvals_ps[1,], "fdr"))
names(fdr_ps) <- "fdr_ps"
pvals_fdr_ps <- as.data.frame(cbind(pvals_ps_t, fdr_ps))
pvals_fdr_ps <- t(pvals_fdr_ps)

# Combine p values
colnames(pvals_fdr_td) <- c("rpall_pval", "lpall_pval", "pall_pval", "racc_pval", "lacc_pval", "acc_pval", "rcau_pval", "lcau_pval", "cau_pval", "rput_pval", "lput_pval", "put_pval")
colnames(pvals_fdr_ps) <- c("rpall_pval", "lpall_pval", "pall_pval", "racc_pval", "lacc_pval", "acc_pval", "rcau_pval", "lcau_pval", "cau_pval", "rput_pval", "lput_pval", "put_pval")
pvals_fdr_td <- as.data.frame(pvals_fdr_td)
pvals_fdr_ps <- as.data.frame(pvals_fdr_ps)

all_pvals <- rbind(pvals_fdr_td, pvals_fdr_ps)

final_pvals <- all_pvals[c("pvals_td_t", "pvals_ps_t", "fdr_td", "fdr_ps"),]
print(final_pvals)
##              rpall_pval   lpall_pval    pall_pval    racc_pval
## pvals_td_t 3.534849e-29 2.586817e-25 4.064404e-32 8.876450e-08
## pvals_ps_t 9.104527e-15 5.521619e-15 2.597176e-18 5.362273e-04
## fdr_td     2.120909e-28 1.034727e-24 4.877285e-31 1.256229e-07
## fdr_ps     3.641811e-14 3.312972e-14 3.116611e-17 7.149697e-04
##               lacc_pval     acc_pval    rcau_pval    lcau_pval
## pvals_td_t 9.421719e-08 8.217211e-10 5.204701e-07 0.0002420983
## pvals_ps_t 1.571650e-02 1.440166e-03 2.426957e-05 0.0071107241
## fdr_td     1.256229e-07 1.408665e-09 6.245641e-07 0.0002420983
## fdr_ps     1.571650e-02 1.728200e-03 4.160497e-05 0.0077571536
##                cau_pval    rput_pval    lput_pval     put_pval
## pvals_td_t 2.079886e-06 2.988180e-13 5.915704e-16 2.043481e-17
## pvals_ps_t 5.521710e-05 4.080813e-07 1.461125e-06 2.434185e-08
## fdr_td     2.268966e-06 5.976360e-13 1.419769e-15 6.130443e-17
## fdr_ps     8.282564e-05 9.793952e-07 2.922251e-06 7.302556e-08