文章目录
- 泊松回归模型
-
- 模型形式
- R语言建模
- 参考资料
泊松回归模型
临床中存在大量计数数据，事件发生次数是典型例子。患者随访次数和疾病发生次数等数据都属于计数数据。统计某种疾病的患病人数，分析一段时间和一定地区中事件发生数量均依赖于离散分布。进一步分析这些计数数据的影响因素，则需要引入解释变量，并对计数数据进行建模。
模型形式
泊松分布是常见的离散分布之一，适用于单位时间或空间内的计数变量。线性回归是基本的统计模型。病人情况各不相同，临床资料表现出异质性，因此样本往往不服从同分布，即相应泊松分布的参数不同。但是，临床上的异质性可由某些变量解释。或者说，给定解释变量取值，响应变量的分布参数可能是确定的。更具体地，对给定解释变量，可以对响应变量的条件期望建模。由于泊松分布的概率分布函数通常为指数形式，可取条件期望后与解释变量建立线性回归模型。
y
i
∣
x
i
∼
P
(
μ
i
)
,
i
=
1
,
2
,
.
.
.
,
n
y_i | x_i \\sim P(\\mu_i), i = 1, 2, …, n
yi∣xi∼P(μi),i=1,2,…,n
l
o
g
(
μ
i
)
=
β
0
+
β
1
x
i
1
+
β
2
x
i
2
+
β
3
x
i
3
,
i
=
1
,
2
,
.
.
.
,
n
log(\\mu_i) = \\beta_0 + \\beta_1 x_{i1} + \\beta_2 x_{i2} + \\beta_3 x_{i3}, i = 1, 2, …, n
log(μi)=β0+β1xi1+β2xi2+β3xi3,i=1,2,…,n
其中的对数函数称为连接函数。
模型的拟合优度检验有两种方式。一种是通过Pearson卡方统计量，另一种是计算偏差统计量（deviance statistic）。
泊松分布的期望和方差相等，称为等离散（equidispersion）。实际场景中，方差可能大于期望，即过离散（overdispersion）。方差也可能小于期望，即欠离散（underdispersion）。过离散较为常见，其检测依赖于Pearson统计量和偏差统计量除以自由度。计算结果接近1，表明等离散。结果远大于1，表示过离散。自由度由样本量n减去回归模型的系数个数p给出。一旦发现过离散，基本的泊松模型便不适用。负二项
分布是适用于过离散数据的模型之一。
R语言建模
先随机生成符合泊松模型的数据：
set.seed(6543) # 保证每次运行结果一致
sample_size = 50000
# 生成符合泊松模型的数据
x1 = runif(sample_size)
x2 = runif(sample_size)
x3 = runif(sample_size)
log_mu = 1 + 0.75 * x1 – 1.25 * x2 + .5 * x3
mu = exp(log_mu)
py = rpois(sample_size, mu)
对计数数据进行建模：
# 基于泊松模型建模数据
glm.pois = glm(py ~ x1 + x2 + x3, family = poisson)
confint(glm.pois) # 通过基于likelihood的标准误计算置信区间
confint.default(glm.pois) # 通过基于模型的标准误计算置信区间
summary(glm.pois) # 概括建模结果
pchi2 = sum(resid(glm.pois, type = "pearson") ^ 2) # 计算Pearson离散统计量
# Pearson统计量除以自由度，值接近1，表明等离散（equidispersed）
(pdisp = pchi2 / glm.pois$df.residual)
(beta = coef(glm.pois)) # 模型系数
exp(beta) # 系数指数化，表示解释变量的单位增加能引起响应变量增加的倍数
模型估计的系数与生成时给定的系数接近。
运行Monte Carlo模拟数据100次迭代，可计算出统计量的均值。
# Monte Carlo模拟100次，计算模型系数的平均值
set.seed(6543)
mysim = function() {
sample_size = 50000
# 生成符合泊松模型的数据
x1 = runif(sample_size)
x2 = runif(sample_size)
x3 = runif(sample_size)
log_mu = 1 + 0.75 * x1 – 1.25 * x2 + .5 * x3
mu = exp(log_mu)
py = rpois(sample_size, mu)
# 基于泊松模型建模数据
glm.pois = glm(py ~ x1 + x2 + x3, family = poisson)
# 计算Pearson离散统计量
pchi2 = sum(resid(glm.pois, type = "pearson") ^ 2)
# Pearson统计量除以自由度，值接近1，表明等离散（equidispersed）
pdisp = pchi2 / glm.pois$df.residual
beta = coef(glm.pois) # 模型系数
list(beta, pdisp)
}
B = replicate(100, mysim())
# 截距和x1, x2, x3的系数均值
apply(matrix(unlist(B[1,]), 4, 100), 1, mean)
# Pearson离散统计量均值
mean(unlist(B[2, ]))A
