Table 2 Summary of methods compared

MCMSeq

mcmseq

Bayesian negative binomial hierarchical model, assuming the following

hyper-parameters: T=7²I₄, where I₄ is a 4 x 4 identity matrix, U=

V=0.01, A_g= centered on empirical prior as described in 2.3, \(B=4\cdot \sigma ^{2}_{re}\)

Chains were run for 30,000 iterations and 10% was discarded as burn-in.

CPGLMM ^‡

cplm

Compound Poisson generalized linear mixed model fit using maximum

[44]

likelihood.

DESeq2*

DESeq2

DESeq2 analysis using default settings. To account for correlation,

[15]

subjects are included as fixed effects in the model by replacing the

random effects in Equation 4 with fixed effects.

DESeq2

DESeq2

DESeq2 analysis using default settings, including only fixed effects in

[15]

Equation 4. Correlation is ignored.

edgeR*

edgeR

edgeR analysis using default settings. To account for correlation,

[14]

subjects are included as fixed effects in the model by replacing the

random effects in Equation 4 with fixed effects.

edgeR

edgeR

edgeR analysis using default settings, including only fixed effects in

[14]

Equation 4. Correlation is ignored.

limma

limma

Count data were transformed for analysis with limma using voom.

[26]

The duplicateCorrelation function using subject id as the blocking

[27]

factor was used to account for correlation.

LMM

lmerTest

Data were first transformed using DESeq2’s variance stabilizing

[45]

transformation. Linear mixed models, which assume a normally

distributed errors, were then fit to the transformed counts.

MACAU

MACAU2

Poisson GLMM with random effects to account for over-dispersion and

[19]

correlation among observations. A known covariance matrix must be supplied;

we used the covariance of the centered and scaled gene expression matrix as

shown in the third example of [19].

NBGLMM

glmmADMB

Frequentist negative binomial generalized linear mixed model fit using

[46]

maximum likelihood.

ShrinkBayes

ShrinkBayes

Bayesian hierarchical model using shrinkage priors to control false

[23]

discovery rates. We used a negative binomial distribution for the count

data and shrinkage priors for all model parameters involved in statistical

tests. All other parameters were set to default values.

VarSeq ^‡

tcgsaseq

Model free method for identifying gene sets whose expression changes

[25]

over time. Counts are transformed to log(CPM) for analysis. The

asymptotic test was used to generate p-values.