0todd0000 / spm1dmatlab

One-Dimensional Statistical Parametric Mapping in Matlab.
GNU General Public License v3.0
28 stars 13 forks source link

Calculate P values at a point below the critical threshold. #168

Closed ArjS21 closed 1 year ago

ArjS21 commented 2 years ago

Hi, I have used a general linear model implemented using SPM and from the SPM(t) plot, there is one supra-threshold cluster, and another region where the curve approaches the critical threshold closely, but does not cross the threshold. I would like to calculate the P value at this point as a descriptive variable and refer to it when qualitatively describing my data. I have looked at the spmi variable and outputs but am unsure how I can calculate the P value at this region using the outputs available. Can this be done with the t* and z outputs?

Thanks, A

0todd0000 commented 2 years ago

Yes, this can be done, but please note:

That being said, assuming two-tailed inference, you can calculate pointwise p-values for the t-statistic like this in Python:

spm = spm1d.stats.glm(y, X, c)
p   = spm1d.rft1d.f.sf(spm.z**2, spm.df, spm.Q, spm.fwhm, withBonf=True)

I think that should work in MATLAB (after changing spm.z**2 to spm.z.^2 and using 'withBonf', true) but let me know if you're unsure.

ArjS21 commented 2 years ago

Thanks Todd. I have changed the variables mentioned to MATLAB syntax, but get a matrix dimensions error, as below. I'm not too sure on where this is stemming from.

Matrix dimensions must agree.

Error in spm_P_RF (line 82) EM = (R./G).*P; % over D dimensions

Error in spm_P (line 51) [P,p,Em,En,EN] = spm_P_RF(c,k,Z,df,STAT,R,n);

Error in spm1d.rft1d.RFTProbability/upcrossing (line 23) p = spm_P(1, 0, u, c.df, c.STAT, c.resels, c.n, c.Q);

Error in spm1d.rft1d.RFTCalculator/sf (line 58) p = self.p.upcrossing(u);

Error in spm1d.rft1d.f.sf (line 38) p = calc.sf(u);

0todd0000 commented 2 years ago

Apologies for the delay! This slipped through my inbox.

It seems that the MATLAB sf function cannot accept arrays so instead try calculating for each time point like this:

Q = spm.Q;
p = zeros(1,Q);
for q = 1:Q
    p(q) = spm1d.rft1d.f.sf(spm.z(q)^2, spm.df, Q, spm.fwhm);
end