|
21 | 21 | MersenneTwister(r), |
22 | 22 | design, |
23 | 23 | signal, |
24 | | - UniformOnset(; offset=5, width=4), |
| 24 | + NoOnset(), |
25 | 25 | RedNoise(noiselevel=1); |
26 | 26 | return_epoched=true, |
27 | 27 | ) |
28 | 28 | data = reshape(data, size(data, 1), :) |
29 | | - data = data[:, evts.condition.=="small"] .- data[:, evts.condition.=="large"] |
| 29 | + data_diff = data[:, evts.condition.=="small"] .- data[:, evts.condition.=="large"] |
30 | 30 |
|
31 | 31 | return data, |
32 | | - clusterdepth(data'; τ=quantile(TDist(n_subjects - 1), 0.975), nperm=1000, show_warnings=false) |
33 | | -end;</code></pre><h2 id="Understanding-the-simulation"><a class="docs-heading-anchor" href="#Understanding-the-simulation">Understanding the simulation</a><a id="Understanding-the-simulation-1"></a><a class="docs-heading-anchor-permalink" href="#Understanding-the-simulation" title="Permalink"></a></h2><p>let's have a look at the actual data by running it once, plotting condition wise trials, the ERP and histograms of uncorrected and corrected p-values</p><pre><code class="language-julia hljs">data, pval = run_fun(5) |
| 32 | + clusterdepth(data_diff; τ=quantile(TDist(n_subjects - 1), 0.975), nperm=4000, show_warnings=false) |
| 33 | +end;</code></pre><h2 id="Understanding-the-simulation"><a class="docs-heading-anchor" href="#Understanding-the-simulation">Understanding the simulation</a><a id="Understanding-the-simulation-1"></a><a class="docs-heading-anchor-permalink" href="#Understanding-the-simulation" title="Permalink"></a></h2><p>let's have a look at the actual data by running it once, plotting condition wise trials, the ERP and histograms of uncorrected and corrected p-values</p><pre><code class="language-julia hljs">data, pval = run_fun(1) |
34 | 34 | conditionSmall = data[:, 1:2:end] |
35 | 35 | conditionLarge = data[:, 2:2:end] |
36 | 36 | pval_uncorrected = |
|
39 | 39 | TDist(n_subjects - 1), |
40 | 40 | abs.(ClusterDepth.studentt(conditionSmall .- conditionLarge)), |
41 | 41 | ) |
42 | | -sig = pval_uncorrected .<= 0.025;</code></pre><p>For the uncorrected p-values based on the t-distribution, we get a type1 error over "time":</p><pre><code class="language-julia hljs">mean(sig)</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">0.168141592920354</code></pre><div class="admonition is-info" id="Note-f01acb8e6031b57e"><header class="admonition-header">Note<a class="admonition-anchor" href="#Note-f01acb8e6031b57e" title="Permalink"></a></header><div class="admonition-body"><p>Type-I error is not the FWER (family wise error rate). FWER is the property of a set of tests (in this case tests per time-point), we can calculate it by repeating such tests, and checking for each repetition whether any sample of a repetition is significant (e.g. <code>any(sig)</code> followed by a <code>mean(repetitions_anysig)</code>).</p></div></div><pre><code class="language-julia hljs">f = Figure(); |
| 42 | +sig = pval_uncorrected .<= 0.025;</code></pre><p>For the uncorrected p-values based on the t-distribution, we get a type1 error over "time":</p><pre><code class="language-julia hljs">mean(sig)</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">0.017699115044247787</code></pre><div class="admonition is-info" id="Note-f01acb8e6031b57e"><header class="admonition-header">Note<a class="admonition-anchor" href="#Note-f01acb8e6031b57e" title="Permalink"></a></header><div class="admonition-body"><p>Type-I error is not the FWER (family wise error rate). FWER is the property of a set of tests (in this case tests per time-point), we can calculate it by repeating such tests, and checking for each repetition whether any sample of a repetition is significant (e.g. <code>any(sig)</code> followed by a <code>mean(repetitions_anysig)</code>).</p></div></div><pre><code class="language-julia hljs">f = Figure(); |
43 | 43 | series!(Axis(f[1, 1], title="condition==small"), conditionSmall', solid_color=:red) |
44 | 44 | series!(Axis(f[1, 2], title="condition==large"), conditionLarge', solid_color=:blue) |
45 | 45 | ax = Axis(f[2, 1:2], title="ERP (mean over trials)") |
|
53 | 53 |
|
54 | 54 | hist!(Axis(f[3, 1], title="uncorrected pvalues"), pval_uncorrected, bins=0:0.01:1.1) |
55 | 55 | hist!(Axis(f[3, 2], title="clusterdepth corrected pvalues"), pval, bins=0:0.01:1.1) |
56 | | -f</code></pre><img src="8c918a6e.png" alt="Example block output"/><h2 id="Run-simulations"><a class="docs-heading-anchor" href="#Run-simulations">Run simulations</a><a id="Run-simulations-1"></a><a class="docs-heading-anchor-permalink" href="#Run-simulations" title="Permalink"></a></h2><p>This takes some seconds (depending on your infrastructure)</p><pre><code class="language-julia hljs">reps = 500 |
| 56 | +f</code></pre><img src="92f9659c.png" alt="Example block output"/><h2 id="Run-simulations"><a class="docs-heading-anchor" href="#Run-simulations">Run simulations</a><a id="Run-simulations-1"></a><a class="docs-heading-anchor-permalink" href="#Run-simulations" title="Permalink"></a></h2><p>This takes some seconds (depending on your infrastructure)</p><pre><code class="language-julia hljs">reps = 500 |
57 | 57 | res = fill(NaN, reps, 2) |
58 | 58 | Threads.@threads for r = 1:reps |
59 | 59 | data, pvals = run_fun(r) |
60 | | - res[r, 1] = mean(pvals .<= 0.05 / 2) |
| 60 | + res[r, 1] = mean(pvals .<= 0.05) |
61 | 61 | res[r, 2] = |
62 | 62 | mean(abs.(ClusterDepth.studentt(data)) .>= quantile(TDist(n_subjects - 1), 0.975)) |
63 | | -end;</code></pre><p>Finally, let's calculate the percentage of simulations where we find a significant effect somewhere</p><pre><code class="language-julia hljs">mean(res .> 0, dims=1) |> x -> (:clusterdepth => x[1], :uncorrected => x[2])</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">(:clusterdepth => 0.1, :uncorrected => 1.0)</code></pre><p>Nice, correction seems to work in principle :) Clusterdepth is not necessarily exactly 5%, but with more repetitions we should get there (e.g. with 5000 repetitions, we got 0.051%).</p><div class="admonition is-info" id="Info-f4cd68306325666c"><header class="admonition-header">Info<a class="admonition-anchor" href="#Info-f4cd68306325666c" title="Permalink"></a></header><div class="admonition-body"><p>if you look closely, the <code>:uncorrected</code> value (can be around 60%) is not as bad as the 99% promised in the introduction. This is due to the correlation between the tests introduced by the noise. Indeed, a good exercise is to repeat everything, but put <code>RedNoise</code> to <code>WhiteNoise</code></p></div></div><hr/><p><em>This page was generated using <a href="https://github.com/fredrikekre/Literate.jl">Literate.jl</a>.</em></p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../../howto/R_betweengroups/">« R: Between groups</a><a class="docs-footer-nextpage" href="../type1_troendle/">Troendle FWER »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option><option value="catppuccin-latte">catppuccin-latte</option><option value="catppuccin-frappe">catppuccin-frappe</option><option value="catppuccin-macchiato">catppuccin-macchiato</option><option value="catppuccin-mocha">catppuccin-mocha</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.14.1 on <span class="colophon-date" title="Wednesday 1 October 2025 08:15">Wednesday 1 October 2025</span>. Using Julia version 1.11.7.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html> |
| 63 | +end;</code></pre><p>Finally, let's calculate the percentage of simulations where we find a significant effect somewhere</p><pre><code class="language-julia hljs">mean(res .> 0, dims=1) |> x -> (:clusterdepth => x[1], :uncorrected => x[2])</code></pre><pre class="documenter-example-output"><code class="nohighlight hljs ansi">(:clusterdepth => 0.044, :uncorrected => 1.0)</code></pre><p>Nice, correction seems to work in principle :) Clusterdepth is not necessarily exactly 5%, but with more repetitions we should get there (e.g. with 5000 repetitions, we got 0.051%).</p><div class="admonition is-info" id="Info-f4cd68306325666c"><header class="admonition-header">Info<a class="admonition-anchor" href="#Info-f4cd68306325666c" title="Permalink"></a></header><div class="admonition-body"><p>if you look closely, the <code>:uncorrected</code> value (can be around 60%) is not as bad as the 99% promised in the introduction. This is due to the correlation between the tests introduced by the noise. Indeed, a good exercise is to repeat everything, but put <code>RedNoise</code> to <code>WhiteNoise</code></p></div></div><hr/><p><em>This page was generated using <a href="https://github.com/fredrikekre/Literate.jl">Literate.jl</a>.</em></p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../../../howto/R_betweengroups/">« R: Between groups</a><a class="docs-footer-nextpage" href="../type1_troendle/">Troendle FWER »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="auto">Automatic (OS)</option><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option><option value="catppuccin-latte">catppuccin-latte</option><option value="catppuccin-frappe">catppuccin-frappe</option><option value="catppuccin-macchiato">catppuccin-macchiato</option><option value="catppuccin-mocha">catppuccin-mocha</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> version 1.14.1 on <span class="colophon-date" title="Friday 3 October 2025 16:08">Friday 3 October 2025</span>. Using Julia version 1.11.7.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html> |
0 commit comments