Cómo probar si la varianza de dos distribuciones es diferente si las distribuciones no son normales

Estoy estudiando dos poblaciones geográficamente aisladas de la misma especie. Al inspeccionar las distribuciones, veo que ambas son bimodales (hay cierta estacionalidad en su ocurrencia), pero los picos en una población son mucho más altos y mucho más estrechos (es decir, la variación de los picos locales es menor).

¿Qué tipo de prueba estadística sería apropiada para determinar si estas diferencias son significativas?

Para aclarar, mi eje y es el número de individuos identificados en una trampa en un día en particular, y el eje x es el día juliano.

¿Son estas distribuciones de algo a lo largo del tiempo? Cuenta, tal vez? (Si es así, es posible que necesite algo bastante diferente de las discusiones aquí hasta ahora)

Lo que describe no suena como si estuviera muy bien recogido como una diferencia en la variación de las distribuciones.

Parece que está describiendo algo vagamente como esto (ignore los números en los ejes, es solo para dar una idea del tipo general de patrón que parece estar describiendo):

Si eso es correcto, entonces considere:

Si bien el ancho de cada pico alrededor de los centros locales es más estrecho para la curva azul, la variación de las distribuciones roja y azul en general apenas difiere.

Si identifica los modos y antimodos de antemano, podría medir la variabilidad local.

En primer lugar, creo que debería ver las distribuciones estacionales por separado, ya que es probable que la distribución bimodal sea el resultado de dos procesos bastante separados. Las dos distribuciones pueden controlarse mediante diferentes mecanismos, de modo que, por ejemplo, las distribuciones de invierno pueden ser más sensibles al clima anual. Si desea ver las diferencias de población y las razones de esto, creo que es más útil estudiar las distribuciones estacionales por separado.

En cuanto a una prueba, puede probar la prueba de Levine (básicamente una prueba de homocedasticidad), que se utiliza para comparar las variaciones entre los grupos. La prueba de Bartlett es una alternativa, pero se supone que la prueba de Levene es más robusta a la no normalidad (especialmente cuando se usa la mediana para la prueba). En R se encuentran las pruebas de Levene y Bartlett library(car).

Estoy de acuerdo con lo que otros han dicho, a saber, que "varianza" es probablemente la palabra incorrecta para usar (ya que la función que está considerando no es una distribución de probabilidad sino una serie de tiempo).

Creo que es posible que desee abordar este problema desde una perspectiva diferente: simplemente ajuste las dos series de tiempo con curvas BAJAS. Puede calcular intervalos de confianza del 95% y comentar cualitativamente sus formas. No estoy seguro de que necesites hacer algo más elegante que esto.

He escrito un código MATLAB a continuación para ilustrar lo que estoy diciendo. Tengo un poco de prisa pero puedo dar aclaraciones pronto. Gran parte de lo que hice se puede tomar directamente desde aquí: 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')

Es posible que desee normalizar las dos series temporales para comparar sus tendencias relativas en lugar de sus niveles absolutos.

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

Ahora haga ajustes BAJOS ...

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

Finalmente, puede crear bandas de confianza del 95% de la siguiente manera:

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

Ahora puede interpretar la cifra final como lo desee, y tiene los ajustes BAJOS para respaldar su hipótesis de que los picos en la curva roja son en realidad más amplios que la curva azul. Si tiene una mejor idea de cuál es la función, podría hacer una regresión no lineal.

Editar: Basado en algunos comentarios útiles a continuación, estoy agregando algunos detalles más sobre la estimación explícita de los anchos de pico. Primero, debe llegar a una definición de lo que está considerando que es un "pico". Quizás cualquier golpe que se eleve por encima de algún umbral (algo así como 0.05 en las parcelas que hice arriba). El principio básico es que debe encontrar una manera de separar los picos "reales" o "notables" del ruido.

Luego, para cada pico, puede medir su ancho de varias maneras. Como mencioné en los comentarios a continuación, creo que es razonable observar el "ancho medio máximo", pero también podría ver el tiempo total en que el pico se encuentra por encima de su umbral. Idealmente, debe usar varias medidas diferentes de ancho de pico e informar qué tan consistentes se les dieron estas opciones a sus resultados.

Cualquiera sea su métrica (s) de elección, puede usar bootstrapping para calcular un intervalo de confianza para cada pico en cada traza.

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

Este código crea 1000 ajustes de arranque para los trazos azules y rojos en las parcelas anteriores. Un detalle que pasaré por alto es la elección del factor de suavizado 0.15: puede elegir este parámetro para minimizar el error de validación cruzada (consulte el enlace que publiqué). Ahora todo lo que tiene que hacer es escribir una función que aísle los picos y calcule su ancho:

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

Luego ejecuta este código en las 1000 curvas para cada conjunto de datos y calcula los percentiles 2.5 y 97.5 para el ancho de cada pico. Ilustraré esto en la serie temporal Y1: usted haría lo mismo para la serie temporal Y2 o cualquier otro conjunto de datos de interés.

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

Si lo desea, puede realizar pruebas de hipótesis en lugar de calcular intervalos de confianza. Tenga en cuenta que el código anterior es simplista: supone que cada curva de lowess bootstrapped tendrá 2 picos. Es posible que esta suposición no siempre sea válida, así que tenga cuidado. Solo estoy tratando de ilustrar el enfoque que adoptaría.

Nota: la función "mylowess" aparece en el enlace que publiqué anteriormente. Esto es lo que parece...

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