Thank you very much for your answer, that was very helpful!

Impossible to say. It depends on your data, particularly on how the donors are associated with batches. Normally I would prefer to keep donor in the model.

I did some more investigation of my data and totally agree with you. For me, SVA was unable to "pick up" the donor effect (i.e. still siginificant association between donor and PC1-5). To account for this, my idea was to include this in the model that I give to sva like so:

mod <- model.matrix(~donor+group, colData(dds))
mod0 <- model.matrix(~donor, colData(dds))
svseq <- svaseq(norm_counts, mod, mod0, n.sv = 5)

Do you by chance know if this is the correct approach or should mod0 look like this: model.matrix(~1, colData(dds)).

Furthermore, I downstream I am using limma::removeBatchEffect for visualizations. Is the following usage in this case correct given my SVA code:

mod <- model.matrix(~group, tmp_mdata)
formula <- as.formula(paste("~", paste0("SV", 1:5, collapse = " + ")))
mod0 <- model.matrix(formula, tmp_mdata)[, -1]
# adjust the counts using limma
adjusted_counts <- limma::removeBatchEffect(norm_data,
        batch = mdata$donor,
        covariates = mod0,
        design = mod
    )

Any help is much appreciated!