Hi all,

I'm trying to analyse some RNAseq data which has been generated using an experimental design not dissimilar to the example in the user guide (section 3.5, p41), but with some differences that are making applying that model difficult.

My study design:

> Pheno_mat
   disease Patient Treatment
1  healthy       1      None
2  healthy       1    Drug_A
3  healthy       1    Drug_B
4  healthy       2      None
5  healthy       2    Drug_A
6  healthy       2    Drug_B
7  disease       3      None
8  disease       3    Drug_A
9  disease       3    Drug_B
10 disease       4      None
11 disease       4    Drug_A
12 disease       4    Drug_B
13 disease       5      None
14 disease       5    Drug_A
15 disease       5    Drug_B

I want to find: The impact of drugA/B on healthy and diseased samples. The difference in the impact of DrugA/B in disease vs healthy.

The reason I can't apply the model near exactly as shown in 3.5 is that I do not have a balanced number of disease/healthy samples, so revising the targets names to 1-3 in both disease groups doesn't work, I still have design matrix that is not full rank:

> Pheno_mat
   disease Patient Treatment Patient_fix
1  healthy       1      None           1
2  healthy       1    Drug_A           1
3  healthy       1    Drug_B           1
4  healthy       2      None           2
5  healthy       2    Drug_A           2
6  healthy       2    Drug_B           2
7  disease       3      None           1
8  disease       3    Drug_A           1
9  disease       3    Drug_B           1
10 disease       4      None           2
11 disease       4    Drug_A           2
12 disease       4    Drug_B           2
13 disease       5      None           3
14 disease       5    Drug_A           3
15 disease       5    Drug_B           3

So, what are my options? 1: I can't follow the example as I have a different number of patients in healthy/disease, so my design matrix isn't full rank 2: Disease state is explained by patient ID, but I can't use patient:disease then handle disease state in the contrasts because I get n=1 in each group and can't estimate dispersions 3: If I simply fit a nested Disease+Drug model I ignore that my patient data are paired...

Any suggestions would be gratefully received!

Thanks,

Dean

The only reason the matrix isn't full rank is because you have coefficients that don't estimate anything.

> phenomat <- data.frame(disease = factor(rep(c("healthy","disease"), c(6,9))), patient = factor(rep(c(1:2,1:3), each = 3)), treatment = factor(rep(c("Drug_A","Drug_B","None"), 5)))
> design <- model.matrix(~disease + disease:patient + disease:treatment, phenomat)
> which(colSums(design) == 0L)
diseasehealthy:patient3 
                      6 

## Note that there isn't a healthy patient 3!

> design <- design[,-6]
> library(limma)
> is.fullrank(design)
[1] TRUE