R語言進行數值模擬：模擬泊松迴歸模型的數據

原文連接：http://tecdat.cn/?p=6751

模擬迴歸模型的數據

驗證迴歸模型的首選方法是模擬來自它們的數據，並查看模擬數據是否捕獲原始數據的相關特徵。感興趣的基本特徵是平均值。我喜歡這種方法，由於它能夠擴展到廣義線性模型（logistic，Poisson，gamma，...）和其餘迴歸模型，好比_t_ -regression。

您的標準迴歸模型假設存在將預測變量與結果相關聯的真實/固定參數。可是，當咱們執行迴歸時，咱們只估計這些參數。所以，迴歸軟件返回表示係數不肯定性的標準偏差。

我將用一個例子來證實個人意思。

示範

我將使用泊松迴歸來證實這一點。我模擬了兩個預測變量,使用50的小樣本。

n <- 50
set.seed(18050518) 
xc的係數爲0.5 ，xb的係數爲1 。我對預測進行取冪，並使用該rpois()函數生成泊松分佈結果。

# Exponentiate prediction and pass to 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

估計的係數與人口模型相距不太遠，.21表明截距而不是0，.46而不是.5，而0.81而不是1。

接下來模擬模型中的數據，我想要10,000個模擬數據集，爲了捕捉迴歸係數的不肯定性，我假設係數來自多元正態分佈，估計係數做爲均值，迴歸係數的方差 - 協方差矩陣做爲多元正態分佈的方差 - 協方差矩陣。

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) # means of simulated coefficients

(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) # standard deviation of simulated coefficients

(Intercept)          xc          xb
 0.16002806  0.09219235  0.20034148

下一步是模擬模型中的數據。咱們經過將模擬係數的每一行乘以原始預測變量來實現。而後咱們傳遞預測：

# One row per case, one column per simulated set of coefficients


 # Cross product of model matrix by coefficients, exponentiate result,
# then use to simulate Poisson-distributed outcome
for (i in 1:nrow(coefs)) {
  sim.dat[, i] <- rpois(n, exp(fit.p.mat %*% coefs[i ]))
}
rm(i, fit.p.mat) # Clean house

如今一個是完成模擬，將模擬數據集與原始數據集至少比較結果的均值和方差：

c(mean(dat$y), var(dat$y)) # Mean and variance of original outcome

[1] 2.020000 3.366939

c(mean(colMeans(sim.dat)), mean(apply(sim.dat, 2, var))) # average of mean and var of 10,000 simulated outcomes

[1] 2.050724 4.167751

模擬結果的平均值略高於原始數據，平均方差更高。平均而言，能夠預期方差比平均值更偏離目標。方差也將與一些極高的值正誤差，同時，它的界限爲零，所以中位數可能更好地反映了數據的中心：

[1] 3.907143

中位數方差更接近原始結果的方差。

這是模擬均值和方差的分佈：

​

 繪製10,000個模擬數據集中的一些數據集的直方圖並將其與原始結果的直方圖進行比較也是有用的。還能夠測試原始數據和模擬數據集中xb = 1和xb = 0之間結果的平均差別。 

---

回到基礎R，它具備simulate()執行相同操做的功能：

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

與忽略係數不肯定性時相比，均值和方差更接近原始結果的均值和方差。與考慮迴歸係數的不肯定性時相比，這種方法老是會致使方差較小。它要快得多，而且須要零編程來實現，但我不習慣忽略迴歸係數的不肯定性，使模型看起來比它更充分。