3

我正在尝试使用 R 生成带有回归量的 Beta-Binomial 数据。我使用以下代码生成 Beta-Binomial 数据。现在我想在方程中添加一个协变量。感谢任何帮助。

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))和形状或过度离散度参数化基础 Beta 分布，该参数决定方差（通常等效于alpha+beta；较大的 theta 意味着较小的方差）。此外，将均值设为协变量的逻辑函数可能是最简单的方法（如果需要，可以替换为不同的反向链接函数）。

包里的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为零，您应该得到与以前相同的结果（logistic(0)=0.5，与您的平均值相同），但我还没有实际检查过。

编辑：要获得此数据集的 300 个副本，

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

似乎工作（给出一个 20 x 300 矩阵）。替换plyr::raply和replicate转置（r*ply比 提供更多的一致性和控制replicate）也是如此。

于 2012-07-29T01:00:14.890 回答