[FieldTrip] Calculating the variance of a group with given variances of each subject

Schoffelen, J.M. (Jan Mathijs) janmathijs.schoffelen at donders.ru.nl
Sat Oct 1 17:21:46 CEST 2022


Hi Sebastian,

I don’t understand why you think that the variance of ft_timelockgrandaverage leads to wrong statistics. As I said in an earlier e-mail, the traditional way to do second-level statistics is to compute the variance of the mean-per-subject, and not to compute the variance of the individual-trials-concatenated-across-subjects.

The computation that you call ‘real mean and variance’ assume the unit of observation to be the trials, while the unit-of-observation that ft_timelockgrandaverage uses for the mean and variance computation is the subject. While the mean is not going to change (assuming the same number of trials per subject), the variance will turn out differently, because that results from a non-linear operation (squaring).

In your numeric example, the 'FT-code' leads to a rejection of the null-hypothesis due to the fact that you use the incorrect n1 and n2 in the numerator of the t-statistic, these should be 2 and 2, rather than 100 and 100.

Best wishes,
Jan-Mathijs





On 1 Oct 2022, at 08:48, Sebastian.Neudek at med.uni-duesseldorf.de<mailto:Sebastian.Neudek at med.uni-duesseldorf.de> wrote:


Hi Jan-Mathijs,

thank you for your suggestion of mixed-effect models. Probably those are the perfect fit for my case.
But I still think the variance implementation in fieldtrip in ft_timelockgrandaverage is not correct (only calculating the variance of means of subjects as group variance). I looked a little bit deeper into the implementations and it seems the degress of freedom are the number of trials. Therefore the implemented mean in ft_timelockgrandaverage already is the mean from a pooled variance.

Here you have an example matlab script of which I think the variance of ft_timelockgrandaverage leads to wrong statistics. There are only four subjects who belong to 2 groups (Subject 1 and 2 in group A, Subject 3 and 4 in group B). The trials from all subjects are drawn from the same distribution. All have same mean and same variance. Therefore the real group means are the same and don't differ. Here a t-test is evaluated, to test whether the means of the groups differ.



rng(0,'twister');

% all people have the same number of trials and the same(!) mean
N = 100;

n1 = N; % Person 1 belongs to group A
mu1 = 0;
sigma1 = 1;

n2 = N; % Person 2 belongs to group A
mu2 = 0;
sigma2 = 1;

n3 = N; % Person 3 belongs to group B
mu3 = 0;
sigma3 = 1;

n4 = N; % Person 4 belongs to group B
mu4 = 0;
sigma4 = 1;

% normal distributed random numbers as trials for each person
trials_1 = mu1 + randn([1 n1])*sigma1;
trials_2 = mu2 + randn([1 n2])*sigma2;
trials_3 = mu3 + randn([1 n3])*sigma3;
trials_4= mu4 + randn([1 n4])*sigma4;

% calculate mean and variance for each person

avg1 = mean(trials_1);
avg2 = mean(trials_2);
avg3 = mean(trials_3);
avg4 = mean(trials_4);

var1 = var(trials_1);
var2 = var(trials_2);
var3 = var(trials_3);
var4 = var(trials_4);

% real mean and variance

avg_A = mean([trials_1 trials_2]); % avg for group A
avg_B = mean([trials_3 trials_4]); % avg for group B
var_A = var([trials_1 trials_2]);
var_B = var([trials_3 trials_4]);

% mean and variance of fieldtrip

avg_A_ft = mean([avg1 avg2]);
avg_B_ft = mean([avg3 avg4]);
var_A_ft = var([avg1 avg2]);
var_B_ft = var([avg3 avg4]);

% now a t-test, to check if the mean of group A: {person 1 and person 2}
% and group B: {person 3 and person 4} are different
% H0: both groups share the same mean
% H1: both groups have different mean

