10
$\begingroup$

I'm studying two geographically-isolated populations of the same species. Inspecting the distributions, I see that both are bimodal (there's some seasonality to their occurrence), but the peaks in one population are much higher and much narrower (i.e., the variance of the local peaks is smaller).

What sort of statistical test would be appropriate to determine whether these differences are significant?

To clarify, my y-axis is the number of individuals identified in a trap on a particular day, and the x-axis is Julian day.

Improve this question
$\endgroup$
10
  • $\begingroup$ You can try to do some outlier detection. en.wikipedia.org/wiki/Outlier. $\endgroup$
    – Yevgeniy.Chernobrivets
    Commented Nov 18, 2013 at 18:56
  • $\begingroup$ Are you able to write down a statistical model? Also, there are many different ways to specify "the variances are not equal" and "the variances are equal" and your conclusion may well depend on which particular choices you make, especially if it is a subtle difference. So it is better to use a model chosen by you, rather than one chosen by someone with no context. $\endgroup$
    – probabilityislogic
    Commented Nov 19, 2013 at 8:43
  • 1
    $\begingroup$ It's both! You have a time series of counts. $\endgroup$
    – whuber ♦
    Commented Dec 3, 2013 at 17:53
  • 1
    $\begingroup$ It would help immensely to have a model, or at least some suggestive theory, that attempts to explain why some peaks would be narrower and others wider. Because you are interested in the widths of these peaks, you must have at least a conceptual model, if not a quantitative one. What mechanisms do you suppose produce such peaks and govern their widths? Do you have independent information that suggests when the peaks ought to occur? (This reduces uncertainty in peak identification.) Do peaks occur contemporaneously or at different times? $\endgroup$
    – whuber ♦
    Commented Dec 4, 2013 at 16:24
  • 2
    $\begingroup$ @whuber, peaks of the two populations are nearly contemporaneous. One is in temperate latitudes, and one is in tropical latitudes. Our hypothesis is that the tropical population has a narrower ecological niche than the temperate population (i.e., a wider host of predators and pathogens pressures the population into a narrow emergence time). Does that help? $\endgroup$
    – Atticus29
    Commented Dec 5, 2013 at 19:34

3 Answers 3

Reset to default
3
$\begingroup$

Are these distributions of something over time? Counts, perhaps? (If so then you might need something quite different from the discussions here so far)

What you describe doesn't sound like it would be very well picked up as a difference in variance of the distributions.

It sounds like you're describing something vaguely like this (ignore the numbers on the axes, it's just to give a sense of the general kind of pattern you seem to be describing):

bimodal peaks

If that's right, then consider:

While the width of each peak about the local centers is narrower for the blue curve, the variance of the red and blue distributions overall hardly differs.

If you identify the modes and antimodes beforehand, you could then measure the local variability.

Improve this answer
$\endgroup$
6
  • $\begingroup$ this is exactly my question. Thanks! So, would restricting my x-axis range to encompass, say, only the first peak, and then conducting... an F-test??...be the best approach? $\endgroup$
    – Atticus29
    Commented Nov 19, 2013 at 17:38
  • $\begingroup$ You probably wouldn't want to do specifically do an F test for variances, I think, for a couple of reasons (If you do test variance in that way, @fileunderwater has mentioned some alternatives to the F test). But before we get that far, could you address the two questions at the top of my post? Is this distribution of counts over time? $\endgroup$
    – Glen_b
    Commented Nov 19, 2013 at 19:04
  • $\begingroup$ they are (see edits to the question). $\endgroup$
    – Atticus29
    Commented Nov 19, 2013 at 20:30
  • $\begingroup$ with the new information and per my comment to fileunderwater's answer above, do you have any suggestions? $\endgroup$
    – Atticus29
    Commented Nov 20, 2013 at 19:59
  • 1
    $\begingroup$ There appears to be considerable confusion in the question and these comments about what a "variance" is. In Glen_b's examples, the blue data have greater variances than the red data around the two apparent peaks (near x=10 and x=17), because the blue data swing more between low and high values (which are plotted on the vertical axis, not the horizontal, which apparently represents time). $\endgroup$
    – whuber ♦
    Commented Dec 3, 2013 at 15:36
3
$\begingroup$

First of all, I think that you should look at the seasonal distributions separately, since the bimodal distribution is likely to be the outcome of two fairly separate processes. The two distributions might be controlled by different mechanisms, so that e.g. winter distributions could be more sensitive to yearly climate. If you want to look at population differences and reasons for these I think it is therefore more useful to study the seasonal distributions separately.

As for a test, you could try Levine's test (basically a test of homoscedasticity), which is used to compare variances between groups. Bartlett's test is an alternative, but Levene's test is supposed to be more robust to non-normality (especially when using the median for testing). In R the Levene's and Bartlett's tests are found in library(car).

