0
$\begingroup$

My question is related to How do you deal with "nested" variables in a regression model? However, my issue is slightly different.

I'm analyzing RNA-seq gene expression data (>10 000 genes) from a weight loss intervention. The participants were measured and the RNA-seq was performed from samples taken before intervention and at two follow-up points (variable visit). The participants lose different amounts of body weight and some even regain weight from follow-up 1 to follow-up 2. We have also data from other cohort where weight loss intervention and the weight loss result were different. The cohorts cannot be directly compared, because the RNA-seq experiments are separate. However, to produce even somewhat comparable results, I would like to model the intervention effect on gene expression so that the results would show the change from baseline to each follow-up by percent body weight loss (wlp).

In this case, visit is the explanatory variable and wlp is nested in it. If I would have only two time points, I would get the result I desire by using only wlp as the explanatory variable (gives the same result as visit:wlp interaction).

My question is: In this case, when having three time points, does the visit:wlp interaction without main effects give me the change in outcomes (gene expression) per wlp at the two different visit? The estimates are very similar if I analyse the two follow-up points separately; however, p-values are smaller.

$\endgroup$
11
  • $\begingroup$ If wlp is an outcome, why are you using it as a predictor in an interaction term? Please edit the question to clarify; comments are easy to overlook and can be deleted. Also, if you are modeling a change from control, make sure not to use the control value as a predictor in the model. See these pages, for example. $\endgroup$
    – EdM
    Apr 10 at 14:24
  • $\begingroup$ Thank you, @EdM! I now clarified my post. My outcomes are gene expressions, not wlp. My aim is to get estimates of gene expression changes (fold change) per unit of weight loss (%) at different follow-up points. I do not have the control as a predictor in this model. $\endgroup$ Apr 11 at 6:46
  • $\begingroup$ Without getting directly into your question, you have modeled weight in a way that will be unlikely to fit the data. Almost never analyze a simple difference. Model new log weight as a function of initial log weight and of log height. People to not change weight by additive increments. Weight patterns are a mixture of absolute and relative changes, so a flexible model will cover both. $\endgroup$ Apr 11 at 12:17
  • $\begingroup$ Thanks for updating. Are the gene expression levels measured at all visits? How many genes are involved? (E.g., a few by qRT-PCR, or thousands by RNAseq?) In what ways can't the "other cohort" be compared to yours? What's the same and what's different? I'd be reluctant to answer your specific question in the last paragraph if I think there might be a problem with the models accomplishing what you think they do. Again, please provide extra information by editing the question. $\endgroup$
    – EdM
    Apr 11 at 17:31
  • $\begingroup$ Thank you, @Frank Harrell, for the important suggestion! Forgive me, I'm not sure if I understood your suggestion correctly. So could a better way to calculate the relative weight loss be: (log(weight) - log(baseline weight))/log(height). I read your study pubmed.ncbi.nlm.nih.gov/35127128 and I wonder if you meant a similar approach? I find myself struggling when trying to explain changes in a continuous outcome with changes in continuous explanatory variable. (Trying to identify potential mediators between the treatment and outcome.) I know that using change scores is not ideal. $\endgroup$ Apr 12 at 3:07

1 Answer 1

0
$\begingroup$

does the visit:wlp interaction without main effects give me the change in outcomes (gene expression) per wlp at the two different visit?

The problem with omitting main effects as you propose is that you aren't allowing for systematic gene-expression differences due to the value of visit if there is no weight change. I've been fooled too often by Mother Nature to be comfortable with that type of assumption.

To deal with that I'd recommend following the recommendations for combining a factor with a continuous predictor from this Bioconductor vignette on constructing design matrices. In your case, that would suggest something like:

model.matrix(~0+visit+visit:wlp)

That would also include a column in the model matrix for a visit:wlp interaction at the first visit, at which there is 0 wlp. Delete that column before you fit the model. Then you will get 3 coefficients for visit at the wlp=0 condition, and the two slopes you want from the remaining visit:wlp coefficients.

A couple of warnings

I think that answers your question, but there some additional matters to consider.

First, the above model assumes that all observations are independent. That's not the case when you have multiple measurements on the same individuals. See Bioinformatics, 37(2), 2021, 192–201 for discussion and a way to handle that in gene-expression studies.

Second, the above model assumes that gene expression is linearly associated with wlp. A strict linear association is seldom the case. Section 8.6 of the above vignette outlines how you can use regression splines to model such continuous predictors. It's explained in the context of time as the continuous variable, but the principle applies to any continuous variable.

$\endgroup$
1
  • $\begingroup$ Thank you, @EdM! In this data, the weight loss is so substantial that the visit:wlp interaction eliminates nearly all associations (for both visit and visit:wlp). Perhaps the best solution is to analyze the gene expression changes first by comparing visits, and then in a separate model, use weight loss as the fixed effect. I've previously included participants as a fixed effect, but thank you for the link to the dream() tool, which allows for more flexible modeling! Also, it's good that you reminded me about the likely nonlinearity. I'll definitely try splines in these datasets. $\endgroup$ Apr 15 at 12:12

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.