0

My goal is to use bootstrapping (1000 reps) to calculate a null distribution, mean, and CI of r (Pearson's correlation coefficient) correlating trait (x) in 20 stimulated random pairs generated from my dataset of 600 unique individuals (ID). I have recently switched to R from SAS where I would use "proc surveyselect" to generate the dataset. Questions:

  1. What would be the most efficient way to generate these results (see my attempt below)?
  2. In my example, how would I use the set.seed command to replicate my results?

Simulated starting dataset with 600 individuals and the associated trait values:

ID <- seq(1, 600, by = 1)
x <- rnorm(600, m = 7, sd = 2)
X <- as.data.frame(cbind(ID, x))

I then generate my 1000 replicates of r and calculate the 95% CI:

for (i in 1:1000) { 
  X.sample <- X[ sample(1:nrow(X), 40, replace = FALSE), ] 
  X.sample.1 <- X.sample[1:20, ]
  X.sample.2 <- X.sample[21:40, ]
  Y <- as.data.frame(cbind(X.sample.1$ID, X.sample.1$x, X.sample.2$ID,  X.sample.2$x))
  cor.results <- cor.test(Y[,2], Y[,4], alternative = c("greater"), method = c("pearson"))
  Z[i] <- cor.results$estimate
}

error <- qt(0.975, df = (length(Z) - 1)) * (sd(Z))/sqrt(length(Z))
1
  • 1
    Just a couple of comments re code brevity... The ID column appears extraneous here, but if you do want it, ID <- 1:600 would do the trick. I can't see any reason to use a data.frame in this case, as your ID and x are the same data type (numeric). matrix operations are in general faster than data.frame operations, to my knowledge. See my solution below for some other time-savers.
    – jbaums
    Mar 15, 2012 at 13:57

2 Answers 2

1

Try this on for size:

# generate dataset
set.seed(1)
X <- rnorm(600, 7, 2)

# Create a function that samples 40 elements from X,
#  and calculates Pearson's r for the first 20 elements 
#  against the last 20 elements.
booties <- function(x) {
  X.samp <- sample(x, 40)
  cor(X.samp[1:20], X.samp[21:40])
}

# Replicate this function 1000 times (spits out a vector of cor estimates)
Z <- replicate(1000, booties(X))
error <- qt(0.975, length(Z)-1 * sd(Z)/sqrt(length(Z)))

1000 replicates take around 0.08 sec to complete at my end (about an order of magnitude faster than the for loop you were experimenting with).

0

In general implicit loops are faster the explicit loops. Try taking the code inside your loop and place it in a function, and then using that function in an lapply, or sapply statement.

myfunction = function(<insert relevant parameters here>)
{ 
  X.sample <- X[ sample(1:nrow(X), 40, replace = FALSE), ] 
  X.sample.1 <- X.sample[1:20, ]
  X.sample.2 <- X.sample[21:40, ]
  Y <- as.data.frame(cbind(X.sample.1$ID, X.sample.1$x, X.sample.2$ID,  X.sample.2$x))
  cor.results <- cor.test(Y[,2], Y[,4], alternative = c("greater"), method = c("pearson"))
  cor.results$estimate
}

Z  = sapply(x, myfunction)
#Here every element of x contains the arguments you want to pass to my function
#You can pass multiple arguments separated by commas after the function name

error <- qt(0.975, df = (length(Z) - 1)) * (sd(Z))/sqrt(length(Z))

You could do this, but I find it's probably better to just use the boot() function in the boot package if you can.

As for the set.seed() You need to set it directly before EVERY time you generate random anything. See below.

> rnorm(6)
[1]  1.0915017 -0.6229437 -0.9074604 -1.5937133  0.3026445  1.6343924
> set.seed(1001)
> rnorm(6)
[1]  2.1886481 -0.1775473 -0.1852753 -2.5065362 -0.5573113 -0.1435595
> set.seed(1001)
> rnorm(6)
[1]  2.1886481 -0.1775473 -0.1852753 -2.5065362 -0.5573113 -0.1435595
> rnorm(6)
[1]  1.0915017 -0.6229437 -0.9074604 -1.5937133  0.3026445  1.6343924


> set.seed(1001)
> sample(1:5,10,replace=T)
 [1] 5 3 3 3 3 5 1 1 2 4
> sample(1:5,10,replace=T)
 [1] 3 1 5 3 2 5 1 2 1 4
> set.seed(1001)
> sample(1:5,10,replace=T)
 [1] 5 3 3 3 3 5 1 1 2 4
> rnorm(6)
[1] -0.1435595  1.0915017 -0.6229437 -0.9074604 -1.5937133  0.3026445
> set.seed(1001)
> rnorm(6)
[1]  2.1886481 -0.1775473 -0.1852753 -2.5065362 -0.5573113 -0.1435595

Hope that helps!

When researching the boot function to give you an example I ran into a snag. It only ever returns one row. Strange! I might start a new question about this. In anycase, I think the bootstrap() function in the bootstrap package will do what you are looking for. Here's my example

set.seed(1001)
X <- rnorm(600, 7, 2)


myStat <- function(x, pairs) {
index = sample(1:length(x),(pairs*2))
Z = cor(X[index[1:(length(index)/2)]], X[index[((length(index)/2)+1):length(index)]])
return(Z)
}

b=bootstrap(X,1000,myStat,pairs=20)
Z <- b$thetastar
error <- qt(0.975, length(Z)-1 * sd(Z)/sqrt(length(Z)))
1
  • Thanks for the good advice. Can you make a suggestion of how to create the dataset (note good advice above the ID field is completely unnecessary) and use the boot command? I tried and that is why I feel back on the less than perfect for-loop. Thanks, Keith Mar 15, 2012 at 21:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.