# Comparison of values by group with covariance

I have 3 groups of different sample sizes with covariance I need to include. I wish to evaluate the value dependence of the groups compared to group 1 (group 2, group 3 values compared to group 1).

I can't do a simple ANCOVA since the data showed to be non-normal (shapiro test).

Sample data:

data <- structure(list(age = c(65.7, 65.7, 68.8, 68.8, 60.9, 60.9, 75,
75, 77, 77, 62.9, 62.9, 69.8, 69.8, 75, 73.3, 59.8, 59.8, 70.6,
70.6, 75.7, 75.7, 61.2, 61.2), value = c(94, 90, 113, 100, 103,
92, 70, 80, 86, 82, 101, 101, 90, 86, 94, 103, 96, 95, 90, 92,
84, 83, 84, 89), Group = structure(c(3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("1", "2", "3"), class = "factor")), row.names = c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 140L, 141L, 142L, 143L,
144L, 145L, 146L, 147L, 215L, 216L, 217L, 218L, 219L, 220L), class = "data.frame")

model1 <- lm(value~age*Group, data=data)
model2 <- lm(value~age+Group, data=data)


anova(model1, model2) results in Pr(>F)= 0.2168, but as age is significant I think model 1 is the right one..

Looking for possible solutions I wonder if I should use

1. Evaluation by linear regression, perhaps by center the age and:
 coeftest(model, vcov = vcovHC(model, type="HC1"))

1. emmeans_test() with covariate=age
2. emmeans and contrast (gives different results than emmeans_test):
emm <-  emmeans(model,~Group|age)

1. The age is very significant. Perhaps I should cut() and create age_groups factors and analyze them by that?

Which option would be the best in my case?

• There is no emmeans_test function, at least not in the emmeans package. So I am not sure what you did. If you ran emmeans, and then ran test on those results, then that tests each estimated marginal mean against zero, it does not test comparisons among the means. Jul 2, 2020 at 0:18
• @RussLenth, emmeans_test is a function from the rstatix package said to: "Perform pairwise comparisons between groups using the estimated marginal means. Pipe-friendly wrapper arround the functions emmans() + contrast from the emmeans package..." So should I stick to the evaluation by coeftest?
– YBB
Jul 2, 2020 at 0:49
• Well, as developer of emmeans, I know what it does and I stand by it. Jul 2, 2020 at 1:32
• I will look at the rstatix oackage, and see what they are trying to do and whether it's right. There are 20,000 R packages, so I suggest you state what package's functions you are asking about. Jul 2, 2020 at 1:35
• I looked at an example in rstatix and it seems consistent with what emmeans does. However, emm_test's output seems (to my eye) cluttered, and they show both unadjusted and Tukey'adjusted P values. In your contrast call, "tukey" as the second argument is just a synonym for "pairwise"; that's not specifying the adjustment method. I'd have used pairs(emm, adjust = "tukey") and pairs(emm, adjust = "none") Jul 2, 2020 at 2:13

First, coeftest() tests the regression coefficients, not the estimated marginal means. So it's apples and oranges.

In model1, you will get different comparisons of groups at different ages, so I suggest (a) plotting the results

> require(emmeans)
> emmip(model1, Group ~ age, cov.reduce = range)


... and (b) doing the comparisons at a few specified ages of interest.

> emm = emmeans(model1, ~ Group | age, at = list(age = c(60,70,80)))
> pairs(emm)
age = 60:
contrast estimate    SE df t.ratio p.value
1 - 2       -8.97  7.47 18 -1.200  0.4683
1 - 3      -14.29  7.72 18 -1.851  0.1818
2 - 3       -5.32  6.55 18 -0.813  0.6998

age = 70:
contrast estimate    SE df t.ratio p.value
1 - 2       -8.18  4.71 18 -1.739  0.2185
1 - 3       -3.48  4.21 18 -0.826  0.6923
2 - 3        4.71  4.20 18  1.121  0.5137

age = 80:
contrast estimate    SE df t.ratio p.value
1 - 2       -7.40 10.01 18 -0.739  0.7437
1 - 3        7.34  8.58 18  0.856  0.6741
2 - 3       14.74  8.99 18  1.639  0.2556

P value adjustment: tukey method for comparing a family of 3 estimates


The default P-value adjustment is Tukey, applied separately for each group.

I'd say the plot is essential. People tend to move way too fast to getting asterisk-decorated numbers, without thinking about what they actually have.

Another analysis worth considering is showing the slopes of those 3 lines, and comparing them:

> emtrends(model1, pairwise ~ Group, var = "age")
$emtrends Group age.trend SE df lower.CL upper.CL 1 -0.122 0.550 18 -1.28 1.034 2 -0.200 0.506 18 -1.26 0.863 3 -1.203 0.432 18 -2.11 -0.296 Confidence level used: 0.95$contrasts
contrast estimate    SE df t.ratio p.value
1 - 2      0.0786 0.747 18 0.105   0.9939
1 - 3      1.0818 0.699 18 1.547   0.2935
2 - 3      1.0032 0.665 18 1.508   0.3107

P value adjustment: tukey method for comparing a family of 3 estimates


This is a pretty small dataset, so it isn't surprising that you have subjectively large differences but no small P values. I hope that your actual dataset is larger.

• Thank you for your explanation! It is much clearer for me now. Unfortunately the same patterns you see in the small data portion show in the larger dataset... I am going back and forth between model1 and model2.. In case of model 2 should I just calculate emm = emmeans(model2, "Group") and 'pairs(emm)`?
– YBB
Jul 2, 2020 at 14:42
• The tests will be more significant with more data. LOOK AT THE PLOT -- are you sure you're comfortable with model2, which fits the same slope to all three? Remember, the goal of statistical analysis is to reasonably explain what's going on in the data. That P value of .2168 comparing the two models is smaller than all but one other P value we have seen so far. Jul 2, 2020 at 16:15