Load Packages

library(mgcv)
library(voxel)
library('knitr')
require(graphics)

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

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

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

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

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

Create models

# pallidum
rpallidumModel <- gam(rpallidum ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
lpallidumModel <- gam(lpallidum ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
pallidumModel <- gam(pallidum ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)

pallidumSubset <- as.data.frame(subset(r2star, 144 < ageAtScan1 & ageAtScan1 < 196)) #this is ages 12-16
pallSubset <- (pallidumSubset$Left.Pallidum + pallidumSubset$Right.Pallidum)/2
pallSex <- pallidumSubset$sex
pallAge <- (pallidumSubset$ageAtScan1)/12
pallidumModel2<- gam(pallSubset ~ s(pallAge,k=4) + pallSex + s(pallAge,by=pallSex,k=4),method="REML", data=pallidumSubset)
summary(pallidumModel2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## pallSubset ~ s(pallAge, k = 4) + pallSex + s(pallAge, by = pallSex, 
##     k = 4)
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.213e-02  8.206e-05  269.70   <2e-16 ***
## pallSex.L   -2.216e-04  1.161e-04   -1.91   0.0569 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df      F  p-value    
## s(pallAge)          1.001  1.002 16.733 5.13e-05 ***
## s(pallAge):pallSex2 1.000  1.000  0.083    0.773    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0696   Deviance explained = 7.67%
## -REML = -1952.3  Scale est. = 2.6156e-06  n = 396
plotGAM(pallidumModel2,"pallAge","pallSex")
## [1] "Working with a GAM Model"

# accumbens
raccumbensModel <- gam(raccumbens ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
laccumbensModel <- gam(laccumbens ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
accumbensModel <- gam(accumbens ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)

accumbensSubset <- as.data.frame(subset(r2star, 192 < ageAtScan1 & ageAtScan1 < 240)) #this is ages 16-20
accSubset <- (accumbensSubset$Left.Accumbens.Area + accumbensSubset$Right.Accumbens.Area)/2
accSex <- accumbensSubset$sex
accAge <- (accumbensSubset$ageAtScan1)/12
accumbensModel2<- gam(accSubset ~ s(accAge,k=4) + accSex + s(accAge,by=accSex,k=4),method="REML", data=accumbensSubset)
summary(accumbensModel2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## accSubset ~ s(accAge, k = 4) + accSex + s(accAge, by = accSex, 
##     k = 4)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0174822  0.0001450 120.558  < 2e-16 ***
## accSex.L    0.0005490  0.0002051   2.677  0.00774 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df     F p-value
## s(accAge)         1.002  1.004 0.518   0.472
## s(accAge):accSex2 1.001  1.002 0.039   0.845
## 
## R-sq.(adj) =  0.014   Deviance explained = 2.16%
## -REML = -1708.8  Scale est. = 8.0865e-06  n = 391
plotGAM(accumbensModel2,"accAge","accSex")
## [1] "Working with a GAM Model"

# caudate
rcaudateModel <- gam(rcaudate ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
lcaudateModel <- gam(lcaudate ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
caudateModel <- gam(caudate ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)

# putamen
rputamenModel <- gam(rputamen ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
lputamenModel <- gam(lputamen ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)
putamenModel <- gam(putamen ~ s(age,k=4) + sex + s(age,by=sex,k=4), method="REML", data = r2star)

Calculate p values and fdr corrections

# get p values
# pallidum
rpall_pval <- as.data.frame(summary.gam(rpallidumModel)$s.table[2,4])
lpall_pval <- as.data.frame(summary.gam(lpallidumModel)$s.table[2,4])
pall_pval <- as.data.frame(summary.gam(pallidumModel)$s.table[2,4])

# accumbens
racc_pval <- as.data.frame(summary.gam(raccumbensModel)$s.table[2,4])
lacc_pval <- as.data.frame(summary.gam(laccumbensModel)$s.table[2,4])
acc_pval <- as.data.frame(summary.gam(accumbensModel)$s.table[2,4])

# caudate
rcau_pval <- as.data.frame(summary.gam(rcaudateModel)$s.table[2,4])
lcau_pval <- as.data.frame(summary.gam(lcaudateModel)$s.table[2,4])
cau_pval <- as.data.frame(summary.gam(caudateModel)$s.table[2,4])

# putamen
rput_pval <- as.data.frame(summary.gam(rputamenModel)$s.table[2,4])
lput_pval <- as.data.frame(summary.gam(lputamenModel)$s.table[2,4])
put_pval <- as.data.frame(summary.gam(putamenModel)$s.table[2,4])
# rename p value data frames
## pallidum
names(rpall_pval) <- "r_pallidum_pval"
names(lpall_pval) <- "l_pallidum_pval"
names(pall_pval) <- "pallidum_pval"

## accumbens
names(lacc_pval) <- "l_accumbens_pval"
names(racc_pval) <- "r_accumbens_pval"
names(acc_pval) <- "accumbens_pval"

## caudate
names(rcau_pval) <- "r_caudate_pval"
names(lcau_pval) <- "l_caudate_pval"
names(cau_pval) <- "caudate_pval"

## putamen
names(lput_pval) <- "l_putamen_pval"
names(rput_pval) <- "r_putamen_pval"
names(put_pval) <- "putamen_pval"
# consolidate 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))))
pvalst <- t(pvals)
# fdr corrections for p values
fdrAgeSex <- as.data.frame(p.adjust(pvals[1,], "fdr"))
names(fdrAgeSex) <- "fdrAgeSex"
fdr <- as.data.frame(cbind(pvalst, fdrAgeSex))
pvalsFdr <- fdr[-c(1), ]

Plot models, p values, and fdr corrections

Right pallidum

plotGAM(rpallidumModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[1,])
##                    pvalst fdrAgeSex
## r_pallidum_pval 0.6609822         1

Left pallidum

plotGAM(lpallidumModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[2,])
##                    pvalst fdrAgeSex
## l_pallidum_pval 0.9836095         1

Pallidum

plotGAM(pallidumModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[3,])
##                  pvalst fdrAgeSex
## pallidum_pval 0.8899524         1

Right accumbens

plotGAM(raccumbensModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[4,])
##                     pvalst fdrAgeSex
## l_accumbens_pval 0.2359971  0.724844

Left accumbens

plotGAM(laccumbensModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[5,])
##                      pvalst fdrAgeSex
## r_accumbens_pval 0.05843617 0.3411771

Accumbens

plotGAM(accumbensModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[6,])
##                    pvalst fdrAgeSex
## accumbens_pval 0.07873317 0.3411771

Right caudate

plotGAM(rcaudateModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[7,])
##                   pvalst fdrAgeSex
## r_caudate_pval 0.6283357         1

Left caudate

plotGAM(lcaudateModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[8,])
##                   pvalst fdrAgeSex
## l_caudate_pval 0.9808758         1

Caudate

plotGAM(caudateModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[9,])
##                 pvalst fdrAgeSex
## caudate_pval 0.8911627         1

Right putamen

plotGAM(rputamenModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[10,])
##                   pvalst fdrAgeSex
## l_putamen_pval 0.0665912 0.3411771

Left putamen

plotGAM(lputamenModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[11,])
##                   pvalst fdrAgeSex
## r_putamen_pval 0.4080206 0.8840446

Putamen

plotGAM(putamenModel,"age","sex")
## [1] "Working with a GAM Model"

print(pvalsFdr[12,])
##                 pvalst fdrAgeSex
## putamen_pval 0.2787861  0.724844

P values and FDR corrections

print(pvalsFdr)
##                      pvalst fdrAgeSex
## r_pallidum_pval   0.6609822 1.0000000
## l_pallidum_pval   0.9836095 1.0000000
## pallidum_pval     0.8899524 1.0000000
## l_accumbens_pval  0.2359971 0.7248440
## r_accumbens_pval 0.05843617 0.3411771
## accumbens_pval   0.07873317 0.3411771
## r_caudate_pval    0.6283357 1.0000000
## l_caudate_pval    0.9808758 1.0000000
## caudate_pval      0.8911627 1.0000000
## l_putamen_pval    0.0665912 0.3411771
## r_putamen_pval    0.4080206 0.8840446
## putamen_pval      0.2787861 0.7248440