[R] Multiple comparisons for a two-factor ANCOVA

[Eric Scott]

I'm trying to do an ANCOVA with two factors (clipping treatment with two
levels, and plot with 4 levels) and a covariate (stem diameter). The
response variable is fruit number. The minimal adequate model looks like
this:

model3<-lm(fruit~clip + plot + st.dia + clip:plot)

I'd like to get some multiple comparisons like the ones from TukeyHSD, but
TukeyHSD doesn't work with the covariate. I've tried using glht() in the
multcomp package, but I'm not sure how to get it to give me the TukeyHSD for
all the interactions (i.e. clipped:plotA vs. unclipped:plotA, etc). It
seems as if I can only specify one factor or the other, and it gives me a
warning about there being interactions.

glht(model3, linfct=mcp(plot="Tukey"))

General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Linear Hypotheses:
Estimate
151Nat - 151Trans == 0 -1.924
HC - 151Trans == 0 -7.942
HE - 151Trans == 0 -4.637
HC - 151Nat == 0 -6.018
HE - 151Nat == 0 -2.712
HE - HC == 0 3.305

Warning message:
In mcp2matrix(model, linfct = linfct) :
covariate interactions found -- default contrast might be inappropriate

How can I get it to give me the same sort of output that TukeyHSD gives with
just a regular two-way ANOVA?

Thanks for your help

-Eric Scott

[[alternative HTML version deleted]]

______________________________________________
R-h...@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Dieter Menne]

Eric Scott-3 wrote:
>
> I'm trying to do an ANCOVA with two factors (clipping treatment with two
> levels, and plot with 4 levels) and a covariate (stem diameter). The
> response variable is fruit number. The minimal adequate model looks like
> this:
>
> model3<-lm(fruit~clip + plot + st.dia + clip:plot)
>
> I'd like to get some multiple comparisons like the ones from TukeyHSD,

> but...
>
>

To quote glht:
The mcp function must be used with care when defining parameters of interest
in two-way ANOVA or ANCOVA models. Here, the definition of treatment
differences (such as Tukey's all-pair comparisons or Dunnett's comparison
with a control) might be problem specific. Because it is impossible to
determine the parameters of interest automatically in this case, mcp in
multcomp version 1.0-0 and higher generates comparisons for the main effects
only, ignoring covariates and interactions (older versions automatically
averaged over interaction terms). A warning is given. We refer to Hsu
(1996), Chapter 7, and Searle (1971), Chapter 7.3, for further discussions
and examples on this issue.
----

Too bad there is no example (Hi, Torsten.... it would be nice).

A thread coming close is in:

http://markmail.org/message/rhdohxlrt3cpsdpx

Dieter

--
View this message in context: http://n4.nabble.com/Multiple-comparisons-for-a-two-factor-ANCOVA-tp1593039p1593082.html
Sent from the R help mailing list archive at Nabble.com.

[RICHARD M. HEIBERGER]

Please see the maiz example in ?MMC in the HH package.

maiz is the last example in the help file. Keep going all the way to the
end of
the help file. See also the

demo("MMC.WoodEnergy-aov", "HH")
These examples show how to use glht in the presence of interactions and
covariates.

Rich

[[alternative HTML version deleted]]

[Eric Scott]

Thanks for the example, but I'm still not sure from this example how to see
the pairwise comparisons for the interaction. For example, if I have two
factors, X and Y; and X has 2 levels, A and B, and Y has 3 levels, 1, 2, and
3, a TukeyHSD would give the following comparisons with p-values for each:

factor X
A-B

factor Y
1-2
1-3
2-3

X:Y
A:1 - A:2
A:1 - A:3
A:2 - A:3
A:1 - B:1
A:1 - B:2
A:1 - B:3
A:2 - B:1
........etc.

