拓端tecdat|R語言輔導Poisson回歸的拟合優度檢驗

R語言Poisson回歸的拟合優度檢驗

在這篇文章中，我們将看一下Poisson回歸的拟合優度測試與個體計數資料。許多軟體包在拟合Poisson回歸模型時在輸出中提供此測試，或者在拟合此類模型（例如Stata）之後執行此測試，這可能導緻研究人員和分析人員依賴它。在這篇文章中，我們将看到測試通常不會按預期執行，是以，我認為，應該謹慎使用。

偏差拟合度檢驗

由于偏差度量衡量了模型預測與觀察結果的接近程度，我們可能會考慮将其作為給定模型拟合度檢驗的基礎。雖然我們希望我們的模型預測接近觀察到的結果，但即使我們的模型被正确指定，它們也不會相同 - 畢竟，模型給出了觀察所遵循的泊松分布的預測平均值。

是以，為了将偏差用作拟合優度檢驗，我們需要弄清楚，假設我們的模型是正确的，在泊松假設下，我們在預測均值周圍觀察到的結果中會有多少變化。由于偏差可以作為将目前模型與飽和模型進行比較的輪廓似然比檢驗得出，是以可能性理論會預測（假設模型被正确指定），偏差遵循卡方分布，自由度等于參數數量的差異。飽和模型可以被視為一個模型，它為每個觀察使用不同的參數，是以它具有參數。如果我們提出的模型具有參數，這意味着将偏差與參數的卡方分布進行比較。

在R中執行拟合優度測試

現在看看如何在R中執行拟合優度測試。首先我們将模拟一些簡單的資料，具有均勻分布的協變量x和泊松結果y：

set.seed（612312）

n < -  1000
x < -  runif（n）
y < -  rpois（n，mean）      

為了使Poisson GLM适合資料，我們隻需使用glm函數：

Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2547  -0.8859  -0.1532   0.6096   3.0254  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.04451    0.05775  -0.771    0.441    
x            1.01568    0.08799  11.543   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1247.7  on 999  degrees of freedom
Residual deviance: 1110.3  on 998  degrees of freedom
AIC: 3140.9

Number of Fisher Scoring iterations: 5      

這裡的偏差被glm函數标記為“剩餘偏差”，這裡是1110.3。有1000個觀測值，我們的模型有兩個參數，是以自由度為998，由R作為殘差df給出。為了計算偏差拟合度檢驗的p值，我們簡單地計算998自由度上卡方分布的偏內插補點右側的機率：

pchisq（mod $ deviance，df = mod $ df.residual，lower.tail = FALSE）
[1] 0.00733294      

零假設是我們的模型被正确指定，我們有強有力的證據來拒絕這個假設。是以，我們有充分的證據表明我們的模型非常适合。 

通過仿真檢驗泊松回歸拟合檢驗的偏差優度

為了研究測試的性能，我們進行了一個小的模拟研究。我們将使用與以前相同的資料生成機制生成10,000個資料集。對于每一個，我們将拟合（正确的）泊松模型，并收集拟合p值的偏差良好性。然後我們将看到它小于0.05的次數：

nSim <- 10000
pvalues <- array(0, dim=nSim)

for (i in 1:nSim) {

n <- 1000
x <- runif(n)
mean <- exp(x)
y <- rpois(n,mean)

mod <- glm(y~x, family=poisson)
pvalues[i] <- pchisq(mod$ , df=mod$df. , lower.tail= )

}

mean(1*(pvalues<0.05))      

最後一行建立一個向量，其中如果p值小于0.05，則每個元素為1，否則為零，然後使用mean（）計算這些元素的比例。當我運作這個時，我得到了0.9437，這意味着偏差測試錯誤地表明我們的模型在94％的情況下被錯誤地指定 

為了在平均值較大時檢視情況是否發生變化，讓我們修改模拟。我們現在将生成具有泊松均值的資料，其結果為20到55：

nSim < -  10000
pvalues < -  array（0，dim = nSim）

for（i in 1：nSim）{

n < -  1000
x < -  runif（n）
 < -  exp（3 + x）
y < -  rpois（n，mean）

mod < -  glm（y~x，family = poisson）
pvalues [i] < -  pchisq（mod $  ，df = mod $ df. ，lower.tail = FALSE）

}      

現在，顯着偏差測試的比例降低到0.0635，更接近标稱的5％1類錯誤率。 

結論