0

私の目標は、ブートストラップ (1000 担当者) を使用して、600 の一意の個人 (ID) のデータセットから生成された 20 の刺激されたランダム ペアで特性 (x) に相関する r (ピアソンの相関係数) の帰無分布、平均、および CI を計算することです。最近、「proc surveyselect」を使用してデータセットを生成する SAS から R に切り替えました。質問:

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

次に、r の 1000 回の複製を生成し、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))
4

2 に答える 2

1

これを試着してサイズを確認してください：

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

for私の最後では、1000回の複製が完了するまでに約0.08秒かかります（実験していたループよりも約1桁高速です）。

于 2012-03-15T13:52:46.797 に答える
0

一般に、暗黙的なループは明示的なループよりも高速です。ループ内のコードを取得して関数に配置し、その関数を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、私は障害に遭遇しました。1行だけを返します。変！私はこれについて新しい質問を始めるかもしれません。いずれにせよ、パッケージ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)))
于 2012-03-15T13:33:00.133 に答える