未能模拟负二项式概率分布的数据

Question

我希望有人可以帮助我。

在我参加的一个初学者研讨会上，在拟合多元回归模型的过程中，讲师最初使用泊松分布对结果进行了预先预测检查。这是分两步完成的。最初，创建了一个函数：

    multiple_regression_poisson_dgp <- function(predictor1, 
          predictor2, alpha_mean, alpha_sd, beta_predictor1_mean, 
          beta_predictor1_sd,  beta_predictor2_mean, 
          beta_predictor2_sd) {
    N <- length(predictor1)
    alpha <- rnorm(1, mean = alpha_mean, sd = alpha_sd);
    beta_predictor1 <- rnorm(1, mean = beta_predictor1_mean, 
            sd = beta_predictor1_sd);
    beta_predictor2 <- rnorm(1, mean = beta_predictor2_mean, 
            sd = beta_super_sd);
    outcome <- rpois(N, lambda = alpha + beta_predictor1 * 
                predictor1 + beta_predictor2 * predictor2)
    return(outcome)
    }

创建此函数后，生成了以下先验：

    multiple_regression_poisson_dgp(dataset$predictor1,
                                dataset$predictor2,
                                alpha_mean = 1,
                                alpha_sd = 0.5,
                                beta_predictor1_mean = -0.25,
                                beta_predictor1_sd = 0.5,
                                beta_predictor2_mean = 0,
                                beta_predictor2_sd = 1)

这很好用。问题是，进一步来看，泊松分布并不是最合适的。下一步建议使用负二项式。不幸的是，当我尝试复制负二项式的过程时，我没有成功。我试图复制上面显示的两个步骤，但对于负二项式。第一步编码为：

    multiple_regression_negative_binomial_dgp <- 
       function(predictor1, predictor2, alpha_mean, alpha_sd, 
       beta_predictor1_mean, beta_predictor1_sd, 
       beta_predictor2_mean, beta_predictor2_sd, phi_mean, 
       phi_sd) {
    N <- length(predictor1)
    alpha <- rnorm(1, mean = alpha_mean, sd = alpha_sd);
    beta_predictor1 <- rnorm(1, mean = beta_predictor1_mean, 
        sd = beta_predictor1_sd); 
    beta_predictor2 <- rnorm(1, mean = beta_predictor2_mean, 
        sd = beta_super_sd);
    phi <- rnorm(1, mean = phi_mean, sd = phi_sd);
    outcome<- rnbinom(N, size = mu + mu^2/phi, mu = alpha + 
          beta_predictor1 * predictor1 + beta_predictor2 * 
               predictor2)
    return(outcome)
    }

因为负二项式中有一个 phi，并且考虑到它将是我要计算其先验的参数，我假设需要将它添加到方程中。此外，鉴于rnbinom() 的文档，我认为我可以像对待泊松生成中的 lambda 一样对待 mu，并将回归方程提供给它。

这个函数可能不够完善，但是在我创建它并进入第二步之后，错误就出现了。第二步我编码为：

    multiple_regression_negative_binomial_dgp(dataset$predictor1,
                                dataset$predictor2,
                                alpha_mean = 1,
                                alpha_sd = 0.5,
                                beta_predictor1_mean = -0.25,
                                beta_predictor1_sd = 0.5,
                                beta_predictor2_mean = 0,
                                beta_predictor2_sd = 1,
                                phi_mean = 0,
                                phi_sd = 1)

但是，一旦我尝试运行此数据生成过程，就会收到警告说明：

    Error in rnbinom(N, size = mu, mu = alpha + beta_predictor1 * predictor1 + beta_predictor2 * predictor2 : object 'mu' not found

任何帮助将不胜感激，我意识到我正在应用一种更机械的思维方式来尝试为负二项式复制泊松数据生成过程，但我一直无法找到任何线索来解决这个问题.我遇到的大多数示例都为 mu 和 size 定义了一个值，而不是“输入”公式。

Answer 1

在您的multiple_regression_negative_binomial_dgp 函数中，您调用rnbinom。该函数需要一个size 参数，并且您将mu + mu^2/phi 分配给它，但mu 没有在函数中定义，也没有传递给它。 rnbinom 包含您提供的 mu 参数 (alpha + beta_predictor1 * predictor1 + beta_predictor2 * predictor2) 的事实并没有处理它，因为 rnbinom 不会将该信息传递给 size。我建议你试试：

multiple_regression_negative_binomial_dgp <- function(predictor1, predictor2,  
                                alpha_mean, alpha_sd, beta_predictor1_mean, 
                                beta_predictor1_sd, beta_predictor2_mean, 
                                beta_predictor2_sd, phi_mean,  phi_sd) {
  N               <- length(predictor1)
  alpha           <- rnorm(1, mean=alpha_mean,           sd=alpha_sd)
  beta_predictor1 <- rnorm(1, mean=beta_predictor1_mean, sd=beta_predictor1_sd) 
  beta_predictor2 <- rnorm(1, mean=beta_predictor2_mean, sd=beta_super_sd);
  phi             <- rnorm(1, mean=phi_mean,             sd=phi_sd)
  Mu              <- alpha + beta_predictor1*predictor1 + beta_predictor2*predictor2
  outcome         <- rnbinom(N, size=Mu + Mu^2/phi, mu=Mu)
  return(outcome)
}