r - 使用协变量生成 Beta-二项式响应

Question

我正在尝试使用 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]))

谢谢阿纳米卡

Answer 1

我的猜测是您想要从一个模型生成数据，其中结果的概率(在您的案例中是恶心)是协变量的函数。最标准(尽管不是唯一)的方法是根据其均值 (alpha/(alpha+beta)) 和形状或过度分散参数来参数化基础 Beta 分布确定方差(通常等同于 alpha+beta；较大的 theta 意味着较小方差)。此外，使均值成为协变量的逻辑函数可能是最简单的(如果需要，您可以替换为不同的反向链接函数)。

emdbook 包中的rbetabinom 函数已经用这种方式参数化了(你可以看看代码——它不是很复杂)。

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< 提供了更多的一致性和控制力)。