Post date: Apr 10, 2017 3:25:14 AM
Patrik ran a desiccation experiment (March/April 2017) to test whether green and dark (melanic) bugs differ in terms of survival times from desiccation.
The data from 32 trials are here: /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_fluct/timemaMelExp/dessication_trials_2017.csv
I first fit a simple linear model to these data, then I used the more appropriate Cox proportional hazards regression model (using the 'survival' package in R, see background here). I ran this with two methods for dealing with ties (death at the same time): the Efron option and the exact partial likelihood method. The latter is preferred when time is discrete (true here, as the data are collected at 20 minute intervals) and when there are only a few observed times for death (also true here, depending on how you define 'a few'), but the Efron method is more common. But, it doesn't really matter as the general conclusion is the same either way.
Efron:
exp(B) = 2.58, p = 0.0163, 95% CIs = 1.19-5.61
Exact likelihood
exp(B) = 3.57, p = 0.0111, 95% CIs = 1.34-9.51
Note that exp(B) > 1 means melanic bugs die faster and that exp(B), that is e^B, is the more natural parameter for the model than B (the raw regression coefficient).
Here is the full code from commands.R:
## comparison desiccation tolerance for melanic vs. green timema morphs
dat<-read.csv("dessication_trials_2017.csv")
## first as a simple linear model
o<-lm(dat$death_time_as_last_time_measure ~ dat$mel_yn)
summary(o)
#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 136.67 11.58 11.800 1.03e-11 ***
#dat$mel_yn -35.33 15.54 -2.274 0.0318 *
#Residual standard error: 40.12 on 25 degrees of freedom
# (5 observations deleted due to missingness)
#Multiple R-squared: 0.1714, Adjusted R-squared: 0.1382
#F-statistic: 5.171 on 1 and 25 DF, p-value: 0.03182
## convert to survival data
tm<-seq(0,180,20)
sdat<-as.matrix(dat[,13:22])
time1<-rep(0,32)
time2<-rep(NA,32)
for(i in 1:32){
if(sum(sdat[i,]==0)>0){ ## death
time2[i]<-min(tm[which(sdat[i,]==0)])
}
else{
time2[i]<-180
}
}
## basic cox proportional hazards model, assumes constant effect of covariate
sob<-Surv(time=time1,time2=time2,event=abs(1-dat$alive_at_end))
o<-coxph(sob ~ dat$mel_yn)
#n= 32, number of events= 27
# coef exp(coef) se(coef) z Pr(>|z|)
#dat$mel_yn 0.9494 2.5842 0.3951 2.403 0.0163 *
# exp(coef) exp(-coef) lower .95 upper .95
#dat$mel_yn 2.584 0.387 1.191 5.606
#Concordance= 0.667 (se = 0.065 )
#Rsquare= 0.165 (max possible= 0.992 )
#Likelihood ratio test= 5.77 on 1 df, p=0.01628
#Wald test = 5.77 on 1 df, p=0.01627
#Score (logrank) test = 6.17 on 1 df, p=0.01299
## exact ml method for tie breaking
o<-coxph(sob ~ dat$mel_yn,method="exact")
# n= 32, number of events= 27
# coef exp(coef) se(coef) z Pr(>|z|)
#dat$mel_yn 1.2713 3.5656 0.5006 2.54 0.0111 *
# exp(coef) exp(-coef) lower .95 upper .95
#dat$mel_yn 3.566 0.2805 1.337 9.511
#Rsquare= 0.196 (max possible= 0.96 )
#Likelihood ratio test= 6.97 on 1 df, p=0.008272
#Wald test = 6.45 on 1 df, p=0.01109
#Score (logrank) test = 6.99 on 1 df, p=0.008204
## plot of fitted survival with kaplan-meier, we used this the the mating trials in the Nature E&E paper
ofit<-survfit(sob ~ dat$mel_yn,type="kaplan-meier")
pdf("survPlot.pdf",width=5,height=5)
par(mar=c(5,5,0.5,0.5))
plot(ofit,col=c("forestgreen","brown"),cex.lab=1.3,lwd=1.5,xlab="Time",ylab="Survival proportion")
legend(7,0.18,c("green morph","melanic morph"),col=c("forestgreen","brown"),bty='n',lty=1,lwd=1.5)
dev.off()