# Questions about the glht function for planned comparison

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?).

lvl.treatments=table(treatments)
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)
# view default contrasts in options
#$contrasts # 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. - 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. – user18555 Jan 10 '13 at 22:18 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. –  MattBagg Jan 10 '13 at 23:27