Стандартные оценки ошибок R glm для SAS PROC GENMOD

Я преобразовываю пример SAS PROC GENMOD в R, используя glm в R. Код SAS был:

Когда я использую summary(parmsg2), я получаю те же оценки коэффициентов, что и в SAS, но мои стандартные ошибки сильно отличаются.

Важность различия в стандартных ошибках состоит в том, что коэффициенты SAS являются статистически значимыми, но коэффициенты RACE1 и WEEKEND в выводе R не являются. Я нашел формулу для вычисления доверительных интервалов Вальда в R, но это бессмысленно, учитывая разницу в стандартных ошибках, так как я не получу одинаковых результатов.

По-видимому, SAS использует стабилизированный хребтом алгоритм Ньютона-Рафсона для своих оценок, которые являются ML. Информация, которую я прочитал о функции glm в R, состоит в том, что результаты должны быть эквивалентны ML. Что я могу сделать, чтобы изменить свою процедуру оценки в R, чтобы получить эквивалентные коэффициенты и стандартные оценки ошибок, которые были созданы в SAS?

Чтобы обновить, благодаря ответу Spacedman, я использовал вес, потому что данные получены от людей в диетическом опросе, а REPLICATE_VAR - сбалансированный повторный вес репликации, который является целым числом (и довольно большим, в порядке 1000 или 10000 с). Веб-сайт, описывающий вес, здесь. Я не знаю, почему в SAS использовалась команда FREQ, а не команда WEIGHT. Теперь я проверю, расширив количество наблюдений, используя REPLICATE_VAR, и повторю анализ.

Спасибо Ben ответ ниже, код, который я использую сейчас:

Ответы

Ответ 1

изменить: чтение документации SAS для FREQ и ваших ответов выше и ниже, вот что, я думаю, вам следует попробуйте использовать weights=REPLICATE_VAR в инструкции glm для настройки относительного взвешивания групп (равенство найденных выше коэффициентов предполагает, что это правильный путь), затем используйте N=sum(REPLICATE_VAR) в настройке, предложенной ниже ( Я также думаю, что вы можете использовать lm вместо glm для этой проблемы... это не будет иметь большого значения, но должно быть немного быстрее и надежнее.) Что-то вроде:

s <- coef(summary(lm(y~x,data=data2, weights=REPLICATE_VAR)))
s[,"Std. Error"] <- s[,"Std. Error"]/sqrt(sum(data2$REPLICATE_VAR))
s[,"t value"] <- s[,"Estimate"]/s[,"Std. Error"]
s[,"Pr(>|t|)"] <- 2*pt(abs(s[,"t value"]),df=g$df.resid)

Ответ 2

FREQ в SAS не совпадает с весами в R glm. В SAS - количество вхождений этого события. Для R его "что каждый ответ y_i является средним для w_i единичных весовых наблюдений". Эти две вещи не совпадают.

Если вы хотите, чтобы R выдавал тот же результат, что и SAS (не могу понять почему), вам, вероятно, потребуется повторить каждую строку в вашем "фрейме" вашего количества раз.

Здесь данные 10 строк со всеми весами = 2, а data2 - 20 строк (2 копии каждой строки данных) со всеми весами = 1:

> summary(glm(y~x,data=data2,weights=weights))$coef
              Estimate Std. Error   t value   Pr(>|t|)
(Intercept) 0.32859847 0.13413683 2.4497259 0.02475748
x           0.01540002 0.02161811 0.7123667 0.48537003
> summary(glm(y~x,data=data,weights=weights))$coef
              Estimate Std. Error   t value  Pr(>|t|)
(Intercept) 0.32859847 0.20120525 1.6331506 0.1410799
x           0.01540002 0.03242716 0.4749111 0.6475449

Сжатие немного, N наблюдений с одним и тем же значением имеют меньшую нечеткость, чем говорят, что это наблюдение является средним из N наблюдений, поэтому SE с повторными наблюдениями будет иметь меньший SE, чем средний.