/
SensitivityAnalyses_rest_noMeds.R
215 lines (152 loc) · 8.03 KB
/
SensitivityAnalyses_rest_noMeds.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
207
208
209
210
211
212
213
214
215
########################################################
#### SENSITIVITY ANALYSES WITH RESTING STATE SAMPLE ####
########################################################
###############################
### Load data and libraries ###
###############################
subjData <- readRDS("/data/jux/BBL/projects/pncHeterogeneity/subjectData/MergedWithHYDRA/n763_DepAnxTd_rest_hydra_noMeds_subjData.rds")
#Load libraries
library(mgcv)
library(dplyr)
#Total sample size
n <- nrow(subjData)
######################
### Rest ALFF ROIs ###
######################
#Get variable names
restAlff_all <- subjData[c(grep("rest_jlf_alff",names(subjData)))]
restAlff_short <- restAlff_all[,-grep("Cerebell",names(restAlff_all))]
restAlffRegions <- names(restAlff_short)
##Run models (all three groups included).
#TD is the comparison group in Hydra_k2. This model will provide estimates and p-values for TD vs. S1 and TD vs. S2.
restAlffGam <- lapply(restAlffRegions, function(x) {
gam(substitute(i ~ s(age) + sex + restRelMeanRMSMotion + Hydra_k2, list(i = as.name(x))), method="REML", data=subjData)
})
#Look at model summaries
restAlffGamSumm <- lapply(restAlffGam, summary)
##Pass the results to the anova() function to fit ANOVAs (only anova() and not aov() can fit GAMs)
#These are your omnibus ANOVA tests (tells you if the three groups are significantly different)
restAlffAnova <- lapply(restAlffGam, anova)
#Pull the p-values
restAlff_p <- sapply(restAlffAnova, function(v) v$pTerms.pv[3])
#Convert to data frame
restAlff_p <- as.data.frame(restAlff_p)
#Add row names for ease of viewing
rownames(restAlff_p) <- restAlffRegions
#Print original p-values to three decimal places
restAlff_p_round <- round(restAlff_p,3)
##Pull the F-values for the table
restAlff_F <- sapply(restAlffAnova, function(v) v$pTerms.table[3,2])
#Convert to data frame
restAlff_F <- as.data.frame(restAlff_F)
#Print F-values to two decimal places
restAlff_F_round <- round(restAlff_F,2)
#Add region names
restAlffNames <- as.data.frame(restAlffRegions)
restAlff_F_names <- cbind(restAlffNames,restAlff_F_round)
##FDR correction across all omnibus anova tests
#FDR correct p-values
restAlff_p_fdr <- p.adjust(restAlff_p[,1],method="fdr")
#Convert to data frame
restAlff_p_fdr <- as.data.frame(restAlff_p_fdr)
#Print fdr-corrected p-values to three decimal places
restAlff_p_fdr_round <- round(restAlff_p_fdr,3)
#Add region names
restAlff_p_names <- cbind(restAlffNames,restAlff_p_fdr_round)
#Merge F values
restAlff_omnibus <- merge(restAlff_F_names,restAlff_p_names, by="restAlffRegions")
#Keep the original order of the regions
restAlff_omnibus <- restAlff_omnibus[restAlff_p_names$restAlffRegions,]
#Trim table to only regions that passed FDR correction
restAlff_omnibus_signif <- restAlff_omnibus[restAlff_p_fdr<0.05,]
#Save omnibus results as a .csv file
write.table(restAlff_omnibus_signif, file = "/data/jux/BBL/projects/pncHeterogeneity/tablesFigures/pfdr_omnibus_restAlff_noMeds.csv", sep = ",", row.names = FALSE, col.names = TRUE, quote = FALSE)
##Run same model with Hydra_k2_reordered to make S2 the comparison group in order to see S1 vs S2 differences in the model.
#This model will provide estimates and p-values for S1-S2.
restAlffGam_reordered <- lapply(restAlffRegions, function(x) {
gam(substitute(i ~ s(age) + sex + restRelMeanRMSMotion + Hydra_k2_reordered, list(i = as.name(x))), method="REML", data=subjData)
})
#Look at model summaries
restAlffGamSumm_reordered <- lapply(restAlffGam_reordered, summary)
##Pairwise comparisons
#Pull uncorrected p-values
restAlff_S1vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[4,4])
restAlff_S2vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[5,4])
restAlff_S1vsS2 <- sapply(restAlffGam_reordered, function(v) summary(v)$p.table[4,4])
#Combine the pairwise p values
restAlff_pairs <- cbind(restAlff_S1vsTd,restAlff_S2vsTd,restAlff_S1vsS2)
#Convert to data frame
restAlff_pairs <- as.data.frame(restAlff_pairs)
#Add row names for ease of viewing
rownames(restAlff_pairs) <- restAlffRegions
#Print original p-values to three decimal places
restAlff_pairs_round <- round(restAlff_pairs,3)
##Pull the B-values for the table
restAlff_B_S1vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[4,1])
restAlff_B_S2vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[5,1])
restAlff_B_S1vsS2 <- sapply(restAlffGam_reordered, function(v) summary(v)$p.table[4,1])
#Convert to data frame
restAlff_B_S1vsTd <- as.data.frame(restAlff_B_S1vsTd)
restAlff_B_S2vsTd <- as.data.frame(restAlff_B_S2vsTd)
restAlff_B_S1vsS2 <- as.data.frame(restAlff_B_S1vsS2)
#Print B-values to two decimal places
restAlff_B_S1vsTd_round <- round(restAlff_B_S1vsTd,2)
restAlff_B_S2vsTd_round <- round(restAlff_B_S2vsTd,2)
restAlff_B_S1vsS2_round <- round(restAlff_B_S1vsS2,2)
#Add tract names
restAlff_B_names <- cbind(restAlffNames,restAlff_B_S1vsTd_round,restAlff_B_S2vsTd_round,restAlff_B_S1vsS2_round)
##Pull the SE values for the table
restAlff_SE_S1vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[4,2])
restAlff_SE_S2vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[5,2])
restAlff_SE_S1vsS2 <- sapply(restAlffGam_reordered, function(v) summary(v)$p.table[4,2])
#Convert to data frame
restAlff_SE_S1vsTd <- as.data.frame(restAlff_SE_S1vsTd)
restAlff_SE_S2vsTd <- as.data.frame(restAlff_SE_S2vsTd)
restAlff_SE_S1vsS2 <- as.data.frame(restAlff_SE_S1vsS2)
#Print SE values to two decimal places
restAlff_SE_S1vsTd_round <- round(restAlff_SE_S1vsTd,2)
restAlff_SE_S2vsTd_round <- round(restAlff_SE_S2vsTd,2)
restAlff_SE_S1vsS2_round <- round(restAlff_SE_S1vsS2,2)
#Add tract names
restAlff_SE_names <- cbind(restAlffNames,restAlff_SE_S1vsTd_round,restAlff_SE_S2vsTd_round,restAlff_SE_S1vsS2_round)
##Pull the t-values for the table
restAlff_t_S1vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[4,3])
restAlff_t_S2vsTd <- sapply(restAlffGam, function(v) summary(v)$p.table[5,3])
restAlff_t_S1vsS2 <- sapply(restAlffGam_reordered, function(v) summary(v)$p.table[4,3])
#Convert to data frame
restAlff_t_S1vsTd <- as.data.frame(restAlff_t_S1vsTd)
restAlff_t_S2vsTd <- as.data.frame(restAlff_t_S2vsTd)
restAlff_t_S1vsS2 <- as.data.frame(restAlff_t_S1vsS2)
#Print t-values to two decimal places
restAlff_t_S1vsTd_round <- round(restAlff_t_S1vsTd,2)
restAlff_t_S2vsTd_round <- round(restAlff_t_S2vsTd,2)
restAlff_t_S1vsS2_round <- round(restAlff_t_S1vsS2,2)
#Add tract names
restAlff_t_names <- cbind(restAlffNames,restAlff_t_S1vsTd_round,restAlff_t_S2vsTd_round,restAlff_t_S1vsS2_round)
##FDR correction across the three pairwise group comparisons (S1vsTd, S2vsTd, and S1vsS2)
#Create an empty table for the fdr results
restAlff_fdrTable<-as.data.frame(matrix(nrow=112,ncol=3))
colnames(restAlff_fdrTable)[1]<-"restAlff_S1vsTd_pfdr"
colnames(restAlff_fdrTable)[2]<-"restAlff_S2vsTd_pfdr"
colnames(restAlff_fdrTable)[3]<-"restAlff_S1vsS2_pfdr"
#FDR correct across rows
for(i in 1:nrow(restAlff_pairs)) {
row <- restAlff_pairs[i,]
restAlff_fdrTable[i,] <- p.adjust(restAlff_pairs[i,],method="fdr")
}
#Print fdr-corrected p-values to three decimal places
restAlff_fdrTable_round <- round(restAlff_fdrTable,3)
#Add region names
restAlff_p_pairwise_names <- cbind(restAlffNames,restAlff_fdrTable_round)
#Merge B, SE, t, and pfdr values
restAlff_B_SE <- merge(restAlff_B_names,restAlff_SE_names, by="restAlffRegions")
restAlff_B_SE_t <- merge(restAlff_B_SE,restAlff_t_names, by="restAlffRegions")
restAlff_B_SE_t_pfdr <- merge(restAlff_B_SE_t,restAlff_p_pairwise_names, by="restAlffRegions")
#Keep the original order of the regions
restAlff_pairwise <- restAlff_B_SE_t_pfdr[restAlff_p_pairwise_names$restAlffRegions,]
#Reorder the columns
restAlff_pairwise <- restAlff_pairwise[,c(1,2,5,8,11,3,6,9,12,4,7,10,13)]
#Trim table to only regions that passed FDR correction for the omnibus test
restAlff_pairwise_signif <- restAlff_pairwise[restAlff_p_fdr<0.05,]
#Save the pairwise results as a .csv file
write.table(restAlff_pairwise_signif, file = "/data/jux/BBL/projects/pncHeterogeneity/tablesFigures/pfdr_pairwise_restAlff_noMeds.csv", sep = ",", row.names = FALSE, col.names = TRUE, quote = FALSE)