comparing mixed binomial model against the same model without random effect

[lincoln]

Hi everybody,

If I am correct, you can compare a model with random effect with the same model without the random effect by using the nlme function, like this:

no.random.model <- gls(Richness ~ NAP * fExp,
              method = "REML", data = RIKZ)
random.model <- lme(Richness ~NAP * fExp, data = RIKZ,
          random = ~1 | fBeach, method = "REML")
anova(no.random.model,random.model)

But, nlme is valid only for the gaussian family, isn't it? In my case I have a mixed model with binomial family, like this:

random.model <- lme(sex ~hwp+hcp, data = mydata,

          random = ~1 | colony, method = "REML")

where "sex" is a binary variable, "hwp" and "hcp" are continuous variable and "colony" is a factor with two levels.
I want to compare this model with another one without the random effect, I have tried with the lme4 but after this I cannot figure out how to build this same model without the random effect in order to make it comparable to the random effect model.

Thanks for any help.
Simone  

     
        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.

[lincoln]

Any answer to this?

I really need to compare a mixed model with binomial error against the same model without the random effect. I would use anova() but I don't know how to specify both models in order to make them comparable.

Thanks for any answer

Simone

[Billy.Requena]

Hello Simone,

Given that your response variable is binary and, consequently, you should
use generalized models, just occurs to me a comparison between a Generalized
Linear Model (the model without the random effect) and a Generalized Linear
Mixed Model (the model with the random effect).

You could write them as follows:

no.random.model <- glm(sex ~ hwp + hcp, data = mydata, family=binomial)
random.model <- glmer(sex ~ hwp + hcp + (1 | colony), data = mydata,
family=binomial)

My only doubt here is if one can directly compare both models, built under
different algorithms and from different packages.
Furthermore, I have no idea if the function 'anova' is able to compare
models produced by the 'glmer' function. One possibility is to compare them
based on Akaike Information Criteria (AIC) or any one of its corrected
versions.

I hope it helps you.
best wishes

--
Gustavo Requena
PhD student - Laboratory of Arthropod Behavior and Evolution
Universidade de São Paulo
Correspondence adress:
a/c Glauco Machado
Departamento de Ecologia - IBUSP
Rua do Matão - Travessa 14 no 321 Cidade Universitária, São Paulo - SP,
Brasil
CEP 05508-900
Phone number: 55 11 3091-7488

http://ecologia.ib.usp.br/opilio/gustavo.html

        [[alternative HTML version deleted]]


______________________________________________
[hidden email] 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.

[lincoln]

Thank you very much for answering,

I have just tried it and these are the results:

> random.model<-glmer(sex~hwp+hcp+(1|colony),family=binomial)
Mensajes de aviso perdidos
glm.fit: fitted probabilities numerically 0 or 1 occurred

> no.random.model<-glm(sex~hwp+hcp,family=binomial)
Mensajes de aviso perdidos
glm.fit: fitted probabilities numerically 0 or 1 occurred

> anova(no.random.model,random.model,test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: sex

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev P(>|Chi|)    
NULL                   425     581.51              
hwp   1   33.578       424     547.93 6.846e-09 ***
hcp   1  231.266       423     316.66 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> AIC(no.random.model,random.model)
                df      AIC
no.random.model  3 322.6621
random.model     4 324.5072


I believe that the warning message arises from the fact that males have almost all the values of "hcp" higher than zero and females tend to have "zero" for that variable.
The anova procedure to compare models doesn't seem to work I would like, in fact it seem that it is giving me the anova(model), i.e. the values of intercept, slopes and their p values.
Even though, AIC() gives me two different values, I guess I could use them to make this comparation.
I am worried about the algorithms beyond these two procedures (glm and glmer) because if they calculate the likelihood in a different way they would not be comparable neither the values of AIC.

Any other commentary/suggestion on this?
Thanks

Simone