Проверка предположения о пропорциональных шансах выполняется в порядковой логистической регрессии с использованием функции polr

Я использовал функцию 'polr' в пакете MASS, чтобы запустить порядковую логистическую регрессию для порядковой категориальной переменной ответа с 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))

Проверка предположения. Таким образом, следующий код оценит значения, которые будут отображены. Сначала он показывает нам логит-преобразования вероятностей того, что они больше или равны каждому значению целевой переменной

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

Извиняюсь за то, что я не эксперт по статистике и, возможно, мне здесь не хватает чего-то очевидного. Однако я потратил много времени, пытаясь выяснить, есть ли проблема в том, как я проверял допущение модели, а также пытался выяснить другие способы запуска модели такого же типа.

Например, я прочитал во многих списках рассылки справки, что другие используют функцию vglm (в пакете VGAM) и функцию lrm (в пакете rms) (например, см. Здесь: Предположение пропорциональных шансов в порядковой логистической регрессии в R с пакетами VGAM и RMS ). Я пытался запустить те же модели, но постоянно сталкиваюсь с предупреждениями и ошибками.

Например, когда я пытаюсь согласовать модель vglm с аргументом «parallel = FALSE» (как упоминалось в предыдущей ссылке, важно для проверки предположения о пропорциональных коэффициентах), я сталкиваюсь со следующей ошибкой:

Ошибка в lm.fit (X.vlm, y = z.vlm, ...): NA / NaN / Inf in 'y'
Дополнительно: предупреждающее сообщение:
In Deviance.categorical.data.vgam (mu = mu, y = y, w = w, остатки = остатки: установленные значения, близкие к 0 или 1

Я хотел бы спросить, пожалуйста, есть ли кто-нибудь, кто мог бы понять и объяснить мне, почему график, который я создал выше, выглядит так, как он. Если это действительно означает, что что-то не так, не могли бы вы помочь мне найти способ проверить предположение о пропорциональных коэффициентах при использовании только функции polr. Или, если это просто невозможно, тогда я прибегну к попытке использовать функцию vglm, но тогда мне понадобится некоторая помощь, чтобы объяснить, почему я продолжаю получать ошибку, приведенную выше.

ПРИМЕЧАНИЕ. В качестве фона здесь есть 1000 точек данных, которые фактически являются точками расположения в изучаемой области. Я смотрю, есть ли какие-либо отношения между категориальной переменной ответа и этими 15 объясняющими переменными. Все эти 15 объясняющих переменных являются пространственными характеристиками (например, высота, координаты xy, близость к лесу и т. Д.). 1000 точек данных были распределены случайным образом с использованием ГИС, но я выбрал стратифицированную выборку. Я удостоверился, что 125 точек были случайно выбраны в каждом из 8 различных категорийных уровней ответа. Я надеюсь, что эта информация также полезна.

Зависимая переменная имеет 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]))

Итак, я нашел это через поиск в Google, и я думаю, что ответ по-прежнему может быть полезным по этой причине. Я думаю, что ошибка в

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. Используя пример из стратегий регрессионного моделирования Харрелла:

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)

против