2

我正在对以下数据框实施逻辑回归并得到合理（与使用 STATA 相同）的结果。但是我从 R 得到的 Pearson 卡方和自由度与 STATA 非常不同，这反过来又给了我一个非常小的 p 值。而且我无法获得 ROC 曲线下的面积。谁能帮我找出为什么残差（）不适用于具有先验权重的 glm（），以及如何处理 ROC 曲线下的面积？

以下是我的代码和输出。

1. 数据

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

 y     x1        x2     freq
 0     No        0      268
 0     No        1       14
 0    Yes        0      109
 0    Yes        1        1
 1     No        0       31
 1     No        1        6
 1    Yes        0       45
 1    Yes        1        6

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

2. GLM 模型

然后我做了逻辑回归：

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

输出显示：

偏差残差：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 ***

意义。代码：0''0.001''0.01''0.05 '。' 0.1''1

（二项式族的分散参数取为 1）

   Null deviance: 457.35  on 7  degrees of freedom

残余偏差：5 个自由度上的 416.96 AIC：422.96

Fisher 评分迭代次数：5

系数与 STATA 相同。

3. 卡方统计

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

pearson_chisq = sum(residuals(logit, type = "pearson", weights=test_data$freq)^2)

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

4. 歧视

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

logit_mf = model.frame(logit)

roc.area(logit_mf $y, 拟合(logit))$A

输出是：

[1] 0.5

Warning message:

在 wilcox.test.default(pred[obs == 1], pred[obs == 0], alternative = "great") 中：无法计算带关系的精确 p 值

谢谢！

4

1 回答 1

0

我最终想出了如何解决这个问题。我上面使用的数据集应该总结为协变量模式。然后使用皮尔逊卡方的定义进行计算。我提供的 R 代码如下：

# 提取协变量模式
库(dplyr)
temp=test_data %>% group_by(x1, x2) %>% summarise(m=sum(freq),y=sum(freq*y)) temp$pred=fitted(p01_logit_j)[ 1:4]

# 计算皮尔逊卡方
temp=mutate(temp, pearson=(ym pred)/sqrt(m pred*(1-pred)))
pearson_chi2 = with(temp, sum(pearson^2)) temp_dof = 4-(2+ 1) #自由度=J-(p+1)

# 计算 p 值
pchisq(pearson_chi2, temp_dof, lower.tail=F)

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

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

# 展开观察
y_expand=c(rep(0,392),rep(1,88))
logit_mf = model.frame(logit)
logit_pred = fit(logit)
logit_mf$freq=test_data$freq

# 扩展预测
yhat_expand=with(logit_mf, rep(pred, freq))
library(verification)
roc.area(y_expand, yhat)$A

AUC=0.6760，与 STATA 相同。

于 2014-11-17T02:26:54.613 回答