Logistic regression fitting problem

Question

I'm trying to make a GLM with binomial distribution to evaluate the effect of the diameter at breast height and the occupation of tree termites in species of tree.

I'm using the following model:

model <- glm(ocupation ~ dap * specie, data = termites, family = binomial(link = logit))

My response variable is the occupation of termites on the tree (data of absence and presence); the independent variable is the dap and the interaction with kind of species.

But when I run the model the following errors appear:

Warning message:
glm.fit: odds adjusted numerically 0 or 1 occurred

Error in seq.default(1, length(levels(eval(as.name(qvar)))) - 1, by = 1) : 
  wrong sign in 'by' argument


Error in solve.default(t(X) %*% W %*% X) : 
  system is computationally singular: reciprocal condition number = 1.93217e-38
> 

Warning messages:
1: not plotting observations with leverage one:
  10, 24, 26, 36, 37, 40, 43, 57, 71, 72, 73, 116, 138 
2: not plotting observations with leverage one:
  10, 24, 26, 36, 37, 40, 43, 57, 71, 72, 73, 116, 138 
3: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
4: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced

The same thing happens when I rotate the contrasts in the analysis

Error in seq.default(1, length(levels(eval(as.name(qvar)))) - 1, by = 1) : 
  wrong sign in 'by' argument

Can someone help me? I understand very little about R.

Answer 1

Most likely you have a problem of perfect separation, meaning there are combinations of dap and specie that have all ocupation = 1 or all ocupation = 0. In the first case, the value of $$ \log \frac{p(\texttt{ocupation}=1)}{p(\texttt{ocupation}=0)}, $$ (the quantity that the linear part of your model is attempting to fit) that makes the data most likely is positive infinity (and in the second case, negative infinity). Since this is silly, the model has objected.

This problem occurs because although you know that the coefficients cannot really be large enough to predict such huge values for the logit, a model fit with maximum likelihood does not.

The second answer here lays out your alternatives.

So much for the statistics. Here are some substantive suggestions:

1. Can you get a more informative dependent variable? Are termites an all or nothing affair, or are some infestations larger than others? Thinking of this quantity as a magnitude rather than presence/absence often opens up other more stable model choices.
2. Perhaps one or other of specie and dap is actually the driving force. How does the model look if you remove the interaction between them? I would guess that it works better.
3. We don't know what dap actually is, but if it is some feature of the environment and the 'species packing' assumption is approximately true, then specie will have a preferred range of dap. And if that is true, an interaction will never work well in a binary data context because there will always be missing combinations of variables.

Answer 2

Just need to add a bit more to the answer by @conjugateprior than can fit into a comment easily.

The usual rule of thumb for logistic regression is that you need about 15 events per predictor variable that you are considering in the model. For this purpose, "events" are the least frequent of occupied/not occupied, so that's 45 events. Also, all except one of your species counts as a separate predictor, so in the interaction model with 28 different species and the continuous dap predictor you are trying to evaluate 28 individual predictors PLUS 27 interactions between dap and species, a total of 55 predictors. So you don't even have 1 event per variable. Regularization, as suggested by @conjugateprior in a comment, might help, but at this point it seems that you really don't have enough information to proceed this way. Think carefully about what you are trying to accomplish; might there be some way to combine species (e.g., into genus) in some way to cut your number of predictors down? (That combination of species into groups, however, would need to be done without reference to the occupancy data.)