Comparison of values by group with covariance

Question

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"))

2. emmeans_test() with covariate=age
3. emmeans and contrast (gives different results than emmeans_test):

emm <-  emmeans(model,~Group|age)
contrast(emm, method="tukey", adjust="none")

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

Answer 1

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.