logo标签列表

如何在R中使用anova()函数比较两个glm模型时提取p值

rmodelregressionglmanova
7

因此，我正在尝试比较两个模型fit1和fit2。

最初，我只是做了anova（fit1，fit2），这产生了我理解的输出（包括p值）。

然而，当我将我的模型从基于lm（）的模型切换到基于glm（）的模型时，anova（fit1，fit2）现在产生了剩余自由度，剩余偏差和Df偏差，我很难解释这些指标（解释这些指标的资源似乎很少）。我希望提取比较两个模型的p值，但由于某种原因anova（fit1，fit2，test ='Chisq'）不起作用。有什么建议吗？

我意识到，根据我的glms中的链接函数，卡方可能不是最合适的测试，但是在适当的情况下也使用了'F'，但结果令人失望。

有没有其他人遇到过这个问题？有什么建议吗？非常感谢！

示例：

make_and_compare_models <- function(fitness_trait_name, data_frame_name, vector_for_multiple_regression, predictor_for_single_regression, fam){
        fit1<-glm(formula=as.formula(paste(fitness_trait_name,"~", paste(vector_for_multiple_regression, sep="+"))), family=fam, data=data_frame_name)
        print ("summary fit 1")
        print(summary(fit1))
        fit2<- glm(data=data_frame_name, formula=as.formula(paste(fitness_trait_name,"~",predictor_for_single_regression)), family=fam)

        print("summary fit 2")
        print(summary(fit2))
        print("model comparison stats:")
        mod_test<-anova(fit2,fit1)

        ##suggestion #1
        print(anova(fit2,fit1, test="Chisq"))

        #suggestion #2
        print ("significance:")
    print (1-pchisq( abs(mod_test$Deviance[2]),df=abs(mod_test$Df[2])))

        }


