I am analysing my own data with a linear mixed model and estimated marginal means and have made a strange observation. I have run an experimental study with three treatment groups, measuring a dependent variable at three points. I also measured a continuous covariate. The structure of the data is as follows:
library(tidyverse)
library(afex)
library(emmeans)
df <- tibble(
id = rep(1:3, 100),
treatment = rep(1:3, 100),
t1 = runif(300, min = 1, max = 4),
t2 = runif(300, min = 1, max = 4),
t3 = runif(300, min = 1, max = 4),
cov = runif(300, min = 1, max = 4)
)
str(df)
I restructured the data..
df_long <- df %>%
pivot_longer(cols = c(t1, t2, t3),
values_to = "y",
names_to = "time")
and computed the lmm using the afex package.
lmm <- mixed(y ~ time*treatment*cov + (1|id), df_long)
There is a significant interaction between time and covariate in my data. So I want to compare the contrasts between measuring points depending on the covariate. I use certain values (-sd, m, +sd) for the covariate
lsmeans(lmm, ~time*cov,at = list(cov = c(1, 2, 3))) %>%
contrast(interaction = "pairwise", adjust = "bonferroni")
time_pairwise cov_pairwise estimate SE df t.ratio p.value
t1 - t2 1 - 2 -0.1555 0.0791 887 -1.967 0.4454
t1 - t3 1 - 2 -0.1251 0.0791 887 -1.582 1.0000
t2 - t3 1 - 2 0.0304 0.0791 887 0.385 1.0000
t1 - t2 1 - 3 -0.3110 0.1581 887 -1.967 0.4454
t1 - t3 1 - 3 -0.2501 0.1581 887 -1.582 1.0000
t2 - t3 1 - 3 0.0609 0.1581 887 0.385 1.0000
t1 - t2 2 - 3 -0.1555 0.0791 887 -1.967 0.4454
t1 - t3 2 - 3 -0.1251 0.0791 887 -1.582 1.0000
t2 - t3 2 - 3 0.0304 0.0791 887 0.385 1.0000
Degrees-of-freedom method: kenward-roger
P value adjustment: bonferroni method for 9 tests
In the output, I see that the contrasts between the values of the covariates all have the same standard error, the same t-values and the same p-values. What is the reason for this?
EDIT: When I plot my data, I can clearly see that the slope is different depending on the covariate. However, I still get the following output:
lsmeans(lmm,~time*cov, at = list(cov = c(-0.88, 0, 0.88))) %>% contrast(interaction = "pairwise", adjust = "bonferroni")
time_pairwise cov_pairwise estimate SE df t.ratio p.value
MZP2 - MZP1 (-0.88) - 0 -0.161 0.0755 201 -2.136 0.1016
MZP2 - MZP1 (-0.88) - 0.88 -0.323 0.1510 201 -2.136 0.1016
MZP2 - MZP1 0 - 0.88 -0.161 0.0755 201 -2.136 0.1016
Results are averaged over the levels of: treatment
Degrees-of-freedom method: kenward-roger
P value adjustment: bonferroni method for 3 tests
Depending on the covariate, is there another way to analyse the contrasts between the measurement points?
