Как я могу вычислить статистику теста Пирсона на отсутствие соответствия модели логистической регрессии в R?

10

Коэффициент отношения правдоподобия (он же отклонение) и критерий несоответствия (или качества соответствия) довольно просто получить для модели логистической регрессии (подгонка с использованием функции) в R. Однако это может быть легко подсчитать количество клеток в конечном итоге достаточно низко, чтобы тест был ненадежным. Один из способов проверить надежность теста отношения правдоподобия при отсутствии соответствия - это сравнить его статистику теста и P- значение с данными теста Пирсона хи-квадрат (или ) на отсутствие соответствия.G2glm(..., family = binomial)χ2

Ни glmобъект, ни его summary()метод не сообщают статистику теста для критерия хи-квадрат Пирсона на отсутствие соответствия. В моем поиске единственной вещью, которую я придумал, является chisq.test()функция (в statsпакете): ее документация гласит: « chisq.testвыполняет тесты таблицы непредвиденных обстоятельств по критерию хи-квадрат и проверки на соответствие». Тем не менее, документации о том, как выполнять такие тесты, немного:

Если xэто матрица с одной строкой или столбцом, или если xэто вектор и yон не задан, то выполняется проверка на соответствие критерия соответствия ( xрассматривается как одномерная таблица сопряженности). Записи xдолжны быть неотрицательными целыми числами. В этом случае проверяется гипотеза о том, равны ли вероятности совокупности с вероятностями pили все они равны, если pне даны.

Я полагаю, что вы можете использовать yкомпонент glmобъекта для xаргумента chisq.test. Однако вы не можете использовать fitted.valuesкомпонент glmобъекта для pаргумента chisq.test, потому что вы получите ошибку: " probabilities must sum to 1."

Как я могу (в R) по крайней мере рассчитать статистику теста Пирсона на отсутствие соответствия без необходимости выполнять шаги вручную?χ2

Firefeather
источник

Ответы:

13

Сумма квадратов остатков Пирсона точно равна статистике теста Пирсона на отсутствие соответствия. Поэтому, если вызывается ваша подобранная модель (т. Е. Объект) , следующий код вернет статистику теста:χ2glmlogistic.fit

sum(residuals(logistic.fit, type = "pearson")^2)

См. Документацию residuals.glmдля получения дополнительной информации, в том числе о том, какие другие остатки доступны. Например, код

sum(residuals(logistic.fit, type = "deviance")^2)

получит статистику теста , точно так же, как предоставляет.G2deviance(logistic.fit)

