So testen Sie, ob die Varianz zweier Verteilungen unterschiedlich ist, wenn die Verteilungen nicht normal sind

Ich untersuche zwei geografisch isolierte Populationen derselben Art. Wenn ich die Verteilungen betrachte, sehe ich, dass beide bimodal sind (es gibt eine gewisse Saisonalität für ihr Auftreten), aber die Peaks in einer Population sind viel höher und viel schmaler (dh die Varianz der lokalen Peaks ist kleiner).

Welche Art von statistischem Test wäre geeignet, um festzustellen, ob diese Unterschiede signifikant sind?

Zur Verdeutlichung ist meine y-Achse die Anzahl der Personen, die an einem bestimmten Tag in einer Falle identifiziert wurden, und die x-Achse ist der julianische Tag.

Sind diese Verteilungen etwas im Laufe der Zeit? Zählt vielleicht? (Wenn ja, dann brauchen Sie vielleicht etwas ganz anderes als die bisherigen Diskussionen hier)

Was Sie beschreiben, hört sich nicht so an, als würde es als Unterschied in der Varianz der Verteilungen sehr gut verstanden.

Es hört sich so an, als würden Sie etwas vage beschreiben (ignorieren Sie die Zahlen auf den Achsen, um nur einen Eindruck von der allgemeinen Art von Muster zu vermitteln, die Sie zu beschreiben scheinen):

Wenn das stimmt, dann überlegen Sie:

Während die Breite jedes Peaks um die lokalen Zentren für die blaue Kurve schmaler ist, unterscheidet sich die Varianz der roten und blauen Verteilung insgesamt kaum.

Wenn Sie die Modi und Antimoden vorher identifizieren, können Sie die lokale Variabilität messen.

Zunächst denke ich, dass Sie die saisonalen Verteilungen separat betrachten sollten, da die bimodale Verteilung wahrscheinlich das Ergebnis von zwei ziemlich getrennten Prozessen ist. Die zwei Verteilungen könnten durch unterschiedliche Mechanismen gesteuert werden, so dass beispielsweise Winterverteilungen empfindlicher auf das jährliche Klima reagieren könnten. Wenn Sie die Bevölkerungsunterschiede und die Gründe dafür betrachten möchten, ist es meiner Meinung nach sinnvoller, die saisonalen Verteilungen separat zu untersuchen.

Für einen Test können Sie den Levine-Test (im Grunde ein Test der Homoskedastizität) ausprobieren, mit dem Varianzen zwischen Gruppen verglichen werden. Der Bartlett-Test ist eine Alternative, aber der Levene-Test soll robuster gegenüber Nicht-Normalität sein (insbesondere wenn der Median zum Testen verwendet wird). In R finden sich die Tests von Levene und Bartlett in library(car).

Ich stimme dem zu, was andere gesagt haben - nämlich, dass "Varianz" wahrscheinlich das falsche Wort ist (da die von Ihnen in Betracht gezogene Funktion keine Wahrscheinlichkeitsverteilung, sondern eine Zeitreihe ist).

Ich denke, Sie möchten dieses Problem möglicherweise aus einer anderen Perspektive betrachten - passen Sie einfach die beiden Zeitreihen mit LOWESS-Kurven an. Sie können 95% -Konfidenzintervalle berechnen und deren Formen qualitativ kommentieren. Ich bin mir nicht sicher, ob Sie etwas ausgefalleneres tun müssen.

Ich habe unten einen MATLAB-Code geschrieben, um zu veranschaulichen, was ich sage. Ich bin in Eile, kann aber bald Klarheit schaffen. Vieles von dem, was ich getan habe, kann direkt von hier übernommen werden: http://blogs.mathworks.com/loren/2011/01/13/data-driven-fitting/

