R2Star Sex Models
Load Packages
library(mgcv)
library(voxel)
library('knitr')
library(plyr)
library(ggplot2)
library(visreg)
Read in r2star data
r2star <- read.csv("/Users/imaging/Desktop/r2star/r2star_finalSample_03152017.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
# Pallidum
pallidum <- (r2star$Left.Pallidum + r2star$Right.Pallidum)/2
# Accumbens
accumbens <- (r2star$Left.Accumbens.Area + r2star$Right.Accumbens.Area)/2
# Caudate
caudate <- (r2star$Left.Caudate + r2star$Right.Caudate)/2
# Putamen
putamen <- (r2star$Left.Putamen + r2star$Right.Putamen)/2
Pallidum
# Set variables
pallidum <- as.data.frame(pallidum)
agePall <- as.data.frame(age)
sexPall <- as.data.frame(sex)
pallidumClean <- cbind(pallidum,agePall,sexPall)
pallidumClean <- pallidumClean[complete.cases(pallidumClean), ]
pallidumClean2 <- pallidumClean$pallidum
agePall <- pallidumClean$age
sexPall <- pallidumClean$sex
# Create model
pallGLM <- glm(pallidumClean2 ~ agePall + sexPall)
# Create plots
#plot(predict(pallGLM), pallidumClean2)
#plot(pallGLM$fitted.values, pallGLM$y)
visreg(pallGLM,"agePall","sexPall",overlay=TRUE)
# Print statistics
cor(pallGLM$fitted.values, pallGLM$y)
## [1] 0.5352662
Accumbens
# Set variables
accumbens <- as.data.frame(accumbens)
ageAcc <- as.data.frame(age)
sexAcc <- as.data.frame(sex)
accumbensClean <- cbind(accumbens,ageAcc,sexAcc)
accumbensClean <- accumbensClean[complete.cases(accumbensClean), ]
accumbensClean2 <- accumbensClean$accumbens
ageAcc <- accumbensClean$age
sexAcc <- accumbensClean$sex
# Create model
accumGLM <- glm(accumbensClean2 ~ ageAcc + sexAcc)
# Create plots
#plot(predict(accumGLM), accumbensClean2)
#plot(accumGLM$fitted.values, accumGLM$y)
visreg(accumGLM,"ageAcc","sexAcc",overlay=TRUE)
# Print statistics
cor(accumGLM$fitted.values, accumGLM$y)
## [1] 0.2850843
Caudate
# Set variables
caudate <- as.data.frame(caudate)
ageCau <- as.data.frame(age)
sexCau <- as.data.frame(sex)
caudateClean <- cbind(caudate,ageCau,sexCau)
caudateClean <- caudateClean[complete.cases(caudateClean), ]
caudateClean2 <- caudateClean$caudate
ageCau <- caudateClean$age
sexCau <- caudateClean$sex
# Create model
cauGLM <- glm(caudateClean2 ~ ageCau + sexCau)
# Create plots
#plot(predict(cauGLM), caudateClean2)
#plot(cauGLM$fitted.values, cauGLM$y)
visreg(cauGLM,"ageCau","sexCau",overlay=TRUE)
# Print correlation
cor(cauGLM$fitted.values, cauGLM$y)
## [1] 0.2228193
Putamen
# Set variables
putamen <- as.data.frame(putamen)
agePut <- as.data.frame(age)
sexPut <- as.data.frame(sex)
putamenClean <- cbind(putamen,agePut,sexPut)
putamenClean <- putamenClean[complete.cases(putamenClean), ]
putamenClean2 <- putamenClean$putamen
agePut <- putamenClean$age
sexPut <- putamenClean$sex
# Create model
putGLM <- glm(putamenClean2 ~ agePut + sexPut)
# Create plots
#plot(predict(putGLM), putamenClean2)
#plot(putGLM$fitted.values, putGLM$y)
visreg(putGLM,"agePut","sexPut",overlay=TRUE)
# Print statistics
cor(putGLM$fitted.values, putGLM$y)
## [1] 0.3740499
Summaries
summary(pallGLM)
##
## Call:
## glm(formula = pallidumClean2 ~ agePall + sexPall)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.0057847 -0.0009973 -0.0000066 0.0010050 0.0083191
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.770e-02 2.168e-04 81.634 <2e-16 ***
## agePall 2.995e-04 1.412e-05 21.213 <2e-16 ***
## sexPall.L -9.665e-05 7.223e-05 -1.338 0.181
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.907064e-06)
##
## Null deviance: 0.0045674 on 1123 degrees of freedom
## Residual deviance: 0.0032588 on 1121 degrees of freedom
## AIC: -11134
##
## Number of Fisher Scoring iterations: 2
summary(accumGLM)
##
## Call:
## glm(formula = accumbensClean2 ~ ageAcc + sexAcc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.0092004 -0.0018228 -0.0004506 0.0012495 0.0163502
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.355e-02 3.555e-04 38.117 <2e-16 ***
## ageAcc 2.206e-04 2.318e-05 9.519 <2e-16 ***
## sexAcc.L 3.004e-04 1.186e-04 2.533 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 7.858847e-06)
##
## Null deviance: 0.0096233 on 1127 degrees of freedom
## Residual deviance: 0.0088412 on 1125 degrees of freedom
## AIC: -10052
##
## Number of Fisher Scoring iterations: 2
summary(cauGLM)
##
## Call:
## glm(formula = caudateClean2 ~ ageCau + sexCau)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.0041665 -0.0007916 0.0000443 0.0008475 0.0048045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.398e-02 1.726e-04 80.974 < 2e-16 ***
## ageCau 8.616e-05 1.127e-05 7.645 4.48e-14 ***
## sexCau.L -1.445e-05 5.759e-05 -0.251 0.802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.843469e-06)
##
## Null deviance: 0.0021706 on 1121 degrees of freedom
## Residual deviance: 0.0020628 on 1119 degrees of freedom
## AIC: -11626
##
## Number of Fisher Scoring iterations: 2
summary(putGLM)
##
## Call:
## glm(formula = putamenClean2 ~ agePut + sexPut)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.0045524 -0.0008061 0.0000591 0.0007715 0.0055447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.465e-02 1.600e-04 91.557 <2e-16 ***
## agePut 1.412e-04 1.044e-05 13.516 <2e-16 ***
## sexPut.L -4.959e-05 5.351e-05 -0.927 0.354
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.599247e-06)
##
## Null deviance: 0.0020900 on 1126 degrees of freedom
## Residual deviance: 0.0017976 on 1124 degrees of freedom
## AIC: -11838
##
## Number of Fisher Scoring iterations: 2