1
votes

I'm trying to replicate the results of SAS's PROC GENMOD with glm in R. The model I'm trying to fit is

log[E(Yij|Yearij,Treati)]=Β1+B2Yearij+B3Treati*Yearij

In SAS, the code and result is:

proc sort data=skin; by id year;
run;
proc genmod data=skin;
class id yearcat;
model y=year trt*year / dist=poisson link=log type3 wald waldci;
repeated subject=id / withinsubject=yearcat type=un;
run;
-----------------------------------------------------------------------------------------------
Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Standard 95% Confidence
Parameter Estimate Error       Limits       Z    Pr > |Z|
Intercept -1.3341  0.0815 -1.4938 -1.1743 -16.37 <.0001
year      -0.0090  0.0271 -0.0622  0.0441 -0.33  0.7392
year*trt   0.0429  0.0319 -0.0195  0.1053  1.35  0.1781

As I want, there are only three coefficients estimated, for intercept, year, and year*treat.

In R, however, four coefficients are estimated, even though my model only specifies three:

> glm1<-glm(Y~year+treat*year,data=skin,family="poisson")
> summary(glm1)

Call:
glm(formula = Y ~ year + treat * year, family = "poisson", data = skin)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.34810    0.07647 -17.629   <2e-16 ***
year        -0.01192    0.02528  -0.472    0.637    
treat1       0.05850    0.10468   0.559    0.576    
year:treat1  0.03113    0.03454   0.901    0.367

Does anyone have a suggestion on how to specify my glm() command in R to obtain estimates of only year and year*treat, and not treatment alone?

3
Try glm1<-glm(Y~year+treat:year,data=skin,family="poisson"). The * notation adds the full and interaction terms, : only adds the interactionuser20650
That's it! Thank you so much!KES

3 Answers

2
votes

R interprets the formula A*B as A+B+A:B

try

glm1<-glm(Y~year+treat:year,data=skin,family="poisson")
summary(glm1)
1
votes

As user20650 pointed out, using a colon rather than an asterisk will estimate the interaction term only. Thus, the R code is

glm1<-glm(Y~year+treat:year,data=skin,family="poisson")
1
votes

I would use gee from library(gee) instead. The R code would be:

gee(Y ~ year + treat:year, data = skin, family = poisson, corstr = "unstructured", id = id)