Имитированный набор данных для анализа обвязки ботинком

Question

Моя цель - использовать загрузку (1000 повторений) для вычисления нулевого распределения, среднего и CI коэффициента корреляции (Pearson) корреляции (x) в 20 стимулированных случайных пар, генерируемых из моего набора данных из 600 уникальных людей (ID). Недавно я переключился на R из SAS, где я бы использовал «proc surveyselect» для создания набора данных. Вопросы:

1. Что было бы самым эффективным способом генерации этих результатов (см. Мою попытку ниже)?
2. В моем примере, как я могу использовать команду set.seed для репликации моих результатов?

смоделированные начиная набора данных с 600 лиц и связанных с ними признака значений:

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

Я тогда генерируют мои 1000 повторов г и вычислить 95% ДИ:

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)) 

Answer 1

В общем неявное циклы быстрее явных циклов. Попробуйте взять код внутри вашего цикла и поместить его в функцию, а затем использовать эту функцию в сжатом виде или в инструкции sapply.

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)) 

Вы можете сделать это, но я считаю, что это, вероятно, лучше всего использовать функцию boot() в boot пакете, если вы можете.

Что касается set.seed() Вам необходимо установить его непосредственно перед КАЖДОМ временем, когда вы произвольно генерируете случайное. Смотри ниже.

> 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 

Надеюсь, что это поможет!

При исследовании функции boot, чтобы дать вам пример, я столкнулся с препятствием. Он только возвращает одну строку. Странный! Я мог бы задать новый вопрос. В любом случае, я думаю, что функция bootstrap() в пакете bootstrap сделает то, что вы ищете. Вот мой пример

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))) 

Answer 2

Попробуйте это для размера:

# 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 размножается займет около 0,08 секунд, чтобы завершить в моем конце (примерно на порядок быстрее, чем for петли вы экспериментировали с).