Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add SDM-PSI Estimator #591

Draft
wants to merge 61 commits into
base: main
Choose a base branch
from
Draft
Changes from 1 commit
Commits
Show all changes
61 commits
Select commit Hold shift + click to select a range
52b0ef8
Start working on SDM anisotropic kernel.
tsalo Nov 18, 2020
02c5595
Merge branch 'neurostuff:main' into sdm
tsalo Oct 27, 2021
721ebcf
Run black and improve docstring.
tsalo Oct 29, 2021
0a9b40a
Merge branch 'main' into sdm
tsalo Dec 19, 2021
26cb5e8
Start looking at the Albajes-Eizagirre et al. (2019) supplement.
tsalo Dec 19, 2021
907ac54
Improve documentation a bit.
tsalo Dec 20, 2021
61e8d7d
Use PyMARE instead of metafor's rma.
tsalo Dec 21, 2021
ca6ddc0
Get the permutation functions working.
tsalo Dec 21, 2021
8c530f6
Draft hedges g function.
tsalo Dec 21, 2021
b5be408
Use formulae from wikipedia for hedges' g.
tsalo Dec 21, 2021
f0eee9c
Add empty functions for calculating Hedges g to transforms.
tsalo Dec 21, 2021
b3cf2a7
Try improving Hedges g functions.
tsalo Dec 21, 2021
2067196
Add simulation code (thanks @JulioAPeraza!).
tsalo Dec 22, 2021
57f5d89
Get the simulations working.
tsalo Jan 14, 2022
910a44a
Support different sample sizes per study.
tsalo Jan 14, 2022
0e9207d
I think the single-voxel imputation makes sense now.
tsalo Jan 15, 2022
05260fc
Lots of progress on the simulations.
tsalo Jan 15, 2022
0f0f49a
Update sdm.py
tsalo Jan 15, 2022
3e0a1a2
Document some stuff.
tsalo Jan 16, 2022
1575455
Import stuff.
tsalo Jan 16, 2022
1996c6f
Update setup.cfg
tsalo Jan 16, 2022
e87eea7
Update sdm.py
tsalo Jan 16, 2022
c1c61ed
Some notes.
tsalo Jan 16, 2022
f8b4dd6
Fix the simulated data.
tsalo Jan 16, 2022
61a2769
Documentation improvements.
tsalo Jan 16, 2022
f600bd3
Update sdm.py
tsalo Jan 16, 2022
0d3feb3
Just use the masker to simulate the subject data.
tsalo Jan 18, 2022
fe22b9c
Merge branch 'main' into sdm
tsalo Mar 3, 2022
cb38ec0
Remove supplement R values.
tsalo Mar 3, 2022
679f051
Update sdm.py
tsalo Mar 3, 2022
afa59ca
Run black.
tsalo Mar 3, 2022
4c315fb
Merge branch 'main' into sdm
tsalo Mar 8, 2022
e74be9c
Update sdm.py
tsalo Mar 8, 2022
1855caf
Test out the simulation method.
tsalo Mar 30, 2022
c6d9371
Keep working on the simulation code.
tsalo Mar 30, 2022
e575b0e
Add some tests.
tsalo Mar 30, 2022
34547dc
Update sdm.py
tsalo Mar 30, 2022
a24104f
docstrings.
tsalo Mar 30, 2022
1e0c9cc
Try writing _calculate_hedges_maps.
tsalo May 11, 2022
6df11b3
Run black.
tsalo May 11, 2022
59777cf
Merge remote-tracking branch 'upstream/main' into sdm
tsalo May 11, 2022
c350dd9
Update.
tsalo May 11, 2022
58d9bb4
Update.
tsalo May 11, 2022
63e878e
Move some functions.
tsalo May 11, 2022
d4c2c65
More.
tsalo May 11, 2022
103f5b9
Update sdm.py
tsalo May 13, 2022
267ceaf
Update setup.cfg
tsalo May 14, 2022
7e7a6ac
Update sdm.py
tsalo May 16, 2022
eb52eb9
Update sdm.py
tsalo May 16, 2022
3410c6b
Update test_sdm.py
tsalo May 17, 2022
177c3cc
Update sdm.py
tsalo May 17, 2022
6119e7c
Update sdm.py
tsalo May 17, 2022
ff933c7
Update test_sdm.py
tsalo May 17, 2022
6f7993c
Minor stuff from last night.
tsalo May 18, 2022
54bc79c
More stuff.
tsalo May 18, 2022
ef2e5ab
Update.
tsalo May 19, 2022
ec769ba
Just use PyMARE for the heterogeneity stats.
tsalo May 23, 2022
57d0ebb
Update sdm.py
tsalo May 24, 2022
8bcb8e3
Update sdm.py
tsalo May 24, 2022
2861950
Merge remote-tracking branch 'upstream/main' into sdm
tsalo Aug 12, 2022
ee28f73
[skip ci] Update sdm.py.
tsalo Aug 12, 2022
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Just use PyMARE for the heterogeneity stats.
  • Loading branch information
