# An independent variable that is correlated with another variable in a regression model

I'm doing a regression analysis to understand the relationship between disease severity (Severity) and viral load (VL). The VL data were collected on different days after symptom onset in different patients, and I know that VL is highly correlated with the number of days between symptom onset and sample collection (DAY) (adjusted R-squared: 0.088, p-value: 1.7e-09), with the following linear relationship:

VL = 4.7 - 0.2 * DAY


However, the variable DAY wouldn't affect Severity.

My question is how I should set up the model for Severity and VL. I came up with 4 options and couldn't figure out which one I should use. If anyone could shed some light on this, it would be much appreciated.

Option 1: include both variables

Severity ~ VL + DAY


Option 2: include the interaction term

Severity ~ VL + DAY + VL:DAY


Option 3: include only the interaction term

Severity ~ VL:DAY


Option 4: create a new variable VL_NEW to adjust for DAY, and include only the new variable VL_NEW

VL_NEW = VL + 0.2 * DAY
Severity ~ VL_NEW


My thought is to use option 4 since you seem quite confident of the relationship between day and severity.

You write "day wouldn't affect severity". This is impossible. If day is highly correlated with VL and VL affects severity, then day has to affect it. I think what you mean is that day wouldn't affect severity outside of its relationship with VL.

One reason not to use model 4 would be if severity changes quickly. Then I would look at options 1 and 2. This varies by disease (e.g. the severity of a cold varies a lot from day to day; the severity of HIV symptoms changes more slowly).

Option 3 is a no-no. It is almost never right to include the interaction and not include the main effects. (David Rindskopf of CUNY has written about some exceptions, but they are very rare).

I would also investigate the relationship between VL and day. It seems to me unlikely that it would be linear.

(Also, I think you have a typo of $$R^2 = 0.088$$. That would be not be strong correlation by any definition).

• Thanks Peter. It really helps! I looked at my data again – I guess VL and DAY are not that highly correlated in my dataset, but it's clear that VL decreases over time. How would you suggest I choose between option 1 and 2? Commented Sep 10, 2023 at 17:25
• Look at both and compare AIC or BIC. Commented Sep 10, 2023 at 18:04

As alluded to above, the relationship between viral load in an individual and the time after infection is highly nonlinear (and you may well be interested in the kinetics of that behaviour scientifically).

Fortunately, this happens to be a very well studied area of epidemiology.

For diseases that you tend to completely recover from, a large number of models have been proposed. These tend to feature an early exponential growth phase, and a later exponential elimination phase representing various forms of lagged immune responses. For a recent example, upper respiratory viral load following COVID infection approximately follows the following set of DEs:

\begin{align} \frac{dI_\text{(infected cells)}}{dt} &= \beta V\left[1-\frac{I}{I_\text{max}}\right] -k_a I \frac{A_{3}^{2}}{A_{3}^{2} + A_{50}^{2}}\\ \frac{dV_\text{(viral load)}}{dt} &=p I -\gamma V\\ \frac{dA_{1\text{(immune response phase 1)}}}{dt} &=\frac{I}{I+I_{50}} -k_c A_1 \\ \frac{dA_{2\text{(immune response phase 2)}}}{dt} &=k_c A_1 - k_c A_2 \\ \frac{dA_{3\text{(immune response phase 3)}}}{dt} &=k_c A_2 \\ \end{align}

where $$I_{50}$$, $$I_\text{max}$$, $$A_{50}$$, $$k_a$$ and $$k_c$$ are parameters related to immune function and $$\beta$$, $$p$$ and $$\gamma$$ are parameters related to viral function. This corresponds to this diagram of immune function, provided by the authors:

This is a nice, abstracted model that is based in the biology of how the innate and adaptive immune system works (at least as I understand it).

Note that in panel B the plotted model has these two phases -- "going up" and "going down". I don't know what data you're dealing with, but there's a very real risk that you confidently draw a straight line bang through the middle, from 0 onwards, with a nonzero offset, whereas a slightly more complex model (not necessarily this complex!) may be far more appropriate -- and ascertaining which is better for your data is a perfect example of where AICs/BICs are appropriate (I think).

For diseases that you don't eliminate, such as HIV or hepatitis B and C, there are still excellent models for viral load as a function of time in the absence of treatment -- have a look at A methodology for deriving the sensitivity of pooled testing, based on viral load progression and pooling dilution, which covers both an overview of the models based on clinical data and a nice little aside about how to infer steady-state dynamics for batched (cheaper) testing. Because (unfortunately) your immune system doesn't eliminate the virus completely, you end up with effectively a straight line with time at longer times – where a linear model may well be appropriate.

• Thank you for these interesting references and the detailed explanation! My data are mostly in the exponential elimination phase where late response reduces infected cells and viral load. Commented Sep 11, 2023 at 21:51
• @Michael thank you for the kind words. If you expect to be in the exponential elimination phase, have you taken logs of the viral load in your regression? (You can use model selection methods to see if that is numerically justified too). Then you really only have one causal relationship going on – symptoms ~ viral load – which, well, is what you have hypothesised above (and given in the answer by Peter, which I agree with) Commented Sep 12, 2023 at 10:13
• yes, the viral load has been log transformed before putting into the regression model. In my dataset, the regression line is pretty similar to the figure you showed above - a linear relationship between log viral load and day since symptom onset. Commented Sep 14, 2023 at 22:02

In R:

model1 = lm(Severity ~ VL + DAY + VL:DAY, data=obs)
model2 = lm(Severity ~ VL + DAY, data=obs)
anova(model1, model2)


If you get P<0.05 in the model comparison then the models are different and you must stick to model1. Otherwise, you choose the simpler model2. Remember to check the residuals:

plot(model1, 1)
plot(model1, 2)


Assess the two plots visually.