Verificar a suposição proporcional de probabilidades é válida em uma regressão logística ordinal usando a função polr

Eu usei a função 'polr' no pacote MASS para executar uma regressão logística ordinal para uma variável de resposta categórica ordinal com 15 variáveis ​​explicativas contínuas.

Utilizei o código (mostrado abaixo) para verificar se meu modelo atende à premissa proporcional de chances, seguindo as orientações fornecidas no guia da UCLA . No entanto, estou um pouco preocupado com o resultado, sugerindo que não apenas os coeficientes em vários pontos de corte são semelhantes, mas eles são exatamente os mesmos (veja o gráfico abaixo).

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)

Veja um resumo do modelo:

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

E agora podemos observar os intervalos de confiança para as estimativas de parâmetros:

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

Mas esses resultados ainda são difíceis de interpretar, então vamos converter os coeficientes em odds ratio

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

Verificando a suposição. Portanto, o código a seguir estimará os valores a serem representados graficamente. Primeiro, ele mostra as transformações logit das probabilidades de ser maior ou igual a cada valor da variável de destino.

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

A tabela acima exibe os valores previstos (lineares) que obteríamos se regredíssemos nossa variável dependente em nossas variáveis ​​preditoras uma de cada vez, sem a hipótese de declives paralelos. Portanto, agora, podemos executar uma série de regressões logísticas binárias com pontos de corte variáveis ​​na variável dependente para verificar a igualdade de coeficientes entre os pontos de corte

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

Desculpas por não ser especialista em estatística e talvez esteja perdendo algo óbvio aqui. No entanto, passei muito tempo tentando descobrir se havia algum problema em como testei a suposição do modelo e também tentando descobrir outras maneiras de executar o mesmo tipo de modelo.

Por exemplo, eu li em muitas listas de ajuda que outros usam a função vglm (no pacote VGAM) e a função lrm (no pacote rms) (por exemplo, veja aqui: Pressuposto de probabilidades proporcionais na regressão logística ordinal em R com os pacotes VGAM e rms ). Eu tentei executar os mesmos modelos, mas estou continuamente enfrentando avisos e erros.

Por exemplo, quando tento ajustar o modelo vglm com o argumento 'parallel = FALSE' (como o link anterior menciona é importante para testar a suposição de chances proporcionais), encontro o seguinte erro:

Erro no lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf em 'y'
Além disso: Mensagem de aviso:
No Deviance.categorical.data.vgam (mu = mu, y = y, w = w, resíduos = resíduos,: valores ajustados próximos de 0 ou 1

Gostaria de perguntar, por favor, se existe alguém que possa entender e possa me explicar por que o gráfico que produzi acima se parece com ele. Se, de fato, significa que algo não está certo, você poderia me ajudar a encontrar uma maneira de testar a suposição de probabilidades proporcionais ao usar a função polr. Ou, se isso não for possível, tentarei usar a função vglm, mas precisarei de ajuda para explicar por que continuo recebendo o erro acima.

NOTA: Como pano de fundo, existem 1000 pontos de dados aqui, que na verdade são pontos de localização em uma área de estudo. Eu estou olhando para ver se há alguma relação entre a variável de resposta categórica e essas 15 variáveis ​​explicativas. Todas essas 15 variáveis ​​explicativas são características espaciais (por exemplo, elevação, coordenadas xy, proximidade com a floresta etc.). Os 1000 pontos de dados foram alocados aleatoriamente usando um SIG, mas adotei uma abordagem de amostragem estratificada. Assegurei-me de que 125 pontos fossem escolhidos aleatoriamente em cada um dos 8 níveis de resposta categórica diferentes. Espero que esta informação também seja útil.

A variável dependente possui 8 níveis ordenados, portanto, no gráfico, para testar a suposição de probabilidades proporcionais, você deverá ver 8 símbolos diferentes para cada variável independente. Você vê apenas 2 símbolos para cada variável independente, provavelmente porque escolheu um intervalo muito curto para os valores do eixo x. Se minha conjectura estiver correta, você só precisará usar um intervalo maior para os valores do eixo x. Tente este código:

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

Então, eu achei isso através do Google e acho que uma resposta ainda pode ser útil por esse motivo. Eu acho que o erro está em

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

onde você usa FG1_val_cate não y. Usando o exemplo das Estratégias de Modelagem de Regressão 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