Load Packages

library(mgcv)
library(voxel)
library('knitr')
library(visreg)

Read in r2star & bifactor data

r2star <- read.csv("r2star_finalSample_03152017.csv")

Set age and sex variables

r2star$sex <- ordered(r2star$sex)
sex <- r2star$sex
age <- (r2star$ageAtScan1)/12

Set putamen

rputamen <- r2star$Right.Putamen
lputamen <- r2star$Left.Putamen
putamen <- (r2star$Right.Putamen + r2star$Left.Putamen)/2

Set bifactor variables

mood <- r2star$mood_4factorv2
psychosis <- r2star$psychosis_4factorv2
externalizing <- r2star$externalizing_4factorv2
phobia <- r2star$phobias_4factorv2
overall <- r2star$overall_psychopathology_4factorv2

Bifactor Models

Mood

rputamenMood <- gam(rputamen ~ s(age) + sex + mood, data=r2star, method="REML")
visreg(rputamenMood,"mood")

rput_mood_pval <- as.data.frame(summary(rputamenMood)$p.table[3,4])
names(rput_mood_pval) <- "rputamen_mood_pval"
print(rput_mood_pval)
##   rputamen_mood_pval
## 1         0.08330902
lputamenMood <- gam(lputamen ~ s(age) + sex + mood, data=r2star, method="REML")
visreg(lputamenMood,"mood")

lput_mood_pval <- as.data.frame(summary(lputamenMood)$p.table[3,4])
names(lput_mood_pval) <- "lputamen_mood_pval"
print(lput_mood_pval)
##   lputamen_mood_pval
## 1          0.7541309
putamenMood <- gam(putamen ~ s(age) + sex + mood, data=r2star, method="REML")
visreg(putamenMood,"mood")

put_mood_pval <- as.data.frame(summary(putamenMood)$p.table[3,4])
names(put_mood_pval) <- "putamen_mood_pval"
print(put_mood_pval)
##   putamen_mood_pval
## 1         0.2840296

Psychosis

rputamenpsychosis <- gam(rputamen ~ s(age) + sex + psychosis, data=r2star, method="REML")
visreg(rputamenpsychosis,"psychosis")

rput_psychosis_pval <- as.data.frame(summary(rputamenpsychosis)$p.table[3,4])
names(rput_psychosis_pval) <- "rputamen_psychosis_pval"
print(rput_psychosis_pval)
##   rputamen_psychosis_pval
## 1               0.5923205
lputamenpsychosis <- gam(lputamen ~ s(age) + sex + psychosis, data=r2star, method="REML")
visreg(lputamenpsychosis,"psychosis")

lput_psychosis_pval <- as.data.frame(summary(lputamenpsychosis)$p.table[3,4])
names(lput_psychosis_pval) <- "lputamen_psychosis_pval"
print(lput_psychosis_pval)
##   lputamen_psychosis_pval
## 1               0.6603105
putamenpsychosis <- gam(putamen ~ s(age) + sex + psychosis, data=r2star, method="REML")
visreg(putamenpsychosis,"psychosis")

put_psychosis_pval <- as.data.frame(summary(putamenpsychosis)$p.table[3,4])
names(put_psychosis_pval) <- "putamen_psychosis_pval"
print(put_psychosis_pval)
##   putamen_psychosis_pval
## 1               0.887016

Externalizing

rputamenexternalizing <- gam(rputamen ~ s(age) + sex + externalizing, data=r2star, method="REML")
visreg(rputamenexternalizing,"externalizing")

rput_externalizing_pval <- as.data.frame(summary(rputamenexternalizing)$p.table[3,4])
names(rput_externalizing_pval) <- "rputamen_externalizing_pval"
print(rput_externalizing_pval)
##   rputamen_externalizing_pval
## 1                   0.3301358
lputamenexternalizing <- gam(lputamen ~ s(age) + sex + externalizing, data=r2star, method="REML")
visreg(lputamenexternalizing,"externalizing")

lput_externalizing_pval <- as.data.frame(summary(lputamenexternalizing)$p.table[3,4])
names(lput_externalizing_pval) <- "lputamen_externalizing_pval"
print(lput_externalizing_pval)
##   lputamen_externalizing_pval
## 1                   0.3871073
putamenexternalizing <- gam(putamen ~ s(age) + sex + externalizing, data=r2star, method="REML")
visreg(putamenexternalizing,"externalizing")

