I have a nested model with the following effects

  • fixed: treatments
  • random: experiment_date

I used lme() to model the data

mod1 <- lme(N_cells ~treatments-1, random=~1|experiment_date, method='ML')

Then I want to compare all the other treatments to the control (included in the "treatments" in mod1). After a fair amount of searching around, I decided to use glht() from the multcomp package (any other suggestions?).

K = contrMat(lvl.treatments,type='Dunnett',base=1)
mc<-glht(mod1, linfct=mcp(treatments=K),alternative='greater')

But I got the following error:

Error in contr.treatment(n = 0L) : not enough degrees of freedom to define contrasts

I tried to extract the df parameter using modelparm(), but the function couldn't be applied to lme

Error in UseMethod("modelparm") : no applicable method for 'modelparm' applied to an object of class "lme"

The degree of freedom of the fixed effect was 194. I tried to specify the number in glht(), but got the same error as "not enough degrees of freedom to define contrasts".

Does anyone know what's happening and how I could possibly solve the problem? Thank you so much.


If your treatments are factors (and not ordered factors), you could add the intercept into the model (i.e. remove the "-1") and just do summary(mod1)

The default contrasts, as set in options, is to use contr.treatment for factors. This sounds like what you want. contr.treatment means that each coefficient represents a comparison of that level with level 1 (omitting level 1 itself).

# view default contrasts in options
#        unordered           ordered 
#"contr.treatment"      "contr.poly" 

When you do summary(mod1), the first level will not be labelled, but all the other levels will be in comparison to it.

If your control condition is not your first level, you need to use factor() with a levels argument or relevel() to make it first.

  • $\begingroup$ Thank you, MattBagg. This makes a lot sense now. I was wondering if there is a way to specify the alternative hypothesis, e.g., my testing treatment is lower than the control. I have looked through the help documents of contr.treatment etc, but I did not find anything useful. $\endgroup$ – user18555 Jan 10 '13 at 22:18
  • $\begingroup$ I would probably just extract the t- or p- value and calculate that myself. For example, summary(mod1)$tTable[,"p-value"] The t-value would remain the same and the one-tailed p-value is either 0.5*(two-tailed p-value) or 1-0.5*(two-tailed p-value) depending on the direction. $\endgroup$ – MattBagg Jan 10 '13 at 23:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.