tsalo committed May 23, 2022
commit ec769bae387412b44823fc6a203d0d654a77642b
61 changes: 7 additions & 54 deletions nimare/meta/cbma/sdm.py
Original file line number Diff line number Diff line change
@@ -544,7 +544,7 @@ def run_variance_meta(study_hedges_imgs, study_hedges_var_imgs, design):

.. todo::

Implement heterogeneity statistics and linear contrast calculation.
Implement linear contrast calculation.

Parameters
----------
@@ -589,63 +589,16 @@ def run_variance_meta(study_hedges_imgs, study_hedges_var_imgs, design):
pymare_dset = pymare.Dataset(y=study_hedges_imgs, v=study_hedges_var_imgs, X=design)
est.fit_dataset(pymare_dset)
est_summary = est.summary()
meta_effect_size_img = est_summary.est.squeeze()
fe_stats = est_summary.get_fe_stats()
hg_stats = est_summary.get_heterogeneity_stats()
meta_tau2_img = est_summary.tau2.squeeze()
meta_effect_size_img = fe_stats["est"].squeeze()
cochrans_q_img = hg_stats["Q"].squeeze()
i2_img = hg_stats["I^2"].squeeze()

# NOTE: Should the linear hypothesis contrast per conducted here?

# Just a guess at degrees of freedom for now
df = design.shape[0] - design.shape[1]

for i_voxel in range(n_voxels):
# voxel_effect_size = meta_effect_size_img[i_voxel]
voxel_tau2 = meta_tau2_img[i_voxel]
voxel_hedges_values = study_hedges_imgs[:, i_voxel]
voxel_hedges_var_values = study_hedges_var_imgs[:, i_voxel]

# NOTE: Figure out W and W_FE
# "They then commonly combine the effect sizes using a weighted linear model, in which each
# element of the diagonal of the weighting matrix W is the inverse of the variance of the
# corresponding study plus the heterogeneity."
W = np.diag((1 / (voxel_hedges_var_values**2)) + voxel_tau2)

# "Some meta-analysts assume that tau2 is zero, conducting a so-called "fixed-effects"
# meta-analysis, but this assumption is generally discouraged."
# I'm guessing this means that W_FE comes from a WLS estimate with tau2 == 0,
# which might mean that diag(W_FE) = diag(W) - tau2.
W_FE = np.diag(1 / voxel_hedges_var_values)

# P_FE = W_FE - W_FE * X * (X.T * W_FE * X)^-1 * X.T * W_FE
# P_FE = W_FE - (W_FE * X * np.inv(X.T * W_FE * X) * X.T * W_FE)
P_FE = W_FE - W_FE.dot(design).dot(np.inv(design.T.dot(W_FE).dot(design))).dot(
design.T
).dot(W_FE)

# "trace" is the sum of the main diagonal elements
trace_P_FE = np.sum(np.diag(P_FE))

# Calculate Cochran's Q
# Q = Y.T * P_FE * Y, Q >= 0, per (2) in :footcite:t:`albajes2019meta2`
cochrans_q_img = np.dot(np.dot(voxel_hedges_values.T, P_FE), voxel_hedges_values)
cochrans_q_img[cochrans_q_img < 0] = 0

# Calculate H^2 (ratio of total variability to sampling variability)
# H2 = 1 + (tau2 / df) * trace(P_FE), per (2) in :footcite:t:`albajes2019meta2`
h2_img = 1 + (meta_tau2_img / df) * trace_P_FE

# Calculate I^2 (percentage of total variability due to heterogeneity)
# I2 = 1 - (1 / H2), I2 >= 0, per (2) in :footcite:t:`albajes2019meta2`
i2_img = 1 - (1 / h2_img)
i2_img[i2_img < 0] = 0

lambda_ = np.inv(design.T.dot(W).dot(design))
beta = lambda_.dot(design.T).dot(W).dot(voxel_hedges_values)

# I don't know what to do with beta. Is it going to be equivalent to the DerSimonian-Laird
# effect size?
meta_effect_size_img[i_voxel] = beta

return meta_effect_size_img, meta_tau2_img, cochrans_q_img, h2_img, i2_img
return meta_effect_size_img, meta_tau2_img, cochrans_q_img, i2_img


def combine_imputation_results(coefficient_maps, variance_maps, cochrans_q_maps, i2_maps):