But if I add a covariate in addition to these two factors, TukeyHSD won't
work (because it doesn't recognize the covariate as a factor).

With the maiz example in ?MMC, it shows the differences between hybrids, but
doesn't show all the differences between nitrogen levels within and between
hybrids, unless I'm not seeing something (I must admit, the tie-breaker plot
was new to me and I'm still not entirely sure I understand it). glht and
glht.mmc both seem to only allow one factor in the
linfct=mcp(factor="Tukey") argument. For example, glht.mmc(maiz2.aov,
linfct=mcp(hibrido:nitrogeno="Tukey")) doesn't work.

I apologize if my confusion stems from a lack of knowledge of the
statistical underpinnings of an ANCOVA, but it seems like there should be a
way to get the same types of comparisons that TukeyHSD gives.

-Eric Scott

[RICHARD M. HEIBERGER]

In addtition to the example I mentioned previously,

demo("MMC.WoodEnergy-aov", "HH")

Please also see
demo("MMC.WoodEnergy", "HH")

In this example, since anova(energy.aov.4),
shows that the Wood factor and Stove:Wood interaction are significant,
all possible pairwise comparisons of the 12 Stove:Wood terms are not
interpretable. Only comparisons of Stoves within each of the Woods is
interpretable. These estimates are shown with both tables and graphs.
Since the covariate is also significant, it is necessary to pick a reference
value for the comparisons.

Here is a simplification of the WoodEnergy example to ignore the covariate.
The 66 pairwise comparisons that TukeyHSD provides for the interaction
effect are not interpretable. The significant interaction and one
significant
main effect together are an indicator that
main effects and interactions are not interpretable.
Only simple effects of one factor within
a constant level of the other factor are interpretable.

> energy.aov.4b <- aov(Energy ~ Stove*Wood + Stove:Wood,
+ data=energy)
> anova(energy.aov.4b)
Analysis of Variance Table
Response: Energy
Df Sum Sq Mean Sq F value Pr(>F)
Stove 2 0.007 0.003 0.0078 0.9923
Wood 3 274.768 91.589 209.0130 < 2.2e-16 ***
Stove:Wood 6 34.570 5.762 13.1483 3.781e-10 ***
Residuals 76 33.303 0.438
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> energy.aov.4b.HSD <- TukeyHSD(energy.aov.4b)
> sapply(energy.aov.4b.HSD, dim)
Stove Wood Stove:Wood
[1,] 3 6 66
[2,] 4 4 4
>

About a year after I wrote this example, Torsten extended glht to permit
an option of averaging over other factors and covariates. I need to revise
the WoodEnergy example to use that option.

[Eric Scott]

Thank you for your reply. The WoodEnergy example helped a lot. I
understand now that it is inappropriate to make all pairwise comparisons
with an interaction present and better to make comparisons between levels of
one factor within a constant level of the second factor. As I understand it,
the solution in the WoodEnergy example is to produce separate ANOVAs for
each type of wood and then perform the multiple comparisons between stove
types within each wood type. This makes a lot of sense. For my data, I'm
using glm.nb and if I run the model separately for each level of "site," it
estimates a different theta for each which I immagine is a problem. Is this
ok, or do I need to find a way to use the coefficients from the full model
with the interaction to compare levels of clipping within fixed levels of
site?

-Eric Scott

[RICHARD M. HEIBERGER]

On Wed, Apr 7, 2010 at 9:25 PM, Eric Scott <ersc...@illinois.edu> wrote:

> Thank you for your reply. The WoodEnergy example helped a lot. I
> understand now that it is inappropriate to make all pairwise comparisons
> with an interaction present and better to make comparisons between levels of
> one factor within a constant level of the second factor. As I understand it,
> the solution in the WoodEnergy example is to produce separate ANOVAs for
> each type of wood and then perform the multiple comparisons between stove
> types within each wood type. This makes a lot of sense. For my data, I'm
> using glm.nb and if I run the model separately for each level of "site," it
> estimates a different theta for each which I immagine is a problem. Is this
> ok, or do I need to find a way to use the coefficients from the full model
> with the interaction to compare levels of clipping within fixed levels of
> site?
>
> -Eric Scott
>
>

The "right" solution is to fit one model and then work with its
coefficients. For this example
the R glht function did not, at the time I wrote the example, have the
option of averaging over
the wood types. It now has "experimental" options for
interaction_average covariate_average
These usually, but not always, do the right thing.
For this example, I prefer the analysis in file HH/demo/MMC.WoodEnergy.s.R
in which one aov model is calculated and the adjustments are made for the
levels of Wood.
That file works in S-Plus, but not in R. As I noted before, I still need to
revise
the WoodEnergy example to use the experimental option in glht to duplicate
the results I
get from S-Plus.