zoukankan      html  css  js  c++  java
  • R语言代写回归中的Hosmer-Lemeshow拟合优度检验

    原文链接:http://tecdat.cn/?p=6166

    在依赖模型得出结论或预测未来结果之前,我们应尽可能检查我们假设的模型是否正确指定。也就是说,数据不会与模型所做的假设冲突。对于二元结果,逻辑回归是最流行的建模方法。在这篇文章中,我们将看一下 Hosmer-Lemeshow逻辑回归的拟合优度检验。

    Hosmer-Lemeshow拟合优度检验


    Hosmer-Lemeshow拟合优度检验是基于根据预测的概率或风险将样本分开。具体而言,基于估计的参数值,对于样本中的每个观察,基于每个观察的协变量值计算概率。

    然后根据样本的预测概率将样本中的观察分成g组(我们回过头来选择g)。假设(通常如此)g = 10。然后第一组由具有最低10%预测概率的观察组成。第二组由预测概率次之小的样本的10%等组成。

     在实践中,只要我们的一些模型协变量是连续的,每个观测将具有不同的预测概率,因此预测的概率将在我们形成的每个组中变化。为了计算我们预期的观察数量,Hosmer-Lemeshow测试取组中预测概率的平均值,并将其乘以组中的观察数。测试也执行相同的计算,然后计算Pearson拟合优度统计量

       

    选择组的数量


    就我所见,关于如何选择组数g的指导很少。Hosmer和Lemeshow的模拟结论是基于使用的,建议如果我们在模型中有10个协变量 。

    直观地说,使用较小的g值可以减少检测错误规范的机会。 

    R 

    首先,我们将使用一个协变量x模拟逻辑回归模型中的一些数据,然后拟合正确的逻辑回归模型。 

    n < -  100
    x < - rnorm(n)
    xb < - x
    pr < - exp(xb)/(1 + exp(xb))
    y < - 1 *(runif(n)<pr)
    mod < - glm(y~x,family = binomial)

    接下来,我们将结果y和模型拟合概率传递给hoslem.test函数,选择g = 10组:

            Hosmer and Lemeshow goodness of fit (GOF) test
    
    data:  mod$y, fitted(mod)
    X-squared = 7.4866, df = 8, p-value = 0.4851

    这给出p = 0.49,表明没有合适的不良证据。 我们还可以从我们的hl对象中获得一个观察到的与预期的表:

    cbind(hl$observed,hl$expected)
    y0 y1 yhat0 yhat1
    [0.0868,0.219] 8 2 8.259898 1.740102
    (0.219,0.287] 7 3 7.485661 2.514339
    (0.287,0.329] 7 3 6.968185 3.031815
    (0.329,0.421] 8 2 6.194245 3.805755
    (0.421,0.469] 5 5 5.510363 4.489637
    (0.469,0.528] 4 6 4.983951 5.016049
    (0.528,0.589] 5 5 4.521086 5.478914
    (0.589,0.644] 2 8 3.833244 6.166756
    (0.644,0.713] 6 4 3.285271 6.714729
    (0.713,0.913] 1 9 1.958095 8.041905

    为了帮助我们理解计算,现在让我们自己手动执行测试。首先,我们计算模型预测概率,然后根据预测概率的十分位数对观测值进行分类:

    pihat <- mod$fitted
    pihatcat <- cut(pihat, brks=c(0,quantile(pi 1,0.9,0.1)),1),  els=FALSE)

    接下来,我们循环通过组1到10,计算观察到的0和1的数量,并计算预期的0和1的数量。为了计算后者,我们找到每组中预测概率的均值,并将其乘以组大小,这里是10:

    meanprobs <- array(0, dim=c(10,2))
    expevents <- array(0, dim=c(10,2))
    obsevents <- array(0, dim=c(10,2))
    
    for (i in 1:10) {
    	meanprobs[i,1] <- mean(pihat[pihatcat==i])
     
    	obsevents[i,2] <- sum(1-y[pihatcat==i])
    }

    最后,我们可以通过表格的10x2单元格中的(观察到的预期)^ 2 /预期的总和来计算Hosmer-Lemeshow检验统计量:

     
    [1] 7.486643

    与hoslem.test函数的测试统计值一致。

    改变组的数量
    接下来,让我们看看测试的p值如何变化,因为我们选择g = 5,g = 6,直到g = 15。我们可以通过一个简单的for循环来完成:

    for(i in 5:15){
    	print(hoslem.test(mod $ y,fits(mod),g = i)$ p.value)
    }
    [1] 0.4683388
    [1] 0.9216374
    [1] 0.996425
    [1] 0.9018581
    [1] 0.933084
    [1] 0.4851488
    [1] 0.9374381
    [1] 0.9717069
    [1] 0.5115724
    [1] 0.4085544
    [1] 0.8686347

    虽然p值有所改变,但它们都显然不重要,所以他们给出了类似的结论,没有证据表明不合适。因此,对于此数据集,选择不同的g值似乎不会影响实质性结论。

    通过模拟检查Hosmer-Lemeshow测试


    要完成,让我们进行一些模拟,以检查Hosmer-Lemeshow测试在重复样本中的表现。首先,我们将从先前使用的相同模型重复采样,拟合相同(正确)模型,并使用g = 10计算Hosmer-Lemeshow p值。我们将这样做1000次,并将测试p值存储在一个数组中:

    pvalues < -  array(0,1000)
    
    for(i in 1:1000){
    	n < -  100
    	x < -  rnorm(n)
     	pr < -  exp(xb)/(1 + exp(xb))
     	mod < -  glm(y~x,family = binomial)
     }

    完成后,我们可以计算出p值小于0.05的比例。由于此处正确指定了模型,因此我们希望这种所谓的类型1错误率不大于5%:

     
    [1] 0.04

    因此,在1,000次模拟中,Hosmer-Lemeshow测试在4%的情况下给出了显着的p值,表明不合适。所以测试错误地表明在我们预期的5%限制内不合适 - 它似乎工作正常。

    现在让我们改变模拟,以便我们适合的模型被错误地指定,并且应该很难适应数据。希望我们会发现Hosmer-Lemeshow测试在5%的时间内正确地找到了不合适的证据。具体来说,我们现在将生成跟随具有协变量的逻辑模型,但我们将继续使用线性协变量拟合模型,以便我们的拟合模型被错误地指定。 

     

    我们发现,计算p值小于0.05的比例


    [1] 0.648

    因此,Hosmer-Lemeshow测试为我们提供了65%的不合适的重要证据。

    如果您有任何疑问,请在下面发表评论。

     
  • 相关阅读:
    POJ 1681 Painter's Problem(高斯消元法)
    HDU 3530 Subsequence(单调队列)
    HDU 4302 Holedox Eating(优先队列或者线段树)
    POJ 2947 Widget Factory(高斯消元法,解模线性方程组)
    HDU 3635 Dragon Balls(并查集)
    HDU 4301 Divide Chocolate(找规律,DP)
    POJ 1753 Flip Game(高斯消元)
    POJ 3185 The Water Bowls(高斯消元)
    克琳:http://liyu.eu5.org
    WinDbg使用
  • 原文地址:https://www.cnblogs.com/tecdat/p/11448145.html
Copyright © 2011-2022 走看看