I am trying to analyze the data with the pre-post-control design in the context of RNA-seq analysis.

I have read Best practice when analyzing pre-post treatment-control designs, but I am still confused because the method recommended by this manual(page 34) is not mentioned in the above post.

To briefly describe my data, the response is gene expression measurement and explanatory variables contain 3 factors(treatment, patient, and time): there are 2 levels(drug & control) in treatment, total 10 patients, and two levels(pre and post) in time.

The goal is to figure out if the drug has any effect on the gene expression while controlling for heterogeneity among patients.

The software manual suggests using interaction terms to control for individual effects, and use the following model:

$$expression \sim treatment + treatment:patient + treatment :time$$

and the code snippet below describes the idea:

dat <- structure(list(treatment = c("control", "control", "control", 
"control", "control", "control", "control", "control", "control", 
"control", "drug", "drug", "drug", "drug", "drug", "drug", "drug", 
"drug", "drug", "drug"), time = c("pre", "post", "pre", "post", 
"pre", "post", "pre", "post", "pre", "post", "pre", "post", "pre", 
"post", "pre", "post", "pre", "post", "pre", "post"), patient = c(1, 
1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -20L), .Names = c("treatment", 
"time", "patient"))

# prepare for the design matrix 
dat <- dat %>% group_by(treatment) %>% 
              mutate(patient=paste0("p",rep(c(1:5),each=2))) %>%
              ungroup() %>%  
              mutate(time=factor(time, levels = c("pre","post")))

   treatment time patient
1    control  pre      p1
2    control post      p1
3    control  pre      p2
4    control post      p2
5    control  pre      p3
6    control post      p3
7    control  pre      p4
8    control post      p4
9    control  pre      p5
10   control post      p5
11      drug  pre      p1
12      drug post      p1
13      drug  pre      p2
14      drug post      p2
15      drug  pre      p3
16      drug post      p3
17      drug  pre      p4
18      drug post      p4
19      drug  pre      p5
20      drug post      p5

# construct a model matrix
mmat <- model.matrix(~treatment + treatment:patient + treatment:time,data = dat)

# terms that will be in my linear model
 [1] "(Intercept)"                "treatmentdrug"             
 [3] "treatmentcontrol:patientp2" "treatmentdrug:patientp2"   
 [5] "treatmentcontrol:patientp3" "treatmentdrug:patientp3"   
 [7] "treatmentcontrol:patientp4" "treatmentdrug:patientp4"   
 [9] "treatmentcontrol:patientp5" "treatmentdrug:patientp5"   
[11] "treatmentcontrol:timepost"  "treatmentdrug:timepost"

To see the effect of the drug, the manual suggests check if the coefficient "treatmentdrug:timepost" is significant by computing log likelihood ratio between the model without that term and the model with that term and see if the nested model is sufficiently different from the full model.

Now my question is :

  • what's weakness of this approach of using interaction terms to control for patient effects? In other words, what assumptions am I breaking(if any)?
  • It seems that random effects models are used in this context, but what are the benefits?
  • If you were me, how would you have tackled this problem?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.