Load Packages

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

Read in r2star, group, and psychosis 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"))

psych <- read.csv("n1601_goassess_itemwise_psychosis_JCPP_factor_scores_20170331.csv")
r2star_group_psychosis <- merge(r2star_groups, psych, by=c("bblid","scanid"))

# include only PS subjects
ps <- subset(r2star_group_psychosis, goassessDxpmr6 == "PS")

Average brain regions

wpallidum <- (ps$Left.Pallidum + ps$Right.Pallidum)/2
waccumbens <- (ps$Left.Accumbens.Area + ps$Right.Accumbens.Area)/2
wcaudate <- (ps$Left.Caudate + ps$Right.Caudate)/2
wputamen <- (ps$Left.Putamen + ps$Right.Putamen)/2

Set age, sex, and psychosis variables

ps$sex <- ordered(ps$sex)
sex <- ps$sex

age <- (ps$ageAtScan1)/12

psychosis <- ps$Overall_Psychosis

Regional

Pallidum

# right pallidum
rpallidum <- gam(Right.Pallidum ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(rpallidum,"psychosis")

rpall_pval <- as.data.frame(summary(rpallidum)$p.table[3,4])
colnames(rpall_pval) <- "rpall_pval"
print(rpall_pval)
##   rpall_pval
## 1  0.9582044
# left pallidum
lpallidum <- gam(Left.Pallidum ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(lpallidum,"psychosis")

lpall_pval <- as.data.frame(summary(lpallidum)$p.table[3,4])
colnames(lpall_pval) <- "lpall_pval"
print(lpall_pval)
##   lpall_pval
## 1  0.4483615
# whole pallidum
pallidum <- gam(wpallidum ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(pallidum,"psychosis")

pall_pval <- as.data.frame(summary(pallidum)$p.table[3,4])
colnames(pall_pval) <- "pall_pval"
print(pall_pval)
##   pall_pval
## 1 0.9997109

Accumbens

# right accumbens
raccumbens <- gam(Right.Accumbens.Area ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(raccumbens,"psychosis")

racc_pval <- as.data.frame(summary(raccumbens)$p.table[3,4])
colnames(racc_pval) <- "racc_pval"
print(racc_pval)
##   racc_pval
## 1 0.7398627
# left accumbens
laccumbens <- gam(Left.Accumbens.Area ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(laccumbens,"psychosis")

lacc_pval <- as.data.frame(summary(laccumbens)$p.table[3,4])
colnames(lacc_pval) <- "lacc_pval"
print(lacc_pval)
##   lacc_pval
## 1 0.7105031
# whole accumbens
accumbens <- gam(waccumbens ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(accumbens,"psychosis")

acc_pval <- as.data.frame(summary(accumbens)$p.table[3,4])
colnames(acc_pval) <- "acc_pval"
print(acc_pval)
##    acc_pval
## 1 0.7153477

Caudate

# right caudate
rcaudate <- gam(Right.Caudate ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(rcaudate,"psychosis")

rcau_pval <- as.data.frame(summary(rcaudate)$p.table[3,4])
colnames(rcau_pval) <- "rcau_pval"
print(rcau_pval)
##   rcau_pval
## 1 0.3145495
# left caudate
lcaudate <- gam(Left.Caudate ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(lcaudate,"psychosis")

lcau_pval <- as.data.frame(summary(lcaudate)$p.table[3,4])
colnames(lcau_pval) <- "lcau_pval"
print(lcau_pval)
##   lcau_pval
## 1 0.1638047
# whole caudate
caudate <- gam(wcaudate ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(caudate,"psychosis")

cau_pval <- as.data.frame(summary(caudate)$p.table[3,4])
colnames(cau_pval) <- "cau_pval"
print(cau_pval)
##    cau_pval
## 1 0.1524965

Putamen

# right putamen
rputamen <- gam(Right.Putamen ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(rputamen,"psychosis")

rput_pval <- as.data.frame(summary(rputamen)$p.table[3,4])
colnames(rput_pval) <- "rput_pval"
print(rput_pval)
##   rput_pval
## 1 0.9696526
# left putamen
lputamen <- gam(Left.Putamen ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(lputamen,"psychosis")

lput_pval <- as.data.frame(summary(lputamen)$p.table[3,4])
colnames(lput_pval) <- "lput_pval"
print(lput_pval)
##   lput_pval
## 1 0.9106338
# whole putamen
putamen <- gam(wputamen ~ s(age,k=4) + sex + psychosis, data=ps, method="REML")
visreg(putamen,"psychosis")

put_pval <- as.data.frame(summary(putamen)$p.table[3,4])
colnames(put_pval) <- "put_pval"
print(put_pval)
##    put_pval
## 1 0.9045056

P Values

pvalues = list(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 <- 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
## rpall_pval 0.9582044   1
## lpall_pval 0.4483615   1
## pall_pval  0.9997109   1
## racc_pval  0.7398627   1
## lacc_pval  0.7105031   1
## acc_pval   0.7153477   1
## rcau_pval  0.3145495   1
## lcau_pval  0.1638047   1
## cau_pval   0.1524965   1
## rput_pval  0.9696526   1
## lput_pval  0.9106338   1
## put_pval   0.9045056   1