data<-structure(list(ID = c(1L, 2L, 4L, 7L, 9L, 10L, 12L, 13L, 14L, 
15L, 16L, 17L, 18L, 20L, 21L, 22L, 23L, 24L, 25L, 27L, 28L, 29L, 
31L, 34L, 37L, 38L, 39L, 40L, 41L, 43L, 44L, 45L, 46L, 47L, 48L, 
49L, 52L, 55L, 56L, 59L, 60L, 61L, 62L, 63L, 65L, 66L, 67L, 68L, 
69L, 71L), QnWeight_initial = c(158L, 165L, 137L, 150L, 153L, 
137L, 158L, 163L, 159L, 151L, 145L, 144L, 157L, 144L, 133L, 148L, 
151L, 151L, 147L, 158L, 178L, 164L, 134L, 151L, 148L, 142L, 127L, 
179L, 162L, 150L, 151L, 153L, 163L, 155L, 163L, 170L, 149L, 165L, 
128L, 134L, 145L, 147L, 148L, 160L, 131L, 155L, 169L, 143L, 123L, 
151L), Survived_eclosion = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), Days_wrkr_eclosion_minus20 = c(NA, 
1L, NA, 3L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 1L, NA, 0L, 7L, 1L, 0L, 
1L, 0L, 1L, 2L, 2L, NA, 2L, 3L, 2L, 2L, NA, 0L, 1L, NA, NA, 0L, 
0L, 0L, 0L, 3L, 3L, 3L, 1L, 0L, 2L, NA, 1L, 0L, 1L, 1L, 3L, 1L, 
2L), MLH = c(0.5, 0.666666667, 0.555555556, 0.25, 1, 0.5, 0.333333333, 
0.7, 0.5, 0.7, 0.5, 0.666666667, 0.375, 0.4, 0.5, 0.333333333, 
0.4, 0.375, 0.3, 0.5, 0.3, 0.2, 0.4, 0.875, 0.6, 0.4, 0.222222222, 
0.222222222, 0.6, 0.6, 0.3, 0.4, 0.714285714, 0.4, 0.3, 0.6, 
0.4, 0.7, 0.625, 0.555555556, 0.25, 0.5, 0.5, 0.6, 0.25, 0.428571429, 
0.3, 0.25, 0.375, 0.555555556), Acon5 = c(0.35387674, 0.35387674, 
0.35387674, 0.35387674, 0.35387674, 0.35387674, 0.35387674, 0, 
0, 1, 0, 1, 0.35387674, 0, 0, 0.35387674, 1, 1, 0, 0, 0, 1, 0, 
0.35387674, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 
0, 0, 1, 0, 0, 0, 1, 0, 0.35387674), Baez = c(1, 1, 1, 0.467836257, 
1, 1, 0, 0, 1, 1, 0, 0.467836257, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 0, 0.467836257, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 
1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1), C294 = c(0, 1, 0, 0, 1, 
0.582542694, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 
0, 1, 1, 0, 0, 0.582542694, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1), C316 = c(1, 1, 0, 0, 0.519685039, 
0.519685039, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0.519685039, 0, 
1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0.519685039, 1, 0, 1, 
1, 0, 0.519685039, 1, 0.519685039, 1, 1, 1, 0.519685039, 0.519685039, 
0, 0.519685039, 0.519685039, 0), i_120_PigTail = c(1, 1, 0, 1, 
0.631236443, 0.631236443, 1, 1, 1, 1, 1, 0, 0.631236443, 1, 1, 
1, 0, 0.631236443, 1, 1, 1, 0, 0, 1, 1, 1, 0.631236443, 0, 1, 
1, 0, 1, 0.631236443, 1, 0, 1, 0, 0, 1, 0.631236443, 0.631236443, 
0, 1, 0, 0.631236443, 0.631236443, 1, 0.631236443, 0.631236443, 
1), i129 = c(0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L), Jackstraw_PigTail = c(0L, 1L, 1L, 0L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 
0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Neil_Young = c(0.529636711, 
0, 1, 0, 0.529636711, 0.529636711, 1, 1, 0, 1, 1, 1, 0, 0, 1, 
1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 
1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1), Ramble = c(0, 0, 0, 
0, 0.215163934, 0.215163934, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 
0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.215163934, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0.215163934, 0, 0, 0, 0), Sol_18 = c(1, 
0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0.404669261, 
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1)), .Names = c("ID", "QnWeight_initial", 
"Survived_eclosion", "Days_wrkr_eclosion_minus20", "MLH", "Acon5", 
"Baez", "C294", "C316", "i_120_PigTail", "i129", "Jackstraw_PigTail", 
"Neil_Young", "Ramble", "Sol_18"), class = "data.frame", row.names = c(NA, 
-50L))

make_and_compare_models("QnWeight_initial", data, c("Acon5","Baez","C294","C316","i_120_PigTail","i129","Jackstraw_PigTail","Neil_Young","Ramble","Sol_18"), "MLH", "gaussian")
3
除非嵌套模型具有相同的拟合效果，否则anova(fit1,fit2,test="Chisq")应该有效。是否能提供更多细节呢？ - Ben Bolker
还有，是否有一个类似于test="F"的模拟？我在手册中找不到关于anova()函数中test参数的任何信息... - Atticus29
(1) 是的，对于固定的卡方检验和估计的F检验。 (2) 您能提供一个可重现的模式示例吗？这将非常有帮助。一个令人困惑的“特性”是，如果偏差差异小于或大约为零（即模型基本相同），则R不会打印p值，而您真正看到的是偏差差异（例如，请参见@DWin的示例）。 (3) 对于“F”，“Chisq”等，请参见?anova.glm，该文档又引用了?stat.anova。 - Ben Bolker
好的@BenBolker，我已经添加了一个例子。希望这能为您提供所需的信息... - Atticus29
1
你的例子表明你正在比较非嵌套模型：df差异（在Df列中显示）为零！所有anova()框架（如下面的答案所讨论的）都是围绕嵌套模型构建的。如果你想比较非嵌套模型的拟合优度，可以使用AIC（谨慎使用）或Vuong测试... - Ben Bolker
显示剩余7条评论
2个回答
8

