4

I am trying to recreate a PROC GLIMMIX command in R using glmer. Here is a link to the data (from SAS product support GLIMMIX documentation): https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_glimmix_a0000001403.htm

The SAS code is

 proc glimmix data=multicenter;

      class center group;

      model sideeffect/n = group / solution;

      random intercept / subject=center;

   run;

The coefficients are:

Intercept: -0.8071

Group A : -0.4896

Group B: 0

Here is the R command, after swapping 1 and 0 in the sideeffect column to align the defaults in R and SAS:

mt <- glmer(sideeffect/n ~ group + (1|center), data = test, family = binomial, 
    weights=n, control = glmerControl(optimizer = "bobyqa"), nAGQ = 10)

The coefficients are:

Intercept: 1.3379

Group B: -0.4966

The different coefficients are not concerning on their face, since I transformed the data. However, when I compute the associated probabilities, I get from SAS:

P(sideeffect | A) = 0.2147

P(sideeffect | B) = 0.3085

and from R:

P(sideeffect | A) = 0.2078556

P(sideeffect | B) = 0.3018929

These estimates have a discrepancy of approximately 2%. I know that R and SAS use slightly different approaches to the generalized linear models - is this enough to explain this discrepancy? If not, what should I do to get R to conform to the SAS code? I have seen an example online using different data where R code exactly replicated the SAS results: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q3/004002.html

Silverfish
  • 20,678
  • 23
  • 92
  • 180
banwarp
  • 41
  • 1
  • 3

1 Answers1

3

The difference is because PROC GLIMMIX does not use maximum likelihood with adaptive quadrature as the default method. You need to specify the method that matches the method used by your R code:

proc glimmix data=multicenter method=quad(qpoints=10);
  class center group;
  model sideeffect/n = group / solution;
  random intercept / subject=center;
  lsmeans group / ilink;
run;

Then you get the estimates 0.20785716897455 and 0.30189193311442 which matches pretty close to what you report from R.

JimB
  • 2,043
  • 8
  • 14