Interpretation of interaction coefficients and p-values with categorical predictors

Question

Let's say that we have gene expression measurements from bacteria after heatshock. The expression changes are collected as log2 changes, and are approximately normal and centered around 0. Ribosomal genes decrease strongly in response to heatshock. Conserved genes also decrease in response to heatshock, but this gene group overlaps a lot with the ribosomal gene group.

My advisor asks:

• Do conserved genes decrease on their own or only because they include ribosomal genes?
• What is the statistical significance of the effect on ribosomal genes? On conserved genes?

The overarching goal is to know which attribute (ribosome or conserved) is most correlated with gene expression changes after heat shock.

My first thought is that people generally use a two-way ANOVA to measure these effects. However, the data, like most biological data in my experience, doesn't meet the assumptions of ANOVA -- specifically both unequal variance between groups and extremely different group sizes.

Let's say that in case A, the effect is only from the ribosome genes and conserved genes have no effect. In case B, conserved genes also have an effect. In both cases, 150/200 ribosome genes are also in the class of conserved genes.

Here is some example synthetic data / test case in Python:

import numpy as np
import pandas as pd

seed = 0
rand = np.random.RandomState(seed)  

# 4000 background genes, 200 ribosome genes, 200 conserved genes
n_bg, n_ribo, n_cons = (4000, 200, 200)

# Create random normal distribution for each group of genes 
y = rand.normal(0, 3, n_bg)
# Ribosome genes mean change = -4 and group has less variance than other genes
y_ribo = rand.normal(-4, 1, n_ribo)
# Conserved genes either do not change (case A) or also decrease (case B)
y_cons_A = rand.normal(0, 2, n_cons)
y_cons_B = rand.normal(-2, 2, n_cons)
# 150 of the ribosome genes are also conserved
num_ol = 150

# Case A: conserved, non-ribo genes don't change
df_A = pd.DataFrame(pd.concat([pd.Series(y), pd.Series(y_ribo), pd.Series(y_cons_A)]), columns=['log_change'])
# Case B: conserved, non-ribo genes also change
df_B = pd.DataFrame(pd.concat([pd.Series(y), pd.Series(y_ribo), pd.Series(y_cons_B)]), columns=['log_change'])

dfs = {'A':df_A, 'B':df_B}

for case in dfs:
    dfs[case]['ribosome'] = [0]*n_bg + [1]*n_ribo + [0]*n_cons
    dfs[case]['conserved'] = [0]*(n_bg+n_ribo-num_ol) + [1]*(num_ol+n_cons)

Here is what the data looks like plotted:

And here are the results of running OLS regression to model the change in expression either with or without interactions:

We can see right away that adding the interaction term was important in the case B, where both conserved and ribosome groups had their own effects and the interaction term was significant. Addition of the interaction term caused the ribosome coefficient to approach its true value (-4).

The two things I am confused about:

1. I do not really understand the interaction coefficient for case B +interaction (1.68). If you think about the effect of conserved on ribosome, then it could be positive because conserved has less extreme changes on average than ribosome (-2 vs. -4). But if you think about the effect of ribosome on conserved, then it should be negative, because ribosome has more extreme changes on average than conserved (-4 vs. -2). So what is this positive interaction coefficient telling us?

2. Is it acceptable to use the more negative interaction coefficient for ribosome vs. conserved in model B +interaction to conclude that ribosome has a larger effect than conserved? They both have very significant p-values. If not, what is the alternative?

This type of problem is really common in exploratory data analysis of biological data, but somehow I am having trouble figuring out what the interpretation is and what statistics are appropriate to report. I am happy to hear suggestions for alternative methods and please feel free to write your response in R.

Answer 1

This situation is best handled with tools that account for unequal variances and differences in gene-set sizes. The Bioconductor limma package has extensive tools for analyzing gene-expression data, including ways to handle RNA-seq data like I assume you have here. It's been developed extensively over a couple of decades to address the specific types of issues that you correctly raise (and many other issues, too). It can use all the gene-expression data to estimate things like how variance changes with expression levels, and then do post-modeling analysis in many different ways.

Its mroast() function allows for comparisons among pre-specified gene sets in terms of responses to experimental manipulation (heat shock in your case), although I haven't used it myself. To answer your scientific question, I think you could directly compare the set of ribosomal "conserved" genes against the set of non-ribosomal "conserved" genes to evaluate whether those 2 groups differ in differential expression following heat shock.

To answer your specific question, interaction models can be tricky to interpret. In your most complicated model (B + interaction), the intercept is the estimated log-2 change for non-conserved and non-ribosomal genes. The coefficient for ribosome is the difference from the intercept for all ribosome genes. The coefficient for conserved is the difference from the intercept for all conserved genes. The interaction coefficient is the extra difference for conserved ribosomal genes, beyond what you would expect from the sum of the individual conserved and ribosome coefficients.

In your case, to get the estimated log-2 change specifically for conserved ribosomal genes, you would use the sum -0.0764 - 3.7988 -1.9199 + 1.6834 = -4.1117.

Comparing individual coefficients can be tricky. You need to consider both the variances of the individual coefficient estimates and their covariances. Given the multiple difficulties of trying to use simple ANOVA for differential gene expression, use the tools in limma or in other established software packages that account for those difficulties.