Load Packages

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

Read in r2star & corrtraits 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 corrtraits variables

mood <- r2star$mood_corrtraitsv2
psychosis <- r2star$psychosis_corrtraitsv2
externalizing <- r2star$externalizing_corrtraitsv2
fear <- r2star$fear_corrtraitsv2

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.03504195
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.03869431
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.01216773

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.001823309
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.00744302
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.000870209

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.01181279
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.0364483
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.005234008

Fear

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

rput_fear_pval <- as.data.frame(summary(rputamenfear)$p.table[3,4])
names(rput_fear_pval) <- "rputamen_fear_pval"
print(rput_fear_pval)
##   rputamen_fear_pval
## 1        0.001718764
lputamenfear <- gam(lputamen ~ s(age) + sex + fear, data=r2star, method="REML")
visreg(lputamenfear,"fear")

lput_fear_pval <- as.data.frame(summary(lputamenfear)$p.table[3,4])
names(lput_fear_pval) <- "lputamen_fear_pval"
print(lput_fear_pval)
##   lputamen_fear_pval
## 1         0.01338387
putamenfear <- gam(putamen ~ s(age) + sex + fear, data=r2star, method="REML")
visreg(putamenfear,"fear")

put_fear_pval <- as.data.frame(summary(putamenfear)$p.table[3,4])
names(put_fear_pval) <- "putamen_fear_pval"
print(put_fear_pval)
##   putamen_fear_pval
## 1      0.0009817216

P values

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_fear_pval, lput_fear_pval, put_fear_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.03504195 0.041918831
## lputamen_mood_pval            0.03869431 0.041918831
## putamen_mood_pval             0.01216773 0.019332255
## rputamen_psychosis_pval      0.001823309 0.005925754
## lputamen_psychosis_pval       0.00744302 0.016126543
## putamen_psychosis_pval       0.000870209 0.005925754
## rputamen_externalizing_pval   0.01181279 0.019332255
## lputamen_externalizing_pval    0.0364483 0.041918831
## putamen_externalizing_pval   0.005234008 0.013608422
## rputamen_fear_pval           0.001718764 0.005925754
## lputamen_fear_pval            0.01338387 0.019332255
## putamen_fear_pval           0.0009817216 0.005925754