## 示范

``````n <- 50
set.seed(18050518)
xc的系数为0.5 ，xb的系数为1 。我对预测进行取幂，并使用该函数生成泊松分布结果。
rpois()

``````
``````# 指数预测传递给 rpois()

summary(dat)

xc                  xb             y
Min.   :-2.903259   Min.   :0.00   Min.   :0.00
1st Qu.:-0.648742   1st Qu.:0.00   1st Qu.:1.00
Median :-0.011887   Median :0.00   Median :2.00
Mean   : 0.006109   Mean   :0.38   Mean   :2.02
3rd Qu.: 0.808587   3rd Qu.:1.00   3rd Qu.:3.00
Max.   : 2.513353   Max.   :1.00   Max.   :7.00
``````

``````
Call:
glm(formula = y ~ xc + xb, family = poisson, data = dat)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.9065  -0.9850  -0.1355   0.5616   2.4264

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.20839    0.15826   1.317    0.188
xc           0.46166    0.09284   4.973 6.61e-07 ***
xb           0.80954    0.20045   4.039 5.38e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 91.087  on 49  degrees of freedom
Residual deviance: 52.552  on 47  degrees of freedom
AIC: 161.84

Number of Fisher Scoring iterations: 5
``````

``````coefs <- mvrnorm(n = 10000, mu = coefficients(fit.p), Sigma = vcov(fit.p))
``````

``````coefficients(fit.p)

(Intercept)          xc          xb
0.2083933   0.4616605   0.8095403

colMeans(coefs) # 模拟系数的均值

(Intercept)          xc          xb
0.2088947   0.4624729   0.8094507
``````

``````sqrt(diag(vcov(fit.p)))

(Intercept)          xc          xb
0.15825667  0.09284108  0.20044809

apply(coefs, 2, sd) #模拟系数的标准差

(Intercept)          xc          xb
0.16002806  0.09219235  0.20034148
``````

``````#每种情况一行，每组模拟系数一行

#模型矩阵与系数的乘积，取幂，
#然后用于模拟泊松分布的结果
for (i in 1:nrow(coefs)) {
sim.dat[, i] <- rpois(n, exp(fit.p.mat %*% coefs[i ]))
}
rm(i, fit.p.mat) # ``````

``````c(mean(dat\$y), var(dat\$y))#原始结果的均值和方差

[1] 2.020000 3.366939

c(mean(colMeans(sim.dat)), mean(apply(sim.dat, 2, var)))#10,000个模拟结果的平均值和变量的平均值

[1] 2.050724 4.167751
``````

``````
[1] 3.907143
``````

绘制10,000个模拟数据集中的一些数据集的直方图并将其与原始结果的直方图进行比较也是有用的。还可以测试原始数据和模拟数据集中xb = 1和xb = 0之间结果的平均差异。

``````sim.default <- simulate(fit.p, 10000)
``````

``````sim.default <- replicate(10000, rpois(n, fitted(fit.p)))
``````

`fitted(fit.p)`是响应尺度的预测，或指数线性预测器，因为这是泊松回归。因此，我们将使用模型中的单组预测值来重复创建模拟结果。

``````c(mean(colMeans(sim.default)), mean(apply(sim.default, 2, var)),

[1] 2.020036 3.931580 3.810612
``````