R2Star Sex Models

Load Packages

library(mgcv)
library(voxel)
library('knitr')
library(plyr)
library(ggplot2)

Read in r2star data

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

Set sex and age variables

r2star$sex <- ordered(r2star$sex)
r2star$sex <- revalue(r2star$sex, c('1'='male', '2'='female'))
sex <- r2star$sex
age <- (r2star$ageAtScan1)/12

Create and plot models

Right pallidum

rpallidum <- r2star$Right.Pallidum
rpallidumModel <- gam(rpallidum ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(rpallidumModel,"age","sex")

## [1] "Working with a GAM Model"

rpall_pval <- as.data.frame(summary(rpallidumModel)$p.pv)
print(rpall_pval[2,])

## [1] 0.4182165

Left pallidum

lpallidum <- r2star$Left.Pallidum
lpallidumModel <- gam(lpallidum ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(lpallidumModel,"age","sex")

## [1] "Working with a GAM Model"

lpall_pval <- as.data.frame(summary(lpallidumModel)$p.pv)
print(lpall_pval[2,])

## [1] 0.06578102

Pallidum

pallidum <- (r2star$Left.Pallidum + r2star$Right.Pallidum)/2
pallidumModel <- gam(pallidum ~ s(age,k=4) + sex, method="REML", data = r2star)

pallidumGraph <- plotGAM(pallidumModel,"age","sex")

## [1] "Working with a GAM Model"

pallidumGraph <- pallidumGraph + ggtitle("Pallidum R2* vs Age") + labs(x="Age (years)", y="R2* (1/sec)",caption = "(p < 0.001)")

plot(pallidumGraph)

pall_pval <- as.data.frame(summary(pallidumModel)$p.pv)
print(pall_pval[2,])

## [1] 0.1720676

Right accumbens

raccumbens <- r2star$Right.Accumbens.Area
raccumbensModel <- gam(raccumbens ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(raccumbensModel,"age","sex")

## [1] "Working with a GAM Model"

racc_pval <- as.data.frame(summary(raccumbensModel)$p.pv)
print(racc_pval[2,])

## [1] 0.037054

Left accumbens

laccumbens <- r2star$Left.Accumbens.Area
laccumbensModel <- gam(laccumbens ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(laccumbensModel,"age","sex")

## [1] "Working with a GAM Model"

lacc_pval <- as.data.frame(summary(laccumbensModel)$p.pv)
print(lacc_pval[2,])

## [1] 0.02918618

Accumbens

accumbens <- (r2star$Left.Accumbens.Area + r2star$Right.Accumbens.Area)/2
accumbensModel <- gam(accumbens ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(accumbensModel,"age","sex")

## [1] "Working with a GAM Model"

accumbensGraph <- plotGAM(accumbensModel,"age","sex")

## [1] "Working with a GAM Model"

accumbensGraph <- accumbensGraph + ggtitle("Accumbens R2* vs Age") + labs(x="Age (years)", y="R2* (1/sec)",caption = "(p < 0.001)")

plot(accumbensGraph)

acc_pval <- as.data.frame(summary(accumbensModel)$p.pv)
print(acc_pval[2,])

## [1] 0.01310872

Right caudate

rcaudate <- r2star$Right.Caudate
rcaudateModel <- gam(rcaudate ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(rcaudateModel,"age","sex")

## [1] "Working with a GAM Model"

rcau_pval <- as.data.frame(summary(rcaudateModel)$p.pv)
print(rcau_pval[2,])

## [1] 0.8370361

Left caudate

lcaudate <- r2star$Left.Caudate
lcaudateModel <- gam(lcaudate ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(lcaudateModel,"age","sex")

## [1] "Working with a GAM Model"

lcau_pval <- as.data.frame(summary(lcaudateModel)$p.pv)
print(lcau_pval[2,])

## [1] 0.4549892

Caudate

caudate <- (r2star$Left.Caudate + r2star$Right.Caudate)/2
caudateModel <- gam(caudate ~ s(age,k=4) + sex, method="REML", data = r2star)

caudate <- (r2star$Left.Caudate + r2star$Right.Caudate)/2
caudateModel <- gam(caudate ~ s(age,k=4) + sex, method="REML", data = r2star)

caudateGraph <- plotGAM(caudateModel,"age","sex")

## [1] "Working with a GAM Model"

caudateGraph <- caudateGraph + ggtitle("Caudate R2* vs Age") + labs(x="Age (years)", y="R2* (1/sec)",caption = "(p < 0.001)")

plot(caudateGraph) 

cau_pval <- as.data.frame(summary(caudateModel)$p.pv)
print(cau_pval[2,])

## [1] 0.7562452

Right putamen

rputamen <- r2star$Right.Putamen
rputamenModel <- gam(rputamen ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(rputamenModel,"age","sex")

## [1] "Working with a GAM Model"

rput_pval <- as.data.frame(summary(rputamenModel)$p.pv)
print(rput_pval[2,])

## [1] 0.9668501

Left putamen

lputamen <- r2star$Left.Putamen
lputamenModel <- gam(lputamen ~ s(age,k=4) + sex, method="REML", data = r2star)
plotGAM(lputamenModel,"age","sex")

## [1] "Working with a GAM Model"

lput_pval <- as.data.frame(summary(lputamenModel)$p.pv)
print(lput_pval[2,])

## [1] 0.1014401

Putamen

putamen <- (r2star$Left.Putamen + r2star$Right.Putamen)/2
putamenModel <- gam(putamen ~ s(age,k=4) + sex, method="REML", data = r2star)

putamenGraph <- plotGAM(putamenModel,"age","sex")

## [1] "Working with a GAM Model"

putamenGraph <- putamenGraph + ggtitle("Putamen R2* vs Age") + labs(x="Age (years)", y="R2* (1/sec)",caption = "(p < 0.001)")

plot(putamenGraph)

put_pval <- as.data.frame(summary(putamenModel)$p.pv)
print(put_pval[2,])

## [1] 0.362697

P values

pvalues = list(rpall_pval,lpall_pval,pall_pval,lacc_pval,racc_pval,acc_pval,rcau_pval,lcau_pval,cau_pval,lput_pval,rput_pval,put_pval)
pvals <- Reduce(merge, lapply(pvalues, function(x) data.frame(x, rn = row.names(x))))
names(pvals) <- c("region","rpall_pval","lpall_pval","pall_pval","lacc_pval","racc_pval","acc_pval","rcau_pval","lcau_pval","cau_pval","lput_pval","rput_pval","put_pval")
print(pvals[2,])

##   region rpall_pval lpall_pval pall_pval  lacc_pval racc_pval   acc_pval
## 2  sex.L  0.4182165 0.06578102 0.1720676 0.02918618  0.037054 0.01310872
##   rcau_pval lcau_pval  cau_pval lput_pval rput_pval put_pval
## 2 0.8370361 0.4549892 0.7562452 0.1014401 0.9668501 0.362697