R2Star Group Factor Models

Load Packages

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

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"))

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

# include only TD & PS subjects
r2star_final <- subset(r2star_groups_psychosis, goassessDxpmr6 == "TD" | goassessDxpmr6 == "PS")

r2star_final <- r2star_final[complete.cases(r2star_final$Overall_Psychosis),]

Average brain regions

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

Set age and sex variables

r2star_final$sex <- ordered(r2star_final$sex)
sex <- r2star_final$sex

age <- (r2star_final$ageAtScan1)/12

psychosis <- r2star_final$Overall_Psychosis

Striatal Models

Pallidum

pallidum <- gam(wpallidum ~ s(age) + sex + psychosis, data=r2star_final, method="REML")
visreg(pallidum,"psychosis")

pall_pval <- as.data.frame(summary(pallidum)$p.table[3,4])
names(pall_pval) <- "pallidum_pval"
print(pall_pval)

##   pallidum_pval
## 1      0.321141

Accumbens

accumbens <- gam(waccumbens ~ s(age) + sex + psychosis, data=r2star_final, method="REML")
visreg(accumbens,"psychosis")

acc_pval <- as.data.frame(summary(accumbens)$p.table[3,4])
names(acc_pval) <- "accumbens_pval"
print(acc_pval)

##   accumbens_pval
## 1      0.3299606

Caudate

caudate <- gam(wcaudate ~ s(age) + sex + psychosis, data=r2star_final, method="REML")
visreg(caudate,"psychosis")

cau_pval <- as.data.frame(summary(caudate)$p.table[3,4])
names(cau_pval) <- "caudate_pval"
print(cau_pval)

##   caudate_pval
## 1    0.5899932

Putamen

putamen <- gam(wputamen ~ s(age) + sex + psychosis, data=r2star_final, method="REML")
visreg(putamen,"psychosis")

put_pval <- as.data.frame(summary(putamen)$p.table[3,4])
names(put_pval) <- "putamen_pval"
print(put_pval)

##   putamen_pval
## 1   0.06777008

P values

pvalues = list(pall_pval, acc_pval, cau_pval, put_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
print(pvals)

##   rn pallidum_pval accumbens_pval caudate_pval putamen_pval
## 1  1      0.321141      0.3299606    0.5899932   0.06777008