Lethal concentrations LCx for each combination of additional stressors are estimated by fitting a Binomial regression model that relates survival to log toxin concentration, adjusting for background mortality and assuming a linear relation between log LCx and descriptors of the additional stressors.
More precisely, the a Binomial regression model parameterized in terms of log LCx is used to relate survival to log toxin concentration [ \begin{align} y_{ijk} &\sim \operatorname{Bin}(N_{ijk},\pi_{ijk})\ \pi_{ijk} &= \begin{cases} p_{ijk}q_{k} & C_{ijk} > 0\ q_{k} & C_{ijk}= 0 \end{cases}\ \ell_{1}(p_{ijk})&= \alpha_{k} (\log C_{ijk} - \log\mathrm{LCx}{k}) + z\ \log \mathrm{LCx}{k} &= \beta_{0} + \beta_{1} x_{1k} + \dots + \beta_{m}x_{mk}\ \ell_{2}(q_{k}) &= \gamma_{k} \end{align} ]
Here (y_{ijk}) is the number of survivors at the end of the trial from (N_{ijk}) individuals in replicate (i) receiving target toxin concentration (j) and combination (k) of additional stressors and actual measured concentration (C_{ijk}), (\pi_{ijk}) is the fraction of individuals expected to survive which is decomposed into a fraction (q_{i}) of individuals expected to survive in the absence of the toxin, and a fraction (p_{ijl}) expected to survive the additional mortality induced by the toxin. The fraction (p_{ij}) is related through a monotonic link function (\ell_{1}) to the log measured toxin concentration (\log C_{ijk}), the log lethal concentration (\log \mathrm{LCx}{i}) and a rate parameter (\alpha{i}) specific to the combination of stressors. Similarly, the baseline survival fraction (q_{i}) is related to a survival parameter (\gamma_{i}) specific to the combination of stressors through a monotonic link function (\ell_{2}). The lethal concentration (\mathrm{LCx}{i}) specific to a stressor combination is related through a log link to a linear predictor constructed from covariates (x{1},x_{2},\ldots,x_{m}) describing the combination of additional stressors. The level of lethality is modelled is determined by the offset (z = \ell_{1}(x/100)), (that is x=50 if LCx=LC50).
The baseline survival parameters (\gamma_{i}) and rate parameters (\alpha_{i}) are viewed as nuisance parameters, and inference focuses on the effects of the additional stressors and the LCx lethal concentrations.
The toxicity
dataset is simulated data showing survival following
exposure to a known toxin in the presence of the additional stressors
temperature and salinity. In the imagined scenario, samples of an
aqueous solution of a toxin, of varying concentrations, are prepared
for each of three salinities and are held at one of three constant
temperatures. Approximately 10 individuals were added to each sample
and the number of survivors recorded at the end of four days exposure.
library(LC50) head(toxicity)
Plotting the survival fraction against concentration for each combination of additional stressors shows that survival does decrease with increasing concentrations of the toxin, and in many samples there is a non-negligible level of control mortality
library(ggplot2) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ facet_grid(temperature~salinity)
The LCx lethal concentrations for each temperature and salinity
combination can be estimated with lcx
. This function fits the
Binomial model described in the previous section by maximum
likelihood. The details of the model are specified by three
arguments:
a formula, the left hand side of which is a two column matrix constructed from the number of survivors and mortalities in each sample, while the right hand side specifies the terms in the linear predictor for log LCx
the concentration
argument names the variable in the dataframe
recording the concentration of the toxin, and
the group
argument names a factor in the dataframe that
distinguishes the combinations of additional stressors. The model
requires that the terms of linear predictor for log LCx are
constant within levels of this factor.
fit <- lcx(cbind(alive,dead)~salinity*temperature,concentration=conc,group=group,data=toxicity) summary(fit)
By default, separate background survival rates are estimated for each group, and the LC50 lethal concentration is modelled using a probit link.
An analysis of deviance table for the fit is constructed with anova
.
anova(fit,test="Chisq")
This table shows no evidence that log LC50 depends on the interaction, and the model can be simplified to a main effects only model for log LC50. Although this fits a restricted model for log LC50, it still fits a separate control parameter (\gamma_{i}) and rate parameter (\alpha_{i}) for each group.
fit <- lcx(cbind(alive,dead)~salinity+temperature,concentration=conc,group=group,data=toxicity) summary(fit)
The table of lethal concentrations can be extracted from the summary and displayed graphically
library(gridExtra) d.lc50 <- with(summary(fit), cbind(lcx, toxicity[match(rownames(lcx),toxicity$group), c("group","salinity","temperature")])) grid.arrange( ggplot(d.lc50,aes(x=group,y=LC50))+ geom_linerange(mapping=aes(x=group,ymin=`Lwr 95%`,ymax=`Upr 95%`),col="dodgerblue1",lwd=3)+ geom_point(), ggplot(d.lc50,aes(x=salinity,y=LC50,group=temperature,colour=temperature))+ geom_point()+ geom_line()+ geom_errorbar(mapping=aes(ymin=`Lwr 95%`,ymax=`Upr 95%`),width=0.05), ggplot(d.lc50,aes(x=temperature,y=LC50,group=salinity,colour=salinity))+ geom_point()+ geom_line()+ geom_errorbar(mapping=aes(ymin=`Lwr 95%`,ymax=`Upr 95%`),width=0.05), layout_matrix=rbind(1,2:3))
Alternately the effects package can be used to construct a graphical summary of the relation between LC50 and the additional stressors.
library(effects) plot(Effect(c("salinity","temperature"),fit))
Predictions from the fitted model can be generated with predict
.
d.pr <- expand.grid(temperature=factor(c(12,14,16)), salinity=factor(c(31,33,35)), conc=0:90) d.pr$group <- with(d.pr,interaction(temperature,salinity,sep="|")) d.pr$survival <- predict(fit,d.pr) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ geom_line(aes(x=conc,y=survival,colour=group),data=d.pr)+ geom_vline(aes(xintercept=LC50,colour=group),data=d.lc50)+ geom_rect(aes(xmin=`Lwr 95%`,xmax=`Upr 95%`,ymin=0,ymax=1,fill=group), inherit.aes=FALSE,data=d.lc50,alpha=0.2,color=NA)+ facet_grid(temperature~salinity)
By default the model is fit with a standard Binomial likelihood. When there is evidence that the data are overdispersed, a quasi-Binomial model can be fitted by specifying quasi=TRUE
fitq <- lcx(cbind(alive,dead)~salinity+temperature,concentration=conc,group=group,data=toxicity,quasi=TRUE) summary(fitq)
Alternately, confidence intervals for the lethal concentrations can be calculated by parametric bootstrap - new data sets are simulated from the fitted model and the model refitted to the simulated data, confidence intervals are then derived from the variability seen in the estimates from the simulated data.
d.boot <- toxicity boot <- lapply(simulate(fit,100),function(s) { ## Update data d.boot[,colnames(s)] <- s; ## Refit model and extract lc50s update(fit,data=d.boot,start=fit)}) ## Summaries of bootstrap estimates lc50 <- sapply(boot,function(fit) exp(fit$loglcx)) t(apply(lc50,1,function(lc50) c(LC50=mean(lc50),quantile(lc50,c(0.025,0.975)))))
Confidence intervals for predictions can be generated similarly
pr <- sapply(boot,function(fit) predict(fit,d.pr)) d.pr <- cbind(d.pr,t(apply(pr,1,quantile,c(0.025,0.975)))) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ geom_ribbon(aes(x=conc,y=survival,ymin=`2.5%`,ymax=`97.5%`,fill=group),data=d.pr,alpha=0.2,color=NA)+ geom_line(aes(x=conc,y=survival,colour=group),data=d.pr)+ facet_grid(temperature~salinity)
Multiplicity adjusted pairwise comparison tests can be constructed with
glht
from the multcomp package.
library(multcomp) mc <- glht(fit,linfct=mcp(salinity="Tukey",temperature="Tukey")) summary(mc)
Multiplicity adjusted confidence intervals for
differences in log LC50 can be obtained with confint
confint(mc)
These intervals can be converted to confidence intervals for ratios of LC50s by taking exponentials
exp(confint(mc)$confint)
If JAGS and the rjags package are installed, Bayesian
estimates of LC50 can be obtained with lcxJAGS
. This function is
an analog of lcx
that constructs an object of class jags
that can
be used to sample from the posterior for the model using the
facilities of the rjags package.
To generate Bayesian estimates for the reduced model, lcxJAGS
is
used to generate a model object
library(rjags) model <- lcxJAGS(cbind(alive,dead)~salinity+temperature,concentration=conc,group=group,data=toxicity)
and an initial burn-in sample is drawn
update(model,1000)
A final sample is then drawn and stored. It is possible to sample for any of the model variables "alpha", "beta", "gamma", "lcx", "loglcx", "p" or "q" as required.
s <- coda.samples(model,c("beta","q","lcx"),n.iter=10000,n.thin=10) summary(s)
In trials with few distinct concentrations, the rate parameters
(\alpha_{i}) can be poorly resolved. When these rates are
overestimated the confidence intervals for the LC50 can be too
short. The rate.shrink
parameter adds a penalty to the likelihood to
shrink the (\alpha_{i}) towards zero.
Without shrinkage, the estimates of some of the rate parameters for the toxicity data are particularly high, and the confidence intervals for the corresponding LC50s are not very wide.
fit <- lcx(cbind(alive,dead)~salinity+temperature,concentration=conc,group=group,data=toxicity) summary(fit,rate=TRUE)
d.pr <- expand.grid(temperature=factor(c(12,14,16)), salinity=factor(c(31,33,35)), conc=0:90) d.pr$group <- with(d.pr,interaction(temperature,salinity,sep="|")) d.pr$survival <- predict(fit,d.pr) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ geom_line(aes(x=conc,y=survival,colour=group),data=d.pr)+ #geom_vline(aes(xintercept=LC50,colour=group),data=d.lc50)+ #geom_rect(aes(x=Estimate,y=0,xmin=`Lwr 95%`,xmax=`Upr 95%`,ymin=0,ymax=1,fill=group),data=d.lc50,alpha=0.2,color=NA)+ facet_grid(temperature~salinity)
With shrinkage,
fit <- lcx(cbind(alive,dead)~salinity+temperature,concentration=conc,group=group,data=toxicity,rate.shrink=0.001) summary(fit,rate=TRUE)
d.pr <- expand.grid(temperature=factor(c(12,14,16)), salinity=factor(c(31,33,35)), conc=0:90) d.pr$group <- with(d.pr,interaction(temperature,salinity,sep="|")) d.pr$survival <- predict(fit,d.pr) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ geom_line(aes(x=conc,y=survival,colour=group),data=d.pr)+ #geom_vline(aes(xintercept=LC50,colour=group),data=d.lc50)+ #geom_rect(aes(x=Estimate,y=0,xmin=`Lwr 95%`,xmax=`Upr 95%`,ymin=0,ymax=1,fill=group),data=d.lc50,alpha=0.2,color=NA)+ facet_grid(temperature~salinity)
By default all calculations are performed based on LC50, but it is possible to perform the analysis based on other lethal concentrations.
To repeat the basic analysis above for LC10
fit <- lcx(cbind(alive,dead)~salinity*temperature,concentration=conc,group=group,data=toxicity,lethal=10) summary(fit)
Changing the reference lethality impacts inferences - for the LC10 there is strong evidence of a salinity by temperature interaction
anova(fit,test="Chisq")
d.lcx <- with(summary(fit), cbind(lcx, toxicity[match(rownames(lcx),toxicity$group), c("group","salinity","temperature")])) ggplot(d.lcx,aes(x=group,y=LC10))+ geom_linerange(mapping=aes(x=group,ymin=`Lwr 95%`,ymax=`Upr 95%`),col="dodgerblue1",lwd=3)+ geom_point()
d.pr <- expand.grid(temperature=factor(c(12,14,16)), salinity=factor(c(31,33,35)), conc=0:90) d.pr$group <- with(d.pr,interaction(temperature,salinity,sep="|")) d.pr$survival <- predict(fit,d.pr) ggplot(toxicity,aes(x=conc,y=alive/total,colour=group))+ geom_point()+ geom_line(aes(x=conc,y=survival,colour=group),data=d.pr)+ geom_vline(aes(xintercept=LC10,colour=group),data=d.lcx)+ geom_rect(aes(x=Estimate,y=0,xmin=`Lwr 95%`,xmax=`Upr 95%`,ymin=0,ymax=1,fill=group),data=d.lcx,alpha=0.2,color=NA)+ facet_grid(temperature~salinity)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.