Firefeather
источник
Я согласен с Макро. Если вы хотели получить ответы от группы, вам следовало бы дождаться, чтобы услышать, что другие скажут в первую очередь. Все, что вы можете получить сейчас, пристрастно, увидев ваш ответ. Кроме того, если вы знаете ответ, что вы пытаетесь доказать, делая это?
Майкл Р. Черник
@Macro - Firefeather разместил на этом сайте четыре вопроса (включая этот) и сам ответил на три из них (включая этот), приняв один из своих ответов один раз. Еще несколько, как это, и я мог бы начать видеть образец!
jbowman
@jbowman, я могу представить себе ответ на свой вопрос, если вы выяснили его самостоятельно, прежде чем кто-либо другой оставил ответ, но firefeather разместил ответ буквально менее чем через 5 минут после публикации вопроса, кажется, ему / ей не нужна помощь Именно это заставило меня спросить, почему ... Я не совсем понимаю мотивацию ...
Макро
3
@Macro, перейдите по этой официальной ссылке: blog.stackoverflow.com/2011/07/… (ссылка на странице «Задать вопрос» в ярлыке флажка внизу: «Ответь на свой вопрос - поделись своими знаниями в стиле Q & A») «). У меня был этот вопрос, когда я выполнял домашнюю работу (я решил использовать R вместо Minitab, хотя Minitab был продемонстрирован в классе), но у меня не было достаточно времени, чтобы напечатать вопрос и дождаться ответа. Я нашел этот обходной путь и решил поделиться им с сообществом.
Firefeather
2
@Macro, пожалуйста! Я хотел бы задать больше вопросов, где я не даю ответ, и ответить на больше вопросов, которые я не задавал. Но jbowman прав в отношении шаблона: мой вклад в сообщество имеет тенденцию говорить с самим собой. :) (По крайней мере, я как-то
помогаю
10

Статистика Пирсона имеет вырожденное распределение, поэтому вообще не рекомендуется для соответствия логистической модели. Я предпочитаю структурированные тесты (линейность, аддитивность). Если вы хотите провести комплексный тест, посмотрите единую степень свободы теста невзвешенной суммы квадратов Ле Сесси - ван Хоуилинген - Копас - Хосмер, реализованную в функции rmsпакета R.residuals.lrm

Фрэнк Харрелл
источник
2
-1: Спасибо за понимание! Однако это не отвечает на мой вопрос. Поскольку это релевантный комментарий / обсуждение заявления, которое я сделал на заднем плане моего вопроса, ваш ответ, вероятно, принадлежит комментарию, а не ответу.
Firefeather
2
Я думаю, что четыре человека, которые проголосовали за мой ответ, не согласны с вами. И вы не имели дело с вырожденным распределением.
Фрэнк Харрелл
@FrankHarrell: Эта мера GOF отличается от теста GOF Хосмера-Лемешоу (HL)? Предположим, так из-за названия, а также сравнили два: Проведенный тест HL GOF, как найдено в ResourceSelectionпакете, и его результат отличается от того, что я получаю при запуске resid(lrm_object, 'gof')после подгонки моей модели логистической регрессии как lrm_object <- lrm(...). Если они действительно разные, можете ли вы прокомментировать, как тест HL соотносится с тем, который вы упомянули здесь? Спасибо!
Мег
1
Два совершенно разные. Статистика HL (в настоящее время устарела) имеет фиксированное значение df и обычно основана на децилях риска. Таким образом, статистика HL не вырождается как . С другой стороны, остерегайтесь любой статистики , где ФР продолжает расширяться с . N χ 2 Nχ2Nχ2N
Фрэнк Харрелл
Я хотел бы увидеть симуляцию, которая показывает это вырождение.
wdkrnls
0

Спасибо, я не осознавал, что это было так просто: сумма (остатки (f1, type = "pearson") ^ 2) Однако, обратите внимание, что остаток Пирсона варьируется в зависимости от того, рассчитывается ли он по ковариантной группе или по отдельным. Простой пример:

m1 - это матрица (это голова более крупной матрицы):

m1 [1: 4,1: 8]

    x1 x2 x3 obs    pi   lev  yhat y
obs  1  1 44   5 0.359 0.131 1.795 2
obs  0  1 43  27 0.176 0.053 4.752 4
obs  0  1 53  15 0.219 0.062 3.285 1
obs  0  1 33  22 0.140 0.069 3.080 3

Где x1-3 - предикторы, obs - нет. наблюдения в каждой группе, pi - вероятность членства в группе (предсказано по уравнению регрессии), lev - плечо, диагональ шляпной матрицы, yhat - предсказанное нет. (из у = 1) в группе и у фактического нет.

Это даст вам Пирсона группой. Обратите внимание, как это отличается, если y == 0: ' 'fun1 <- function(j){        if (m1[j,"y"] ==0){ # y=0 for this covariate pattern     Pr1 <- sqrt( m1[i,"pi"] / (1-m1[i,"pi"]))     Pr2 <- -sqrt (m1[i,"obs"]) res <- round( Pr1 * Pr2, 3) return(res) } else {  Pr1 <- m1[j,"y"] - m1[j,"yhat"] Pr2 <- sqrt(   m1[j,"yhat"] * ( 1-(m1[j,"pi"]) )   ) res <- round( Pr1/Pr2, 3) return(res) }    }

таким образом

nr <-nrow(m1)
nr1 <- seq(1,nr)
Pr <- sapply(1:nrow[m1], FUN=fun1)
PrSj <- sum(Pr^2)

Если имеется большое количество субъектов с ковариатными шаблонами y = 0, то остаток Пирона будет намного больше при расчете с использованием метода «по группам», а не «по индивидуальным».

См., Например, Hosmer & Lemeshow "Прикладная логистическая регрессия", Wiley, 200.

dardisco
источник
0

Вы также можете использовать, c_hat(mod)что даст тот же результат, что и sum(residuals(mod, type = "pearson")^2).

user54098
источник
1
В какой упаковке находится c_hat?
Firefeather