Solved – glmer with binary response variable: how to select fixed effects

glmmlme4-nlme

I have measured nest building (building / not building) 5 times over the breeding season and want to see if there is an effect of my treatment (treated / control).
After selecting random effect with likelihood ratio tests (anova(model1,model2) ), I have a model like this:

glmer(nest ~ treatment*time + (1|ID), family=binomial)

Now I want to select my fixed effects. According to Bolker et al. TREE 2009 I should use Wald Z or Chi squared tests to select my fixed effects. How do I do this? Anova{car} won't test 2 glmer's against each other. Can I use wald.test {aod}? If so, how?

Thanks for your help!

Laura

Best Answer

First of all: do you really need to do a model selection? MS with subsequent naive regression tests inflates type I error, so if you want to test for effects of your explanatory variables and you have enough data (no issue with power), there is no reason to run a MS.

If you must do a MS, I don't think anova(model1,model2) is appropriate to select on the random effects - likely this function will use naive degrees of freedom (1 or k, not sure) for the random effects. I would recommend fixing the random structure a priori and stick with it. If you must select, simulated LRTs would probably be safe.

Once the random effect structure is fixed, degrees of freedom are less problematic, so it should usually be OK to run a normal AICc basesd model selection on the fixed effects, e.g. with the MuMIn package - all the caveats that I noted first remain though.

Related Solutions

Model Selection in Longitudinal Data – Testing the Need for Random-Effects Terms in Longitudinal Data Analysis

The likelihood ratio test is slightly incorrect (in general, conservative) for testing the significance of a random effect, because the null value ($\sigma^2=0$) is at the boundary of the feasible space, but in this case there is overwhelmingly strong evidence against the null hypothesis. The model with random effects of individual is 15713-6772=8941 log-likelihood units better; twice the log-likelihood value is $\chi^2$ distributed, so the direct p-value calculation would give you ...

pchisq(2*8941,df=1,lower.tail=FALSE,log.p=TRUE)/log(10)
## -3885.251

... a p-value of approximately $10^{-3885}$.

You should really consider a random-slope model (random = ~time|id) as well.

Update: relative to the random-intercept model, the random-slopes model is again much better. The improvement is now 935 log-likelihood units, which doing the equivalent calculation as above corresponds to a rejection of the null hypothesis (among-individual variation in slope is equal to zero) with a p-value of "only" $10^{-408}$.

R – Addressing gls with ML Not Working for Best Model Selection

I initially used a GLM (because LMA data are neither normal nor heteroskedastic [maybe you meant "homoskedastic"?])

That might not be necessary. The distributional issues have to do with errors around the model predictions, not the raw data.

we supplemented LMA measurements with average values for some species

That's not typically a good way to deal with missing data. See Stef van Buuren's Flexible Imputation of Missing Data for the advantages of creating multiple imputed data sets so that you have the best chance of avoiding bias while including the variance arising from filling in missing data.

Example of anova results for 3 models. The problem is that depending on the models, the significant effects are different ...

Check the documentation for how anova() works with those model objects. It's possible that it uses a sequential Type I method so the order of variables in the model might affect the apparent "significance" of the individual predictors.

It's not clear why you are doing separate models for each of cat_TCT, specimen_TCT and combined_TCT. It's generally better to start with as many predictors in a model as reasonable without overfitting, even if some predictors are somewhat correlated. (If predictors are linearly dependent or close to that, you should choose a linearly-independent set.)

I don't know why you are having troubles with ML versus REML. You can't compare AIC values among models having different fixed effects that are fit via REML. Work through one step at a time to see the source of the error, instead of combining updates with anova().