% https://en.wikipedia.org/wiki/Student%27s_t-test<https://urldefense.com/v3/__https://en.wikipedia.org/wiki/Student*27s_t-test__;JQ!!HJOPV4FYYWzcc1jazlU!6Mx3xotEc6n0WD6ucT1uJXuGoGRb3nBniS8bvTkgu5bBo2C7RGs6aXqagsx4CtsBeIlIV_8l_5FZv7IoifABgNu85rFnFfmYfPKd5WAnGaP2sxaQRB88G4Q$>
alpha = 0.05;

t = avg_A-avg_B/(sqrt((var_A+var_B)/2)*sqrt(2/(n1+n2))); % this formula is only correct for same group sizes
%p = 1-tcdf(t,n1+n2+n3+n4-4) % right tailed t-test
%p = tcdf(t,n1+n2+n3+n4-4) % left tailed t-test
p = 2*tcdf(-abs(t),n1+n2+n3+n4-4); % two tailed t-test

fprintf('Result:\n');
if p < alpha
    fprintf('H0 rejected. Significant difference of means between Group A and B (P = %.4f, T = %.4f)\n\n',p ,t);
else
    fprintf('H0 not rejected. No significant difference of means between Group A and B (P = %.4f, T = %.4f)\n\n',p ,t);
end


% t value with fieldtrips variance

t_ft = avg_A_ft-avg_B_ft/(sqrt((var_A_ft+var_B_ft)/2)*sqrt(2/(n1+n2))); % this formula is only correct for same group sizes
% p_ft = 1-tcdf(t_ft,n1+n2+n3+n4-4) % right tailed t-test
% p_ft = tcdf(t_ft,n1+n2+n3+n4-4) % left tailed t-test
p_ft = 2*tcdf(-abs(t_ft),n1+n2+n3+n4-4); % two tailed t-test



fprintf('Result using variance of ft_timelockgrandaverage:\n');
if p_ft < alpha
    fprintf('H0 rejected. Significant difference of means between Group A and B (P = %.4f, T = %.4f)\n\n',p_ft ,t_ft);
else
    fprintf('H0 not rejected. No significant difference of means between Group A and B (P = %.4f, T = %.4f)\n\n',p_ft ,t_ft);
end


Best,
Sebastian





________________________________
Von: fieldtrip <fieldtrip-bounces at science.ru.nl<mailto:fieldtrip-bounces at science.ru.nl>> im Auftrag von Schoffelen, J.M. (Jan Mathijs) via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>>
Gesendet: Freitag, 30. September 2022 13:46:40
An: FieldTrip discussion list
Cc: Schoffelen, J.M. (Jan Mathijs)
Betreff: [Extern] - Re: [FieldTrip] Calculating the variance of a group with given variances of each subject

Hi Sebastian,