%% Generate Example data
npts = 200;
x = linspace(1,100,npts)';
y1 = (1e3*exp(-(x-25).^2/20) + 5e2*exp(-(x-65).^2/40));
y1_noisy = 50*randn(npts,1) + y1;
y2 = (1e3*exp(-(x-25).^2/60) + 5e2*exp(-(x-65).^2/100));
y2_noisy = 50*randn(npts,1) + y2;

figure; hold on
plot(x,y1_noisy,'ob')
plot(x,y2_noisy,'or')
title('raw data'); ylabel('count'); xlabel('time')
legend('y1','y2')

Möglicherweise möchten Sie die beiden Zeitreihen normalisieren, um ihre relativen Trends und nicht ihre absoluten Werte zu vergleichen.

%% Normalize data sets
figure; hold on
Y1 = y1_noisy./norm(y1_noisy);
Y2 = y2_noisy./norm(y2_noisy);
plot(x,Y1,'ob')
plot(x,Y2,'or')
title('normalized data'); ylabel('normalized count'); xlabel('time')
legend('Y1','Y2')

Jetzt machen LOWESS passt ...

%% Make figure with lowess fits
figure; hold on
plot(x,Y1,'o','Color',[0.5 0.5 1])
plot(x,Y2,'o','Color',[1 0.5 0.5])
plot(x,mylowess([x,Y1],x,0.15),'-b','LineWidth',2)
plot(x,mylowess([x,Y2],x,0.15),'-r','LineWidth',2)
title('fit data'); ylabel('normalized count'); xlabel('time')

Schließlich können Sie 95% -Konfidenzbänder wie folgt erstellen:

%% Use Bootstrapping to determine 95% confidence bands
figure; hold on
plot(x,Y1,'o','Color',[0.75 0.75 1])
plot(x,Y2,'o','Color',[1 0.75 0.75])

f = @(xy) mylowess(xy,x,0.15);
yboot_1 = bootstrp(1000,f,[x,Y1])';
yboot_2 = bootstrp(1000,f,[x,Y2])';
meanloess(:,1) = mean(yboot_1,2);
meanloess(:,2) = mean(yboot_2,2);
upper(:,1) = quantile(yboot_1,0.975,2);
upper(:,2) = quantile(yboot_2,0.975,2);
lower(:,1) = quantile(yboot_1,0.025,2);
lower(:,2) = quantile(yboot_2,0.025,2);

plot(x,meanloess(:,1),'-b','LineWidth',2);
plot(x,meanloess(:,2),'-r','LineWidth',2);
plot(x,upper(:,1),':b');
plot(x,upper(:,2),':r');
plot(x,lower(:,1),':b');
plot(x,lower(:,2),':r');
title('fit data -- with confidence bands'); ylabel('normalized count'); xlabel('time')

Jetzt können Sie die endgültige Zahl nach Ihren Wünschen interpretieren und haben die LOWESS-Passungen, um Ihre Hypothese zu untermauern, dass die Peaks in der roten Kurve tatsächlich breiter als die blaue Kurve sind. Wenn Sie eine bessere Vorstellung von der Funktion haben, können Sie stattdessen eine nichtlineare Regression durchführen.

Bearbeiten: Basierend auf einigen hilfreichen Kommentaren unten füge ich einige weitere Details zum expliziten Schätzen von Peakbreiten hinzu. Zunächst müssen Sie eine Definition für das finden, was Sie als "Peak" betrachten. Vielleicht jede Beule, die über eine Schwelle steigt (so etwas wie 0,05 in den Plots, die ich oben gemacht habe). Das Grundprinzip ist, dass Sie einen Weg finden sollten, "echte" oder "bemerkenswerte" Spitzen von Rauschen zu trennen.

Dann können Sie für jeden Peak seine Breite auf verschiedene Arten messen. Wie ich in den Kommentaren unten erwähnt habe, halte ich es für vernünftig, die "halbe maximale Breite" zu betrachten, aber Sie können auch die Gesamtzeit betrachten, in der der Peak über Ihrer Schwelle steht. Im Idealfall sollten Sie verschiedene Maße für die Peakbreite verwenden und angeben, wie konsistent Ihre Ergebnisse bei diesen Auswahlmöglichkeiten waren.