Improve this answer
$\endgroup$
5
  • $\begingroup$ I'm looking into Levene's test in R (I found it in the library "car"). It looks like it only takes a linear model object as an argument. This doesn't really make sense in my case, because I just want to compare the variance of two distributions (not analyze them with linear models and validate those assumptions). Any advice? $\endgroup$
    – Atticus29
    Commented Nov 20, 2013 at 19:58
  • 1
    $\begingroup$ @Atticus29 Yes, its in car - my mistake. However, it is not based on a strict linear model - you can use leveneTest(y ~ as.factor(group), data= datafile) for a test of difference in variance between groups, and if you use the option `center="median" it is more robust to non-normality. Strictly, I think its called Brown–Forsythe test if based on the median. $\endgroup$
    – fileunderwater
    Commented Nov 20, 2013 at 22:19
  • $\begingroup$ Ok, so, dumb question, but I have two columns of data which are counts of individuals from a particular species caught in traps. These two columns represent counts of the same species on the same days from different locations. I'm not sure how to group them using location without losing date information using the above format, if that makes sense.... $\endgroup$
    – Atticus29
    Commented Nov 21, 2013 at 20:13
  • $\begingroup$ @ Atticus Can you include some sample data to your question (including all columns and classification variables)? This would help to clarify some of the confusion on exactly what type of data you have (see e.g. comments by @whuber). My feeling was that you had pooled all species records from two seasons, but now when I re-read your Q this seems not to be the case, and I'm not sure that my solution is suitable. Do you only have traps at two locations and have daily(?) counts in these over time (for a single year)? $\endgroup$
    – fileunderwater
    Commented Dec 6, 2013 at 11:53
  • $\begingroup$ [cnd]... What do you think the late-season peak is caused by; a second generation within the same year (what taxa are you studying?) or two different phenotypes? @Atticus29 $\endgroup$
    – fileunderwater
    Commented Dec 6, 2013 at 11:54
3
$\begingroup$

I agree with what others have said -- namely that "variance" is probably the wrong word to use (seeing as the function you are considering isn't a probability distribution but a time-series).

I think you may want to approach this problem from a different perspective -- just fit the two time series with LOWESS curves. You can calculate 95% confidence intervals and qualitatively comment on their shapes. I'm not sure you need to do anything more fancy than this.

I've written some MATLAB code below to illustrate what I'm saying. I'm in a bit of a rush but can provide clarifications soon. Much of what I did can be taken directly from here: 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')

You may want to normalize the two time-series to compare their relative trends rather than their absolute levels.

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

Now make LOWESS fits...

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

enter image description here

Finally, you can create 95% confidence bands as follows:

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

Now you can interpret the final figure as you wish, and you have the LOWESS fits to back up your hypothesis that the peaks in the red curve are actually broader than the blue curve. If you have a better idea of what the function is you could do non-linear regression instead.

Edit: Based on some helpful comments below, I am adding some more details about estimating peak widths explicitly. First, you need to come up with some definition for what you are considering a "peak" to be in the first place. Perhaps any bump that rises above some threshold (something like 0.05 in the plots I made above). The basic principle is that you should find a way from separating "real" or "notable" peaks from noise.

Then, for each peak, you can measure its width in a couple of ways. As I mentioned in the comments below, I think it is reasonable to look at the "half-max-width" but you could also look at the total time the peak stands above your threshold. Ideally, you should use several different measures of peak width and report how consistent your results were given these choices.

Whatever your metric(s) of choice, you can use bootstrapping to calculate a confidence interval for each peak in each trace.

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

This code creates 1000 bootstrapped fits for the blue and red traces in the plots above. One detail that I will gloss over is the choice of the smoothing factor 0.15 -- you can choose this parameter such that it minimizes cross validation error (see the link I posted). Now all you have to do is write a function that isolates the peaks and estimates their width:

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

Then you run this code on the 1000 curves for each dataset and calculate the 2.5th and 97.5th percentiles for the width of each peak. I'll illustrate this on the Y1 time series - you would do the same for the the Y2 time series or any other data set of interest.

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

If you desire, you can perform hypothesis tests rather than calculating confidence intervals. Note that the code above is simplistic - it assumes each bootstrapped lowess curve will have 2 peaks. This assumption may not always hold, so be careful. I'm just trying to illustrate the approach I would take.

Note: the "mylowess" function is given in the link I posted above. This is what it looks like...

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
Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ Welcome to our site and thank you for posting a clear, well-illustated answer. This looks like a good approach and a promising technique. It appears to fall short of answering the question, however: just how would you go about (a) identifying "peaks" and (b) formally testing their widths? $\endgroup$
    – whuber ♦
    Commented Dec 9, 2013 at 19:22
  • $\begingroup$ My inclination would be to show the above plots and provide interpretation: "The red and blue populations each show two peaks around t=25 and t=65. The red population, however, approaches these peaks at a slower rate (e.g. for the first peak, beginning around t=10 vs. t=15 for the blue population)..." The 95% confidence bands give the reader a sense of what bends in the curves are noise vs real effects. I think this should be sufficient to explain the original data set described for publication (if that is the end goal). $\endgroup$
    – digbyterrell
    Commented Dec 9, 2013 at 19:56
  • $\begingroup$ Many peer reviewers would point out that (a) these CIs are not CIs for the peak widths and (b) even if they were, direct comparison of CIs is not a legitimate statistical procedure with known Type I and Type II error rates. Whence the original question: how does one formally test the visually apparent differences? $\endgroup$
    – whuber ♦
    Commented Dec 9, 2013 at 20:01
  • $\begingroup$ If you really wanted to do some "formal" calculations... I suppose you could find all local min/max in the lowess fit (points where the first derivative is zero), then calculate the amplitude for each peak (you may have to ignore small-amplitude peaks), and finally calculate the "half-max-width" of each peak (the time between when the curve is halfway up to when the curve is halfway down). Then you could do a similar bootstrapping procedure as described in my answer above to see whether the red "half-max-width" is consistently larger. I can provide more details if there is interest. $\endgroup$
    – digbyterrell
    Commented Dec 9, 2013 at 20:03
  • $\begingroup$ Bootstrapping is appealing, but it is not at all clear how it should be conducted, given that no specific statistical model has been proposed in the question. Some kind of model appropriate for the data is essential because (at a minimum) these time series will likely exhibit strong serial correlation. Other details are almost as important: how does one determine which peaks are "small" and which are not? Should the peak widths be measured at half-height or at some other point? What degree of smoothing should be used for the lowess fit? (There is at least one arbitrary parameter to set.) $\endgroup$
    – whuber ♦
    Commented Dec 9, 2013 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.