There is a similar, if not duplicate, question from 9 years ago on this but it never got answered so I'm hoping to get some clarity on this.
Let's say we have N
treatments, such as medicines or fertilizers, that can be applied at varying intensities (continuous) and we want to fit an ANCOVA model to make inference about the difference in response to increasing intensity among the N
treatments. So we want to know if the slope is different between treatments. The issue is that we have some individuals with intensity 0, where it is not meaningful to assign those individuals to any one of the treatments.
I just want to confirm that the following approach in R is correct. I assigned the individuals with intensity 0 to a treatment of "none"
or something similar, so we now have N+1
unique values in the treatment column. Then fit the ANCOVA model as normal and also do post-hoc contrasts comparing the slopes between all pairs of treatments, but exclude the "none"
treatment from those contrasts. Does this properly take account of the information from the individuals where intensity is 0 and treatment is "none"
?
Example with two treatments
exdat <- structure(list(treatment = c("none", "none", "none", "none",
"none", "none", "trtA", "trtA", "trtA", "trtB", "trtB", "trtB",
"trtA", "trtA", "trtA", "trtB", "trtB", "trtB", "trtA", "trtA",
"trtA", "trtA", "trtB", "trtB", "trtB", "trtA", "trtA", "trtA",
"trtA", "trtB", "trtB", "trtB"), y = c(1.887069649, 2.2721258855,
1.9459101491, 1.8405496334, 2.2300144002, 2.2192034841, 2.9123506646,
2.8094026954, 1.974081026, 1.9169226122, 2.1041341543, 1.960094784,
2.7472709143, 2.9391619221, 2.6532419646, 2.4680995315, 2.0149030205,
2.7850112422, 2.2925347571, 2.3513752572, 3.1945831323, 2.6532419646,
3.3068867022, 2.9285235239, 2.5649493575, 2.9231615807, 3.0155349009,
2.9391619221, 2.7536607124, 2.9444389792, 2.8033603809, 2.9069010598
), intensity = c(0L, 0L, 0L, 0L, 0L, 0L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L,
8L, 9L, 11L, 12L, 13L, 15L, 16L, 17L, 19L, 20L, 21L, 23L, 24L,
25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L), class = "data.frame")
library(emmeans)
model <- lm(y ~ intensity + treatment + intensity:treatment, data = exdat)
summary(model)
emt <- emtrends(model, ~ treatment, var = 'intensity', at = list(treatment = c('trtA', 'trtB')))
contrast(emt, 'pairwise')
In this case because there are only two treatments other than the "none"
treatment, the contrast gives us the same t-ratio and p-value as we get from looking at the interaction coefficient in the model summary:
contrast estimate SE df t.ratio p.value
trtA - trtB -0.203 0.079 27 -2.567 0.0161