1

我在 R 中运行了 Hosmer Lemeshow 统计数据，但我得到的 p 值为 1。这对我来说似乎很奇怪。我知道高 p 值意味着我们不会拒绝观察到的和预期的相同的原假设，但是我有可能在某处有错误吗？

我如何解释这样的 p 值？

下面是我用来运行测试的代码。我还附上了我的模型的外观。响应变量是计数变量，而所有回归变量都是连续的。由于在我的初始泊松模型中检测到过度分散，我运行了一个负二项式模型。

> hosmerlem <- function(y, yhat, g=10)
+     {cutyhat <- cut(yhat, breaks = quantile(yhat, probs=seq(0,1, 1/g)), include.lowest=TRUE)
+     obs <- xtabs(cbind(1 - y, y) ~ cutyhat)
+     expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
+     chisq <- sum((obs - expect)^2/expect)
+     P <- 1 - pchisq(chisq, g - 2)
+     return(list(chisq=chisq,p.value=P))}
>   hosmerlem(y=TOT.N, yhat=fitted(final.model))
$chisq
[1] -2.529054

$p.value
[1] 1

> final.model <-glm.nb(TOT.N ~ D.PARK + OPEN.L + L.WAT.C + sqrt(L.P.ROAD))
> summary(final.model)

Call:
glm.nb(formula = TOT.N ~ D.PARK + OPEN.L + L.WAT.C + sqrt(L.P.ROAD), 
    init.theta = 4.979895131, link = log)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-3.08218  -0.70494  -0.09268   0.55575   1.67860  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     4.032e+00  3.363e-01  11.989  < 2e-16 ***
D.PARK         -1.154e-04  1.061e-05 -10.878  < 2e-16 ***
OPEN.L         -1.085e-02  3.122e-03  -3.475  0.00051 ***
L.WAT.C         1.597e-01  7.852e-02   2.034  0.04195 *  
sqrt(L.P.ROAD)  4.924e-01  3.101e-01   1.588  0.11231    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4.9799) family taken to be 1)

    Null deviance: 197.574  on 51  degrees of freedom
Residual deviance:  51.329  on 47  degrees of freedom
AIC: 383.54

Number of Fisher Scoring iterations: 1


              Theta:  4.98 
          Std. Err.:  1.22 

 2 x log-likelihood:  -371.542
4

1 回答 1

0

正如@BenBolker 正确指出的那样，Hosmer-Lemeshow 是逻辑回归的测试，而不是负二项式广义线性模型的测试。

如果我们考虑将测试应用于逻辑回归，函数的输入（包中函数hosmerlem的副本）应该是： - = 观察的数值向量，二进制 (0/1) - = 期望值（概率)hoslem.testResourceSelection
y
yhat

这是一个说明性示例，显示了如何获得正确的输入：

set.seed(123)
n <- 500
x <- rnorm(n)
y <- rbinom(n, 1, plogis(0.1 + 0.5*x))
logmod <- glm(y ~ x, family=binomial)

# Important: use the type="response" option
yhat <- predict(logmod, type="response")
hosmerlem(y, yhat)

########
$chisq
[1] 4.522719

$p.value
[1] 0.8071559

函数给出相同的结果hoslem.test：

library(ResourceSelection)
hoslem.test(y, yhat)

########
Hosmer and Lemeshow goodness of fit (GOF) test

data:  y, yhat
X-squared = 4.5227, df = 8, p-value = 0.8072
于 2017-06-05T12:47:00.997 回答