Hello,

I wanted to ask if using the following model design was appropriate for the question I'm trying to answer. I have samples that are paired (BG and AG status). The samples come from 8 individuals and each individual has a BG and AG sample. Additionally, there are 4 treatment groups (A, C, T, E) and a control (W).

Sample ID | Treatment | Status |

A01BG | A | BG |

A01AG | A | AG |

C01BG | C | BG |

C01AG | C | AG |

I wanted to compare the difference between the two statuses across the treatments and control. I was able to use this model in DESeq2, but my local statistician suggested that I compare the results to a package that uses a different distribution.

DESeq2 model: ~treatment + treatment:status

results(dds, contrast=list("TreatmentC.StatusBG","TreatmentA.StatusBG"))

Based on this site, I created a similar model, but received this output:

```
# Define the metadata categories
TB.treatment <- pData(n_paired_merged_MRE)$Treatment
TB.ID <- pData(n_paired_merged_MRE)$id
TB.bs <- pData(n_paired_merged_MRE)$BrushStatus
# Define the normalisation factor
normfactor_MRE = normFactors(n_paired_merged_MRE)
normfactor_MRE = log2(normfactor_MRE/median(normfactor_MRE)+1)
# Create the model
TB.mod <- model.matrix(~ TB.treatment + TB.bs + TB.ID + normfactor_MRE)
settings <- zigControl(maxit = 10, verbose = TRUE)
TB.fit <- fitZig(obj = n_paired_merged_MRE, mod = TB.mod, useCSSoffset = FALSE, control = settings)
# output
Coefficients not estimable: TB.IDW21AG TB.IDW21BG normfactor_MRE TB.IDC21BG TB.IDE21BG TB.IDT21BG
it= 0, nll=NaN, log10(eps+1)=Inf, stillActive=1141
```

Coefficients not estimable: TB.IDW21AG TB.IDW21BG normfactor_MRE TB.IDC21BG TB.IDE21BG TB.IDT21BG Error in lm.fit(mmZero, qlogis(pi)) : NA/NaN/Inf in 'y' In addition: Warning messages: 1: Partial NA coefficients for 1141 probe(s) 2: Partial NA coefficients for 1137 probe(s)

Thanks for the help!