[FieldTrip] [Extern] - Re: Calculating the variance of a group with given variances of each subject

Tue Oct 4 11:01:25 CEST 2022

Hi Jan-Mathijs,

thanks for your response. Apparently my example is bad. It will also lead to different p-values if the bugs are corrected, because it violates the independece of samples. That is the reason, why the pooled variance can't be used for statistics.

I totally agree with you with the general procedure, but I had some thoughts about it, which I want to share. I will call the outcomes of this method group-mean, group-variance and the groups standard error of the mean. Also will call the average of one person a 'measurement'. Therefore the groups mean is just the mean of our measurements.
But I think there is an uncertainty to each of these measurements, which in my opinion cannot be neglegted. We can calculate these uncertainties of each measurement by calculating the standard error for each person. The uncertainty of these measurements 'propagates' if they are averaged. In our case the uncertainty for the end measurement is the average of all subjects uncertainties (I found this propagation only on the german wiki-page for propagation of uncertainty. Maybe this should be further checked).

Following the group-mean itself isn't just a value, it is smeared around the value, because we couldn't measure it any more precise. If we plot our ERF, we could change the thickness of our line according to this uncertainty. We now have calculated two errors: Once the groups standard error of the mean and also the uncertainty of our measurement (which is also an error). Therefore we can calculate a total error of our mean if we add both errors together (adding together errors is allowed according to the german wiki-page for propagation of uncertainty). This total error is the error of our ERFs. If we want to compare the ERFs of two groups, we therefore need to calculate back the variance.

With the standard procedure of only calculating the group-variance you don't really get a benefit of measuring multiple trials. With this procedure you reduce the uncertainty of your measurement by using multiple trials per subject and following the total variance of the ERF.

A suggestion would be to rename the variance field from the output ft_timelockgrandaverage to group variance and add a field for the measurement uncertainty. The old variance field would then represent the total variance of the ERF.

Is it already implemented in fieldtrip and I missed it? If not, what are your thoughts about it?

Best,
Sebastian

________________________________
Von: fieldtrip <fieldtrip-bounces at science.ru.nl> im Auftrag von Schoffelen, J.M. (Jan Mathijs) via fieldtrip <fieldtrip at science.ru.nl>
Gesendet: Samstag, 1. Oktober 2022 17:21:46
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 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/20221004/20428657/attachment.htm>