I would say that for visualization purposes (and for traditional parametric second-level statistics) you are interested in the variance across your units-of-observation. In this case across subjects. This is indeed what ft_timelockgrandaverage gives you. If you are into more sophisticated estimates, it smells a bit as if you want to be moving into the direction of mixed-effect models, which are gaining quite a bit of momentum these days. Since I am a bit old-school I have not eaten any cheese from that (http://www.dwotd.nl/2008/02/360-daar-heb-ik-geen-kaas-van-gegeten.html<https://urldefense.com/v3/__http://www.dwotd.nl/2008/02/360-daar-heb-ik-geen-kaas-van-gegeten.html__;!!HJOPV4FYYWzcc1jazlU!6Mx3xotEc6n0WD6ucT1uJXuGoGRb3nBniS8bvTkgu5bBo2C7RGs6aXqagsx4CtsBeIlIV_8l_5FZv7IoifABgNu85rFnFfmYfPKd5WAnGaP2sxaQG6s86NA$>), so perhaps somebody else wants to chime in here.

Best wishes,
Jan-Mathijs


On 26 Sep 2022, at 21:36, Sebastian Neudek via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>> wrote:


Hi all,

I am currently writing my master's thesis using FieldTrip and might need a bit of your expertise for a little hurdle I encountered.

First a few words about what I want to do:
I have data from a MEG-go/nogo study: There are multiple subjects and multiple trials for each condition. I now want to plot the ERF for each condition averaged across all subjects. But together with the ERF I also want to plot the  standard error or variance or confidence interval or standard deviation. Something like this(Figure 3a):
https://www.researchgate.net/figure/ERP-results-including-P3-amplitudes-a-Grand-average-and-standard-errors-of-ERP_fig2_327651209<https://urldefense.com/v3/__https://www.researchgate.net/figure/ERP-results-including-P3-amplitudes-a-Grand-average-and-standard-errors-of-ERP_fig2_327651209__;!!HJOPV4FYYWzcc1jazlU!4AKDCxWgDNOVGlWwimPV3Sc-5suVMj59sPhwC-eHwI7h0FQr10-YxE5fzliZsQft3R0I3s5Swpei12UfNPYgMdNHwshEvUe7t8p0Gw$>
Independent which of these measures I finally use, I need to calculate the variance of the ERF at each timepoint.
To do this I first used ft_timelockanalysis to calculate the averages of each subject. ft_timelockanalysis returns not only the average, but also the variance of the ERF. Afterwards I use ft_timelockgrandaverage to calculate the group average of the ERF.

And this is where my troubles begin. I get how ft_timelockgrandaverage calulates the new average: It averages over the mean values of the subjects (it also accounts for the degrees of freedom). But for the variance it just calculates the variance of the means. I thought, it will also account for the variances of each subject.
In the extreme case it leads to the following:
Imagine you have two subjects and you measure very close averages avg(sub1)=1 and avg(sub2)=1.0001. But their variances are extremly high: var(sub1)=var(sub2)=100. Then the variance calculated by ft_timelockgrandaverage is almost zero, when in reality it should be something like 100.
Maybe this should be noted in ft_timelockgrandaverage, because it can lead to wrong statistics if the user isn't aware of this effect.

Back to my task:
This is surely not the variance I want to use for my ERF plots and therefore I searched for a different calculation.
I found on wikipedia a formula for so called 'pooled variance' (at the bottom of the page for sample-based statistics):
https://en.wikipedia.org/wiki/Pooled_variance<https://urldefense.com/v3/__https://en.wikipedia.org/wiki/Pooled_variance__;!!HJOPV4FYYWzcc1jazlU!4AKDCxWgDNOVGlWwimPV3Sc-5suVMj59sPhwC-eHwI7h0FQr10-YxE5fzliZsQft3R0I3s5Swpei12UfNPYgMdNHwshEvUeBLHOpZA$>

But I am not 100% happy with this pooled variance.
First of all the deegrees of freedom is not accounted for. Then there is a different mean calculated than in ft_timelockgrandaverage: The mean is accounting for the number of trials for each subject. I don't know if the different mean is a pro or a con (On the one hand, means with a higher precisision because of a higher number of trials are weighted higher than means with a lower trial number. On the other hand, each subject on the group should be weighted same).

I also found a second formula on stackexchange, which also calculates a variances of variances:
https://stats.stackexchange.com/questions/300392/calculate-the-variance-from-variances<https://urldefense.com/v3/__https://stats.stackexchange.com/questions/300392/calculate-the-variance-from-variances__;!!HJOPV4FYYWzcc1jazlU!4AKDCxWgDNOVGlWwimPV3Sc-5suVMj59sPhwC-eHwI7h0FQr10-YxE5fzliZsQft3R0I3s5Swpei12UfNPYgMdNHwshEvUdcUgTl6A$>

Therefore I am uncertain which formula to use, because as it seems, they differ.


I hope somebody of you is a little more experienced in statistics than I am and can help me out finding the correct or best calculation for this variance. Best case would be if it is already implemented in fieldtrip and I only missed it, but I can also implement it myself.

Best,
Sebastian







_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://urldefense.com/v3/__https://doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!4AKDCxWgDNOVGlWwimPV3Sc-5suVMj59sPhwC-eHwI7h0FQr10-YxE5fzliZsQft3R0I3s5Swpei12UfNPYgMdNHwshEvUdau-41SA$

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20221001/baa15b7c/attachment.htm>


More information about the fieldtrip mailing list