Unabhängig von der Metrik Ihrer Wahl können Sie mithilfe von Bootstrapping ein Konfidenzintervall für jeden Peak in jedem Trace berechnen.

f = @(xy) mylowess(xy,x,0.15);
N_boot = 1000;
yboot_1 = bootstrp(N_boot,f,[x,Y1])';
yboot_2 = bootstrp(N_boot,f,[x,Y2])';

Dieser Code erstellt 1000 Bootstrap-Anpassungen für die blauen und roten Spuren in den obigen Darstellungen. Ein Detail, das ich beschönigen werde, ist die Wahl des Glättungsfaktors 0,15 - Sie können diesen Parameter so wählen, dass er den Kreuzvalidierungsfehler minimiert (siehe den von mir geposteten Link). Jetzt müssen Sie nur noch eine Funktion schreiben, die die Peaks isoliert und ihre Breite schätzt:

function [t_peaks,heights,widths] = getPeaks(t,Y)
%% Computes a list of times, heights, and widths, for each peak in a time series Y
%% (column vector) with associated time points t (column vector).

% The implementation of this function will be problem-specific...

Anschließend führen Sie diesen Code auf den 1000 Kurven für jeden Datensatz aus und berechnen die 2,5- und 97,5-Perzentile für die Breite jedes Peaks. Ich werde dies anhand der Y1-Zeitreihe veranschaulichen - Sie würden dasselbe für die Y2-Zeitreihe oder einen anderen interessierenden Datensatz tun.

N_peaks = 2;  % two peaks in example data
t_peaks = nan(N_boot,N_peaks);
heights = nan(N_boot,N_peaks);
widths = nan(N_boot,N_peaks);
for aa = 1:N_boot
  [t_peaks(aa,:),heights(aa,:),widths(aa,:)] = getPeaks(x,yboot_1(:,aa));
end

quantile(widths(:,1),[0.025 0.975]) % confidence interval for the width of first peak
quantile(widths(:,2),[0.025 0.975]) % same for second peak width

Wenn Sie möchten, können Sie Hypothesentests durchführen, anstatt Konfidenzintervalle zu berechnen. Beachten Sie, dass der obige Code vereinfacht ist - es wird davon ausgegangen, dass jede Bootstrap-Lowess-Kurve 2 Peaks aufweist. Diese Annahme gilt möglicherweise nicht immer. Seien Sie also vorsichtig. Ich versuche nur, den Ansatz zu veranschaulichen, den ich verfolgen würde.

Hinweis: Die Funktion "mylowess" ist in dem oben angegebenen Link angegeben. So sieht es aus ...

function ys=mylowess(xy,xs,span)
%MYLOWESS Lowess smoothing, preserving x values
%   YS=MYLOWESS(XY,XS) returns the smoothed version of the x/y data in the
%   two-column matrix XY, but evaluates the smooth at XS and returns the
%   smoothed values in YS.  Any values outside the range of XY are taken to
%   be equal to the closest values.

if nargin<3 || isempty(span)
  span = .3;
end

% Sort and get smoothed version of xy data
xy = sortrows(xy);
x1 = xy(:,1);
y1 = xy(:,2);
ys1 = smooth(x1,y1,span,'loess');

% Remove repeats so we can interpolate
t = diff(x1)==0;
x1(t)=[]; ys1(t) = [];

% Interpolate to evaluate this at the xs values
ys = interp1(x1,ys1,xs,'linear',NaN);

% Some of the original points may have x values outside the range of the
% resampled data.  Those are now NaN because we could not interpolate them.
% Replace NaN by the closest smoothed value.  This amounts to extending the
% smooth curve using a horizontal line.
if any(isnan(ys))
  ys(xs<x1(1)) = ys1(1);
  ys(xs>x1(end)) = ys1(end);
end