Отличия оценок стандартной ошибки R glm от SAS PROC GENMOD
Я конвертирую пример SAS PROC GENMOD в R, используя glm в R. Код SAS был:
proc genmod data=data0 namelen=30;
model boxcoxy=boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3 + WEEKEND +
SEQ/dist=normal;
FREQ REPLICATE_VAR;
run;
Мой код R это:
parmsg2 <- glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3 + WEEKEND +
SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR)
Когда я использую summary (parmsg2) , я получаю те же оценки коэффициентов, что и в SAS, но мои стандартные ошибки сильно отличаются.
Итоговый результат SAS:
Name df Estimate StdErr LowerWaldCL UpperWaldCL ChiSq ProbChiSq
Intercept 1 6.5007436 .00078884 6.4991975 6.5022897 67911982 0
agegrp4 1 .64607262 .00105425 .64400633 .64813891 375556.79 0
agegrp5 1 .4191395 .00089722 .41738099 .42089802 218233.76 0
agegrp6 1 -.22518765 .00083118 -.22681672 -.22355857 73401.113 0
agegrp7 1 -1.7445189 .00087569 -1.7462352 -1.7428026 3968762.2 0
agegrp8 1 -2.2908855 .00109766 -2.2930369 -2.2887342 4355849.4 0
race1 1 -.13454883 .00080672 -.13612997 -.13296769 27817.29 0
race3 1 -.20607036 .00070966 -.20746127 -.20467944 84319.131 0
weekend 1 .0327884 .00044731 .0319117 .03366511 5373.1931 0
seq2 1 -.47509583 .00047337 -.47602363 -.47416804 1007291.3 0
Scale 1 2.9328613 .00015586 2.9325559 2.9331668 -127
Итоговый результат R:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.50074 0.10354 62.785 < 2e-16
AGEGRP4 0.64607 0.13838 4.669 3.07e-06
AGEGRP5 0.41914 0.11776 3.559 0.000374
AGEGRP6 -0.22519 0.10910 -2.064 0.039031
AGEGRP7 -1.74452 0.11494 -15.178 < 2e-16
AGEGRP8 -2.29089 0.14407 -15.901 < 2e-16
RACE1 -0.13455 0.10589 -1.271 0.203865
RACE3 -0.20607 0.09315 -2.212 0.026967
WEEKEND 0.03279 0.05871 0.558 0.576535
SEQ -0.47510 0.06213 -7.646 2.25e-14
Важность разницы в стандартных ошибках состоит в том, что все коэффициенты SAS статистически значимы, но RACE1 и WEEKEND коэффициентов в выходе R нет. Я нашел формулу для вычисления доверительных интервалов Вальда в R, но это бессмысленно, учитывая разницу в стандартных ошибках, так как я не получу тех же результатов.
Очевидно, SAS использует алгоритм Ньютона-Рафсона со стабилизацией гребня для своих оценок, которые являются ML. Информация, которую я прочитал о функции glm в R, заключается в том, что результаты должны быть эквивалентны ML.Что я могу сделать, чтобы изменить свою процедуру оценки в R, чтобы получить эквивалентные коэффициенты и оценки стандартной ошибки, которые были получены в SAS?
Чтобы обновить, благодаря ответу Спейсдмана, я использовал веса, потому что данные получены от людей в диетическое обследование, а REPLICATE_VAR - это сбалансированный вес повторной репликации, то есть целое число (и довольно большое, порядка 1000 или 10000). Веб-сайт с описанием веса находится здесь . Я не знаю, почему в SAS использовалась команда FREQ , а не WEIGHT . Теперь я протестирую, увеличив количество наблюдений с помощью REPLICATE_VAR и перезапустив анализ.
Благодаря ответу Бена ниже, я сейчас использую следующий код:
parmsg2 <- coef(summary(glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 + RACE3
+ WEEKEND + SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR)))
#clean up the standard errors
parmsg2[,"Std. Error"] <- parmsg2[,"Std. Error"]/sqrt(mean(data0$REPLICATE_VAR))
parmsg2[,"t value"] <- parmsg2[,"Estimate"]/parmsg2[,"Std. Error"]
#note: using the t-distribution for p-values, correct the t-values
allsummary <- summary.glm(glm(boxcoxxy ~ AGEGRP4 + AGEGRP5 + AGEGRP6 + AGEGRP7 + AGEGRP8 + RACE1 +
RACE3 + WEEKEND + SEQ , data=data0, family=gaussian, weights = REPLICATE_VAR))
parmsg2[,"Pr(>|t|)"] <- 2*pt(-abs(parmsg2[,"t value"]),df=allsummary$df.resid)