| title | Computing and using estimates of effect size | ||
|---|---|---|---|
| tags |
|
The following code demonstrates how you can compute and plot the effect size.
% find the interesting segments of data
cfg = [];
cfg.dataset = 'Subject01.ds';
cfg.trialdef.eventtype = 'backpanel trigger';
cfg.trialdef.prestim = 1;
cfg.trialdef.poststim = 2;
cfg.trialdef.eventvalue = [3 5 9];
% 3 = FIC
% 5 = FC
% 9 = IC
cfg = ft_definetrial(cfg);
% preprocess the data
cfg.channel = {'MEG', '-MLP31', '-MLO12'};
cfg.demean = 'yes';
cfg.baselinewindow = [-0.2 0];
% cfg.lpfilter = 'yes';
% cfg.lpfreq = 35;
data = ft_preprocessing(cfg);
%%
cfg = [];
cfg.keeptrials = 'yes';
cfg.trials = (data.trialinfo==3);
timelock_FIC = ft_timelockanalysis(cfg, data);
cfg.trials = (data.trialinfo==5);
timelock_FC = ft_timelockanalysis(cfg, data);
cfg.trials = (data.trialinfo==9);
timelock_IC = ft_timelockanalysis(cfg, data);
%%
cfg = [];
cfg.parameter = 'trial';
cfg.method = 'analytic';
cfg.statistic = 'indepsamplesT';
cfg.design = [1*ones(1,size(timelock_FC.trial,1)) 2*ones(1,size(timelock_FIC.trial,1))];
cfg.ivar = 1;
stat_FCvsFIC = ft_timelockstatistics(cfg, timelock_FC, timelock_FIC);
%%
cfg = [];
cfg.layout = 'CTF151_helmet';
cfg.parameter = 'stat';
ft_multiplotER(cfg, stat_FCvsFIC);
%%
% MLF32 shows a large positive p-value that peaks around 600ms
chansel = match_str(timelock_FC.label, 'MLF32');
timesel = nearest(timelock_FC.time, 0.6);
%%
x1 = timelock_FC.trial(:, chansel, timesel)*1e12;
x2 = timelock_FIC.trial(:, chansel, timesel)*1e12;
x1 = mean(mean(x1,3),2);
x2 = mean(mean(x2,3),2);
n1 = length(x1);
n2 = length(x2);
if n1==n2
% this fails if x1 and x2 are of different length
figure
hist([x1 x2], 50); legend({'FC', 'FIC'})
end
pooled_sd = sqrt( ((n1-1)*std(x1)^2 + (n2-1)*std(x2)^2) / (n1+n2-1) );
cohensd = (mean(x1)-mean(x2)) / pooled_sd
% see https://en.wikipedia.org/wiki/Effect_size#Cohen.27s_d
% Very small 0.01
% Small 0.20
% Medium 0.50
% Large 0.80
% Very large 1.20
% Huge 2.00
It is interesting to see how the effect size increases by taking the average over more channels and time points.
%%
chansel = match_str(timelock_FC.label, {'MLF22', 'MLF23', 'MLF32', 'MLF33', 'MLF42', 'MLF43', 'MLF52'});
timesel = nearest(timelock_FC.time, 0.496607) : nearest(timelock_FC.time, 0.765893);
% now repeat the computation of Cohen's d above
Proper preprocessing of the data also increases the effect size.
%%
cfg = [];
cfg.method = 'summary';
data = ft_rejectvisual(cfg, data);
% now repeat the computation of Cohen's d above
% this requires the following code for plotting
edges = linspace(-0.4, 0.4, 30);
h1 = histcounts(x1,edges);
h2 = histcounts(x2,edges);
bar(edges(1:end-1),[h1; h2]'); legend({'FC', 'FIC'})