<div dir="ltr"><div><div>Hello, I cannot help you with the inference on subject level, however we do have nice example for the group level inference. It is what we used in pandora data paper <a href="http://f1000research.com/articles/4-174/v1">http://f1000research.com/articles/4-174/v1</a> (review pending) to produce figure 4 and table3. Code for the replication of the analysis is available at <a href="https://github.com/psychoinformatics-de/paper-f1000_pandora_data">https://github.com/psychoinformatics-de/paper-f1000_pandora_data</a><br><br></div>Best wishes,<br></div>Richard<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Aug 11, 2015 at 11:13 PM, Yaroslav Halchenko <span dir="ltr"><<a href="mailto:debian@onerussian.com" target="_blank">debian@onerussian.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
On Tue, 11 Aug 2015, Roni Maimon wrote:<br>
<br>
> Hi all,<br>
<br>
Hi Roni,<br>
<br>
> I started by implementing the inference at the subject level (attaching the<br>
> code). Is this how I'm supposed to evaluate the p values of the<br>
> classifications for a single subject? What is the differences between<br>
> adding the null_dist to the sl level and the cross validation level?<br>
<br>
if you add it at the CV level (which is also legit) you would need some<br>
way to "collect" all of the stats across all the searchlights, e.g.<br>
replace output of CV accuracy/error with the p value, or make two<br>
samples (if you care to see also effects -- i.e. CV<br>
error/accuracy), one of which would be accuracy, another p-value. Also<br>
it might be "noisier" since different searchlights then might permute<br>
differently.  Also I think it will take a bit longer<br>
<br>
So the cleaner way  -- at the searchlight level, where the entire<br>
searchlight then estimated with the same permutation at a time.<br>
<br>
> My code:<br>
> clf  = LinearCSVMC()<br>
> splt = NFoldPartitioner(attr='chunks')<br>
<br>
> repeater = Repeater(count=100)<br>
> permutator = AttributePermutator('targets', limit={'partitions': 1},<br>
> count=1)<br>
> null_cv = CrossValidation(clf, ChainNode([splt,<br>
> permutator],space=splt.get_space()),<br>
>                           postproc=mean_sample())<br>
> null_sl = sphere_searchlight(null_cv, radius=3, space='voxel_indices',<br>
>                              enable_ca=['roi_sizes'])<br>
> distr_est = MCNullDist(repeater,tail='left', measure=null_sl,<br>
>                        enable_ca=['dist_samples'])<br>
<br>
I see.. you are trying to maintain the same assignment in testing dataset...<br>
IMHO (especially since you seems to just use betas from GLM, so I guess a beta<br>
per run) it is not necessary.  Just make a straightforward permutator in all<br>
chunks (but within each chunk) at the level of the searchlight without trying<br>
to do that fancy dance (I know -- the one our tutorial tutors atm):<br>
<br>
        permutator = AttributePermutator('targets', count=100,<br>
                                         limit='chunks')<br>
        distr_est = MCNullDist(permutator, tail='left',<br>
                               enable_ca=['dist_samples'])<br>
<br>
> cv = CrossValidation(clf,splt,<br>
>                      enable_ca=['stats'], postproc=mean_sample() )<br>
> sl = sphere_searchlight(cv, radius=3, space='voxel_indices',<br>
>                         null_dist=distr_est,<br>
>                         enable_ca=['roi_sizes'])<br>
> ds = glm_dataset.copy(deep=False,<br>
>                        sa=['targets','chunks'],<br>
>                        fa=['voxel_indices'],<br>
>                        a=['mapper'])<br>
> sl_map = sl(ds)<br>
> p_values = distr_est.cdf(sl_map.samples) # IS THIS THE RIGHT WAY??<br>
<br>
just access   sl.ca.null_prob   for the p value (you could also use .ca.null_t<br>
(if you enable it for searchlight) to get a corresponding t ...  actually<br>
z-score (i.e. it is not t-score since assumption of the distribution is<br>
normal), which would be easier to visualize and comprehend (2.3 must correspond<br>
to p=0.01 at a one-sided test ;-))<br>
<br>
<br>
> Is there a way to make sure the permutations are exhaustive?<br>
<br>
if they are exhaustive -- you might need more data ;)  especially with<br>
the searchlight, 100 permutations would just give you p no smaller than ~0.01<br>
(well 1/101 to be precise).  with any correction for multiple comparisons (you<br>
have many searchlights) you would need "stronger" p's , thus more permutations.<br>
<br>
> In order to make an inference on the group level I understand I can<br>
> use GroupClusterThreshold.<br>
<br>
Correct<br>
<br>
> Does anyone have a code sample for that? Do I use the MCNullDist's created<br>
> at the subject level?<br>
<br>
Unfortunately we don't have a full nice example yet, and Michael whom we could<br>
blame (d'oh -- thank!) for this functionality is on vacation (CCing though the<br>
very original author -- Richard, please correct me if I am wrong<br>
anywhere) BUT your code above, with my addition (note me enabling<br>
'dist_samples') is pretty much ready.<br>
<br>
1. Just run your searchlights per subject e.g. 20-50 times, collect per each<br>
subject  distr_est.ca.dist_samples, and place them all into a single dataset<br>
where each permutation will be a "sample", while you set sa.chunks  to group<br>
subjects (i.e. all permutations from subject 1 should have chunks == 1).<br>
Call it e.g.   perms<br>
<br>
2. Create bootstrapping of those permutation results: e.g.<br>
<br>
clthr = gct.GroupClusterThreshold() # defaults might be adequate ;), otherwise adjust<br>
<br>
3. Train it on your collection of<br>
<br>
clthr.train(perms)<br>
<br>
which will do all the bootstrapping (would take awhile)<br>
<br>
4. Estimate significance/treshold your original map (mean of errors across<br>
subjects without permutations)<br>
<br>
res = clthr(mean_map)<br>
<br>
<br>
5. look into res.a  for all kinds of stats (# of clusters, their locations, significance etc)<br>
and then<br>
<br>
res.fa.clusters_fwe_thresh<br>
<br>
will contain actual map of clusters indices which passed fwe thresholding.<br>
<br>
<br>
Hope this helps!<br>
<span class="HOEnZb"><font color="#888888">--<br>
Yaroslav O. Halchenko, Ph.D.<br>
<a href="http://neuro.debian.net" rel="noreferrer" target="_blank">http://neuro.debian.net</a> <a href="http://www.pymvpa.org" rel="noreferrer" target="_blank">http://www.pymvpa.org</a> <a href="http://www.fail2ban.org" rel="noreferrer" target="_blank">http://www.fail2ban.org</a><br>
Research Scientist,            Psychological and Brain Sciences Dept.<br>
Dartmouth College, 419 Moore Hall, Hinman Box 6207, Hanover, NH 03755<br>
Phone: <a href="tel:%2B1%20%28603%29%20646-9834" value="+16036469834">+1 (603) 646-9834</a>                       Fax: <a href="tel:%2B1%20%28603%29%20646-1419" value="+16036461419">+1 (603) 646-1419</a><br>
WWW:   <a href="http://www.linkedin.com/in/yarik" rel="noreferrer" target="_blank">http://www.linkedin.com/in/yarik</a><br>
</font></span></blockquote></div><br></div>