0

I am running a mixed model on something akin to this data:

df<-data.frame(stage=c("a","a","a","a","b","b","b","b","c","c","c","c"),
nematode=c("fn","fn","bn","bn","fn","fn","bn","bn","fn","fn","bn","bn"),
id2=c(1,2,3,4,1,2,3,4,1,2,3,4),
value=c(1,0,0,2,3,1,1,2,0,0,0,2))

The model I am trying to fit is:

stage.id <- function(x) round(summary(glmer(value ~ stage + (1 | id2),family="poisson", data = x))$coefficients[2, c(1, 2, 4)], 3)
models.id0  <- ddply(tree2, .(stage, nematode), stage.id)

However, when I run this, I continually get an error:

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

which doesn't make sense to me given that I've used the nlevels() command on each of the factors (df$stage and df$nematode) and they are 3 and 2, respectively. Any sense what might be awry?

elduderino260
  • 223
  • 3
  • 12

1 Answers1

3

You have stage as a fixed effect in your model, but you're trying to fit the model for every combination of stage and nematode (ddply(tree2, .(stage, nematode), ...)). Thus in each chunk of the data there's only a single stage value, which causes the error.

You could:

  • apply only over nematode values, i.e. ddply(tree2,.(nematode), ...)
  • leave stage out of your model, i.e. fit the model value ~ 1 + (1 | id2)

According to your comment ("my goal is to compare the stages for each of the types of nematodes. I.e., for each type of nematode (e.g., bn) do the stages differ"), you want the former solution (apply only over nematode values).

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • I should say that my goal is to compare the stages for each of the types of nematodes. I.e., for each type of nematode (e.g., bn) do the stages differ. Either one of your suggestions seems to yield the same results. I believe they will allow me to determine whether the stages are significantly different for each type of nematode, but is there any way to do any sort of post-hoc comparison test to figure out which stages differ for each type of nematode? Thanks. – elduderino260 Dec 18 '16 at 00:19
  • Actually, I think I got it: bn<-subset(df,nematode=="bn") bn.mod<-glmer(value ~ stage + (1 | id2), family="poisson", data = bn) bn.tuk<-glht(bn.mod,mcp(stage="Tukey")) summary(bn.tuk) tuk.cld <- cld(bn.tuk) # letter-based display opar <- par(mai=c(1,1,1.5,1)) plot(tuk.cld) par(opar) – elduderino260 Dec 18 '16 at 00:54
  • But, if I do that then my comparison tests find significant differences in nematodes between stages for nematodes the mixed model did not find significantly different. I.e., the mixed model found that bn did not differ between the stages (p=0.34) whereas if I run the post hoc Tukey test, it finds significant differences (p<0.01)...That said, this is more of a question for Cross Validated... – elduderino260 Dec 18 '16 at 00:57
  • Maybe I am misinterpreting the mixed model results? – elduderino260 Dec 18 '16 at 01:04
  • 1
    maybe. We'd kind of need to see the full (even reproducible) example to know ... – Ben Bolker Dec 18 '16 at 01:09