Имитация набора данных для анализа загрузки багажника - PullRequest
Новая
       12

Имитация набора данных для анализа загрузки багажника

0 голосов
/ 15 марта 2012

Моя цель - использовать начальную загрузку (1000 повторений) для вычисления нулевого распределения, среднего значения и КИ r (коэффициент корреляции Пирсона), коррелирующего признак (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 повторений r и вычислю 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))

Ответы [ 2 ]

1 голос
/ 15 марта 2012

Примерьте размер: 

# 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, с которым вы экспериментировали).

0 голосов
/ 15 марта 2012

Вообще неявные циклы быстрее явных циклов.Попробуйте взять код внутри вашего цикла и поместить его в функцию, а затем использовать эту функцию в выражении lapply или 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)))
...