I have this experimental design:

Two groups A and B. Individuals from group B were genetically manipulated such that when they are given a certain drug the drug turns on a gene that was inserted through this genetic manipulation and that gene is supposed to curb a reaction to a stimulus.

In order to test whether this gene really curbs the reaction to a stimulus I've set the following scheme:

My individuals are mice and the stimulus is placing a piece of their food in their cage, and the reaction time to that stimulus, which is the response that I'm measuring, is how long it takes them to go and eat that piece of food.

Each group consists of 7 individuals. At trial 1 I just placed the food inside the cage of each mouse and measured the reaction time. At trial 2, I treated all mice with the drug and again introduced the stimulus and measured the reaction time. The trial was ended after 5 minutes at which not all mice went and ate that piece of food. Hence there is censoring. Finally, at trial 3, which proceeded trial 2 by enough time for the effect of the drug to disappear, again I introduced the stimulus and measured the reaction time, again with the 5 minute censoring.

Here are example data:

df <- data.frame(mouse_ID = rep(c(paste0("A",1:7),paste0("B",1:7)),3),
                 group = rep(c(rep("A",7),rep("B",7)),3),
                 trial = c(rep("t1",14),rep("t2",14),rep("t3",14)),
                 fetch = c(rep(1,14),c(rep(1,7),0,rep(1,3),0,1,0),rep(1,6),rep(0,4),1,rep(0,3)),
                 time = c(runif(14,50,150),runif(7,50,150),300,runif(3,120,280),300,runif(1,120,280),300,runif(6,50,150),rep(300,4),runif(1,120,280),rep(300,3)))

df$group <- factor(df$group, levels = c("A","B"))
df$trial <- factor(df$trial, levels = c("t1","t2","t3"))

So given this design, my hypothesis for testing if the gene curbs the reaction to the stimulus means that I want to compare the reaction times of the two groups accounting for the censoring issue that I have and for the animal-specific variation that exists across the trials. More specifically, in trial 2 are the response times different between the groups given what I observed in trial 1. And in trial 3, again, are the response times different between the groups given what I observed in trials 1 and 2.

I thought that using an accelerated failure time with a frailty function is a good way to go, but I'm not exactly sure how to set it up for the three trials and I'm also not sure if I'm encoding the fetch status properly.

This is what I'm doing so far:

aft.fit <- survreg(Surv(df$time, df$fetch) ~ df$group * df$trial + frailty(df$mouse_ID))

Which gives:

> summary(aft.fit)

survreg(formula = Surv(df$time, df$fetch) ~ df$group * df$trial + 
                   Value Std. Error      z       p
(Intercept)           4.8264     0.2789  17.31 < 2e-16
df$groupB            -0.2172     0.3969  -0.55    0.58
df$trialt2            0.1395     0.1002   1.39    0.16
df$trialt3            0.0129     0.1209   0.11    0.92
df$groupB:df$trialt2  0.8533     0.1876   4.55 5.4e-06
df$groupB:df$trialt3  1.4697     0.2382   6.17 6.8e-10
Log(scale)           -1.7579     0.1559 -11.27 < 2e-16

Scale= 0.172 

Weibull distribution
Loglik(model)= -152.2   Loglik(intercept only)= -199.3
    Chisq= 94.32 on 14.7 degrees of freedom, p= 1.1e-13 
Number of Newton-Raphson Iterations: 10 35 
n= 42 

Which seems to make sense, looking at the df$groupB:df$trialt2 and df$groupB:df$trialt3 lines.

But I was under the impression that in censored data 0 encodes alive and 1 dead, opposite from my fetch encoding.

If I try flipping my fetch encoding (i.e., df$fetch <- -1*(df$fetch-1)), trying to fit the mode above throws this error:

Error in solve.default(matrix(fit0$var, 1 + nstrat2)) : 
  Lapack routine dgesv: system is exactly singular: U[1,1] = 0

So my questions are:

  1. Is the group * trial interaction term the correct way to answer my question?

  2. Am I encoding the fetch censoring correctly?


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.