Simple randomisation
As an illustration of both randomisation and simple usage of the permute package we consider a small data set of mandible length measurements on specimens of the golden jackal (Canis aureus) from the British Museum of Natural History, London, UK. These data were collected as part of a study comparing prehistoric and modern canids (Higham, Kijngam, and Manly 1980), and were analysed by Manly (2007). There are ten measurements of mandible length on both male and female specimens. The data are available in the jackal data frame supplied with permute.
## Length Sex
## 1 120 Male
## 2 107 Male
## 3 110 Male
## 4 116 Male
## 5 114 Male
## 6 111 Male
## 7 113 Male
## 8 117 Male
## 9 114 Male
## 10 112 Male
## 11 110 Female
## 12 111 Female
## 13 107 Female
## 14 108 Female
## 15 110 Female
## 16 105 Female
## 17 107 Female
## 18 106 Female
## 19 111 Female
## 20 111 FemaleThe interest is whether there is a difference in the mean mandible length between male and female golden jackals. The null hypothesis is that there is zero difference in mandible length between the two sexes or that females have larger mandibles. The alternative hypothesis is that males have larger mandibles. The usual statistical test of this hypothesis is a one-sided \(t\) test, which can be applied using t.test()
##
## Two Sample t-test
##
## data: Length by Sex
## t = 3.4843, df = 18, p-value = 0.001324
## alternative hypothesis: true difference in means between group Male and group Female is greater than 0
## 95 percent confidence interval:
## 2.411156 Inf
## sample estimates:
## mean in group Male mean in group Female
## 113.4 108.6The observed \(t\) is 3.484 with 18 df. The probability of observing a value this large or larger if the null hypothesis were true is 0.0013. Several assumptions have been made in deriving this \(p\)-value, namely
- random sampling of individuals from the populations of interest,
- equal population standard deviations for males and females, and
- that the mandible lengths are normally distributed within the sexes.
Assumption 1 is unlikely to be valid for museum specimens such as these, that have been collected in some unknown manner. Assumption 2 may be valid, Fisher’s \(F\)-test and a Fligner-Killeen test both suggest that the standard deviations of the two populations do not differ significantly
var.test(Length ~ Sex, data = jackal)
fligner.test(Length ~ Sex, data = jackal)This assumption may be relaxed using var.equal = FALSE (the default) in the call to t.test(), to employ Welch’s modification for unequal variances. Assumption 3 may be valid, but with such a small sample we are unable to reliably test this.
A randomisation test of the same hypothesis can be performed by randomly allocating ten of the mandible lengths to the male group and the remaining lengths to the female group. This randomisation is justified under the null hypothesis because the observed difference in mean mandible length between the two sexes is just a typical value for the difference in a sample if there were no difference in the population. An appropriate test statistic needs to be selected. We could use the \(t\) statistic as derived in the \(t\)-test. Alternatively, we could base our randomisation test on the difference of means \(D_i\) (male - female).
The main function in permute for providing random permutations is shuffle(). We can write our own randomisation test for the jackal data by first creating a function to compute the difference of means for two groups
which can be used in a simple for() loop to generate the null distribution for the difference of means. First, we allocate some storage to hold the null difference of means; here we use 4999 random permutations so allocate a vector of length 5000. Then we iterate, randomly generating an ordering of the Sex vector and computing the difference of means for that permutation.
Djackal <- numeric(length = 5000)
N <- nrow(jackal)
set.seed(42)
for(i in seq_len(length(Djackal) - 1)) {
perm <- shuffle(N)
Djackal[i] <- with(jackal, meanDif(Length, Sex[perm]))
}
Djackal[5000] <- with(jackal, meanDif(Length, Sex))The observed difference of means was added to the null distribution, because under the null hypothesis the observed allocation of mandible lengths to male and female jackals is just one of the possible random allocations.
The null distribution of \(D_i\) can be visualised using a histogram, as shown in Figure~@ref{draw_hist_jackal}. The observed difference of means (4.8) is indicated by the red tick mark.
hist(Djackal, main = "",
xlab = expression("Mean difference (Male - Female) in mm"))
rug(Djackal[5000], col = "red", lwd = 2)The number of values in the randomisation distribution equal to or larger than the observed difference is
## [1] 12giving a permutational \(p\)-value of
## [1] 0.0024which is comparable with that determined from the frequentist \(t\)-test, and indicates strong evidence against the null hypothesis of no difference.
(#fig:draw_hist_jackal)Distribution of the difference of mean mandible length in random allocations, ten to each sex.
In total there \(^{20}C_{10} = 184,756\) possible allocations of the 20 observations to two groups of ten
## [1] 184756so we have only evaluated a small proportion of these in the randomisation test.
The main workhorse function we used above was shuffle(). In this example, we could have used the base R function sample() to generate the randomised indices perm that were used to permute the Sex factor. Where shuffle() comes into it’s own is for generating permutation indices from restricted permutation designs.