[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