I have been testing how different design matrices and different ways of specifying a linear model via R's formula interface affect modeling results, within the context of differential expression analysis with DESeq2. I noticed that when I specified two models that from what I understand should decompose the data variance in equivalent ways, and then checked the corresponding coefficients from the two models that should have the same interpretation, the numerical values of the coefficients (log2FoldChange) and P values have some minor differences. However I could be wrong in my statistical claims above (I'm not a statistician). So I decided to seek expert help here to double check whether it's indeed just a numerical issue or maybe I missed something.

As a toy example, I'm using the same "pasilla" dataset used in the DESeq2 tutorial. The sample meta data is given below. The group variable was manually created simply by combining condition and type:

             condition   type           group
  treated1     treated single   treatedsingle
  treated2     treated paired   treatedpaired
  treated3     treated paired   treatedpaired
  untreated1 untreated single untreatedsingle
  untreated2 untreated single untreatedsingle
  untreated3 untreated paired untreatedpaired
  untreated4 untreated paired untreatedpaired

In the code snippet below I fit two models, one with interacting condition and type terms, the other simply with group, which I think should be equivalent. The base levels of the factors were also specified manually, which allowed me to then check for coefficients from each model that should both correspond to the difference between the "treated-single" group and the "untreated-single" group (at least that's from what I understand). However, the resulting log2FoldChange values and P values have some minor differences, although the correlation across all genes are essentially 1. My question is: Is this simply a numerical issue caused by loss of precision during computation, or is there another technical cause, or was I wrong in my understanding of the models and their interpretations?

Thanks a lot in advance!

library(DESeq2)
library(pasilla)

# count matrix
cts <- data.matrix(read.csv(system.file("extdata", "pasilla_gene_counts.tsv", package="pasilla", mustWork=TRUE), sep="\t", row.names="gene_id"))
cts <- cts[, c(5:7, 1:4)]
# meta data
coldata <- read.csv(system.file("extdata", "pasilla_sample_annotation.csv", package="pasilla", mustWork=TRUE), row.names=1)
coldata <- coldata[, c("condition", "type")]
rownames(coldata) <- sub("fb", "", rownames(coldata))
coldata$condition <- factor(coldata$condition, levels=c("untreated", "treated"))
coldata$type <- factor(substr(coldata$type, 1, 6), levels=c("single", "paired"))
coldata$group <- with(coldata, factor(paste0(condition, type), levels=c("untreatedsingle", "treatedsingle", "untreatedpaired", "treatedpaired")))

# model 1
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition * type)
dds <- DESeq(dds)
res1 <- results(dds, contrast=c("condition", "treated", "untreated"))

# model 2
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ group)
dds <- DESeq(dds)
res2 <- results(dds, contrast=c("group", "treatedsingle", "untreatedsingle"))

# compare the results from the two models:
all.equal(res1$log2FoldChange, res2$log2FoldChange)
# [1] "Mean relative difference: 0.008926539"
all.equal(res1$pvalue, res2$pvalue)
# [1] "Mean relative difference: 0.001391931"

Yes there can be small numerical differences in the estimated quantities because of how the algorithms decide to stop iterating.