模拟回歸模型的資料
驗證回歸模型的首選方法是模拟來自它們的資料，并檢視模拟資料是否捕獲原始資料的相關特征。感興趣的基本特征是平均值。我喜歡這種方法，因為它可以擴充到廣義線性模型（logistic，Poisson，gamma，...）和其他回歸模型，比如t -regression。
您的标準回歸模型假設存在将預測變量與結果相關聯的真實/固定參數。但是，當我們執行回歸時，我們隻估計這些參數。是以，回歸軟體傳回表示系數不确定性的标準誤差。
我将用一個例子來證明我的意思。
示範
我将使用泊松回歸來證明這一點。我模拟了兩個預測變量,使用50的小樣本。
n <- 50
set.seed(18050518)
xc的系數為0.5 ，xb的系數為1 。我對預測進行取幂，并使用該rpois()函數生成泊松分布結果。
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 與忽略系數不确定性時相比，均值和方差更接近原始結果的均值和方差。與考慮回歸系數的不确定性時相比，這種方法總是會導緻方差較小。它要快得多，并且需要零程式設計來實作，但我不習慣忽略回歸系數的不确定性，使模型看起來比它更充分。