[R] first post and bootstarpping problems

Tim Hesterberg timhesterberg at gmail.com
Fri Oct 8 15:30:18 CEST 2010


>I use R for a year now and am dealing with geometric morphometrics of deer skulls. Yes, I am a biologist and my math skills are just beginning to brush up. To cut to the chase...
>
>I have two groups in my data (males and females) and my data is in a simple vector form. Now I need a bootstrap test for this value
>
>szc1 <- ((mean(maleCent)-mean(femaCent))^ 2)/(var(maleCent)+var(femaCent))
>
>which concerns two aforementioned groups. I have 39 males and 11 females totaling to 50 individuals. Now I don`t know how to assign this to a bootstrap boot() function. Any ideas?

You should use a two-sample permutation test rather than doing a bootstrap.

The basic idea is:
compute the statistic for the original data
compute the statistic for many permutations of the data
compute upper and lower P-values using 
  (1+number of permutation statistics that exceed the original) / (r+1)
two-sided P-value = 2 * smaller of one-sided P-values

Here is example code:

maleCent <- rnorm(39) # example data
femaCent <- rnorm(11)
xBoth <- c(maleCent, femaCent)
statistic <- function(x){
  x1 <- x[1:39]
  x2 <- x[40:50]
  (mean(x1)-mean(x2))^2 / (var(x1) + var(x2))
}
theta <- statistic(xBoth)

r <- 10^5-1
permutationDistribution <- numeric(r)
for(i in 1:r) {
  permutationDistribution[i] <- statistic(sample(xBoth, size=50, replace=TRUE))
}
  
pValueLower <- (1 + sum(permutationDistribution <= theta)) / (r+1)
pValueUpper <- (1 + sum(permutationDistribution >= theta)) / (r+1)
pValueTwosided <- 2 * min(pValueLower, pValueUpper)

Tim Hesterberg



More information about the R-help mailing list