一个“较大”或更复杂的模型与嵌套或“简化”的模型之间的差异在两个模型的自由度差异下（渐近地）分布为卡方变量。因此，您需要提取偏差估计值和自由度差异，并将其与pchisq（偏差，diff（df））进行比较。 “p值”只是1减去该值。

> 1-pchisq(3.84,1)
[1] 0.05004352

如果您在glm帮助页面中运行第一个示例，然后添加一个没有“treatment”变量的简化模型，您将获得以下结果：

glm.D93.o <- glm(counts ~ outcome, family=poisson())
 anova.res <-anova(glm.D93, glm.D93.o)
 anova.res
#------------
Analysis of Deviance Table

Model 1: counts ~ outcome + treatment
Model 2: counts ~ outcome
  Resid. Df Resid. Dev Df    Deviance
1         4     5.1291               
2         6     5.1291 -2 -2.6645e-15
#---------------
 str(anova.res)
Classes ‘anova’ and 'data.frame':   2 obs. of  4 variables:
 $ Resid. Df : num  4 6
 $ Resid. Dev: num  5.13 5.13
 $ Df        : num  NA -2
 $ Deviance  : num  NA -2.66e-15
 - attr(*, "heading")= chr  "Analysis of Deviance Table\n" "Model 1: counts ~ outcome + treatment\nModel 2: counts ~ outcome"

所以在查看对象本身存储方式后，这为“outcome”给出了P值：

 1-pchisq( abs(anova.res$Deviance[2]), abs(anova.res$Df[2]))
[1] 1

这是在治疗+结果模型与仅治疗模型上的相应过程：

> glm.D93.t <- glm(counts ~ treatment, family=poisson())
> anova.res2 <-anova(glm.D93, glm.D93.t)
> 1-pchisq( abs(anova.res2$Deviance[2]), abs(anova.res2$Df[2]))
[1] 0.06547071
谢谢，DWin！这回答了我的问题！ - Atticus29
1-pchisq()可能不正确。我已经运行了完全混淆的数据模拟（即，两个模型之间不应该有显著差异，因为两个模型都没有成功地预测响应），并且报告的p值始终为“0”。您确定在这种情况下不只是pchisq()吗？ - Atticus29
我非常确定 1-pchisq(3.84,1) 返回 0.05。你需要确保在第一个参数中放入正确的偏差差值的绝对值，并在第二个参数中放入正确的自由度。模型参数的顺序将翻转 anova $Deviance 结果的符号，但是 abs() 应该可以解决这个问题。 - IRTFM
明白了。绝对值在那里。嗯...好的，我刚刚将第二个参数特别指定为“df=model$DF[2]”，问题就解决了。有趣... - Atticus29
这是一个有点病态的例子！由于我尚未理解的原因，treatment变量是多余的（具有零预测能力），所以即使请求了chisq p-value，R也不会打印它。 glm.D93.i <- glm（counts〜1，family = poisson）; anova（glm.D93.i，glm.D93.o，test =“Chisq”）更容易理解一些。 - Ben Bolker
1
如果您的2个模型是嵌套的，那么您可以使用这2个模型的偏差变化来查看包含额外参数的模型是否提供了更好的拟合。如果模型1包含k个参数，而模型2包含相同的k个参数加上额外的m个参数，则偏差变化遵循（近似）自由度为m的卡方分布。您可以使用此检验统计量来查看模型2是否比模型1更好。
如果您对此领域不熟悉，我强烈建议您阅读有关广义线性模型的入门文本。
那很完美，只是我不太确定如何实际实现它。也就是说，你知道它的R语法吗？ - Atticus29
很不幸，我已经好几年没有使用R了。据我回忆，glm.summary输出曾经提供了这个计算所需的一切信息。希望你能得到一个针对R的具体答案，而不仅仅是理论上的。 - mathematician1975
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接