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
Best Answer
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:
Then you get the estimates 0.20785716897455 and 0.30189193311442 which matches pretty close to what you report from R.