Two-sample permutation tests for earthquake magnitudes before and after a damaging quake

Background and data. An earthquake of magnitude 5.17 stuck near Yountville, California in the early morning hours of September 9, 2000, injuring about 25 people and doing about $50 million damage. The magnitudes of 24 smaller earthquakes nearby are given in the vector Mag below. The first 12 occurred within the few hours before the big quake and the last 12 occurred within a few hours after. The issue is whether the Before and After quakes differ significantly from each other in magnitude.  Mag=c(1.00, 1.14, 1.05, 1.79, 1.64, 2.76, 1.42, 1.30, 1.14, 0.96, 1.62, 1.40, 2.06, 1.64, 1.87, 1.41, 1.18, 1.59, 0.95, 1.54, 2.00, 1.16, 0.93, 1.90)  The two boxplots below show magnitudes of these two groups of earthquakes. Within the boxes, vertical lines show the sample medians and dots show their sample means. The outlier at 2.76 among those before the main quake might be called a foreshock. The quakes that occurred after the main one include some mild aftershocks. Traditional tests may be inappropriate. We wish to test$H_0: \mu_B = \mu_A$against$H_a: \mu_B \ne \mu_A.$Because the general population of California earthquake magnitudes is known to be distinctly right-skewed and our samples are small, one may be reluctant to use a two-sample t test to determine whether the two population means differ. (This test, which assumes normal data, gives the p-value 0.66, indicating no significant difference.) Ordinarily, the substitute would be a Mann-Whitney-Wilcoxon nonparmatric rank sum test, which tests whether population medians are equal. However, tied magnitudes (with values 1.14 and 1.64) would require adjustments in order test whether the population medians are equal (it shows an approximate p-value of 0.42). (And there are other quibbles with the applicability of the MWW test.) By contrast, there are no problematic assumptions for permutation tests. Our purpose here is use appropriate$permutation\: tests$to see whether the Before and After populations differ, and to compare the results of permutation tests with those for the t and MWW tests. 1 Answer Under the null hypothesis$H_0: \mu_B = \mu_A,$observations Before and After the big earthquake are effectively from the same population. Thus the observations in Mag may be permuted. The first question is what 'metric' we will use to measure the distance between the two samples. We begin with the difference$D$between the two sample means. There are${24, 12} = 2,704,156$possible permutations. To find the exact permutation distribution of the differences$D$in sample means, we would find the difference in sample means between the first and last 12 observations of each permuted sample. Then to test$H_0$against$H_a: \mu_B \ne \mu_A,$we would judge whether the observed$D$from the original (unpermutated) data is surprisingly large. But it seems there are too many permutations to construct and consider individually and sometimes two different permutations give the same$D$, so tallying results would be intricate. Instead, we look at 100,000 random permutations. For each permutation, we find the difference in sample means between the first 12 and the last 12 observations. Taken together, these results give a reasonably accurate imitation of the permutation distribution. We can find the p-value of the test from the simulated permutation distribution. The observed absolute difference is -0.08417, which nowhere near the tails of the permutation distribution. The p-value of this permutation test is . 0.6546, almost the same as for the traditional t test mentioned in the Question. The R code for approximating the permutation distribution and finding the p-value is shown below, followed by a graphic display of the results.  Mag = c(1.00, 1.14, 1.05, 1.79, 1.64, 2.76, 1.42, 1.30, 1.14, 0.96, 1.62, 1.40, 2.06, 1.64, 1.87, 1.41, 1.18, 1.59, 0.95, 1.54, 2.00, 1.16, 0.93, 1.90) BefAft = rep(1:2, each=12) mean.bef = mean(Mag[BefAft==1]); mean.aft = mean(Mag[BefAft==2]) d.obs = mean.bef - mean.aft # t.obs = t.test(Mag ~ BefAft, var.eq=T)$statistic
m = 10000;  d.perm = t.perm  = numeric(m)
for(i in 1:m) {
perm = sample(Mag, 24)
d.perm[i] = mean(perm[BefAft==1]) - mean(perm[BefAft==2])
# t.perm[i] = t.test(perm ~ BefAft, var.eq=T)$statistic } mean(abs(d.perm) > abs(d.obs)) # mean(abs(t.perm) > abs(t.obs)) mean(abs(d.perm) > abs(d.obs)) ## 0.6546 # mean(abs(t.perm) > abs(t.obs)) ## 0.6546  The p-value for the metric$D$is 0.65, which is roughly the same as we obtained in the Question using the traditional t test. This illustrates the 'robustness' (insensitivity to departure from assumptions) of the traditional t test. Another valid metric is the$T$statistic (pooled t test). By removing the single # symbols in the program above, we see that the p-value for the$T$metric is the same as for$D$. The histograms below show the simulated permutation distributions for the$D$and$T$metrics, respectively. The density of Student's t distribution with 22 degrees of freedom is shown in the right-hand panel. In both panels, the total area outside of the vertical lines is the p-value. The observed value of the metric is at the dashed line. A third useful metric for a permutation test might be to look at the differences between the two sample medians. The corresponding permutation distribution and p-value can be found by making minor changes in the program above, replacing mean by median at appropriate places, mostly inside the loop. The p-value is about 0.40. This would be a test for the difference in population medians. It is not surprising that the p-value is smaller here because the figure in the Question shows that the two sample medians are more widely separated than the two sample means. (Also note that this p-value is similar to that from the Mann-Whitney-Wilcoxon rank sum test, which also tests for differences in population medians.) • How would you go about determining the power of such a test? I suppose it would require more simulations under a variety of distributional assumptions? – heropup Dec 27 '15 at 2:36 • @heropup. Usually yes. More simulations are needed to get power. Notice my code for$T\$ metric uses built-in R function t.test. This is quick for beginners to grasp. But this is also a lazy and slow way to simulate because t.test has lots of built-in overhead. Important to write fast, clean code if you're doing a more intricate power study. – BruceET Dec 27 '15 at 6:45
• @heropup: Less vaguely, you might want to look at page 4 of https://inst.eecs.berkeley.edu/~cs174/sp08/lecs/lec11/lec11.pdf – BruceET Dec 27 '15 at 17:15