Load Packages

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

Read in r2star data

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

Set age and sex variables

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

Set striatal regions

rcaudate <- r2star$Right.Caudate
lcaudate <- r2star$Left.Caudate
caudate <- (r2star$Right.Caudate + r2star$Left.Caudate)/2

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

Set cognitive factor variables

f1ExecAcc <- r2star$F1_Exec_Comp_Res_Accuracy   
f2SocialAcc <- r2star$F2_Social_Cog_Accuracy    
f3MemoryAcc <- r2star$F3_Memory_Accuracy

f1CompEff <- r2star$F1_Complex_Reasoning_Efficiency 
f2MemoryEff <- r2star$F2_Memory.Efficiency  
f3ExecEff <- r2star$F3_Executive_Efficiency 
f4SocialEff <- r2star$F4_Social_Cognition_Efficiency

Striatal Models (for those found significant)

Caudate Accuracy

caudateF1ExecAcc <- gam(caudate ~ s(age) + sex + f1ExecAcc, data=r2star, method="REML")
visreg(caudateF1ExecAcc,"f1ExecAcc")

cau_F1ExecAcc_pval <- as.data.frame(summary(caudateF1ExecAcc)$p.table[3,4])
names(cau_F1ExecAcc_pval) <- "caudate_F1ExecAcc_pval"
print(cau_F1ExecAcc_pval)
##   caudate_F1ExecAcc_pval
## 1             0.02422407
caudateF2SocialAcc <- gam(caudate ~ s(age) + sex + f2SocialAcc, data=r2star, method="REML")
visreg(caudateF2SocialAcc,"f2SocialAcc")

cau_F2SocialAcc_pval <- as.data.frame(summary(caudateF2SocialAcc)$p.table[3,4])
names(cau_F2SocialAcc_pval) <- "caudate_F2SocialAcc_pval"
print(cau_F2SocialAcc_pval)
##   caudate_F2SocialAcc_pval
## 1               0.02803628
caudateF3MemoryAcc <- gam(caudate ~ s(age) + sex + f3MemoryAcc, data=r2star, method="REML")
visreg(caudateF3MemoryAcc,"f3MemoryAcc")

cau_F3MemoryAcc_pval <- as.data.frame(summary(caudateF3MemoryAcc)$p.table[3,4])
names(cau_F3MemoryAcc_pval) <- "caudate_F3MemoryAcc_pval"
print(cau_F3MemoryAcc_pval)
##   caudate_F3MemoryAcc_pval
## 1               0.03301555

Caudate Efficiency

caudatef1CompEff <- gam(caudate ~ s(age) + sex + f1CompEff, data=r2star, method="REML")
visreg(caudatef1CompEff,"f1CompEff")

cau_F1CompEff_pval <- as.data.frame(summary(caudatef1CompEff)$p.table[3,4])
names(cau_F1CompEff_pval) <- "caudate_F1CompEff_pval"
print(cau_F1CompEff_pval)
##   caudate_F1CompEff_pval
## 1              0.1691573
caudatef2MemoryEff <- gam(caudate ~ s(age) + sex + f2MemoryEff, data=r2star, method="REML")
visreg(caudatef2MemoryEff,"f2MemoryEff")

cau_F2MemoryEff_pval <- as.data.frame(summary(caudatef2MemoryEff)$p.table[3,4])
names(cau_F2MemoryEff_pval) <- "caudate_F2MemoryEff_pval"
print(cau_F2MemoryEff_pval)
##   caudate_F2MemoryEff_pval
## 1                0.0826403
caudatef3ExecEff <- gam(caudate ~ s(age) + sex + f3ExecEff, data=r2star, method="REML")
visreg(caudatef3ExecEff,"f3ExecEff")

cau_F3ExecEff_pval <- as.data.frame(summary(caudatef3ExecEff)$p.table[3,4])
names(cau_F3ExecEff_pval) <- "caudate_F3ExecEff_pval"
print(cau_F3ExecEff_pval)
##   caudate_F3ExecEff_pval
## 1             0.02708993
caudatef4SocialEff <- gam(caudate ~ s(age) + sex + f4SocialEff, data=r2star, method="REML")
visreg(caudatef4SocialEff,"f4SocialEff")

cau_F4SocialEff_pval <- as.data.frame(summary(caudatef4SocialEff)$p.table[3,4])
names(cau_F4SocialEff_pval) <- "caudate_F4SocialEff_pval"
print(cau_F4SocialEff_pval)
##   caudate_F4SocialEff_pval
## 1               0.04409771

Putamen Efficiency

putamenf1CompEff <- gam(putamen ~ s(age) + sex + f1CompEff, data=r2star, method="REML")
visreg(putamenf1CompEff,"f1CompEff")

put_F1CompEff_pval <- as.data.frame(summary(putamenf1CompEff)$p.table[3,4])
names(put_F1CompEff_pval) <- "putamen_F1CompEff_pval"
print(put_F1CompEff_pval)
##   putamen_F1CompEff_pval
## 1             0.04960271
putamenf2MemoryEff <- gam(putamen ~ s(age) + sex + f2MemoryEff, data=r2star, method="REML")
visreg(putamenf2MemoryEff,"f2MemoryEff")

put_F2MemoryEff_pval <- as.data.frame(summary(putamenf2MemoryEff)$p.table[3,4])
names(put_F2MemoryEff_pval) <- "putamen_F2MemoryEff_pval"
print(put_F2MemoryEff_pval)
##   putamen_F2MemoryEff_pval
## 1                0.1050839
putamenf3ExecEff <- gam(putamen ~ s(age) + sex + f3ExecEff, data=r2star, method="REML")
visreg(putamenf3ExecEff,"f3ExecEff")

put_F3ExecEff_pval <- as.data.frame(summary(putamenf3ExecEff)$p.table[3,4])
names(put_F3ExecEff_pval) <- "putamen_F3ExecEff_pval"
print(put_F3ExecEff_pval)
##   putamen_F3ExecEff_pval
## 1             0.06071977
putamenf4SocialEff <- gam(putamen ~ s(age) + sex + f4SocialEff, data=r2star, method="REML")
visreg(putamenf4SocialEff,"f4SocialEff")

put_F4SocialEff_pval <- as.data.frame(summary(putamenf4SocialEff)$p.table[3,4])
names(put_F4SocialEff_pval) <- "putamen_F4SocialEff_pval"
print(put_F4SocialEff_pval)
##   putamen_F4SocialEff_pval
## 1                0.1997631

P values

pvalues = list(cau_F1ExecAcc_pval, cau_F2SocialAcc_pval, cau_F3MemoryAcc_pval, cau_F1CompEff_pval, cau_F2MemoryEff_pval, cau_F3ExecEff_pval, cau_F4SocialEff_pval, put_F1CompEff_pval, put_F2MemoryEff_pval, put_F3ExecEff_pval, put_F4SocialEff_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
print(pvals)
##   rn caudate_F1ExecAcc_pval caudate_F2SocialAcc_pval
## 1  1             0.02422407               0.02803628
##   caudate_F3MemoryAcc_pval caudate_F1CompEff_pval caudate_F2MemoryEff_pval
## 1               0.03301555              0.1691573                0.0826403
##   caudate_F3ExecEff_pval caudate_F4SocialEff_pval putamen_F1CompEff_pval
## 1             0.02708993               0.04409771             0.04960271
##   putamen_F2MemoryEff_pval putamen_F3ExecEff_pval putamen_F4SocialEff_pval
## 1                0.1050839             0.06071977                0.1997631