-
Notifications
You must be signed in to change notification settings - Fork 0
/
natalie_comparison.Rmd
115 lines (95 loc) · 5.69 KB
/
natalie_comparison.Rmd
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
---
title: "Natalie Randomise comparison"
author: "John Flournoy"
date: "`r Sys.Date()`"
output:
html_document:
toc: true
toc_float: true
code_folding: hide
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
library(RNifti)
library(ggplot2)
colors <- c('#5C2A70',
'#C09A9E',
'#F3EDE8',
'#F3D0D0',
'#44ADF1')
```
This is a comparison of the test statics generated from the original parametric group FEAT model, and the same model re-estimated using non-parametric permutation testing. We can also compare the original *p*-value (in FSL land, it's actually always 1-*p*) with the TFCE corrected *p*-value. Below, find the two sets of plots. _Note that TFCE is a single-sided procedure_.
```{r}
orig_tstat <- dir('/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/dep_threat_pipeline/Group/ThreatReactivity/Natalie/20200430_ThreatReactivity_n151_Tanner_Age_cope2.gfeat/cope1.feat/stats/', pattern = '^zstat[0-9].nii.gz', full.names = TRUE)
randomise_tstat <- dir('/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/dep_threat_pipeline/Group/ThreatReactivity/Natalie/20200430_ThreatReactivity_n151_Tanner_Age_cope2.gfeat/randomise/', pattern = '.*gfeat_tstat[0-9].nii.gz', full.names = TRUE)
randomise_p <- dir('/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/dep_threat_pipeline/Group/ThreatReactivity/Natalie/20200430_ThreatReactivity_n151_Tanner_Age_cope2.gfeat/randomise/', pattern = '.*gfeat_tfce_corrp_tstat[0-9].nii.gz', full.names = TRUE)
mask_file <- dir('/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/dep_threat_pipeline/Group/ThreatReactivity/Natalie/20200430_ThreatReactivity_n151_Tanner_Age_cope2.gfeat/cope1.feat/', pattern = 'mask.nii.gz', full.names = TRUE)
mask <- RNifti::readNifti(mask_file)
```
# Comparison Plots (*Z* scores) {.tabset}
```{r results='asis'}
nvx <- 5e4
j <- sample(1:sum(mask != 0), size = nvx)
if(length(orig_tstat) == length(randomise_tstat)){
for(i in 1:length(orig_tstat)){
cat(paste0('\n\n## ', basename(orig_tstat[[i]]), '\n\n'))
orig <- RNifti::readNifti(orig_tstat[[i]])
randomise <- RNifti::readNifti(randomise_tstat[[i]])
if(!length(randomise) == length(orig)){
stop('Mismatch in dimensions of comparison images')
}
voxel_df <- data.frame(o = as.vector(orig[mask != 0])[j],
r = as.vector(randomise[mask != 0])[j])
p <- ggplot2::ggplot(voxel_df, aes(x = o, y = r)) +
ggplot2::geom_hex(aes(fill = log(..count..), color = log(..count..)), bins = 25, size = .25) +
ggplot2::scale_fill_gradient(breaks = log(c(10, 100, 1000, 10000)),
labels = c(10, 100, 1000, 10000), name = 'Count',
low = colors[[1]], high = colors[[5]]) +
ggplot2::scale_color_gradient(breaks = log(c(10, 100, 1000, 10000)),
labels = c(10, 100, 1000, 10000), name = 'Count',
low = colors[[1]], high = colors[[5]]) +
ggplot2::geom_line(stat = 'smooth', method = 'gam', color = colors[[2]], size = 1) +
ggplot2::labs(x = 'Original parametric Z', y = 'Randomise Z') +
theme(panel.background = element_rect(fill = colors[[3]], size = 0),
panel.grid = element_line(color = colors[[4]], linetype = 'dotted'),
panel.grid.minor = element_blank())
print(p)
}
} else {
stop('Length of original and randomise stats is not the same')
}
```
# Comparison Plots (*p*-values) {.tabset}
The horizontal line is at the threshold p = .025 to show what survives a two-sided test at alpha = .05. Note that since this is just a sample of all voxels in the mask, this may not agree perfectly with what you would see on thresholded maps. This is just to check that higher uncorrected *p*-values tend to go with higher TFCE corrected *p*-values.
```{r, results='asis'}
if(length(orig_tstat) == length(randomise_tstat)){
for(i in 1:length(orig_tstat)){
cat(paste0('\n\n## ', basename(orig_tstat[[i]]), '\n\n'))
orig <- RNifti::readNifti(orig_tstat[[i]])
randomise <- RNifti::readNifti(randomise_p[[i]])
if(!length(randomise) == length(orig)){
stop('Mismatch in dimensions of comparison images')
}
voxel_df <- data.frame(o = pnorm(as.vector(orig[mask != 0])[j]),
r = as.vector(randomise[mask != 0])[j])
p <- ggplot2::ggplot(voxel_df, aes(x = o, y = r)) +
ggplot2::geom_hex(aes(fill = log(..count..), color = log(..count..)), bins = 30, size = .25) +
ggplot2::scale_fill_gradient(breaks = log(c(10, 100, 1000, 10000)),
labels = c(10, 100, 1000, 10000), name = 'Count',
low = colors[[1]], high = colors[[5]]) +
ggplot2::scale_color_gradient(breaks = log(c(10, 100, 1000, 10000)),
labels = c(10, 100, 1000, 10000), name = 'Count',
low = colors[[1]], high = colors[[5]]) +
ggplot2::geom_line(stat = 'smooth', method = 'gam', color = colors[[2]], size = 1) +
ggplot2::labs(x = expression(paste('Original parametric ', P(Z),' (uncorrected)')),
y = expression(paste('Randomise TFCE ', italic(p),'-value (corrected)'))) +
theme(panel.background = element_rect(fill = colors[[3]], size = 0),
panel.grid = element_line(color = colors[[4]], linetype = 'dotted'),
panel.grid.minor = element_blank()) +
ggplot2::geom_hline(yintercept = .975, color = 'black', alpha = .5)
print(p)
}
} else {
stop('Length of original and randomise stats is not the same')
}
```