模拟速率数据（负二项分布）？

Question

我正在尝试模拟接近速率数据的数据 - 即：计算通常符合负二项分布但也适合调查工作的偏移项的数据。

我认为我可以使用负二项式函数 (rnbinom()) 很好地模拟计数，但是没有办法解释偏移项 - 这每次调查都是随机的。换句话说：

模拟非整数速率数据的最佳方法是什么？

是否有模拟偏移项值的实际方法？

我是否需要使用除负二项分布之外的其他分布来生成实际范围内的非整数值？

背景：我们的调查衡量每单位调查工作量（时间）的个人数量，得出的比率是正非整数 ( >= 0)。调查计数数据似乎很好地模拟了负二项分布，并且在 GLM 框架中，我会使用调查时间作为抵消项来计算工作量。在下面的模拟数据中，我生成了一个负二项分布来表示我实际调查中的数据。偏移项被模拟为 2-10（搜索时间范围，以分钟为单位）之间的随机均匀变量。然后计算速率为 counts/time.

我绘制了计数和比率的直方图，以帮助证明此处的比率在整数之间采用许多小数值。由于调查计数通常与调查工作量相关，因此关键我最终使用比率数据进行分析（即下图 'B'）。


library(tidyverse)
theme_set(theme_classic())


d = data.frame(counts = rnbinom(n = 500,mu = 5,size = 1),  ## dispersion parameter 'theta' set to 1
               time = runif(500,2,10)) %>% 
  mutate(rate = counts/time)


## Count data histogram
 ggplot(d, aes(x = counts))+
  geom_histogram(fill = 'peachpuff',color = 'black')+
  ylab('frequency')+
  scale_y_continuous(expand = c(0,0))+
  ggtitle('A: Histogram of Counts')

## Rate data histogram
ggplot(d, aes(x = rate))+
  geom_histogram(binwidth = .1, fill = 'dodgerblue1',color = 'black')+
  scale_x_continuous(breaks = seq(0,10,1))+
  scale_y_continuous(expand = c(0,0))+
  ylab('frequency')+
  ggtitle('B: Histogram of Rate')

下面我可以很容易地从我的原始调查数据中模拟计数，但不知道如何正确模拟回溯数据作为比率。例如，如果我适合仅截距 nbinom GLM，我可以使用该系数来模拟与原始数据非常相似的计数的新负二项式分布（即对 'mu' 使用相似的值）

[我意识到在这个例子中这似乎是循环的，但这是我处理真实数据的方法。首先用 GLM 描述平均值和离散度 'theta'，然后模拟模仿我的原始数据集的反向数据集]

我在下面使用这种方法来生成回溯计数数据，而且还通过使用偏移项拟合模型来模拟回图 'B' 中具有均值 'rate' 的分布.


### Simulate back count data from the original survey data:

##describe mean value 'mu' by finding intercept
## 'theta'  could also be calculated
m1 = MASS::glm.nb(counts ~ 1, data = d)

# summary(m1)
# mean(d$counts)
# exp(m1$coefficients[1])


## simulated negative-binomial distribution using calculated 'mu'
d.sim = data.frame(new.counts = rnbinom(500,
                                        mu = as.numeric(exp(m1$coefficients[1])),  ## coef on log-scale, exponentiate to use
                                        size = 1))  ## holding dispersion parameter 'theta'  constant at 1


## Plot and compare with plot 'A' above
ggplot(d.sim, aes(x = new.counts))+
  geom_histogram(fill = 'peachpuff3',color = 'black')+
  ylab('frequency')+
  scale_y_continuous(expand = c(0,0))+
  ggtitle('C: Simulated Counts')


###########################################
###########################################

### Simulate back 'rate' data by including an offset term for effort in the GLM model
  ## the exponentiated coefficient should equal the mean of the raw rate data

m2 = MASS::glm.nb(counts ~ 1 + offset(log(time)), data = d)

# summary(m2) 
# mean(d$rate)
# exp(m2$coefficients[1])


d.sim.2 = data.frame(new.counts = rnbinom(500,
                                        mu = as.numeric(exp(m2$coefficients[1])),  ## coef on log-scale, exponentiate to use
                                        size = 1))    ## holding dispersion parameter 'theta'  constant at 1


## compare these simulated 'rate' data with the non-integer 'true rate' data in figure D
ggplot(d.sim.2, aes(x = new.counts))+
  geom_histogram(binwidth = .1, fill = 'dodgerblue3',color = 'black')+
  scale_x_continuous(breaks = seq(0,10,1))+
  ylab('frequency')+
  scale_y_continuous(expand = c(0,0))+
  ggtitle('D: Simulated Rate')

所以在这一点上，我生成了图 'C' 作为模拟数据集，代表我在现实生活中观察到的计数，它与图 'A' 中的原始数据非常匹配。图'D'中的'rate'数据是（必然）从rnbinom()中提取的所有整数值，而图 'D' 的均值近似于图 'B' 的均值，我的感觉是这两个分布并不完全等价。

所以我的问题又来了：

有没有一种方法可以代替模拟数据来匹配数字 'B'（非整数率数据）？
您认为图 'D' 中的数据与 'B' 的数据近似，因为平均值（和离差）相似吗？

对于其他上下文，我将使用模拟数据集（其中很多）来运行其他蒙特卡洛类型的模拟分析（例如功率分析）。我担心如果我使用图 'D' 中生成的数据，它不会真正代表我的实际调查数据（图 'B'）。

Answer 1

您生成样本数据的方式（代替您的经验数据）与您描述的数据生成过程不一致。来自 rnbinom(n = 500, mu = 5, size = 1) 的计数数据不依赖于时间。 mu 需要是时间变量的函数，否则计数与时间无关。

另外，设置size = 1意味着没有过度离散（也没有欠离散），因此它应该被称为泊松分布，这是负二项分布的特例。但是根据您对 DGP 的描述，听起来经验数据中会出现过度分散。

要回答您的第一个问题，请参阅下面的代码示例。关于你的第二个问题，不，我认为这不是个好主意。

library(tidyverse)
library(rstanarm)
options(mc.cores = parallel::detectCores())

n <- 1000

empirical <-
  tibble(
    time = runif(n, 2, 10),
    count = rnbinom(n = n, mu = time, size = 1) # Generate count data that actually depends on time
    ) |> 
  mutate(rate = count/time)

m_stan <- stan_glm.nb(count ~ time, data = empirical)
 
simulated <- 
  tibble(
    time = runif(n, 2,10),
  ) %>%
  mutate(
    count = posterior_predict(m_stan, ., draws = 1) |> 
      as.vector(),
    rate = count/time
    )

d <- lst(simulated, empirical) |> 
  bind_rows(.id = "data")

d |>
  select(data, count, rate) |> 
  pivot_longer(c(count, rate)) |> 
  ggplot() +
  geom_histogram(aes(value), binwidth = .2) +
  facet_grid(data ~ name, scales = "free")

^{由 reprex package (v2.0.1)}

创建于 2022-02-03

模拟速率数据（负二项分布）？

Simulating rate data (negative-binomial distribution)?

r

simulation

distribution

montecarlo

poisson