R2Star K Models

Load Packages

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

Read in data

r2star <- read.csv("~/Desktop/r2star/r2star_finalSample_03152017.csv")
kvals <- read.csv("~/Desktop/r2star/n453_Kvalues_08162016.csv")
r2star <- merge(r2star,kvals,by=c("bblid"))

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 k variables

k <- r2star$k
logk <- r2star$logk

Striatal Models

Pallidum

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

pall_k_pval <- as.data.frame(summary(pallidum_k)$p.table[3,4])
names(pall_k_pval) <- "pallidum_k_pval"
print(pall_k_pval)

##   pallidum_k_pval
## 1       0.9011253

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

pall_logk_pval <- as.data.frame(summary(pallidum_logk)$p.table[3,4])
names(pall_logk_pval) <- "pallidum_logk_pval"
print(pall_logk_pval)

##   pallidum_logk_pval
## 1          0.7421413

Accumbens

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

acc_k_pval <- as.data.frame(summary(accumbens_k)$p.table[3,4])
names(acc_k_pval) <- "accumbens_k_pval"
print(acc_k_pval)

##   accumbens_k_pval
## 1        0.1721679

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

acc_logk_pval <- as.data.frame(summary(accumbens_logk)$p.table[3,4])
names(acc_logk_pval) <- "accumbens_logk_pval"
print(acc_logk_pval)

##   accumbens_logk_pval
## 1           0.5949585

Caudate

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

cau_k_pval <- as.data.frame(summary(caudate_k)$p.table[3,4])
names(cau_k_pval) <- "caudate_k_pval"
print(cau_k_pval)

##   caudate_k_pval
## 1      0.9755573

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

cau_logk_pval <- as.data.frame(summary(caudate_logk)$p.table[3,4])
names(cau_logk_pval) <- "caudate_logk_pval"
print(cau_logk_pval)

##   caudate_logk_pval
## 1         0.2682225

Putamen

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

put_k_pval <- as.data.frame(summary(putamen_k)$p.table[3,4])
names(put_k_pval) <- "putamen_k_pval"
print(put_k_pval)

##   putamen_k_pval
## 1      0.9286563

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

put_logk_pval <- as.data.frame(summary(putamen_logk)$p.table[3,4])
names(put_logk_pval) <- "putamen_logk_pval"
print(put_logk_pval)

##   putamen_logk_pval
## 1          0.519523

P values

pvalues = list(pall_k_pval, pall_logk_pval, acc_k_pval, acc_logk_pval, cau_k_pval, cau_logk_pval, put_k_pval, put_logk_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
print(pvals)

##   rn pallidum_k_pval pallidum_logk_pval accumbens_k_pval
## 1  1       0.9011253          0.7421413        0.1721679
##   accumbens_logk_pval caudate_k_pval caudate_logk_pval putamen_k_pval
## 1           0.5949585      0.9755573         0.2682225      0.9286563
##   putamen_logk_pval
## 1          0.519523