Effect of wound on HR and PLUD
Jean-Baptiste Ferdy
25/01/2022
Survival analysis
Plotting survival curves for each DGRP line, wounded or not wounded.
require(survival)## Loading required package: survivalsymbol <- c(3,21:24)
names(symbol) <- genotype_wPBS
msurv <- survfit(Surv(t,Censor)~factor(Line,levels=genotype_woPBS)+Treatment,
data=subset(d_surv,Line %in% genotype_woPBS))
dmsurv <- data.frame(Line=rep(rep(genotype_woPBS,each=2),msurv$strata),
Treatment=rep(rep(c("no_wound","wound"),length(genotype_woPBS)),msurv$strata),
time=msurv$time,proportion=msurv$surv)
dmsurv <- rbind(data.frame(Line=rep(genotype_woPBS,each=2),
Treatment=rep(c("no_wound","wound"),length(genotype_woPBS)),
time=0,proportion=1),dmsurv)
par(ann=FALSE)
plot(proportion~time,data=dmsurv,type="n",las=1)
foo <- sapply(split(dmsurv,with(dmsurv,paste(Line,Treatment))),function(xx) lines(proportion~time,
data=xx,type="b",
cex=0.75,pch=symbol[as.character(Line)],
lty=ifelse(xx$Line=="PBS",2,1),
col=ifelse(xx$Treatment=="wound","indianred","steelblue"),
bg=adjustcolor(ifelse(xx$Treatment=="wound","indianred","steelblue"),alpha.f=0.75)))
points(0,1,pch=21,col="gray",bg="gray",cex=1.25)
nw <- with(subset(d_surv,Treatment=="wound"),tapply(Censor,Line,length))
nnw <- with(subset(d_surv,Treatment=="No_wound"),tapply(Censor,Line,length))
l1 <- legend("topright",legend=paste(genotype_woPBS," (N=",nnw[genotype_woPBS],"/",nw[genotype_woPBS],")",sep=""),
pch=symbol[genotype_woPBS],bty="n")
legend(100,l1$rect$top,legend=c("No wound","Wound"),
fill=adjustcolor(c("steelblue","indianred"),alpha.f=0.75),bty="n",xjust=1)
title(xlab="Hours post injection")
title(ylab="Proportion surviving")
Estimating the effect of wound on survival
require(coxme)
m.hr.m <- list()
for(i in 1:4) {
m.hr.m[[i]] <- coxme(Surv(t,Censor)~Treatment+(1|Day_of_injection),
data=subset(d_surv,Line == genotype_woPBS[i] & Dose=="P. rettgeri" & Sex=="Male"))
print(m.hr.m[[i]])
}## Cox mixed-effects model fit by maximum likelihood
## Data: subset(d_surv, Line == genotype_woPBS[i] & Dose == "P. rettgeri" & Sex == "Male")
## events, n = 241, 299
## Iterations= 6 34
## NULL Integrated Fitted
## Log-likelihood -1227.375 -1089.561 -1085.118
##
## Chisq df p AIC BIC
## Integrated loglik 275.63 2.00 0 271.63 264.66
## Penalized loglik 284.51 2.93 0 278.65 268.44
##
## Model: Surv(t, Censor) ~ Treatment + (1 | Day_of_injection)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Treatmentwound 2.928119 18.69244 0.2006103 14.6 0
##
## Random effects
## Group Variable Std Dev Variance
## Day_of_injection Intercept 0.8307245 0.6901031
## Cox mixed-effects model fit by maximum likelihood
## Data: subset(d_surv, Line == genotype_woPBS[i] & Dose == "P. rettgeri" & Sex == "Male")
## events, n = 241, 266
## Iterations= 7 39
## NULL Integrated Fitted
## Log-likelihood -1161.342 -1049.04 -1046.615
##
## Chisq df p AIC BIC
## Integrated loglik 224.61 2.00 0 220.61 213.64
## Penalized loglik 229.46 1.98 0 225.50 218.61
##
## Model: Surv(t, Censor) ~ Treatment + (1 | Day_of_injection)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Treatmentwound 2.778515 16.0951 0.2128559 13.05 0
##
## Random effects
## Group Variable Std Dev Variance
## Day_of_injection Intercept 0.8250847 0.6807648
## Cox mixed-effects model fit by maximum likelihood
## Data: subset(d_surv, Line == genotype_woPBS[i] & Dose == "P. rettgeri" & Sex == "Male")
## events, n = 271, 289
## Iterations= 13 69
## NULL Integrated Fitted
## Log-likelihood -1310.727 -1198.527 -1193.785
##
## Chisq df p AIC BIC
## Integrated loglik 224.40 2.00 0 220.40 213.20
## Penalized loglik 233.88 2.95 0 227.99 217.38
##
## Model: Surv(t, Censor) ~ Treatment + (1 | Day_of_injection)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Treatmentwound 2.233361 9.331174 0.1685292 13.25 0
##
## Random effects
## Group Variable Std Dev Variance
## Day_of_injection Intercept 0.9435444 0.8902760
## Cox mixed-effects model fit by maximum likelihood
## Data: subset(d_surv, Line == genotype_woPBS[i] & Dose == "P. rettgeri" & Sex == "Male")
## events, n = 254, 328
## Iterations= 12 64
## NULL Integrated Fitted
## Log-likelihood -1321.697 -1206.79 -1202.571
##
## Chisq df p AIC BIC
## Integrated loglik 229.81 2.00 0 225.81 218.74
## Penalized loglik 238.25 2.91 0 232.42 222.11
##
## Model: Surv(t, Censor) ~ Treatment + (1 | Day_of_injection)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Treatmentwound 2.138708 8.488464 0.1551419 13.79 0
##
## Random effects
## Group Variable Std Dev Variance
## Day_of_injection Intercept 0.7270316 0.5285750m.hr <- coxme(Surv(t,Censor)~Treatment*Line+(1|Day_of_injection),
data=subset(d_surv,Line %in% genotype_woPBS & Dose=="P. rettgeri" & Sex=="Male"))
print(m.hr)## Cox mixed-effects model fit by maximum likelihood
## Data: subset(d_surv, Line %in% genotype_woPBS & Dose == "P. rettgeri" & Sex == "Male")
## events, n = 1007, 1182
## Iterations= 5 29
## NULL Integrated Fitted
## Log-likelihood -6431.988 -5954.866 -5946.318
##
## Chisq df p AIC BIC
## Integrated loglik 954.25 8.00 0 938.25 898.93
## Penalized loglik 971.34 10.79 0 949.76 896.73
##
## Model: Surv(t, Censor) ~ Treatment * Line + (1 | Day_of_injection)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Treatmentwound 2.5078069 12.2779736 0.1444134 17.37 0.0e+00
## Line584 0.7953847 2.2152931 0.1527533 5.21 1.9e-07
## Line630 0.4485291 1.5660070 0.1556473 2.88 4.0e-03
## Line818 0.2307055 1.2594882 0.1591877 1.45 1.5e-01
## Treatmentwound:Line584 -0.5793245 0.5602767 0.1875672 -3.09 2.0e-03
## Treatmentwound:Line630 -0.3007019 0.7402985 0.1922737 -1.56 1.2e-01
## Treatmentwound:Line818 0.4433085 1.5578528 0.1988239 2.23 2.6e-02
##
## Random effects
## Group Variable Std Dev Variance
## Day_of_injection Intercept 0.5784215 0.3345714The wound significantly decreases survival in all lines, but some are more impacted than others.
PLUD estimation from count data
Load data and estimate PLUDs
d_cfu <- read.table(file="data_CFU_wound.csv",header=TRUE)
require(CFUfit)## Loading required package: CFUfit## estimates load using a Poisson/binomial model
dPLUDest <- with(subset(d_cfu,BACTERIA>=0),loadEstimate(BACTERIA,1/(4*Dilution)
,sampleID=sample_ID
,volume=1 ## volume is 5 microliters but dividing dilution by four accounts for this
,maxCFU=50 ## counts above 50 should not be trusted
,ignoreAboveMaxCFU=FALSE
,logbase = 10))
rownames(dPLUDest) <- dPLUDest$sampleID
pos <- match(dPLUDest$sampleID,d_surv$sample_ID)
dPLUDest$Sex <- d_surv[pos,"Sex"]
dPLUDest$Line <- d_surv[pos,"Line"]; dPLUDest$Line <- factor(dPLUDest$Line)
dPLUDest$Plate <- d_surv[pos,"Plate"]
dPLUDest$Treatment <- d_surv[pos,"Treatment"]
dPLUDest$Dose <- d_surv[pos,"Dose"]
dPLUDest$Replicate <- d_surv[pos,"Replicate"]
dPLUDest$t <- d_surv[pos,"t"]
dPLUDest$Day_of_injection <- d_surv[pos,"Day_of_injection"]We use here the home-made package CFUfit to perform the estimation. This is because the experimental method we use to estimate load produces right truncated Poisson data: CFUare counted in 5 microliters droplet, and counts above 50 should be considered as unreliable. The function loadEstimate takes this specificity into account.
Outlier detection, using Rosner test
require(EnvStats)
dPLUDest$outlier <- FALSE
dPLUDest$grpe <- with(dPLUDest,factor(paste(Line,Treatment)))
for(g in levels(dPLUDest$grpe)) {
l <- dPLUDest$grpe == g & !is.na(dPLUDest$load)
test <- try(rosnerTest(dPLUDest$logload[l],k = floor(sum(l,na.rm=TRUE)/2)))
if("try-error" %in% class(test)) {
dPLUDest$outlier[l] <- NA
} else {
dPLUDest$outlier[which(l)[test$all.stats$Obs.Num[test$all.stats$Outlier]]] <- TRUE
}
}
with(dPLUDest,table(grpe,outlier))## outlier
## grpe FALSE TRUE
## 559 No_wound 29 7
## 559 wound 48 4
## 584 No_wound 74 0
## 584 wound 65 2
## 630 No_wound 44 1
## 630 wound 47 0
## 818 No_wound 40 0
## 818 wound 28 0(list_outliers <- rownames(dPLUDest)[dPLUDest$outlier])## [1] "B1_Plate11" "B2_Plate11" "B2_Plate63" "C2_Plate11" "D2_Plate11"
## [6] "D7_Plate63" "D9_Plate65" "E1_Plate11" "E11_Plate64" "F1_Plate11"
## [11] "F5_Plate62" "G1_Plate11" "G2_Plate11" "H1_Plate11"require(beanplot)## Loading required package: beanplotrequire(beeswarm)## Loading required package: beeswarmgr <- beanplot(logload~Treatment+factor(Line,levels=genotype_woPBS),
ylab="",#ylab="log10(PLUD)",
col=list(adjustcolor("steelblue",alpha.f=0.75),adjustcolor("indianred",alpha.f=0.75)),
what=c(F,T,T,F),
axes=FALSE,
ylim=c(2,8),
data=subset(dPLUDest,Sex=="Male" & Line %in% genotype_woPBS & !outlier))
axis(1,2*(1:4)-0.5,gsub("wound.","",gr$names[c(F,T)]))
axis(2,with(dPLUDest,pretty(logload)),las=1)
box()
pos <- 1:8; names(pos) <- gr$names
with(subset(dPLUDest, Line %in% genotype_woPBS),
beeswarm(logload~pos[paste(Treatment,Line,sep=".")],
pch=21,
cex=0.5,
pwbg=ifelse(outlier,"red","white"),
add=T))
n <- with(subset(dPLUDest,!is.na(logload) & Sex=="Male" & Line %in% genotype_woPBS),tapply(logload,factor(Line,levels=genotype_woPBS),length))
mtext(text=paste0("(N=",as.numeric(n),")"),side=3,line=-1,at=2*(1:4)-0.5)
title(ylab="log PLUD",srt=90)
title(xlab="D. melanogaster DGRP line")
legend("bottomleft",legend=c("No wound","Wound"),fill=adjustcolor(c("steelblue","indianred"),alpha.f=0.75),bty="n",xjust=1)
Outliers with particularly low values of PLUD might correspond to flies which died but not from the infection.
Testing the wound effect on PLUD
A first (crude) approach, using a Wilcoxon test
dspl <- split(dPLUDest,dPLUDest$Line)
## outliers included
sapply(dspl,function(dtmp) wilcox.test(logload~Treatment,data=dtmp,alternative="less"))## Warning in wilcox.test.default(x = c(7.00858977779907, 6.51850354591504, :
## cannot compute exact p-value with ties## Warning in wilcox.test.default(x = c(6.67208746397287, 6.16434246182159, :
## cannot compute exact p-value with ties## 559
## statistic 595
## parameter NULL
## p.value 0.001927058
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 584
## statistic 2205
## parameter NULL
## p.value 0.1294069
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 630
## statistic 534
## parameter NULL
## p.value 2.201878e-05
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 818
## statistic 495
## parameter NULL
## p.value 0.2107596
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"## outliers removed
sapply(dspl,function(dtmp) wilcox.test(logload~Treatment,data=subset(dtmp,!outlier),alternative="less"))## Warning in wilcox.test.default(x = c(6.35697100703029, 6.61277346275689, :
## cannot compute exact p-value with ties## Warning in wilcox.test.default(x = c(7.00858977779907, 6.51850354591504, :
## cannot compute exact p-value with ties## Warning in wilcox.test.default(x = c(6.67208746397287, 6.16434246182159, :
## cannot compute exact p-value with ties## 559
## statistic 470
## parameter NULL
## p.value 0.008874682
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 584
## statistic 2057
## parameter NULL
## p.value 0.0711882
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 630
## statistic 534
## parameter NULL
## p.value 3.637232e-05
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"
## 818
## statistic 495
## parameter NULL
## p.value 0.2107596
## null.value 0
## alternative "less"
## method "Wilcoxon rank sum test with continuity correction"
## data.name "logload by Treatment"The PLUD is always greater in wounded flies, but the difference is significant only for Line 630 and 559. This approach is crude first because it ignores differences in precision among load estimates, and second because it does not consider random fluctuations among experiments. Excluding the outliers does not change much to this conclusion.
A global analysis using spaMM
require(spaMM)## Loading required package: spaMM## Registered S3 methods overwritten by 'registry':
## method from
## print.registry_field proxy
## print.registry_entry proxy## spaMM (Rousset & Ferdy, 2014, version 3.9.62) is loaded.
## Type 'help(spaMM)' for a short introduction,
## 'news(package='spaMM')' for news,
## and 'citation('spaMM')' for proper citation.
## Further infos, slides, etc. at https://gitlab.mbb.univ-montp2.fr/francois/spamm-ref.require(DHARMa)## Loading required package: DHARMa## This is DHARMa 0.4.5. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')m_r <- fitme(logload~Treatment*Line+(1|Day_of_injection),
prior.weights=1/(se.logload^2), ## load estimates with large se count less in model fitting
#resid.model = grpe, ## model heteroskedasticity, not usefull here
data=subset(dPLUDest,!outlier))
sim <- simulateResiduals(fittedModel = m_r,refit = TRUE) ## refit is advised for random models, if I got it wright## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.plot(sim)
with(subset(dPLUDest,!outlier),plotResiduals(sim, form = grpe))
Here we use a model with a random effect that reproduces fluctuations among experiment replicates, and with a weights computed as te inverse of the variance in load estimates. This weight dicreases the influence of samples for which load is hard to estimate.
The package DHARMa (https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html) is used to test gooness of fit. As the model includes a random effect, we use the option refit=TRUE which refits the model for each new simulated dataset. In the absence of random effects, the authors advise to use the default refit=FALSE see here for details on this issue.
m_r_noint <- update(m_r,.~.-Treatment:Line)
m_r_noT <- update(m_r_noint,.~.-Treatment)
## test the interaction
anova(m_r,m_r_noint)## chi2_LR df p_value
## p_v 1.257257 3 0.739307## test the Treatment effect
anova(m_r_noint,m_r_noT)## chi2_LR df p_value
## p_v 11.29718 1 0.0007762495# > anova(m_r,m_r_noint)
# chi2_LR df p_value
# p_v 1.257257 3 0.739307
# > anova(m_r_noint,m_r_noT)
# chi2_LR df p_value
# p_v 11.29718 1 0.0007762495The wound has a significant effect on PLUD, with no significant interaction between wound and Line.
We now analyze the dataset with outliers included, in order to test for their effect.
m_r_wo <- update(m_r,data=dPLUDest)
sim <- simulateResiduals(fittedModel = m_r_wo,refit = TRUE)## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.plot(sim) ## DHARMa nicely detects the presence of outliers !
with(dPLUDest,plotResiduals(sim, form = grpe))
m_r_wo_noint <- update(m_r_wo,.~.-Treatment:Line)
m_r_wo_noT <- update(m_r_wo_noint,.~.-Treatment)
## test the interaction
anova(m_r_wo,m_r_wo_noint)## chi2_LR df p_value
## p_v 1.902321 3 0.5929256## test the Treatment effect
anova(m_r_wo_noint,m_r_wo_noT)## chi2_LR df p_value
## p_v 11.04414 1 0.0008896792The results are almost unchanged.
PLUD as a function of HR
Extraction of HR estimates and computation of CI. The coxme class has no simulate function, so that confidence interval must be computed from estimates and standard error.
est <- fixef(m.hr)
pos <- grep("Treatmentwound",names(est))
est <- est[pos]
est[-1] <- est[-1]+est[1] ## the first estimate is that of the reference line; the following are interaction terms
vest <- diag(vcov(m.hr))[pos]
vest[-1] <- vest[-1]+vest[1] ## variance of the sum is sum of variances... not quite, because estimates are not independent!
vest <- sqrt(vest)
hr_est <- data.frame(est=est,
inf=est-qnorm(0.975)*vest,
sup=est+qnorm(0.975)*vest)
hr_est$Line <- factor(levels(d_surv$Line))
plot(est~Line,data=hr_est,ylim=c(1.5,3.5),
medcol=NA,boxcol=NA,staplecol=NA,
xlab="Line",ylab="log HR (wound relative to noWound)")
with(hr_est,segments(1:4,inf,1:4,sup))
with(hr_est,points(1:4,est,pch=21,bg="white"))
We shall now extract log PLUD ratios (and corresponding CI) from the global model. The class spaMM has a simulate function, so that CI can be estimated by boostraping. The estimates are in decimal log scale while the HR is in plain log scale. This comes from the option logbase=10 in loadEstimate. Rather than changing this option (which gives violin plots that are easier to interpret, I think) I change the scale of estimates in the computation of CIs.
ci <- list()
ci[[1]] <- confint(m_r,parm=function(v) fixef(v)[2]*log(10), boot_args=list(nsim=99, seed=123))## Bootstrap replicates: bootstrap took 3.27 s.
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 99 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = list(99, "parametric"), conf = 0.95, type = c("basic",
## "perc", "norm"), t0 = 0.282564537061954, `str(long numeric 't')` = " num [1:99] 0.135 0.282 0.292 0.558 0.438 ...")
##
## Intervals :
## Level Normal Basic Percentile
## 95% (-0.0925, 0.5965 ) (-0.0409, 0.6425 ) (-0.0773, 0.6060 )
## Calculations and Intervals on Original Scale
## Some basic intervals may be unstable
## Some percentile intervals may be unstableci[[2]] <- confint(m_r,parm=function(v) (fixef(v)[2]+fixef(v)[6])*log(10), boot_args=list(nsim=99, seed=123))## Bootstrap replicates: bootstrap took 3.92 s.
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 99 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = list(99, "parametric"), conf = 0.95, type = c("basic",
## "perc", "norm"), t0 = 0.216920087369519, `str(long numeric 't')` = " num [1:99] 0.2814 -0.0806 0.0676 0.1815 0.2524 ...")
##
## Intervals :
## Level Normal Basic Percentile
## 95% (-0.0147, 0.4639 ) (-0.0272, 0.4944 ) (-0.0606, 0.4610 )
## Calculations and Intervals on Original Scale
## Some basic intervals may be unstable
## Some percentile intervals may be unstableci[[3]] <- confint(m_r,parm=function(v) (fixef(v)[2]+fixef(v)[7])*log(10), boot_args=list(nsim=99, seed=123))## Bootstrap replicates: bootstrap took 3.02 s.
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 99 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = list(99, "parametric"), conf = 0.95, type = c("basic",
## "perc", "norm"), t0 = 0.438744786950566, `str(long numeric 't')` = " num [1:99] 0.643 0.388 0.621 0.149 0.657 ...")
##
## Intervals :
## Level Normal Basic Percentile
## 95% ( 0.0841, 0.7577 ) ( 0.0556, 0.7666 ) ( 0.1109, 0.8219 )
## Calculations and Intervals on Original Scale
## Some basic intervals may be unstable
## Some percentile intervals may be unstableci[[4]] <- confint(m_r,parm=function(v) (fixef(v)[2]+fixef(v)[8])*log(10), boot_args=list(nsim=99, seed=123))## Bootstrap replicates: bootstrap took 3.4 s.
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 99 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = list(99, "parametric"), conf = 0.95, type = c("basic",
## "perc", "norm"), t0 = 0.20909083659694, `str(long numeric 't')` = " num [1:99] -0.0465 -0.2548 0.2384 -0.1641 0.719 ...")
##
## Intervals :
## Level Normal Basic Percentile
## 95% (-0.2393, 0.6833 ) (-0.2804, 0.6925 ) (-0.2743, 0.6986 )
## Calculations and Intervals on Original Scale
## Some basic intervals may be unstable
## Some percentile intervals may be unstableplud_est <- data.frame(est=sapply(ci,function(z) z$t0),
inf=sapply(ci,function(z) z$normal[2]),
sup=sapply(ci,function(z) z$normal[3]))
plud_est$Line <- levels(dPLUDest$Line)
plud_est <- plud_est[match(plud_est$Line,hr_est$Line),]plot(hr_est$est,plud_est$est,
xlim=c(0,3.5),ylim=c(-0.1,0.75),type="n",las=1,
xlab="log hazard ratio",ylab="log PLUD ratio")
abline(h=0,lty=2)
abline(v=0,lty=2)
segments(hr_est$est,plud_est$inf,hr_est$est,plud_est$sup,col="indianred")
segments(hr_est$inf,plud_est$est,hr_est$sup,plud_est$est,col="indianred")
points(hr_est$est,plud_est$est,pch=21,col="white",bg="white",cex=2)
points(hr_est$est,plud_est$est,pch=21,
col="indianred",bg=adjustcolor("indianred",alpha.f=0.75))
points(0,0,pch=21,col="white",bg="white",cex=2)
points(0,0,pch=21,bg=adjustcolor("steelblue",alpha.f=0.75),col="steelblue")
for(i in 1:4) {
text(hr_est$est[i],plud_est$est[i],plud_est$Line[i],adj=c(1.5,1.5))
}
text(0,0,"no wound",adj=c(-0.25,-0.5))