Vérification de l'hypothèse de cotes proportionnelles dans une régression logistique ordinale à l'aide de la fonction polr

J'ai utilisé la fonction «polr» dans le package MASS pour exécuter une régression logistique ordinale pour une variable de réponse catégorielle ordinale avec 15 variables explicatives continues.

J'ai utilisé le code (illustré ci-dessous) pour vérifier que mon modèle répond à l'hypothèse de cotes proportionnelles en suivant les conseils fournis dans le guide de l'UCLA . Cependant, je suis un peu inquiet à propos de la sortie, ce qui implique que non seulement les coefficients à travers différents points de coupure sont similaires, mais qu'ils sont exactement les mêmes (voir le graphique ci-dessous).

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)

Voir un résumé du modèle:

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

Et maintenant, nous pouvons regarder les intervalles de confiance pour les estimations des paramètres:

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

Mais ces résultats sont encore assez difficiles à interpréter, alors convertissons les coefficients en odds ratios

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

Vérification de l'hypothèse. Ainsi, le code suivant estimera les valeurs à représenter graphiquement. Elle nous montre tout d'abord les transformations logit des probabilités d'être supérieures ou égales à chaque valeur de la variable cible

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

Le tableau ci-dessus affiche les valeurs prédites (linéaires) que nous obtiendrions si nous régressions notre variable dépendante sur nos variables prédictives une à la fois, sans l'hypothèse de pentes parallèles. Alors maintenant, nous pouvons exécuter une série de régressions logistiques binaires avec différents points de coupure sur la variable dépendante pour vérifier l'égalité des coefficients entre les points de coupure

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

Je m'excuse de ne pas être un expert en statistiques et peut-être qu'il me manque quelque chose d'évident ici. Cependant, j'ai passé beaucoup de temps à essayer de déterminer s'il y a un problème dans la façon dont j'ai testé l'hypothèse du modèle et à essayer de trouver d'autres façons d'exécuter le même type de modèle.

Par exemple, j'ai lu dans de nombreuses listes de diffusion d'aide que d'autres utilisent la fonction vglm (dans le package VGAM) et la fonction lrm (dans le package rms) (par exemple, voir ici: Hypothèse de cotes proportionnelles dans la régression logistique ordinale en R avec les packages VGAM et rms ). J'ai essayé d'exécuter les mêmes modèles, mais je suis constamment confronté à des avertissements et à des erreurs.

Par exemple, lorsque j'essaie d'adapter le modèle vglm avec l'argument 'parallel = FALSE' (comme le lien précédent le mentionne est important pour tester l'hypothèse de cotes proportionnelles), je rencontre l'erreur suivante:

Erreur dans lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf dans 'y'
En outre: Message d'avertissement:
dans Deviance.categorical.data.vgam (mu = mu, y = y, w = w, résidus = résidus,: valeurs ajustées proches de 0 ou 1

Je voudrais demander s'il y a quelqu'un qui pourrait comprendre et être en mesure de m'expliquer pourquoi le graphique que j'ai produit ci-dessus ressemble à cela. Si en effet cela signifie que quelque chose ne va pas, pourriez-vous s'il vous plaît m'aider à trouver un moyen de tester l'hypothèse de cotes proportionnelles lorsque vous utilisez simplement la fonction polr. Ou si cela n'est tout simplement pas possible, je vais essayer d'utiliser la fonction vglm, mais j'aurais alors besoin d'aide pour expliquer pourquoi je continue à obtenir l'erreur donnée ci-dessus.

REMARQUE: En arrière-plan, il y a 1000 points de données ici, qui sont en fait des points de localisation dans une zone d'étude. Je cherche à voir s'il existe des relations entre la variable de réponse catégorielle et ces 15 variables explicatives. Toutes ces 15 variables explicatives sont des caractéristiques spatiales (par exemple, l'altitude, les coordonnées xy, la proximité de la forêt, etc.). Les 1000 points de données ont été attribués au hasard à l'aide d'un SIG, mais j'ai pris une approche d'échantillonnage stratifié. Je me suis assuré que 125 points ont été choisis au hasard dans chacun des 8 différents niveaux de réponse catégorique. J'espère que ces informations sont également utiles.

La variable dépendante a 8 niveaux ordonnés donc dans le graphique pour tester l'hypothèse de cotes proportionnelles, vous devriez voir 8 symboles différents pour chaque variable indépendante. Vous ne voyez que 2 symboles pour chaque variable indépendante, probablement parce que vous avez choisi un intervalle trop court pour les valeurs de l'axe des x. Si ma conjecture est correcte, il vous suffit d'utiliser un intervalle plus large pour les valeurs de l'axe des x. Essayez ce code:

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

Donc, j'ai trouvé cela via googler et je pense qu'une réponse pourrait encore être utile pour cette raison. Je pense que l'erreur est

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

où vous utilisez FG1_val_catplutôt que y. En utilisant l'exemple des stratégies de modélisation de la régression 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)

contre