7
$\begingroup$

Patients received two different vaccines twice and we collected their blood after each vaccination to measure antibody levels in the blood. There are other variables but here I omitted them and just want to know how different vaccines affect the antibody levels. It seems that the longitudinal data mixed effect model is a proper analysis.

Here is the simulated data. 6 patients received Vaccines A or B for the first vaccine and all received Vaccine A as the second vaccine; after each vaccination they visited a hospital to measure blood antibody levels (2 data points at 2 Visits). It is already confirmed with studies that Vaccine A induced higher antibody levels than Vaccine B. So if using only Visit 1 data to do a linear regression, we can see Vaccine A induced significantly higher antibody levels compared to that with Vaccine B. As all patients received A as the second vaccine, if using only Visit 2 data, there is no significance. My question is, if we include all the data in a linear mixed effect model, how to extract coefficients corresponding to each visit data? I only know how to use the summary() to get the coefficients for both visits data, however, I could not find the coefficients to show similar results to linear regression of Visit 2 data that there were no differences in the antibody levels if everyone got Vaccine A,

data_raw = data.frame(ID=c(1,2,3,4,5,6,1,2,3,4,5,6), 
                   Antibody=c(50,60,70,30,40,35,101,102,102,102,101,103), 
                   Visit=c(1,1,1,1,1,1,2,2,2,2,2,2),
                   Vaccine=c("A","A","A","B","B","B","A","A","A","A","A","A"),
                   VaccineChange=c(0,0,0,0,0,0,0,0,0,1,1,1))

Linear regression of Visit 1 data confirmed that Vaccine A induced higher antibody levels.

lm_Visit1=subset(data_raw, Visit ==1)
summary(lm(Antibody~Vaccine,lm_Visit1))


Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   60.000      4.564  13.145 0.000193 ***
VaccineB     -25.000      6.455  -3.873 0.017948 *  

Linear regression of Visit 2 data indicated no significance as everyone got Vaccine A.

lm_Visit2=subset(data_raw, Visit ==2)
summary(lm(Antibody~VaccineChange,lm_Visit2))       


Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   101.6667     0.4714   215.7 2.77e-09 ***
VaccineChange   0.3333     0.6667     0.5    0.643 

If using mixed model with both visits data, the vaccine variable is significant. Is it possible to extract coefficients corresponding to each visit? Just to show similar results to linear regression that Visit 1 showed significance but not Visit 2.

summary(lmer(Antibody~Vaccine+Visit+(1|ID),data_raw))       


Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   18.179      6.438   7.602   2.824 0.023544 *  
VaccineB     -25.012      4.325   8.987  -5.784 0.000266 ***
Visit         41.827      3.702   6.380  11.299 1.87e-05 ***

summary(lmer(Antibody~VaccineChange+Visit+(1|ID),data_raw))


Fixed effects:
              Estimate Std. Error     df t value Pr(>|t|)   
(Intercept)      3.632      7.982  8.657   0.455  0.66031   
VaccineChange   20.930      6.357  4.008   3.292  0.03006 * 
Visit           43.868      4.602  3.587   9.532  0.00113 **

Originally, I performed two linear regression with each visit data. But I was told that for visit 2 data I should use linear mixed model as I got both visit data to account for repeated measurement within patient differences. Other covariates are not included in this simulated model just for simplicity.

Many thanks,

Jordan

$\endgroup$

1 Answer 1

7
$\begingroup$

Context: Since it is known that vaccine A induces higher antibody levels than vaccine B, I interpret the question to mean: How do we use visit 1 and visit 2 data to test that A has the same efficacy on the second visit no matter which vaccine a subject receives on the first visit? I describe two ways to do it, with and without random subject effects.


It's preferable to fit one model to all the data, rather than separate models to subsets of the data. This way information is appropriately pooled across multiple subsets/visits.

Your mixed models treat Visit as a continuous variable because it's of type "numeric". You get a coefficient estimate for the linear effect of number of visits (and it's number, not order).

This is not what you intend. So make Visit a categorical variable by casting it to type "factor".

data_raw$Visit <- factor(data_raw$Visit)

Let's specify a model for the full data by interacting Visit and Vaccine.

m1 <- lmer(Antibody ~ Vaccine * Visit + (1 | ID), data_raw)
#> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

We get a warning message. Since all subjects receive the same vaccine on their second visit, lmer drops the interaction. We can compare the two vaccines on the first visit or receiving A on both visits. But we cannot compare the follow-up efficacy of A between those who receive A first and those who receive B first.

Analysis I: mixed model for repeated measurements

Instead let's introduce a new variable FirstVaccine which indicates the vaccine received on the first visit. (Acknowledgement: The OP mentioned this approach in a comment.)

data_raw$FirstVaccine <- rep(data_raw$Vaccine[data_raw$Visit == 1], 2)

