-
Notifications
You must be signed in to change notification settings - Fork 1
/
fit_ageGams.R
208 lines (170 loc) · 15.8 KB
/
fit_ageGams.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
#### Fit Development GAMs ####
library(dplyr)
library(tidyverse)
library(cifti)
source("/cbica/projects/thalamocortical_development/code/thalamocortical_development/gam_functions/GAM_functions_thalamocortical.R")
############################################################################################################
#### Prepare Data ####
#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -
#S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by="orig_parcelname", sort = F)
#Dataset-specific tract lists for analysis
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R
#Dataset-specific tract diffusion measures
FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv")
MD.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_MD_finalsample_primary.csv")
FA.glasser.pnc$sex <- as.factor(FA.glasser.pnc$sex)
MD.glasser.pnc$sex <- as.factor(MD.glasser.pnc$sex)
FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv")
MD.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_MD_finalsample_primary.csv")
FA.glasser.hcpd$sex <- as.factor(FA.glasser.hcpd$sex)
MD.glasser.hcpd$sex <- as.factor(MD.glasser.hcpd$sex)
############################################################################################################
#### Region-wise GAM Statistics and Derivative-based Temporal Developmental Properties ####
run_gam.fit.smooth <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, covs, k, fixed_edf){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
smooth.characteristics <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.fit.smooth(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, covariates = covs,
knots = k, set_fx = fixed_edf))})
write.csv(smooth.characteristics, sprintf("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_gameffects_%s_%s_%s.csv", dwi.measure, dwi.atlas, dwi.dataset), quote = F, row.names =F)
}
#PNC
run_gam.fit.smooth(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE)
run_gam.fit.smooth(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE)
#HCPD
run_gam.fit.smooth(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE)
run_gam.fit.smooth(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE)
############################################################################################################
#### Region-wise GAM Fitted Value Predictions ####
run_gam.smooth.predict <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, covs, k, fixed_edf, num.pred){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
smooth.fittedvals <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.smooth.predict(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, covariates = covs,
knots = k, set_fx = fixed_edf, increments = num.pred)) %>% mutate(tract = x)})
write.csv(smooth.fittedvals, sprintf("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_fittedvalues_%s_%s_%s.csv", dwi.measure, dwi.atlas, dwi.dataset), quote = F, row.names =F)
}
#PNC
run_gam.smooth.predict(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
run_gam.smooth.predict(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
#HCPD
run_gam.smooth.predict(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
run_gam.smooth.predict(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
############################################################################################################
#### Region-wise GAM Smooth Estimates ####
run_gam.estimate.smooth <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, covs, k, fixed_edf, num.pred){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
smooth.centered <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.estimate.smooth(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, covariates = covs,
knots = k, set_fx = fixed_edf, increments = num.pred)) %>% mutate(tract = x)})
write.csv(smooth.centered, sprintf("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_smoothcentered_%s_%s_%s.csv", dwi.measure, dwi.atlas, dwi.dataset), quote = F, row.names =F)
}
#PNC
run_gam.estimate.smooth(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
run_gam.estimate.smooth(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
#HCPD
run_gam.estimate.smooth(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
run_gam.estimate.smooth(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.pred = 200)
############################################################################################################
#### Region-wise Developmental Derivatives ####
run_gam.derivatives <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, covs, k, fixed_edf, num.draws, num.pred, credible.interval){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
smooth.derivatives <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.derivatives(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, covariates = covs,
knots = k, set_fx = fixed_edf, draws = num.draws,
increments = num.pred, return_posterior_derivatives = credible.interval)) %>% mutate(tract = x)})
write.csv(smooth.derivatives, sprintf("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_derivatives_%s_%s_%s.csv", dwi.measure, dwi.atlas, dwi.dataset), quote = F, row.names =F)
}
#PNC
run_gam.derivatives(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1, num.pred = 200, credible.interval = FALSE)
run_gam.derivatives(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1, num.pred = 200, credible.interval = FALSE)
#HCPD
run_gam.derivatives(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1, num.pred = 200, credible.interval = FALSE)
run_gam.derivatives(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1, num.pred = 200, credible.interval = FALSE)
############################################################################################################
#### Region-wise Age by Sex Interactions ####
FA.glasser.pnc$sex.ordered <- ordered(FA.glasser.pnc$sex, levels = c(1,2))
MD.glasser.pnc$sex.ordered <- ordered(MD.glasser.pnc$sex, levels = c(1,2))
FA.glasser.hcpd$sex.ordered <- ordered(FA.glasser.hcpd$sex, levels = c("F","M"))
MD.glasser.hcpd$sex.ordered <- ordered(MD.glasser.hcpd$sex, levels = c("F","M"))
run_gam.fit.interaction <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, int.var, covs, k, fixed_edf){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
agebysex.interactioneffects <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.factorsmooth.interaction(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, int_var = int.var, covariates = covs,
knots = k, set_fx = fixed_edf))})
write.csv(agebysex.interactioneffects, sprintf("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_sexeffects_%s_%s_%s.csv", dwi.measure, dwi.atlas, dwi.dataset), quote = F, row.names =F)
}
#PNC
run_gam.fit.interaction(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", int.var = "sex.ordered", covs = "sex.ordered + mean_fd", k = 3, fixed_edf = TRUE)
run_gam.fit.interaction(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", int.var = "sex.ordered", covs = "sex.ordered + mean_fd", k = 3, fixed_edf = TRUE)
#HCPD
run_gam.fit.interaction(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", int.var = "sex.ordered", covs = "sex.ordered + mean_fd", k = 3, fixed_edf = TRUE)
run_gam.fit.interaction(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", int.var = "sex.ordered", covs = "sex.ordered + mean_fd", k = 3, fixed_edf = TRUE)
############################################################################################################
#### Region-wise Posterior Smooth Derivatives Correlation with Sensorimotor-Association Axis ####
run_gam.derivatives.SAaxis <- function(dwi.measure, dwi.atlas, dwi.dataset, smooth.var, covs, k, fixed_edf, num.draws, num.pred, credible.interval){
if(dwi.dataset == "hcpd"){
tractlist <- tractlist.hcpd
}
if(dwi.dataset == "pnc"){
tractlist <- tractlist.pnc
}
smooth.derivatives.draws <- map_dfr(tractlist$tract,
function(x){as.data.frame(gam.derivatives(measure = dwi.measure, atlas = dwi.atlas, dataset = dwi.dataset,
region = as.character(x), smooth_var = smooth.var, covariates = covs,
knots = k, set_fx = fixed_edf, draws = num.draws,
increments = num.pred, return_posterior_derivatives = credible.interval)) %>% mutate(tract = x) %>% select(-label)})
# Compute correlations with S-A axis
smooth.derivatives.draws <- left_join(smooth.derivatives.draws, S.A.axis, by = "tract") #assign axis rank to each tract
corr_values <- smooth.derivatives.draws %>%
group_by(draw, age) %>%
do(SAcorrelation = cor(as.numeric(.$SA.axis), as.numeric(.$posterior.derivative), method = c("spearman"))) %>%
unnest(cols = c(SAcorrelation))
corr_values.wide <- corr_values %>% pivot_wider(names_from = "draw", values_from = "SAcorrelation", names_sort = FALSE)
corr_values.wide <- corr_values.wide %>% select(contains("draw"))
return(corr_values.wide)
}
#PNC
ageresolved.FA.pnc <- map(1:10, ~run_gam.derivatives.SAaxis(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1000, num.pred = 100, credible.interval = TRUE)) %>% bind_cols()
write.table(ageresolved.FA.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_temporalSAcorr_FA_glasser_pnc.csv", quote = FALSE, col.names = FALSE, row.names = FALSE, sep = ",")
ageresolved.MD.pnc <- map(1:10, ~run_gam.derivatives.SAaxis(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "pnc", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1000, num.pred = 100, credible.interval = TRUE)) %>% bind_cols()
write.table(ageresolved.MD.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_temporalSAcorr_MD_glasser_pnc.csv", quote = FALSE, col.names = FALSE, row.names = FALSE, sep = ",")
#HCPDs
ageresolved.FA.hcpd <- map(1:10, ~run_gam.derivatives.SAaxis(dwi.measure = "FA", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1000, num.pred = 100, credible.interval = TRUE)) %>% bind_cols()
write.table(ageresolved.FA.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_temporalSAcorr_FA_glasser_hcpd.csv", quote = FALSE, col.names = FALSE, row.names = FALSE, sep = ",")
ageresolved.MD.hcpd <- map(1:10, ~run_gam.derivatives.SAaxis(dwi.measure = "MD", dwi.atlas = "glasser", dwi.dataset = "hcpd", smooth.var = "age", covs = "sex + mean_fd", k = 3, fixed_edf = TRUE, num.draws = 1000, num.pred = 100, credible.interval = TRUE)) %>% bind_cols()
write.table(ageresolved.MD.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_temporalSAcorr_MD_glasser_hcpd.csv", quote = FALSE, col.names = FALSE, row.names = FALSE, sep = ",")