La comprobación del supuesto de probabilidades proporcionales se cumple en una regresión logística ordinal utilizando la función polr

He usado la función 'polr' en el paquete MASS para ejecutar una regresión logística ordinal para una variable de respuesta categórica ordinal con 15 variables explicativas continuas.

He utilizado el código (que se muestra a continuación) para verificar que mi modelo cumpla con el supuesto de probabilidades proporcionales siguiendo los consejos proporcionados en la guía de UCLA . Sin embargo, estoy un poco preocupado porque la salida implica que no solo los coeficientes en varios puntos de corte son similares, sino que son exactamente los mismos (vea el gráfico a continuación).

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)

Ver un resumen del modelo:

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

Y ahora podemos ver los intervalos de confianza para las estimaciones de parámetros:

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

Pero estos resultados aún son bastante difíciles de interpretar, así que vamos a convertir los coeficientes en odds ratios

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

Comprobando el supuesto. Entonces, el siguiente código estimará los valores a graficar. Primero nos muestra las transformaciones logit de las probabilidades de ser mayor o igual a cada valor de la variable objetivo

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)))

La tabla anterior muestra los valores pronosticados (lineales) que obtendríamos si retrocediéramos nuestra variable dependiente en nuestras variables predictoras de una en una, sin el supuesto de pendientes paralelas. Entonces, ahora podemos ejecutar una serie de regresiones logísticas binarias con puntos de corte variables en la variable dependiente para verificar la igualdad de coeficientes entre los puntos de corte

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

Disculpas por no ser un experto en estadísticas y tal vez me estoy perdiendo algo obvio aquí. Sin embargo, he pasado mucho tiempo tratando de descubrir si hay un problema en cómo probé el supuesto del modelo y también tratando de descubrir otras formas de ejecutar el mismo tipo de modelo.

Por ejemplo, leí en muchas listas de correo de ayuda que otros usan la función vglm (en el paquete VGAM) y la función lrm (en el paquete rms) (por ejemplo, vea aquí: Suposición de probabilidades proporcionales en regresión logística ordinal en R con los paquetes VGAM y rms ). He intentado ejecutar los mismos modelos, pero continuamente me encuentro con advertencias y errores.

Por ejemplo, cuando trato de ajustar el modelo vglm con el argumento 'paralelo = FALSO' (como menciona el enlace anterior es importante para probar el supuesto de probabilidades proporcionales), me encuentro con el siguiente error:

Error en lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf en 'y'
Además: Mensaje de advertencia:
En Deviance.categorical.data.vgam (mu = mu, y = y, w = w, residuales = residuales,: valores ajustados cercanos a 0 o 1

Me gustaría preguntar por favor si hay alguien que pueda entender y poder explicarme por qué el gráfico que produje arriba se ve como se ve. Si de hecho significa que algo no está bien, ¿podría ayudarme a encontrar una manera de probar la suposición de probabilidades proporcionales cuando solo use la función polr? O si eso simplemente no es posible, entonces recurriré a tratar de usar la función vglm, pero luego necesitaré ayuda para explicar por qué sigo recibiendo el error anterior.

NOTA: Como fondo, hay 1000 puntos de datos aquí, que en realidad son puntos de ubicación en un área de estudio. Estoy buscando ver si hay alguna relación entre la variable de respuesta categórica y estas 15 variables explicativas. Todas esas 15 variables explicativas son características espaciales (por ejemplo, elevación, coordenadas xy, proximidad al bosque, etc.). Los 1000 puntos de datos se asignaron al azar utilizando un SIG, pero tomé un enfoque de muestreo estratificado. Me aseguré de que se eligieran aleatoriamente 125 puntos dentro de cada uno de los 8 niveles de respuesta categórica diferentes. Espero que esta información también sea útil.

La variable dependiente tiene 8 niveles ordenados, por lo que en el gráfico para probar el supuesto de probabilidades proporcionales debería ver 8 símbolos diferentes para cada variable independiente. Verá solo 2 símbolos para cada variable independiente probablemente porque ha elegido un intervalo demasiado corto para los valores del eje x. Si mi conjetura es correcta, solo necesita usar un intervalo más amplio para los valores del eje x. Prueba este código:

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

Entonces, encontré esto buscando en Google y creo que una respuesta aún podría ser útil por esa razón. Creo que el error está en

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)))
}

donde usas en FG1_val_catlugar de y. Usando el ejemplo de las estrategias de modelado de regresión de 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)

vs