R: логистическая регрессия с использованием таблицы частот, не удается найти правильную статистику Пирсона хи-квадрат.

Я реализовал логистическую регрессию к следующему фрейму данных и получил разумные (такие же, как при использовании STATA) результаты. Но хи-квадрат Пирсона и степень свободы, которую я получил от R, сильно отличаются от STATA, который, в свою очередь, дал мне очень маленькое p-значение. И я не могу получить площадь под кривой ROC. Может ли кто-нибудь помочь мне выяснить, почему остаток () не работает в glm () с априорными весами и как работать с площадью под кривой ROC?

Ниже приведены мой код и результат.

1. Данные

Вот мой фрейм данных test_data, y - результат, x1 и x2 - ковариаты:

 y     x1        x2     freq
 0     No        0      268
 0     No        1       14
 0    Yes        0      109
 0    Yes        1        1
 1     No        0       31
 1     No        1        6
 1    Yes        0       45
 1    Yes        1        6

Я сгенерировал этот фрейм данных из исходных данных, подсчитав появление каждого шаблона ковариации, и сохранил число в новой переменной freq.

2. Модель GLM

Затем я сделал логистическую регрессию как:

logit=glm(y~x1+x2, data=test_data, family=binomial, weights=freq)

Вывод показывает:

Остаточные отклонения: 1 2 3 4 5 6 7 8
-7,501 -3,536 -8,818 -1,521 11,957 3,501 10,409 2,129

Коэффициенты: Оценка Станд. Ошибка значения z Pr (> | z |)
(Перехват) -2.2010 0,1892 -11,632 ‹2e-16 ***

x1 1.3538 0.2516 5.381 7.39e-08 ***

x2 1.6261 0.4313 3.770 0.000163 ***

Сигниф. коды: 0 '' 0,001 '' 0,01 '' 0,05 '.' 0,1 дюйма 1

(Параметр дисперсии для биномиального семейства принят равным 1)

   Null deviance: 457.35  on 7  degrees of freedom

Остаточное отклонение: 416,96 на 5 степенях свободы AIC: 422,96

Количество итераций подсчета очков Фишера: 5

Коэффициенты такие же, как у STATA.

3. Статистика хи-квадрат

когда я попытался вычислить хи-квадрат Пирсона:

pearson_chisq = сумма (остатки (logit, type = "pearson", weights = test_data $ freq) ^ 2)

Получил 488 вместо 1.3, выданных STATA. Также степень резкости, сгенерированная R, равна chisq_dof = df.residuals(logit) = 5 вместо 1. Таким образом, я получил чрезвычайно маленькое значение p ~ e ^ -100.

4. Дискриминация

Затем я рассчитал площадь под кривой ROC как: library(verification)

logit_mf = model.frame (логит)

roc.area (logit_mf $ y, подогнано (logit)) $ A

Результат:

[1] 0.5

Warning message:

В wilcox.test.default (pred [obs == 1], pred [obs == 0], alternate = "great"): невозможно вычислить точное значение p со связями

Спасибо!


r diagnostics glm roc pearson
person J.Liu    schedule 16.11.2014    source источник
comment
Я думаю, что это не столько вопрос кодирования, сколько просьба о статистическом обучении и интерпретации, и поэтому он принадлежит CrossValidated.com. Предупреждающее сообщение является стандартным предупреждением для тестов Вилкоксона. P-значения обычно асимптотически верны, и их обычно можно игнорировать. Можно ожидать большого количества связей с такими сгруппированными данными.   -  person IRTFM    schedule 16.11.2014
comment
Привет, BondedDush, спасибо за ответ. Я ожидал такого количества связей, но R, похоже, не способен обрабатывать логистическую регрессию с использованием таблицы частот в качестве STATA. Есть ли способ позволить R выполнять такую ​​диагностику, как STATA?   -  person J.Liu    schedule 17.11.2014
comment
Я понятия не имею, что вы имеете в виду под R, похоже, не способен обрабатывать логистическую регрессию с использованием таблицы частот. Вы уже сказали, что коэффициенты такие же. Если вы не знаете, как справляться с задачами вывода, это не вина R, а отражает ваши ограничения.   -  person IRTFM    schedule 17.11.2014
comment
Вы просто неверно указали свою модель glm. Вы можете указать это двумя разными способами. Здесь следует: total<-with(test_data,freq[1:4]+freq[5:8]); logit=glm(freq/total~x1+x2, data=test_data[1:4,], family=binomial, weights=total); chi_test<-sum(residuals(logit, type = "pearson")^2);   -  person J.R.    schedule 20.11.2014
comment
@ J.R. Это то, что я искал. Я не знал значения весов и должен был раньше извлекать ковариантные паттерны. Спасибо!   -  person J.Liu    schedule 21.11.2014
comment
Добро пожаловать! Меня удивило, что вы получили те же результаты, что и stata.   -  person J.R.    schedule 23.11.2014

Ответы (1)


arrow_upward
0
arrow_downward

Я со временем понял, как решить эту проблему. Набор данных, который я использовал выше, должен быть обобщен для ковариации паттернов. Затем используйте определение хи-квадрат Пирсона для вычисления. Я предоставляю код R следующим образом:

# извлечение ковариантных шаблонов
библиотека (dplyr)
temp = test_data%>% group_by (x1, x2)%>% summarize (m = sum (freq), y = sum (freq * y)) temp $ pred = подогнано (p01_logit_j) [1: 4]

# вычислить хи-квадрат Пирсона
temp = mutate (temp, pearson = (ym pred) / sqrt (m pred * (1-pred)))
pearson_chi2 = with (temp, sum (Пирсон ^ 2)) temp_dof = 4- (2 + 1) # dof = J- (p + 1)

# вычислить p-значение
pchisq (pearson_chi2, temp_dof, lower.tail = F)

Результат p-value равен 0,241941, что совпадает с STATA.

Чтобы вычислить AUC, мы должны сначала расширить шаблон ковариации до «исходных» данных, а затем использовать «расширенные» данные для получения AUC. Заметили, что в таблице частот 392 «0» и 88 «1». Мой код следующий:

# развернуть наблюдение
y_expand = c (rep (0,392), rep (1,88))
logit_mf = model.frame (logit)
logit_pred = fit (logit)
logit_mf $ freq = test_data $ freq

# расширение прогноза
yhat_expand = with (logit_mf, rep (pred, freq))
библиотека (проверка)
roc.area (y_expand, yhat) $ A

AUC = 0,6760, то же, что и для STATA.

person J.Liu    schedule 17.11.2014