m2 <- lmer(Antibody ~ FirstVaccine * Visit + (1 | ID), data_raw)

Now we can estimate the main effects as well as the interaction. We can evaluate the follow-up efficacy of A as the second-visit contrast between those who receive A and those who receive B on their first visit. This is easy to do with the emmeans package.

pairs(emmeans(m2, ~ FirstVaccine | Visit))
#> Visit = 1:
#>  contrast estimate   SE   df t.ratio p.value
#>  A - B      25.000 4.59 7.99   5.448  0.0006
#> 
#> Visit = 2:
#>  contrast estimate   SE   df t.ratio p.value
#>  A - B      -0.333 4.59 7.99  -0.073  0.9439

Analysis II: pre-post treatment comparison

The mixed model accounts for the correlation between observations from the same subject by introducing random subject effects.

Alternatively, we can consider the two visits as paired (before and after) observations and analyze this study as a pre-post treatment design; the first vaccine is the "treatment".

First we rearrange the data so that two observations of the same subject are on the same row. In the "wide" table the rows are independent, so there is no need for random effects.

visit1 <- subset(data_raw, Visit == 1, c("ID", "Antibody", "Vaccine"))
visit2 <- subset(data_raw, Visit == 2, c("ID", "Antibody", "Vaccine"))

data_wide <- merge(visit1, visit2, by = "ID", suffixes = c("_1", "_2"))
data_wide
#>   ID Antibody_1 Vaccine_1 Antibody_2 Vaccine_2
#> 1  1         50         A        101         A
#> 2  2         60         A        102         A
#> 3  3         70         A        102         A
#> 4  4         30         B        102         A
#> 5  5         40         B        101         A
#> 6  6         35         B        103         A

m3 <- lm(Antibody_2 ~ Antibody_1 + Vaccine_1, data = data_wide)

The estimate of the second-visit difference in antibodies between the two groups is bigger (in absolute value), with a smaller standard error, as there is less variability in antibody levels on the second visit than on the first.

pairs(emmeans(m3, ~ Vaccine_1))
#>  contrast estimate   SE df t.ratio p.value
#>  A - B      -0.833 1.65  3  -0.506  0.6475

Keep in mind this distinction between the two analyses: the mixed model assumes the error variance is the same on both visits; the pre-post comparison doesn't and it estimates only the second-visit error variance.

$\endgroup$
6
  • 1
    $\begingroup$ Thank you so much for your time and detailed answer. It really helped me out. One more question, if we did not change the one participant to receive Vaccine B, pairs(emmeans(model, ~ Vaccine | Visit)) can not generate contrast A-B for Visit 2. Is there a way to interpret that Vaccine change does not make a difference to antibody levels at Visit 2? $\endgroup$
    – Jordan Lau
    Commented Aug 18, 2022 at 20:42
  • $\begingroup$ Thanks. Using 2nd visit data only would be my original linear regression test (not favored by a statistician :( ). I think you are right when there is at least one Vaccine B in the 2nd visit. Somebody suggests that if everyone got the same vaccine at Visit 2, we could use the InitialVaccine variable rather than the Actual Vaccine in the model. InitialVaccine=c("A","A","A","B","B","B","A","A","A","B","B","B")). model <- lmer(Antibody ~ InitialVaccine * Visit + (1 | ID), data_raw). Any comments on this? Thanks again! $\endgroup$
    – Jordan Lau
    Commented Aug 18, 2022 at 21:11
  • $\begingroup$ Yes, I think Antibody ~ InitialVaccine * Visit + (1 | ID) works. You can compare the visit2 antibody of those who got A to those who got B on the first visit. You can also estimate the random subject effects because you use all the data. All good. $\endgroup$
    – dipetkov
    Commented Aug 18, 2022 at 22:06
  • $\begingroup$ And here is another option, if you are strictly interested in whether both groups have the same antibodies on the second visit: Consider this a pre-treatment post-treatment design: Antibody_Visit2 ~ Antibody_Visit1 + InitialVaccineB. No need for random effects and you use all the data at once. $\endgroup$
    – dipetkov
    Commented Aug 18, 2022 at 22:41
  • $\begingroup$ I really appreciate your time! Thank you so much! I like the pre-post design. In the real data, we have more variables: eg, drug treatment, age, infection, vaccine, etc. And we want to account for these variables in the 1st visit model and 2nd visit model. For the 1st visit data, the option is a linear regression model. For the 2nd visit data, a linear mixed-effect model or your pre-post design to include visit 1 data, are both fine. If not everyone got same vaccine on the 2nd visit, we might need to stick with the mixed-effect model. Thanks again! $\endgroup$
    – Jordan Lau
    Commented Aug 19, 2022 at 19:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.