put_externalizing_pval <- as.data.frame(summary(putamenexternalizing)$p.table[3,4])
names(put_externalizing_pval) <- "putamen_externalizing_pval"
print(put_externalizing_pval)
##   putamen_externalizing_pval
## 1                  0.2309453

Phobia

rputamenphobia <- gam(rputamen ~ s(age) + sex + phobia, data=r2star, method="REML")
visreg(rputamenphobia,"phobia")

rput_phobia_pval <- as.data.frame(summary(rputamenphobia)$p.table[3,4])
names(rput_phobia_pval) <- "rputamen_phobia_pval"
print(rput_phobia_pval)
##   rputamen_phobia_pval
## 1           0.04490671
lputamenphobia <- gam(lputamen ~ s(age) + sex + phobia, data=r2star, method="REML")
visreg(lputamenphobia,"phobia")

lput_phobia_pval <- as.data.frame(summary(lputamenphobia)$p.table[3,4])
names(lput_phobia_pval) <- "lputamen_phobia_pval"
print(lput_phobia_pval)
##   lputamen_phobia_pval
## 1            0.1754955
putamenphobia <- gam(putamen ~ s(age) + sex + phobia, data=r2star, method="REML")
visreg(putamenphobia,"phobia")

put_phobia_pval <- as.data.frame(summary(putamenphobia)$p.table[3,4])
names(put_phobia_pval) <- "putamen_phobia_pval"
print(put_phobia_pval)
##   putamen_phobia_pval
## 1          0.05190132

Overall Psychopathology

rputamenoverall <- gam(rputamen ~ s(age) + sex + overall, data=r2star, method="REML")
visreg(rputamenoverall,"overall")

rput_overall_pval <- as.data.frame(summary(rputamenoverall)$p.table[3,4])
names(rput_overall_pval) <- "rputamen_overall_pval"
print(rput_overall_pval)
##   rputamen_overall_pval
## 1           0.001989048
lputamenoverall <- gam(lputamen ~ s(age) + sex + overall, data=r2star, method="REML")
visreg(lputamenoverall,"overall")

lput_overall_pval <- as.data.frame(summary(lputamenoverall)$p.table[3,4])
names(lput_overall_pval) <- "lputamen_overall_pval"
print(lput_overall_pval)
##   lputamen_overall_pval
## 1            0.01281675
putamenoverall <- gam(putamen ~ s(age) + sex + overall, data=r2star, method="REML")
visreg(putamenoverall,"overall")

put_overall_pval <- as.data.frame(summary(putamenoverall)$p.table[3,4])
names(put_overall_pval) <- "putamen_overall_pval"
print(put_overall_pval)
##   putamen_overall_pval
## 1          0.001123163

P values & FDR Corrections

pvalues = list(rput_mood_pval, lput_mood_pval, put_mood_pval, rput_psychosis_pval, lput_psychosis_pval, put_psychosis_pval, rput_externalizing_pval, lput_externalizing_pval, put_externalizing_pval, rput_phobia_pval, lput_phobia_pval, put_phobia_pval, rput_overall_pval, lput_overall_pval, put_overall_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
pvalst <- t(pvals)

# fdr corrections
fdr <- as.data.frame(p.adjust(pvals[1,], "fdr"))
names(fdr) <- "fdr"
pvals_fdr <- as.data.frame(cbind(pvalst, fdr))
pvals_fdr <- pvals_fdr[-c(1), ]
print(pvals_fdr)
##                                  pvalst        fdr
## rputamen_mood_pval           0.08330902 0.22215740
## lputamen_mood_pval            0.7541309 0.86186385
## putamen_mood_pval             0.2840296 0.50494148
## rputamen_psychosis_pval       0.5923205 0.78976066
## lputamen_psychosis_pval       0.6603105 0.81268982
## putamen_psychosis_pval         0.887016 0.94615039
## rputamen_externalizing_pval   0.3301358 0.52821723
## lputamen_externalizing_pval   0.3871073 0.56306511
## putamen_externalizing_pval    0.2309453 0.46189060
## rputamen_phobia_pval         0.04490671 0.16608424
## lputamen_phobia_pval          0.1754955 0.40113251
## putamen_phobia_pval          0.05190132 0.16608424
## rputamen_overall_pval       0.001989048 0.01591238
## lputamen_overall_pval        0.01281675 0.06835601
## putamen_overall_pval        0.001123163 0.01591238