使用 polr 函数检查序数逻辑回归中的比例优势假设是否成立

机器算法验证 r 物流 假设 有序的logit 波尔
2022-03-27 23:33:28

我使用 MASS 包中的“polr”函数对具有 15 个连续解释变量的有序分类响应变量运行有序逻辑回归。

我已经使用代码(如下所示)来检查我的模型是否符合UCLA 指南中提供的建议的比例赔率假设。但是,我有点担心输出意味着不仅各个切点之间的系数相似,而且它们完全相同(见下图)。

FGV1b <- data.frame(FG1_val_cat=factor(FGV1b[,"FG1_val_cat"]), 
                    scale(FGV1[,c("X","Y","Slope","Ele","Aspect","Prox_to_for_FG", 
                          "Prox_to_for_mL", "Prox_to_nat_border", "Prox_to_village", 
                          "Prox_to_roads", "Prox_to_rivers", "Prox_to_waterFG", 
                          "Prox_to_watermL", "Prox_to_core", "Prox_to_NR", "PCA1", 
                          "PCA2", "PCA3")]))
b     <- polr(FG1_val_cat ~ X + Y + Slope + Ele + Aspect + Prox_to_for_FG + 
                            Prox_to_for_mL + Prox_to_nat_border + Prox_to_village + 
                            Prox_to_roads + Prox_to_rivers + Prox_to_waterFG + 
                            Prox_to_watermL + Prox_to_core + Prox_to_NR, 
              data=FGV1b, Hess=TRUE)

查看模型摘要:

summary(b)
(ctableb <- coef(summary(b)))
q        <- pnorm(abs(ctableb[, "t value"]), lower.tail=FALSE) * 2
(ctableb <- cbind(ctableb, "p value"=q))

现在我们可以查看参数估计的置信区间:

(cib <- confint(b)) 
confint.default(b)

但是这些结果仍然很难解释,所以让我们将系数转换为优势比

exp(cbind(OR=coef(b), cib))

检查假设。所以下面的代码将估计要绘制的值。首先,它向我们展示了大于或等于目标变量每个值的概率的 logit 变换

FG1_val_cat <- as.numeric(FG1_val_cat)
sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}
(t <- with(FGV1b, summary(as.numeric(FG1_val_cat) ~ X + Y + Slope + Ele + Aspect + 
                             Prox_to_for_FG + Prox_to_for_mL + Prox_to_nat_border + 
                             Prox_to_village + Prox_to_roads + Prox_to_rivers + 
                             Prox_to_waterFG + Prox_to_watermL + Prox_to_core + 
                             Prox_to_NR, fun=sf)))

上表显示了如果我们在没有平行斜率假设的情况下一次将因变量回归到我们的预测变量上,我们将获得的(线性)预测值。所以现在,我们可以在因变量上运行一系列具有不同切点的二元逻辑回归,以检查切点之间系数的相等性

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,7:8]))

polr假设检查

抱歉,我不是统计专家,也许我在这里遗漏了一些明显的东西。但是,我花了很长时间试图弄清楚我如何测试模型假设是否存在问题,并试图找出运行相同类型模型的其他方法。

例如,我在许多帮助邮件列表中读到其他人使用 vglm 函数(在 VGAM 包中)和 lrm 函数(在 rms 包中)(例如,请参见此处: R 中序数逻辑回归中的比例赔率假设与包VGAM 和有效值)。我尝试运行相同的模型,但不断遇到警告和错误。

例如,当我尝试使用 'parallel=FALSE' 参数拟合 vglm 模型时(正如前面的链接提到的对于测试比例赔率假设很重要),我遇到以下错误:

lm.fit(X.vlm, y = z.vlm, ...) 中的错误:'y' 中的 NA/NaN/Inf
另外:警告消息:
在 Deviance.categorical.data.vgam(mu = mu, y = y,w = w,残差 = 残差,:拟合值接近 0 或 1

我想请问是否有人可能理解并能够向我解释为什么我上面制作的图表看起来如此。如果确实意味着某些事情不正确,请您帮我找到一种方法来测试仅使用 polr 函数时的比例赔率假设。或者,如果这是不可能的,那么我将尝试使用 vglm 函数,但是需要一些帮助来解释为什么我不断收到上面给出的错误。

注意:作为背景,这里有 1000 个数据点,它们实际上是整个研究区域的位置点。我正在寻找分类响应变量和这 15 个解释变量之间是否存在任何关系。所有这 15 个解释变量都是空间特征(例如,海拔、xy 坐标、与森林的接近度等)。这 1000 个数据点是使用 GIS 随机分配的,但我采用了分层抽样方法。我确保在 8 个不同的分类响应水平中随机选择 125 个点。我希望这些信息也有帮助。

2个回答

因变量有 8 个有序水平,因此在图表中测试比例优势假设您应该看到每个自变量有 8 个不同的符号。对于每个自变量,您只能看到 2 个符号,这可能是因为您为 x 轴的值选择了太短的间隔。如果我的猜想是正确的,您只需要为 x 轴的值使用更宽的间隔。试试这个代码:

par(mfrow=c(1,1))
plot(t, which=1:8, pch=1:8, xlab='logit', main=' ', xlim=range(s[,3:9]))

因此,我通过谷歌搜索发现了这一点,我认为出于这个原因,答案仍然可能有用。我认为错误在于

sf <- function(y) {
  c('VC>=1' = qlogis(mean(FG1_val_cat >= 1)),
    'VC>=2' = qlogis(mean(FG1_val_cat >= 2)),
    'VC>=3' = qlogis(mean(FG1_val_cat >= 3)),
    'VC>=4' = qlogis(mean(FG1_val_cat >= 4)),
    'VC>=5' = qlogis(mean(FG1_val_cat >= 5)),
    'VC>=6' = qlogis(mean(FG1_val_cat >= 6)),
    'VC>=7' = qlogis(mean(FG1_val_cat >= 7)),
    'VC>=8' = qlogis(mean(FG1_val_cat >= 8)))
}

你在哪里使用FG1_val_cat而不是y. 使用 Harrell 的回归建模策略中的示例:

library(Hmisc)
getHdata(support)
support <- support[complete.cases(support[, c("sfdm2", "adlsc", "sex", "age", "meanbp")]), ]
sfdm <- as.integer (support$sfdm2 ) - 1

sf1 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(sfdm >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(sfdm >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(sfdm >= 3))
  )
}

sf2 <- function (y) {
  c(' Y ≥ 1 ' = qlogis (mean(y >= 1)), 
    ' Y ≥ 2 ' = qlogis (mean(y >= 2)),
    ' Y ≥ 3 ' = qlogis (mean(y >= 3))
  )
}

s1 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf1,
              data = support)
s2 <- summary(sfdm ~ adlsc + sex + age + meanbp, fun=sf2,
              data = support)  

plot(s1, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

plot(s2, which =1:3, pch =1:3, xlab = ' logit ', main = ' ', width.factor = 1.4, cex.lab = 0.75)

在此处输入图像描述对比

在此处输入图像描述