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

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 定义了一个值，而不是 'feeding' 它是公式。

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

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

Failing to simulate data for a negative binomial probability distribution

simulation

r

bayesian

linear-regression