Thanks for your responses. I knew that when you include an interaction term in a model you must include the main effects of each of the factors. Therefore, I assumed that SAS will do that by default. In most statistical packages, as in R, the main effects are automatically added when you include an interaction. But after reading your response I became suspicious and, it turns out, that SAS was indeed not including the main effects. To make sure that the main effects are part of the model in SAS you need to specify: sex eli sex*eli or use sex|eli. When I did this the estimates for the interaction terms now match closely those from R and, yes, the interaction between sex and infection status is now non significant.
Here is the relevant part of the new SAS output: Model Information Variance Matrix Blocked By Site Estimation Technique: Residual PL Degrees of Freedom Method: Containment Fit Statistics -2 Res Log Pseudo-Likelihood: 17868.73 Pseudo-AIC: 17890.73 Pseudo-BIC: 17957.14 Covariance Parameter Estimates Cov Parm Subject Estimate Std Error Intercept Site 0.2975 0.1799 Solutions for Fixed Effects Effect SEX ELI DW DIST SEA Estimate Std. Error DF t Value Pr > |t| Intercept -1.1427 0.4611 17 -2.48 0.0240 SEX 0 -0.6064 0.1673 3077 -3.62 0.0003 SEX 1 0 . . . . ELI 0 -0.1905 0.2197 3077 -0.87 0.3858 ELI 1 0 . . . . SEX*ELI 0 0 -0.4664 0.4617 3077 -1.01 0.3125 SEX*ELI 0 1 0 . . . . SEX*ELI 1 0 0 . . . . SEX*ELI 1 1 0 . . . . DW 0 -0.3308 0.1760 3077 -1.88 0.0603 DW 1 0 . . . . DIST 0 -0.1185 0.3816 3077 -0.31 0.7562 DIST 1 0 . . . . DW*DIST 0 0 -1.0148 0.4064 3077 -2.50 0.0126 DW*DIST 0 1 0 . . . . DW*DIST 1 0 0 . . . . DW*DIST 1 1 0 . . . . SEAS 0 -0.7839 0.1588 3077 -4.94 <.0001 SEAS 1 0 . . . . DEN -0.01343 0.002588 3077 -5.19 <.0001 WT 0.00758 0.01912 3077 0.40 0.6918 And here is again the R output using lmer fro easy comparison: Generalized linear mixed model fit using PQL Formula: SURV ~ SEX * ELI + DW * DIST + SEAS + DEN + WT + (1 | SITE) Family: binomial(logit link) AIC BIC logLik deviance 1539 1606 -758.7 1517 Random effects: Groups Name Variance Std.Dev. SITE (Intercept) 0.27816 0.52741 number of obs: 3104, groups: SITE, 19 Estimated scale (compare to 1 ) 0.9458749 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.144259 0.458672 -2.495 0.012606 SEX -0.606026 0.167289 -3.623 0.000292 *** ELI -0.190757 0.219599 -0.869 0.385034 DW -0.328796 0.175882 -1.869 0.061565 . DIST -0.117745 0.374148 -0.315 0.752989 SEAS -0.784971 0.158748 -4.945 7.62e-07 *** DEN -0.013381 0.002585 -5.176 2.27e-07 *** WT 0.007735 0.019115 0.405 0.685732 SEX:ELI -0.466425 0.461596 -1.010 0.312274 DW:DIST -1.015454 0.404683 -2.509 0.012099 * Now it is strange that after I fitted the model correctly the estimate for the random factor did not change, as well as the fit statistics. How can this be? Regarding the question of What is the Pseudo-Likelihood that is part of the output in SAS and not in R? The manual for the glimmix procedure (which can be found here http://support.sas.com/rnd/app/da/glimmix.html http://support.sas.com/rnd/app/da/glimmix.html ) says that "For a model containing random effects, the GLIMMIX procedure, by default, estimates the parameters by applying pseudo-likelihood techniques as in Wolfinger and O’Connell (1993) and Breslow and Clayton (1993)." That is, by default SAS uses the PQL method that as David Buffy said, it is just an approximation. The procedure involves a series of optimizations obtained via iterative estimation methods based on linearizations (using Taylor series expansions). After each optimization, a new pseudo-model is constructed for the mean response. All the fit statistics (AIC, BIC, etc) that SAS reports are calculated from the likelihood of the final "pseudo" model, thus the term "pseudo-likelihood." The problem is that then these criteria cannot be used to compare models; which is too bad because, at least in my field, it is preferable to use information theory than p-values based on arbitrary set significance levels (see Anderson Burnham and Thompson, 2000). What are the likelihood values that we obtain when using lmer in R? Can they be used to compare models? -- View this message in context: http://www.nabble.com/GLMMs-fitted-with-lmer-%28R%29---glimmix-%28SAS%29-tp14623190p14764018.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.