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:
set.seed(1) 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))
> summary(aft.fit) Call: survreg(formula = Surv(df$time, df$fetch) ~ df$group * df$trial + frailty(df$mouse_ID)) 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
But I was under the impression that in censored data
dead, opposite from my
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:
group * trialinteraction term the correct way to answer my question?
Am I encoding the fetch censoring correctly?