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

rpallidum <- r2star$Right.Pallidum
lpallidum <- r2star$Left.Pallidum
pallidum <- (r2star$Right.Pallidum + r2star$Left.Pallidum)/2

raccumbems <- r2star$Right.Accumbens.Area
laccumbens <- r2star$Left.Accumbens.Area
accumbens <- (r2star$Right.Accumbens.Area + r2star$Left.Accumbens.Area)/2

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 variables

accuracy <- r2star$Overall_Accuracy
efficiency <- r2star$Overall_Efficiency
speed <- r2star$Overall_Speed

Striatal Models

Pallidum

pallidumAcc <- gam(pallidum ~ s(age) + sex + accuracy, data=r2star, method="REML")
visreg(pallidumAcc,"accuracy")

pall_acc_pval <- as.data.frame(summary(pallidumAcc)$p.table[3,4])
names(pall_acc_pval) <- "pallidum_acc_pval"
print(pall_acc_pval)
##   pallidum_acc_pval
## 1         0.2385601
pallidumEff <- gam(pallidum ~ s(age) + sex + efficiency, data=r2star, method="REML")
visreg(pallidumEff,"efficiency")

pall_eff_pval <- as.data.frame(summary(pallidumEff)$p.table[3,4])
names(pall_eff_pval) <- "pallidum_eff_pval"
print(pall_eff_pval)
##   pallidum_eff_pval
## 1         0.1213092
pallidumSpeed <- gam(pallidum ~ s(age) + sex + speed, data=r2star, method="REML")
visreg(pallidumSpeed,"speed")

pall_speed_pval <- as.data.frame(summary(pallidumSpeed)$p.table[3,4])
names(pall_speed_pval) <- "pallidum_speed_pval"
print(pall_speed_pval)
##   pallidum_speed_pval
## 1           0.1889191

Accumbens

accumbensAcc <- gam(accumbens ~ s(age) + sex + accuracy, data=r2star, method="REML")
visreg(accumbensAcc,"accuracy")

acc_acc_pval <- as.data.frame(summary(accumbensAcc)$p.table[3,4])
names(acc_acc_pval) <- "accumbens_acc_pval"
print(acc_acc_pval)
##   accumbens_acc_pval
## 1          0.8317156
accumbensEff <- gam(accumbens ~ s(age) + sex + efficiency, data=r2star, method="REML")
visreg(accumbensEff,"efficiency")

acc_eff_pval <- as.data.frame(summary(accumbensEff)$p.table[3,4])
names(acc_eff_pval) <- "accumbens_eff_pval"
print(acc_eff_pval)
##   accumbens_eff_pval
## 1           0.895044
accumbensSpeed <- gam(accumbens ~ s(age) + sex + speed, data=r2star, method="REML")
visreg(accumbensSpeed,"speed")

acc_speed_pval <- as.data.frame(summary(accumbensSpeed)$p.table[3,4])
names(acc_speed_pval) <- "accumbens_speed_pval"
print(acc_speed_pval)
##   accumbens_speed_pval
## 1            0.6954426

Caudate

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

cau_acc_pval <- as.data.frame(summary(caudateAcc)$p.table[3,4])
names(cau_acc_pval) <- "caudate_acc_pval"
print(cau_acc_pval)
##   caudate_acc_pval
## 1      0.006608071
caudateEff <- gam(caudate ~ s(age) + sex + efficiency, data=r2star, method="REML")
visreg(caudateEff,"efficiency")

cau_eff_pval <- as.data.frame(summary(caudateEff)$p.table[3,4])
names(cau_eff_pval) <- "caudate_eff_pval"
print(cau_eff_pval)
##   caudate_eff_pval
## 1       0.02509784
caudateSpeed <- gam(caudate ~ s(age) + sex + speed, data=r2star, method="REML")
visreg(caudateSpeed,"speed")

cau_speed_pval <- as.data.frame(summary(caudateSpeed)$p.table[3,4])
names(cau_speed_pval) <- "caudate_speed_pval"
print(cau_speed_pval)
##   caudate_speed_pval
## 1          0.5452894

Putamen

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

put_acc_pval <- as.data.frame(summary(putamenAcc)$p.table[3,4])
names(put_acc_pval) <- "putamen_acc_pval"
print(put_acc_pval)
##   putamen_acc_pval
## 1        0.1147318
putamenEff <- gam(putamen ~ s(age) + sex + efficiency, data=r2star, method="REML")
visreg(putamenEff,"efficiency")

put_eff_pval <- as.data.frame(summary(putamenEff)$p.table[3,4])
names(put_eff_pval) <- "putamen_eff_pval"
print(put_eff_pval)
##   putamen_eff_pval
## 1       0.03727838
putamenSpeed <- gam(putamen ~ s(age) + sex + speed, data=r2star, method="REML")
visreg(putamenSpeed,"speed")

put_speed_pval <- as.data.frame(summary(putamenSpeed)$p.table[3,4])
names(put_speed_pval) <- "putamen_speed_pval"
print(put_speed_pval)
##   putamen_speed_pval
## 1         0.08635342

P values

pvalues = list(pall_acc_pval, pall_eff_pval, pall_speed_pval, acc_acc_pval, acc_eff_pval, acc_speed_pval, cau_acc_pval, cau_eff_pval, cau_speed_pval, put_acc_pval, put_eff_pval, put_speed_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
print(pvals)
##   rn pallidum_acc_pval pallidum_eff_pval pallidum_speed_pval
## 1  1         0.2385601         0.1213092           0.1889191
##   accumbens_acc_pval accumbens_eff_pval accumbens_speed_pval
## 1          0.8317156           0.895044            0.6954426
##   caudate_acc_pval caudate_eff_pval caudate_speed_pval putamen_acc_pval
## 1      0.006608071       0.02509784          0.5452894        0.1147318
##   putamen_eff_pval putamen_speed_pval
## 1       0.03727838         0.08635342