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().)