Como testar se a variação de duas distribuições é diferente se as distribuições não são normais

Estou estudando duas populações geograficamente isoladas da mesma espécie. Inspecionando as distribuições, vejo que ambas são bimodais (há certa sazonalidade em sua ocorrência), mas os picos em uma população são muito mais altos e mais estreitos (ou seja, a variação dos picos locais é menor).

Que tipo de teste estatístico seria apropriado para determinar se essas diferenças são significativas?

Para esclarecer, meu eixo y é o número de indivíduos identificados em uma armadilha em um dia específico, e o eixo x é o dia juliano.

Essas distribuições são de alguma coisa ao longo do tempo? Conta, talvez? (Se sim, então você pode precisar de algo bem diferente das discussões aqui até agora)

O que você descreve não parece ser muito bem entendido como uma diferença na variação das distribuições.

Parece que você está descrevendo algo vagamente assim (ignore os números nos eixos, é apenas para dar uma idéia do tipo geral de padrão que você parece estar descrevendo):

Se estiver certo, considere:

Embora a largura de cada pico sobre os centros locais seja mais estreita para a curva azul, a variação geral das distribuições de vermelho e azul dificilmente difere.

Se você identificar os modos e antimodos de antemão, poderá medir a variabilidade local.

Antes de tudo, acho que você deve examinar as distribuições sazonais separadamente, já que a distribuição bimodal provavelmente será o resultado de dois processos bastante separados. As duas distribuições podem ser controladas por mecanismos diferentes, de modo que, por exemplo, as distribuições de inverno sejam mais sensíveis ao clima anual. Se você quiser examinar as diferenças populacionais e as razões para isso, acho que é mais útil estudar as distribuições sazonais separadamente.

Quanto a um teste, você pode tentar o teste de Levine (basicamente um teste de homocedasticidade), usado para comparar variações entre os grupos. O teste de Bartlett é uma alternativa, mas o teste de Levene deve ser mais robusto à não normalidade (especialmente ao usar a mediana para o teste). Em R, os testes de Levene e Bartlett são encontrados em library(car).

Concordo com o que os outros disseram - ou seja, que "variação" é provavelmente a palavra errada a ser usada (visto que a função que você está considerando não é uma distribuição de probabilidade, mas uma série temporal).

Eu acho que você pode querer abordar esse problema de uma perspectiva diferente - basta ajustar as duas séries temporais com curvas LOWESS. Você pode calcular intervalos de confiança de 95% e comentar qualitativamente suas formas. Não sei se você precisa fazer algo mais sofisticado do que isso.

Eu escrevi alguns códigos do MATLAB abaixo para ilustrar o que estou dizendo. Estou com pressa, mas posso fornecer esclarecimentos em breve. Muito do que fiz pode ser retirado diretamente daqui: 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')

Você pode normalizar as duas séries temporais para comparar suas tendências relativas em vez de seus níveis 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')

Agora faça LOWESS se encaixa ...

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

Por fim, você pode criar faixas de confiança de 95% da seguinte maneira:

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

Agora você pode interpretar a figura final como desejar e possui os ajustes LOWESS para apoiar sua hipótese de que os picos na curva vermelha são realmente mais amplos que a curva azul. Se você tiver uma idéia melhor de qual é a função, poderá fazer uma regressão não linear.

Edit: Com base em alguns comentários úteis abaixo, estou adicionando mais alguns detalhes sobre a estimativa de largura de pico explicitamente. Primeiro, você precisa criar uma definição para o que você considera um "pico" em primeiro lugar. Talvez qualquer aumento que ultrapasse algum limite (algo como 0,05 nas parcelas que fiz acima). O princípio básico é que você deve encontrar uma maneira de separar picos "reais" ou "notáveis" do ruído.

Então, para cada pico, você pode medir sua largura de duas maneiras. Como mencionei nos comentários abaixo, acho razoável examinar a "meia largura máxima", mas você também pode observar o tempo total em que o pico está acima do seu limite. Idealmente, você deve usar várias medidas diferentes de largura de pico e relatar a consistência de seus resultados com essas opções.

Quaisquer que sejam suas métricas de escolha, você pode usar a inicialização para calcular um intervalo de confiança para cada pico em cada rastreamento.

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

Esse código cria 1000 ajustes autoinicializados para os traços azul e vermelho nos gráficos acima. Um detalhe que abordarei é a escolha do fator de suavização 0,15 - você pode escolher este parâmetro para minimizar o erro de validação cruzada (consulte o link que eu publiquei). Agora tudo o que você precisa fazer é escrever uma função que isola os picos e estima sua largura:

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

Em seguida, você executa esse código nas 1000 curvas de cada conjunto de dados e calcula os percentis 2,5 e 97,5 para a largura de cada pico. Ilustrarei isso na série temporal Y1 - você faria o mesmo na série temporal Y2 ou em qualquer outro conjunto de dados de interesse.

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

Se desejar, você pode executar testes de hipóteses em vez de calcular intervalos de confiança. Observe que o código acima é simplista - ele assume que cada curva de bootess com bootstrap terá 2 picos. Essa suposição nem sempre é válida, portanto, tenha cuidado. Eu só estou tentando ilustrar a abordagem que eu adotaria.

Nota: a função "mylowess" é fornecida no link que eu postei acima. Isto é o 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