test

library(mgcv)

## Loading required package: nlme

## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.

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

r2star <- read.csv("~/Desktop/r2star/r2star_finalSample_03152017.csv")
factors <- read.csv("~/Desktop/r2star/n1601_goassess_itemwise_psychosis_JCPP_factor_scores_20161219.csv")
r2star <- merge(r2star,factors,by=c("bblid","scanid"))

r2star$Right.Pallidum <- (r2star$Right.Pallidum)*1000
r2star$Left.Pallidum <- (r2star$Left.Pallidum)*1000
r2star$Right.Accumbens.Area <- (r2star$Right.Accumbens.Area)*1000
r2star$Left.Accumbens.Area <- (r2star$Left.Accumbens.Area)*1000
r2star$Right.Caudate <- (r2star$Right.Caudate)*1000
r2star$Left.Caudate <- (r2star$Left.Caudate)*1000
r2star$Right.Putamen <- (r2star$Right.Putamen)*1000
r2star$Left.Putamen <- (r2star$Left.Putamen)*1000

r2star$sex <- ordered(r2star$sex)
sex <- r2star$sex
age <- (r2star$ageAtScan1)/12

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

overallPsychopath <- r2star$overall_psychopathology_4factor
overallPsychosis <- r2star$Overall_Psychosis
pos2fac <- r2star$F1_Positive_2Fac
pos3fac <- r2star$F1_Positive_3Fac
posDelusional <- r2star$F2_Positive_Delusional_3Fac
neg2fac <- r2star$F2_Negative_2Fac
neg3fac <- r2star$F3_Negative_3Fac

rputamen_overallPsychopath <- gam(rputamen ~ s(age) + sex + overallPsychopath, data=r2star, method="REML")

putamen_overallPsychopath <- gam(putamen ~ s(age) + sex + overallPsychopath, data=r2star, method="REML")
put_overallPsychopath_pval <- as.data.frame(summary(putamen_overallPsychopath)$p.table[3,4])
names(put_overallPsychopath_pval) <- "putamen_overallPsychopath_pval"

putamen_neg2fac <- gam(putamen ~ s(age) + sex + neg2fac, data=r2star, method="REML")

put_neg2fac_pval <- as.data.frame(summary(putamen_neg2fac)$p.table[3,4])
names(put_neg2fac_pval) <- "putamen_neg2fac_pval"

putamen_neg3fac <- gam(putamen ~ s(age) + sex + neg3fac, data=r2star, method="REML")

put_neg3fac_pval <- as.data.frame(summary(putamen_neg3fac)$p.table[3,4])
names(put_neg3fac_pval) <- "putamen_neg3fac_pval"

# whole putamen p values
wpvalues = list(put_overallPsychopath_pval, put_neg2fac_pval, put_neg3fac_pval)
wpvals <- as.data.frame(wpvalues)
wpvalst <- t(wpvals)
# whole fdr correctiosn
wfdr <- as.data.frame(p.adjust(wpvals[1,], "fdr"))
names(wfdr) <- "wfdr"
wpvals_fdr <- as.data.frame(cbind(wpvalst, wfdr))
print(wpvals_fdr)

##                                    wpvalst       wfdr
## putamen_overallPsychopath_pval 0.005457402 0.01430206
## putamen_neg2fac_pval           0.012882150 0.01430206
## putamen_neg3fac_pval           0.014302059 0.01430206

visreg(putamen_neg3fac,"neg3fac",main="R2* vs Negative 3 Factor",xlab="Negative: 3 Factor",ylab="R2* (1/sec)",partial=FALSE)
legend("bottomright", bty="n", legend=paste("R2 =", format(summary(putamen_neg3fac)$r.sq, digits=4)))
legend("topright", bty="n", legend=paste("p =", format(wpvals_fdr[3,1], digits=4)))

visreg(putamen_neg2fac,"neg2fac",main="R2* vs Negative 2 Factor",xlab="Negative: 2 Factor",ylab="R2* (1/sec)",partial=FALSE)
legend("bottomright", bty="n", legend=paste("R2 =", format(summary(putamen_neg2fac)$r.sq, digits=4)))
legend("topright", bty="n", legend=paste("p =", format(wpvals_fdr[2,1], digits=4)))

visreg(putamen_overallPsychopath,"overallPsychopath",main="R2* vs Overall Psychopathology",xlab="Overall Psychopathology",ylab="R2* (1/sec)",partial=FALSE)
legend("bottomright", bty="n", legend=paste("R2 =", format(summary(putamen_overallPsychopath)$r.sq, digits=4)))
legend("topright", bty="n", legend=paste("p =", format(wpvals_fdr[1,1], digits=4)))