Absent a probability model, you can consider a permutation test. This requires two things:
A measure of overlap.
A way to re-use the data to create alternative hypothetical datasets for comparison.
A test statistic
One possible measure of overlap is obtained by considering each dataset to represent a union of intervals--that is, a subset of the number line--and to compute the total amount of their set-theoretic intersection. Let's start our programming by implementing this. It will be convenient to represent the data as $2\times n$ arrays, with the first row giving the start and the second row giving the end of each interval, thus:
d.1 <- matrix(c(125,500, 900,1300, 2220,2500), nrow=2)
d.2 <- matrix(c(600,800, 1200,1400, 3020,3500), nrow=2)
It should not matter that the intervals are sorted within each matrix, because conceptually these are just collections of the intervals that were observed. (I won't bother to program any checks to enforce the requirement that intervals within a dataset have no mutual overlaps: the user will be responsible for assuring this. The presence of overlaps will not break any of the code that is later written.)
The amount of overlap can be computed as
overlap <- function(x,y) {
o <- function(i,j) max(0, min(c(i[2]-i[1],i[2]-j[1],j[2]-i[1],j[2]-j[1])))
sum(apply(x, 2, function(u) apply(y, 2, function(v) o(u,v))))
}
stat <- overlap(d.1, d.2)
This is a little crude--it compares all possible ordered pairs of intervals--but it's fast enough to illustrate the ideas. If the permutation test turns out to take too long, this is the function to optimize (by sorting the intervals within each dataset, its timing can be reduced to a linear function of the number of intervals rather than a quadratic function).
I have saved the overlap of the actual data in stat
for future reference. For these data, it equals $100$, which is the length of $[1200, 1300]$, the intersection of the two datasets.
The permutation test
One way to modify the data, while preserving whatever structure is evident--I don't know whether this is appropriate for such an experiment but I have nothing else to go on--is to view the interval lengths as independent realizations of one random variable $X$ and the gaps between the intervals as independent realizations of another variable $Y$. The null hypothesis (that the two datasets do not really differ) is that each is obtained by a suitable number of draws from $(X,Y)$. The permutation test draws without replacement from the combined realizations $(x_i,y_i)$ obtained from pooling the data.
For instance, the within-interval gaps, $x_i$, in the first dataset are $375=500-125$, $400=1300-900$, and $280=2500-2220$. The between-interval gaps, $y_i$, in the first dataset are $125=125-0$, $400=900-500$, and $920=2220-1300$. This function computes them:
gaps <- function(x) {
y <- as.vector(x[, order(x[1,])])
matrix(diff(c(0, y)), nrow=2, dimnames=list(c("Between", "Within"), NULL))
}
(Notice how it protects itself against unordered input by sorting by the starts of each interval.)
Later, after permuting the between-interval gaps among themselves and the within-interval gaps among themselves, we will need to reassemble the gaps into a facsimile of a dataset:
assemble <- function(x) matrix(cumsum(x), nrow=2)
Using these preliminaries, it is now easy to construct a permutation test or even for bootstrapping the distribution (where sampling is done with replacement): an optional parameter, mysteriously indicated by ...
here, will determine what happens. The parameter n
requests the number of iterations. The two datasets are matrices x
and y
:
boot <- function(x, y, n=1, ...) {
m.x <- dim(x)[2]; m.y <- dim(y)[2]
g <- cbind(gaps(x), gaps(y))
trial <- function() {
z <- t(apply(g, 1, sample, ...))
overlap(assemble(z[, 1:m.x]), assemble(z[, 1:m.y + m.x]))
}
replicate(n, trial())
}
The optional parameter ...
is used by sample
. To sample with replacement, set that parameter to replace=TRUE
.
Example
Let's perform both a bootstrap and a permutation test, since we're capable of doing both with the same code. To make the results directly comparable, they each use exactly the same sequence of random numbers:
set.seed(17); sim.1 <- boot(d.1, d.2, 1000) # Permutation test
set.seed(17); sim.2 <- boot(d.1, d.2, 1000, replace=TRUE) # Bootstrap
The p-value depends on the alternative hypothesis, but in any event will be found in terms of the number of trial results having an overlap less than or equal to what was observed:
weight <- function(x) ifelse(x < 0, 1, ifelse(x == 0, 1/2, 0))
p.1 <- mean(weight(sim.1)); p.2 <- mean(weight(sim.2))
For this calculation, these values turn out to be $0.047$ and $0.0465$, respectively. That's just small enough to give a whiff of "significance" to a one-sided test of whether there is too little overlap.
A single number, like p.1
or p.2
, does not fully reveal the simulation results. Let's plot them:
breaks <- hist(c(sim.1,sim.2), plot=FALSE)$breaks # Get common bins for the plots
par(mfrow=c(1,2))
hist(sim.1, probability=TRUE, breaks=breaks, main="Permutation distribution")
abline(v = stat, lwd=2, col="Red")
hist(sim.2, probability=TRUE, breaks=breaks, main="Bootstrap distribution")
abline(v = stat, lwd=2, col="Red")
Comments
To use a different measure of overlap, recode overlap
and proceed as above.
Other probability models can be accommodated within this framework, but require more coding. For instance, we might postulate (perhaps based on scientific grounds) that the width of an interval might be correlated with the gaps to its neighboring intervals. If that is the case, we would not draw the between-interval widths and within-interval widths independently of each other. But now we're starting to fit an actual probability model to the data and we might want to start thinking about a parametric test or a parametric bootstrap.