Sprawdzenie założenia proporcjonalności szans zachodzi w regresji logistycznej porządkowej za pomocą funkcji polr

Użyłem funkcji „polr” w pakiecie MASS do uruchomienia porządkowej regresji logistycznej dla porządkowej zmiennej jakościowej z 15 ciągłymi zmiennymi objaśniającymi.

Użyłem kodu (pokazanego poniżej), aby sprawdzić, czy mój model spełnia założenia proporcjonalnego prawdopodobieństwa zgodnie z poradami zawartymi w przewodniku UCLA . Jednak trochę martwię się o wynik sugerujący, że nie tylko współczynniki w różnych punktach odcięcia są podobne, ale są one dokładnie takie same (patrz rysunek poniżej).

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)

Zobacz podsumowanie modelu:

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

A teraz możemy spojrzeć na przedziały ufności dla oszacowań parametrów:

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

Ale te wyniki są nadal dość trudne do interpretacji, więc przekonwertujmy współczynniki na iloraz szans

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

Sprawdzanie założenia. Tak więc poniższy kod oszacuje wartości do wykreślenia. Najpierw pokazuje nam transformacje logitów prawdopodobieństwa, że ​​będą większe lub równe każdej wartości zmiennej docelowej

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

Tabela powyżej pokazuje (liniowe) przewidywane wartości, które otrzymalibyśmy, gdybyśmy regresowali naszą zmienną zależną od naszych zmiennych predykcyjnych pojedynczo, bez założenia równoległych nachyleń. Teraz możemy uruchomić serię binarnych regresji logistycznych z różnymi punktami odcięcia dla zmiennej zależnej, aby sprawdzić równość współczynników między punktami odcięcia

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

Przepraszam, że nie jestem ekspertem od statystyki i być może brakuje mi czegoś oczywistego. Spędziłem jednak dużo czasu próbując dowiedzieć się, czy jest problem z testowaniem założenia modelu, a także próbując znaleźć inne sposoby uruchamiania tego samego rodzaju modelu.

Na przykład czytałem w wielu listach pomocy, że inni używają funkcji vglm (w pakiecie VGAM) i funkcji lrm (w pakiecie rms) (na przykład patrz tutaj: Proporcjonalne założenie szans w regresji logistycznej reginalnej w R z pakietami VGAM i rms ). Próbowałem uruchomić te same modele, ale ciągle napotykam ostrzeżenia i błędy.

Na przykład, gdy próbuję dopasować model vglm do argumentu „parallel = FALSE” (ponieważ poprzednie odniesienie jest ważne przy testowaniu założenia proporcjonalności szans), napotykam następujący błąd:

Błąd w lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf in 'y'
Dodatkowo: Komunikat ostrzegawczy:
In Deviance.categorical.data.vgam (mu = mu, y = y, w = w, reszty = reszty,: dopasowane wartości bliskie 0 lub 1

Chciałbym zapytać, czy jest ktoś, kto mógłby zrozumieć i być w stanie wyjaśnić mi, dlaczego wykres, który utworzyłem powyżej, wygląda tak jak on. Jeśli rzeczywiście oznacza to, że coś jest nie tak, czy mógłbyś mi pomóc w znalezieniu sposobu na przetestowanie założenia proporcjonalności prawdopodobieństwa podczas korzystania z funkcji polr. A jeśli to po prostu niemożliwe, skorzystam z funkcji vglm, ale potrzebuję pomocy w wyjaśnieniu, dlaczego wciąż pojawia się błąd podany powyżej.

UWAGA: W tle znajduje się tutaj 1000 punktów danych, które w rzeczywistości są punktami lokalizacji w obszarze badań. Szukam, czy istnieją jakiekolwiek powiązania między kategoryczną zmienną odpowiedzi a tymi 15 zmiennymi objaśniającymi. Wszystkie te 15 zmiennych objaśniających to cechy przestrzenne (na przykład wysokość, współrzędne xy, bliskość lasu itp.). 1000 punktów danych przydzielono losowo za pomocą GIS, ale zastosowałem podejście do próbkowania warstwowego. Upewniłem się, że 125 punktów zostało losowo wybranych w ramach każdego z 8 różnych kategorycznych poziomów reakcji. Mam nadzieję, że te informacje również będą pomocne.

Zmienna zależna ma 8 uporządkowanych poziomów, dlatego na wykresie w celu przetestowania założenia proporcjonalności szans należy zobaczyć 8 różnych symboli dla każdej zmiennej niezależnej. Widoczne są tylko 2 symbole dla każdej zmiennej niezależnej, prawdopodobnie dlatego, że wybrałeś zbyt krótki przedział dla wartości osi x. Jeśli moje przypuszczenie jest prawidłowe, wystarczy użyć szerszego interwału dla wartości osi x. Wypróbuj ten kod:

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

Znalazłem to, przeglądając go i myślę, że odpowiedź z tego powodu może być przydatna. Myślę, że pomyłka jest

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

gdzie używasz FG1_val_catzamiast y. Korzystając z przykładu ze strategii modelowania regresji Harrella:

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