R glm стандартная ошибка оценки различий с SAS PROC GENMOD - PullRequest
Новая
       10

R glm стандартная ошибка оценки различий с SAS PROC GENMOD

10 голосов
/ 28 ноября 2011

Я преобразовываю пример 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?

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

Ответы [ 2 ]

12 голосов
/ 28 ноября 2011

ЧАСТОТА в 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

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

1 голос
/ 28 ноября 2011

edit : читая документацию 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)
...