Generalized mixed effects model in R produces different results to SAS GLIMMIX

Viewed 23

I am trying to fit a mixed effects model with a logit link with the following variables:

response(Y) = ratio between (0, 1)

fixed effect(X) = categorical {"red", "blue"} (nominal, 2 levels)

random effect(M) = categorical {"a", ..., "g"} (nominal, 7 levels)

I want to produce R outputs that resemble SAS GLIMMIX results.

On SAS, I used GLIMMIX procedure to fit the model:

## SAS code ###
proc glimmix data = data; 
    class X M;
    model Y = X / link=logit D=binomial ddfm=kr ;
    random M;
    lsmeans X / pdiff ilink;
run;

I do not include an intercept in the model. The output shows F-statistic with DF = 1 for Type 3 test of fixed effect X.

Now, I'm trying to obtain a similar result using R. The following is my code:

### R code ###
library(MASS)     # for glmmPQL
library(nlme)     # for anova.lme

data$X <- as.factor(data$X)
data$M <- as.factor(data$M)

glmm_model <- glmmPQL(fixed = Y ~ X - 1,
                      random =~ 1|M, 
                      data = data,
                      family = binomial, 
                      na.action=na.exclude)

Then, the below code for summary gives me the following output:

summary(glmm_model)

## Linear mixed-effects model fit by maximum likelihood  
##  Data: data 
##  AIC BIC logLik 
##   NA  NA     NA 
##
## Random effects:  
##  Formula: ~1 | M  
##         (Intercept)  Residual  
## StdDev:   0.6305211 0.2050787  
##
##
## Variance function:  
##  Structure: fixed weights  
##  Formula: ~invwt  
## Fixed effects:  Y ~ X - 1  
##           Value  Std.Error DF   t-value p-value  
## XRed   -2.809994 0.3461666 30 -8.117461       0  
## XBlue  -3.277712 0.3528129 30 -9.290227       0  
##  Correlation:   
##                XRed  
## XBlue          0.489      
##
## Standardized Within-Group Residuals:  
##        Min         Q1        Med         Q3        Max   
## -1.2354719 -0.6104660 -0.2999117  0.2181443  2.8621013   
## 
## Number of Observations: 38  
## Number of Groups: 7  

Finally, type 3 test of the fixed effect using anova.lme() from nlme package gives me:

## Type-3 test of the fixed effect ##
anova.lme(glmm_model, type = "marginal", adjustSigma = FALSE)

##        numDF denDF  F-value   p-value   
##    X     2    30    54.40448  <.0001

As you can see, the above code gives me an F-statistic with numDF=2 instead of numDF=1 for X.

I tried using glmer() in lme4 package as well, but I found that glmmPQL() gives the outputs that's closest to what SAS gives me.

What modification should there be in my R code to get a Type-3 test result to SAS glimmix?

Thank you, and I hope you have a wonderful day :)

(edit: I put my R outputs for summary() and anova.lme().)

0 Answers
Related