首页 文章

R：使用频率表进行逻辑回归，找不到正确的Pearson Chi Square统计量

提问于
浏览 次
1

我对以下数据框实现了逻辑回归，得到了一个合理的（与使用STATA相同）结果 . 但是我从R得到的Pearson卡方和自由度与STATA非常不同，后者反过来给了我一个非常小的p值 . 而且我无法将该区域置于ROC曲线下 . 任何人都可以帮我找出为什么residual（）不能用于带有先验权重的glm（），以及如何处理ROC曲线下的面积？

Following is my code and output.

1. Data

这是我的数据框 test_data ，y是结果，x1和x2是协变量：

y x1 x2 freq
0否0 268
0否1 14
0是0 109
0是1 1
1否0 31
1否1 6
1是0 45
1是1 6

我通过计算每个协变量模式的出现次数从原始数据生成此数据帧，并将该数字存储在新变量 freq 中 .

2. GLM Model

然后我做了逻辑回归：

logit=glm(y~x1+x2, data=test_data, family=binomial, weights=freq)

输出显示：

Deviance Residuals：1 2 3 4 5 6 7 8 -7.501 -3.536 -8.818 -1.521 11.957 3.501 10.409 2.129系数：估计标准 . 误差z值Pr（> | z |）（截距）-2.2010 0.1892 -11.632 <2e-16 *** x1 1.3538 0.2516 5.381 7.39e-08 *** x2 1.6261 0.4313 3.770 0.000163 *** Signif . 代码：0''0.001''0.01''0.05' . ' 0.1''1（二项式族的色散参数为1）空偏差：7自由度为457.35
剩余偏差：5自由度上的416.96 AIC：422.96 Fisher评分迭代次数：5

系数与STATA相同 .

3. Chi Square Statistics

当我试图计算皮尔逊卡方时：

pearson_chisq = sum（残差（logit，type =“pearson”，weights = test_data $ freq）^ 2）

我得到488，而不是STATA给出的1.3 . 由R生成的DOF也是 chisq_dof = df.residuals(logit) = 5，而不是1.所以我得到了一个非常小的p值~e ^ -100 .

4. Discrimination

然后我计算了ROC曲线下的面积： library(verification)

logit_mf = model.frame（logit）roc.area（logit_mf $ y，fit（logit））$ A

输出是：

[1] 0.5

警告信息：
在wilcox.test.default（pred [obs == 1]，pred [obs == 0]，alternative =“great”）：无法用关系计算精确的p值

谢谢！

r glm diagnostics roc pearson

1 回答

  • 0

    我想出了最终如何解决这个问题 . 我上面使用的数据集应该归纳为协变量模式 . 然后使用Pearson chi square的定义进行计算 . 我提供的R代码如下：

    ＃extract covariate patterns library（dplyr）temp = test_data％>％group_by（x1，x2）％>％summary（m = sum（freq），y = sum（freq * y））temp $ pred = fits（p01_logit_j） [1：4]＃calculate Pearson chi square temp = mutate（temp，pearson =（y-mpred）/ sqrt（mpred *（1-pred）））pearson_chi2 = with（temp，sum（pearson ^ 2））temp_dof = 4-（2 1）#dof = J-（p 1）＃计算p值pchisq（pearson_chi2，temp_dof，lower.tail = F）

    p值的结果为0.241941，与STATA相同 .

    为了计算AUC，我们首先应该将协变量模式扩展为“原始”数据，然后使用“扩展”数据来获得AUC . 注意，我们在频率表中有392“0”和88“1” . 我的代码如下：

    #expansion observation y_expand = c（rep（0,392），rep（1,88））logit_mf = model.frame（logit）logit_pred = fits（logit）logit_mf $ freq = test_data $ freq #expansion forecast yhat_expand = with（logit_mf， rep（pred，freq））library（verification）roc.area（y_expand，yhat）$ A

    AUC = 0.6760，与STATA相同 .

    回复于

相关问题