Cara menguji apakah varians dua distribusi berbeda jika distribusi tidak normal

Saya sedang mempelajari dua populasi yang secara geografis terisolasi dari spesies yang sama. Memeriksa distribusi, saya melihat bahwa keduanya bimodal (ada beberapa musim yang terjadi), tetapi puncak dalam satu populasi jauh lebih tinggi dan lebih sempit (yaitu, varian dari puncak lokal lebih kecil).

Apa jenis uji statistik yang sesuai untuk menentukan apakah perbedaan ini signifikan?

Untuk memperjelas, sumbu y saya adalah jumlah individu yang diidentifikasi dalam perangkap pada hari tertentu, dan sumbu x adalah hari Julian.

Apakah ini distribusi sesuatu dari waktu ke waktu? Hitungan, mungkin? (Jika demikian maka Anda mungkin perlu sesuatu yang sangat berbeda dari diskusi di sini sejauh ini)

Apa yang Anda gambarkan tidak terdengar seperti itu akan sangat baik diambil sebagai perbedaan dalam varian distribusi.

Sepertinya Anda menggambarkan sesuatu yang samar-samar seperti ini (abaikan angka pada sumbu, itu hanya untuk memberi gambaran tentang jenis pola umum yang tampaknya Anda uraikan):

Jika itu benar, maka pertimbangkan:

Sementara lebar setiap puncak tentang pusat-pusat lokal lebih sempit untuk kurva biru, varians dari distribusi merah dan biru secara keseluruhan hampir tidak berbeda.

Jika Anda mengidentifikasi mode dan antimode sebelumnya, Anda dapat mengukur variabilitas lokal.

Pertama-tama, saya pikir Anda harus melihat distribusi musiman secara terpisah, karena distribusi bimodal kemungkinan merupakan hasil dari dua proses yang cukup terpisah. Kedua distribusi mungkin dikendalikan oleh mekanisme yang berbeda, sehingga distribusi musim dingin misalnya bisa lebih sensitif terhadap iklim tahunan. Jika Anda ingin melihat perbedaan populasi dan alasannya, saya pikir akan lebih berguna untuk mempelajari distribusi musiman secara terpisah.

Sedangkan untuk tes, Anda bisa mencoba tes Levine (pada dasarnya tes homoscedasticity), yang digunakan untuk membandingkan perbedaan antar kelompok. Tes Bartlett adalah alternatif, tetapi tes Levene seharusnya lebih kuat untuk non-normalitas (terutama ketika menggunakan median untuk pengujian). Dalam R, tes Levene dan Bartlett ditemukan di library(car).

Saya setuju dengan apa yang dikatakan orang lain - yaitu bahwa "varians" mungkin kata yang salah untuk digunakan (mengingat fungsi yang Anda pertimbangkan bukan distribusi probabilitas tetapi serangkaian waktu).

Saya pikir Anda mungkin ingin mendekati masalah ini dari perspektif yang berbeda - cukup paskan dua seri waktu dengan kurva LOWESS. Anda dapat menghitung interval kepercayaan 95% dan mengomentari bentuknya secara kualitatif. Saya tidak yakin Anda perlu melakukan sesuatu yang lebih mewah dari ini.

Saya telah menulis beberapa kode MATLAB di bawah ini untuk menggambarkan apa yang saya katakan. Saya agak terburu-buru tetapi dapat segera memberikan klarifikasi. Banyak dari apa yang saya lakukan dapat diambil langsung dari sini: 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')

Anda mungkin ingin menormalkan dua seri waktu untuk membandingkan tren relatif mereka daripada tingkat absolutnya.

%% 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')

Sekarang buat LOWESS cocok ...

%% 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')

Terakhir, Anda dapat membuat band kepercayaan 95% sebagai berikut:

%% 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')

Sekarang Anda dapat menafsirkan angka akhir seperti yang Anda inginkan, dan Anda memiliki LOWESS yang cocok untuk mendukung hipotesis Anda bahwa puncak dalam kurva merah sebenarnya lebih luas daripada kurva biru. Jika Anda memiliki gagasan yang lebih baik tentang fungsi tersebut, Anda bisa melakukan regresi non-linear.

Sunting: Berdasarkan beberapa komentar bermanfaat di bawah, saya menambahkan beberapa detail lebih lanjut tentang memperkirakan lebar puncak secara eksplisit. Pertama, Anda perlu membuat beberapa definisi untuk apa yang Anda pertimbangkan sebagai "puncak" di tempat pertama. Mungkin ada benjolan yang naik di atas ambang tertentu (sekitar 0,05 di plot yang saya buat di atas). Prinsip dasarnya adalah Anda harus menemukan cara untuk memisahkan puncak "nyata" atau "penting" dari kebisingan.

Kemudian, untuk setiap puncak, Anda dapat mengukur lebarnya dalam beberapa cara. Seperti yang saya sebutkan di komentar di bawah, saya pikir masuk akal untuk melihat "setengah-max-lebar" tetapi Anda juga bisa melihat total waktu puncaknya berdiri di atas ambang batas Anda. Idealnya, Anda harus menggunakan beberapa ukuran lebar puncak yang berbeda dan melaporkan seberapa konsisten hasil Anda diberikan pilihan ini.

Apa pun metrik pilihan Anda, Anda dapat menggunakan bootstrap untuk menghitung interval kepercayaan untuk setiap puncak dalam setiap jejak.

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])';

Kode ini menciptakan 1000 bootstrapped cocok untuk jejak biru dan merah di plot di atas. Satu detail yang akan saya bahas adalah pilihan faktor smoothing 0,15 - Anda dapat memilih parameter ini sehingga meminimalkan kesalahan validasi silang (lihat tautan yang saya posting). Sekarang yang harus Anda lakukan adalah menulis fungsi yang mengisolasi puncak dan memperkirakan lebarnya:

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...

Kemudian Anda menjalankan kode ini pada 1000 kurva untuk setiap dataset dan menghitung persentil ke-2,5 dan ke-97,5 untuk lebar setiap puncak. Saya akan mengilustrasikan ini pada seri waktu Y1 - Anda akan melakukan hal yang sama untuk seri waktu Y2 atau set data menarik lainnya.

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

Jika diinginkan, Anda dapat melakukan tes hipotesis daripada menghitung interval kepercayaan. Perhatikan bahwa kode di atas adalah sederhana - ini mengasumsikan setiap kurva lowess bootstrap akan memiliki 2 puncak. Asumsi ini mungkin tidak selalu berlaku, jadi berhati-hatilah. Saya hanya mencoba menggambarkan pendekatan yang akan saya ambil.

Catatan: fungsi "mylowess" diberikan dalam tautan yang saya posting di atas. Ini seperti apa ...

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