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). 

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!