linear mixed model: visualising fixed effect and interaction with continuous variable (with confidence interval)

Question

My aim is to analyse the blood concentration of a biomarker and its annual increase in four different subject groups using longitudinal measurements, i.e. biomarker concentrations of each subject were measured at multiple points of time. The dataset consists of (1) the biomarker concentration as the dependent variable and (2) the independent variables age (in years, centred at the mean age of all subject measurements to obtain meaningful intercepts) and group (categorical factor with four levels, i.e. there were four subject groups A, B, C and D, with A being the reference group). The grouping variable is subject (i.e. the subject pseudonym). The number of observations is roughly 320 (4 groups, 20 subjects per group, 4 visits per subject).

A linear mixed model was run to predict the biomarker concentration by the fixed effects of “group” and “age”, their interaction and random intercepts for the variable “subject”, using the lme4 package in R with the following code:

model <- lmer(biomarker ~ group * age + (1 | subject), 
              data = biomarker_data) 

Now I would like to ask for your help to visualise the model estimates with their confidence intervals in the following manner:

(a) Visualising the estimated biomarker concentration of each group (A, B, C, D) at the mean age (i.e. age = 0) as a bar plot with error bars of the 95% confidence interval: The program output only provides me with the group estimates (of group B, C, D) and confidence intervals relative to the reference group A (i.e. intercepts relative to group A). What I would like to display for each group, however, is the estimate and confidence interval of the absolute biomarker concentration. The confidence interval should relate to the fixed effect without the random effect of the variable subject.

(b) Visualising the estimated annual biomarker increase of each group (A, B, C, D) as a bar plot with error bars of the 95% confidence interval: Again, the program output only provides me with the interaction values of the groups B, C, and D relative to the reference group A. What I would like to display for each group, however, is the estimate and confidence of the absolute annual increase. Again, the confidence interval should relate to the interaction of the fixed effect without the variation of the random variable subject.

(c) Visualising the biomarker concentration (y-axis) as a function of age (x-axis) for each group with the confidence interval of the model estimate. My problem is not to plot the four regression lines, but the confidence interval around each regression line. As in the above, I am interested in the fixed effect without the random subject variation. (My question (a) is a special case of question (c)).

Answer 1

I will answer your questions (a) and (b) more explicitly, since Dimitris already suggested some ways for you to move forward. My reasoning in answering these questions is that it would be helpful for you to see the mechanism involved in answering these questions so that you can understand what is going on behind the hood.

Your current model looks like this:

biomarker_ij = (beta0 + u_i) + beta1*groupB_i + beta2*groupC_i + beta3*groupD_i + beta4*age_i + 
                beta5*groupB_i*age_i + beta6*groupC_i*age_i + \beta7*groupD_i*age_i + \epsilon_ij 

where i indexes the subject, j indexes the visit within subject and u_i is the random intercept effect associated with the i-th subject. This notation assumes that both the group dummies and the age are between-subject variables (their value is the same across visits for the same subject).

For your question (a), when age is 0, then:

1. The mean value of the biomarker concentration for group A is given by beta0 (fixed effect intercept);

2. The mean value of the biomarker concentration for group B is given by beta0 + beta1;

3. The mean value of the biomarker concentration for group C is given by beta0 + beta2;

4. The mean value of the biomarker concentration for group D is given by beta0 + beta3.

So to answer your first question, you can set up 4 linear combinations of beta coefficients and compute 95% confidence intervals for them using the multcomp package. (Your model has a total of 8 beta coefficients, with each of these coefficients representing a fixed effect in the model.) The combinations are:

C1 = beta0 

C2 = beta0 + beta1

C3 = beta0 + beta2 

C4 = beta0 + beta3

For example, given that your fixed effects portion of the model includes 8 beta coefficients (i.e., beta0 through beta7), then the weights of the beta coefficients involved in specifying the combination C1 can be specified in R language as:

W1 <- c(1, 0, 0, 0, 0, 0, 0, 0)

The weights corresponding to the remaining combinations can be specified in R language like this:

W2 <- c(1, 1, 0, 0, 0, 0, 0, 0)

W3 <- c(1, 0, 1, 0, 0, 0, 0, 0)

W4 <- c(1, 0, 0, 1, 0, 0, 0, 0)

These weights can then be aggregated into a matrix of weights with the rbind command:

W <- rbind(W1, W2, W3, W4) 

The rows of W can be conveniently named to indicate what linear combination of betas you are interested in:

rownames(W) <- c("beta0", "beta0 + beta1", "beta0 + beta2", "beta0 + beta3")

The columns of W can be named to indicate what weight goes with what beta:

colnames(W) <- c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "beta7")

The multicomp package can now be invoked as follows:

model.glht <- glht(model, linfct = W)

where glht stands for general linear hypothesis testing. Then the R commands below will

summary(model.glht)

confint(model.glht) 

will give you what you want: p-values for testing the null hypothesis that each combination is zero versus the alternative hypothesis that it is different from zero, as well as simultaneous confidence intervals for all combinations. (See https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf.)

For your question (b), it sounds like you are interested in reporting the effect of age separately for each group. The model written above can be re-expressed to showcase the effect of age:

biomarker_ij = (beta0 + u_i) + beta1*groupB_i + beta2*groupC_i + beta3*groupD_i + 
               (beta4 + beta5*groupB_i + beta6*groupC_i + \beta7*groupD_i)*age_i + \epsilon_ij.

As seen here, the effect of age is quantified by the slope beta4 + beta5*groupB_i + beta6*groupC_i + \beta7*groupD_i. Using the definition of the dummy group variables, that means that:

1. The effect of age on biomarker concentration in group A is given by beta4;

2. The effect of age on biomarker concentration in group B is given by beta4 + beta5;

3. The effect of age on biomarker concentration in group C is given by beta4 + beta6;

4. The effect of age on biomarker concentration in group A is given by beta4 + beta7.

You can set contrasts as before and then derive the associated weights and feed them to the multcomp package.