3

Rを使用してリグレッサを使用してベータ二項データを生成しようとしています。次のコードを使用してベータ二項データを生成しました。ここで、方程式に共変量を追加します。助けに感謝します。

set.seed(111) 
k<-20
n<-60
x<-NULL

p<-rbeta(k,3,3)# so that the mean nausea rate is alpha/(alpha+beta)
for(i in 1:k)
x<-cbind(x,rbinom(300,n,p[i]))

ありがとうアナミカ

4

1 に答える 1

4

私の推測では、結果の確率 (あなたの場合は吐き気) が共変量の関数であるモデルからデータを生成したいと考えています。これを行うための最も標準的な (唯一ではありませんが) 方法は、基になるベータ分布をその平均(alpha/(alpha+beta))と形状、または分散を決定するパラメーター (通常は に相当しalpha+betaます。シータが大きいほど分散が小さいことを意味します) の観点からパラメーター化することです。さらに、平均を共変量のロジスティック関数にするのがおそらく最も簡単です (必要に応じて、別の逆リンク関数に置き換えることもできます)。

パッケージ内のrbetabinom関数は、emdbookこの方法で既にパラメーター化されています (コードを見ることができます。それほど複雑ではありません)。

set.seed(111)
beta0 <- 0 ## logit rate at x=0
beta1 <- 2 ## increase in logit-prob(nausea) per unit x
k <- 20
n <- 60
theta <- 6  ## shape parameter, equivalent to alpha+beta
x <- runif(k)  ## distribution of covariates 
               ## (you might want something different)

library(emdbook)
eta <- beta0+beta1*x  ## linear predictor
prob <- plogis(eta)   ## logistic transform
y <- rbetabinom(k, prob=prob, size=n, theta=6)

beta1ゼロに設定すると、以前と同じ結果が得られるはずです（ロジスティック（0）= 0.5、平均と同じ）が、実際には確認していません。

edit : このデータセットの 300 個の複製を取得するには、

Y <- replicate(300,rbetabinom(k, prob=prob, size=n, theta=6))

動作しているようです (20 x 300 の行列が得られます)。置換と転置も同様です (plyr::raplyより一貫性と制御性が少し向上します)。replicater*plyreplicate

于 2012-07-29T01:00:14.890 に答える