[R-sig-ME] contrasts and factors with 2 or more levels

[Gonçalo M. Rosa]

Hi team,

I hope you’re doing well. I’ve been struggling with fitting a model for my
experiment over the past few weeks, and I was hoping to get your advice.

We’re testing the effect of NE over time using a complex experimental
design with several treatment combinations. Here’s a summary of the setup:

   1. *Phase 1*: Each subject was randomly submitted to one of three
   initial treatments involving exposure to NE: NH, NE+, or NE-.
   2. *Phase 2*: After the initial treatment, subjects were randomly
   re-assigned to one of the 3 treatments. So we ended up with the following
   six possible treatment sequences:
      - NE- > NE-
      - NE- > NE+
      - NE- > NH
      - NE+ > NE+
      - NE+ > NE-
      - NE+ > NH

As response variables, we measure:

- body.index (once per phase)
- feeding.rate (3 times each phase)
- cover.use (7 times each phase)

Here a few summaries:

   1. > sapply(DF_Am[c("id", "treatment", "phase", "day", "body.index",
   "cover.use", "feeding.rate")], class)
   2.           id    treatment        phase          day   body.index
   cover.use feeding.rate
   3.     "factor"     "factor"     "factor"    "numeric"    "numeric"
   "numeric"    "numeric"
   4.
   5.
   6. > table(DF_Am$treatment, DF_Am$phase)
   7.
   8.                 phase1     phase2
   9.   NE- > NH         195        195
   10.   NE+ > NH         182        182
   11.   NE- > NE-        195        195
   12.   NE- > NE+        195        195
   13.   NE+ > NE-        195        195
   14.   NE+ > NE+        195        195
   15.
   16.
   17. > sapply(lapply(DF_Am, unique), length)
   18.                    id             treatment
   phase                  day
   19.                   105                     6
   2                   14
   20.           body.weight                   svl
   body.index             cover.use
   21.                   108                    12
   155                     3
   22.               feeding          feeding.rate
   23.                     7                     7

(all variables have at least 2 levels).

I wanted to model the effect of treatment and phase (i.e., treatment *
phase) on feeding.rate (proportion between 0 and 1) and cover.use (binary:
0 or 1), while considering the effect of the body.index as well. So, I’ve
started with the following model:

   1. model1 <- glmmTMB(cover.use ~ treatment * phase + body.index +
   2.                              (1 | id) + (1 | phase/day),
   3.                            family = binomial(link = "logit"),
   4.                            data = DF_Am)

The model runs smoothly and overall not showing any major issues (no
collinearity, no issues of dispersion, residuals, etc). The only red flag
is the binned_residuals(model1), where it says "Only about 79% of the
residuals are inside the error bounds".

I then tried to replicate the same structure for the feeding rate:

   1. model2 <- glmmTMB(feeding.rate ~ treatment * phase + body.index +
   2.                     (1 | id) + (1 | phase/day),
   3.                   family = binomial(link = "logit"),
   4.                   data = DF_Am)

But I always get the same error:

   1. Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
   2.   contrasts can be applied only to factors with 2 or more levels

I tried using counts instead of a rate, so changing the model to a Poisson
family, but made no difference.

I tried simplifying the random factors or even having the body.index as a
random factor... but that error is always there. The only way I found to
run the model is without the inclusion of body.index at all. I would
discharge it, but it seems too important of a variable given that we expect
a higher body.index to influence the feeding rate positively.


I thought that the fact that we only have one body.index data point per
subject (id) per phase would be the driver of this error, but why is that
not an issue in the first model? What am I missing here?


I’d really appreciate any guidance you can offer on how to address this
issue :) I've spent two weeks running in circles and I'm out of options...

Thank you so much for your time and help!

Gonçalo M. Rosa
IMIB Biodiversity Research Institute (CSIC, UO, PA), Spain
Institute of Zoology, Zoological Society of London, UK
Centre for Ecology, Evolution and Environmental Changes (CE3C), Faculdade
de Ciências UL, Portugal
IUCN SSC Amphibian Specialist Group
ᐧ

	[[alternative HTML version deleted]]