In ahproctor/LC50: Estimate LC50 in the presence of additional stressors
Intro
Model
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.
Example
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)
Graphical Exploration
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)
Model fitting
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)
Multiple Comparisons
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)
Bayesian Estimates
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)
Shrinkage
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)
LCx
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)
ahproctor/LC50 documentation built on May 12, 2019, 2:31 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Intro
Model
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.
Example
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)
Graphical Exploration
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)
Model fitting
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
concentrationargument names the variable in the dataframe recording the concentration of the toxin, and -
the
groupargument 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)
Multiple Comparisons
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)
Bayesian Estimates
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)
Shrinkage
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)
LCx
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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.