search

0todd0000 / spm1dmatlab

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

SPM1D ttest_paired limit number subject ? #94

Closed andriacapai closed 5 years ago

andriacapai commented 5 years ago

Hello,

I encountered a problem after using the function spm1d.stats.normality.ttest_paired with the function spm.inference(0.05); .

Here is the error : "Error using ~ Operands must be real."

Here is my two matrix YA and YB :

I tried to use only 75 strides for each subject (Matrix YA and YB = 150x200) and it works ... Do you have any idea about how to unlock the limit number ? Or in case the problem is not about the number could you advise me what is the best test according to my data set ?

Thank you.

Best,

Andria

0todd0000 commented 5 years ago

Hi Andria, there is no limit on N, so I suspect that there is a problem in data.

Can you please copy and send the full error? The full error will give information about what function and what line of code is causing the problem, and this will help me figure out what's going on.

Also, can you please copy and send the output of:

t = spm1d.stats.ttest_paired(YA, YB);
disp(t)

Thanks!

andriacapai commented 5 years ago

Hello Todd,

Here is the entire error message : """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" Error using ~ Operands must be real.

Error in spm_P_RF (line 94) if ~k || ~D

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/cluster (line 15) p = spm_P(1, k, u, c.df, c.STAT, c.resels, c.n, c.Q);

Error in spm1d.rft1d.chi2.p_cluster_resels (line 50) p = calc.p.cluster(k, u);

Error in spm1d.geom.Cluster/inference (line 55) p = spm1d.rft1d.chi2.p_cluster_resels(k, u, df(2), resels, 'withBonf',withBonf, 'nNodes',nNodes);

Error in spm1d.stats.spm.SPM/cluster_inference (line 124) clusters{i} = clusters{i}.inference(self.STAT, self.df, self.fwhm, self.resels, two_tailed, withBonf, self.nNodes);

Error in spm1d.stats.spm.SPM/inference (line 92) [clusters,p] = self.cluster_inference(clusters, two_tailed, withBonf);

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

And here is the the output of:

t = spm1d.stats.ttest_paired(YA, YB); disp(t)

SPM{t} z: [1x200 double] df: [1 199] fwhm: 78.5959 resels: [1 2.5319]

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

And the last picture that I send you (attached to this email) is the result about :

spm1d.stats.normality.ttest_paired(YA,YB);

Thank you so much for your responsiveness.

Best regards,

Andria Capai

De : Todd Pataky notifications@github.com Envoyé : mardi 11 décembre 2018 14:53:37 À : 0todd0000/spm1dmatlab Cc : CAPAI Andria; Author Objet : Re: [0todd0000/spm1dmatlab] SPM1D ttest_paired limit number subject ? (#94)

Hi Andria, there is no limit on N, so I suspect that there is a problem in data.

Can you please copy and send the full error? The full error will give information about what function and what line of code is causing the problem, and this will help me figure out what's going on.

Also, can you please copy and send the output of:

t = spm1d.stats.ttest_paired(YA, YB); disp(t)

Thanks!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/0todd0000/spm1dmatlab/issues/94#issuecomment-446209733, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ArpSN7gMizvwdgxcrPlIcCeYknkz1AZwks5u37jhgaJpZM4ZLpGq.

0todd0000 commented 5 years ago

Thank you for the details. I'm not yet sure what the problem is, but I think if we can figure out what the values are in for k and u as inputs to spm1d.rft1d.chi2.p_cluster_resels, then we will solve the problem. The p_clister_resels function accepts the following arguments:

From the disp(t) output you copied above, there is no problem with df or resels, so the problem must be with k or u. The easiest way to get these values will be to edit the source code in the file:

./+spm1d/+rft1d/chi2.m

Please open this chi2.m file in Matlab. On Line 47 you should see: "function [p] = p_cluster_resels(k, u, df, resels, varargin)"

Immediately after this line, please add two new lines:

disp(k)
disp(u)

Then re-run your analysis. At the top of the error message you should see two printed values (one for k and one for u). Please tell me what these values are, then please delete the two lines you added above.

andriacapai commented 5 years ago

Hello Todd,

I followed your instructions and here is the output of :

I guess it's about the complex number ?

Thank you,

Andria CAPAI

De : Todd Pataky notifications@github.com Envoyé : mercredi 12 décembre 2018 19:19:55 À : 0todd0000/spm1dmatlab Cc : CAPAI Andria; Author Objet : Re: [0todd0000/spm1dmatlab] SPM1D ttest_paired limit number subject ? (#94)

Thank you for the details. I'm not yet sure what the problem is, but I think if we can figure out what the values are in for k and u as inputs to spm1d.rft1d.chi2.p_cluster_resels, then we will solve the problem. The p_clister_resels function accepts the following arguments:

From the disp(t) output you copied above, there is no problem with df or resels, so the problem must be with k or u. The easiest way to get these values will be to edit the source code in the file:

./+spm1d/+rft1d/chi2.m

Please open this chi2.m file in Matlab. On Line 47 you should see: "function [p] = p_cluster_resels(k, u, df, resels, varargin)"

Immediately after this line, please add two new lines:

disp(k) disp(u)

Then re-run your analysis. At the top of the error message you should see two printed values (one for k and one for u). Please tell me what these values are, then please delete the two lines you added above.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/0todd0000/spm1dmatlab/issues/94#issuecomment-446689854, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ArpSN_dVSQkImnTi6jqTD12TCCL2xBw7ks5u4UjLgaJpZM4ZLpGq.

0todd0000 commented 5 years ago

Yes, the complex number is causing the problem. This is a curious bug! I have no idea how a complex number can be generated. There must be some problem in YA and/or YB, but I'm not sure what it is. To figure out the problem, here are some questions / requests:

spm1d.plot.plot_meanSD(YA, 'color','k');
hold on
spm1d.plot.plot_meanSD(YB, 'color','r');
t = spm1d.stats.ttest_paired(YA, YB);
t.plot()
andriacapai commented 5 years ago

Hi Todd,

I followed your instruction and I did 2 other things (see below) , then here are results and figure attached in this message :

-spm1d.plot.plot_meanSD(YA, 'color','k'); hold on spm1d.plot.plot_meanSD(YB, 'color','r'); (See Figure : meanstd_ya_yb)

-t = spm1d.stats.ttest_paired(YA, YB); t.plot() (See figure results_ttest_pairesyayb)

Thanks

Best, Andria

0todd0000 commented 5 years ago

Thank you for these details. Everything looks OK in these figures, so I'm still not sure what the problem is... Here are two more questions:

0todd0000 commented 5 years ago

One more request: Are you able to share the YA and YB data? If yes, please either attach a data file to this issue, or send the file to me privately at pataky.todd.2m at kyoto-u.ac.jp (Without the data it might be difficult to find the source of the problem.) Todd

andriacapai commented 5 years ago

Hello Todd,

Here are the answers :

And I sent you data by email.

Andria

0todd0000 commented 5 years ago

Hi Andria, I tried analyses in Python, and analysis executes without any errors. I don't have access to Matlab right now, but Mark or Jos will test things out in Matlab then reply to this forum issue. Todd

m-a-robinson commented 5 years ago

Hello Andria, Todd,

Well...to deepen the mystery further I am able to reproduce this error in Matlab (2018a). I have checked the data, plotted it and when running the normality test it stops working when the following data is included

spm = spm1d.stats.normality.ttest_paired(yB(1:164,:), yA(1:164,:));

Specifically the spm.z variable is a 1 x 200 complex double Further investigation seems to indicate it is the last number that is the issue -4.690928480235349e+04 + 6.849996087265469e+04i, which is confirmed as the below code executes fine.

spm = spm1d.stats.normality.ttest_paired(yB(1:164,1:199), yA(1:164,1:199));

Further inspection of the spm object gives no other clues to the origin of this number in any of the other data so I had a look in the ttest_paired function in the +normality folder and then ran the code in the k2residuals.m file.

Unsurprisingly it was this code that generated the complex number

    for i = 1:Q
        k2(i) = spm1d.stats.normality.k2_single_node( x(:,i) );
    end

A few data points from k2 are copied below. Which seems to show complex numbers across the array and not just the last value.

825.629173732941 + 0.00000000000000i 565.388411199511 + 0.00000000000000i 2296.44269873964 + 0.00000000000000i -4042.17882702022 + 3566.56394302302i 3974.46963820205 + 0.00000000000000i 1354.62841700969 + 0.00000000000000i 1040.73472716385 + 0.00000000000000i 783.079057243054 + 0.00000000000000i 597.873608650287 + 0.00000000000000i 388.363503781647 + 0.00000000000000i 234.540107676453 + 0.00000000000000i 159.773702440789 + 0.00000000000000i 136.943013109852 + 0.00000000000000i 143.051717524228 + 0.00000000000000i 167.034382931555 + 0.00000000000000i 194.088406195736 + 0.00000000000000i 223.130863954024 + 0.00000000000000i 283.234806765852 + 0.00000000000000i 456.478039071591 + 0.00000000000000i 1143.97818683881 + 0.00000000000000i -2958.90847559263 + 2282.17818717165i -1387.85415154322 + 638.339209455237i -1201.10289464189 + 472.686282882775i -1219.92830418537 + 489.011330343929i -1337.57043336625 + 593.010723950451i -1517.04206911690 + 758.417078752748i -1913.63297595630 + 1147.65801262948i -4900.53929221056 + 4636.37574504929i 1426.78845666040 + 0.00000000000000i 580.690356908692 + 0.00000000000000i

So it looks like something in the k2_single_node function but I'm stuck here for the time being.

Regards Mark

0todd0000 commented 5 years ago

How bizarre. I think I've found the source of the problem, but not yet a solution.

I think the problem is a coding error in DagosPtest.m, which is available on the MathWorks File Exchange at the link below. spm1d uses a modified version of this code in its normality calculations. https://www.mathworks.com/matlabcentral/fileexchange/3954-dagosptest

First, to test this bug, I suggest we use the sample data file at the link below, which is the fourth column of Andria's data, and which produces numerical problems. If you plot a histogram of these data you'll find a bimodal distribution. Clearly these data are non-normal and should yield a large value for the K2 statistic. sample_data.txt

hist

I think the coding problem stems from the kurtosis “Zg2” calculation on Line 125 in DagosPtest.m:

Zg2 = (1-(2/(9*K))-L^(1/3))/sqrt(2/(9*K));

The numerical problem is that a negative L value produces a complex number or NaN. I verified that the sample data above produces a negative L value, and thus disrupts this kurtosis calculation.

Thus the L value (computed on Line 124) must contain an error, but the error could be propagated from previous lines.

Consider Eqn. 5 in Anscombe & Glynn (1983), which represents “Zg2”. It would appear that the DagosPtest.m implementation of this equation is correct, where L represents the contents of the [square brackets]. The symbol translations are:

Anscombe & Glynn (1983) DagosPtest.m
Eqn (5) Zg2
[square brackets in Eqn (5)] L
A K
x H

I think that Zg2, L and K calculations are correct, but I'm not sure about H. It's difficult for me to go much farther than this without access to Matlab.

If we can't figure it out relatively quickly I suggest we write to the author of DagosPtest.m to report the bug.

andriacapai commented 5 years ago

Hello Mark and Todd,

Thank you for yours answers and details you gave us !

Hope we will find the problem

Best,

[1505403877901_ISM][1505404771538_Aix]

Andria CAPAI - Ingénieur d'études UMR 7287.Sc. du Mouvement (EJ Marey) Institut des Sciences du Mouvement - 163 Avenue de Luminy - 13009 Marseille Tél: +33(0)6 38 23 96 09<tel:+33%206%2038%2023%2096%2009>

De : Mark Robinson notifications@github.com Envoyé : mercredi 19 décembre 2018 23:47:42 À : 0todd0000/spm1dmatlab Cc : CAPAI Andria; Author Objet : Re: [0todd0000/spm1dmatlab] SPM1D ttest_paired limit number subject ? (#94)

Hello Andria, Todd,

Well...to deepen the mystery further I am able to reproduce this error in Matlab (2018a). I have checked the data, plotted it and when running the normality test it stops working when the following data is included

spm = spm1d.stats.normality.ttest_paired(yB(1:164,:), yA(1:164,:));

Specifically the spm.z variable is a 1 x 200 complex double Further investigation seems to indicate it is the last number that is the issue -4.690928480235349e+04 + 6.849996087265469e+04i, which is confirmed as the below code executes fine.

spm = spm1d.stats.normality.ttest_paired(yB(1:164,1:199), yA(1:164,1:199));

Further inspection of the spm object gives no other clues to the origin of this number in any of the other data so I had a look in the ttest_paired function in the +normality folder and then ran the code in the k2residuals.m file.

Unsurprisingly it was this code that generated the complex number

for i = 1:Q
    k2(i) = spm1d.stats.normality.k2_single_node( x(:,i) );
end

A few data points from k2 are copied below. Which seems to show complex numbers across the array and not just the last value.

825.629173732941 + 0.00000000000000i 565.388411199511 + 0.00000000000000i 2296.44269873964 + 0.00000000000000i -4042.17882702022 + 3566.56394302302i 3974.46963820205 + 0.00000000000000i 1354.62841700969 + 0.00000000000000i 1040.73472716385 + 0.00000000000000i 783.079057243054 + 0.00000000000000i 597.873608650287 + 0.00000000000000i 388.363503781647 + 0.00000000000000i 234.540107676453 + 0.00000000000000i 159.773702440789 + 0.00000000000000i 136.943013109852 + 0.00000000000000i 143.051717524228 + 0.00000000000000i 167.034382931555 + 0.00000000000000i 194.088406195736 + 0.00000000000000i 223.130863954024 + 0.00000000000000i 283.234806765852 + 0.00000000000000i 456.478039071591 + 0.00000000000000i 1143.97818683881 + 0.00000000000000i -2958.90847559263 + 2282.17818717165i -1387.85415154322 + 638.339209455237i -1201.10289464189 + 472.686282882775i -1219.92830418537 + 489.011330343929i -1337.57043336625 + 593.010723950451i -1517.04206911690 + 758.417078752748i -1913.63297595630 + 1147.65801262948i -4900.53929221056 + 4636.37574504929i 1426.78845666040 + 0.00000000000000i 580.690356908692 + 0.00000000000000i

So it looks like something in the k2_single_node function but I'm stuck here for the time being.

Regards Mark

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/0todd0000/spm1dmatlab/issues/94#issuecomment-448771951, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ArpSN2rjacVOvuhSNNZRAFsZrxJTOaDOks5u6sIOgaJpZM4ZLpGq.

m-a-robinson commented 5 years ago

Hi both,

Thanks @0todd0000 for continuing down the Rabbit hole. I ran the code in k2_single_node.m and was able to replicate the negative L.

H was -4.3844 and the final k2 statistic was a complex number -4.04..e+03 + 3.56..e+03i

I had a look at the calculation and the original paper but couldn't get to the bottom of this quickly so it might be worth emailing the author as you suggest.

Cheers Mark

0todd0000 commented 5 years ago

I just posted a bug report to the MathWorks File Exchange: https://www.mathworks.com/matlabcentral/fileexchange/3954-dagosptest

andriacapai commented 5 years ago

Thanks a lot for that

Andria

[1505403877901_ISM][1505404771538_Aix]

Andria CAPAI - Ingénieur d'études UMR 7287.Sc. du Mouvement (EJ Marey) Institut des Sciences du Mouvement - 163 Avenue de Luminy - 13009 Marseille Tél: +33(0)6 38 23 96 09<tel:+33%206%2038%2023%2096%2009>

De : Todd Pataky notifications@github.com Envoyé : vendredi 21 décembre 2018 00:59:44 À : 0todd0000/spm1dmatlab Cc : CAPAI Andria; Author Objet : Re: [0todd0000/spm1dmatlab] SPM1D ttest_paired limit number subject ? (#94)

I just posted a bug report to the MathWorks File Exchange: https://www.mathworks.com/matlabcentral/fileexchange/3954-dagosptest

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/0todd0000/spm1dmatlab/issues/94#issuecomment-449181882, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ArpSNyHhxQdnYJ1ZojAap_z4aqnYvQ5Dks5u7CRwgaJpZM4ZLpGq.