Adjusting control mortality in LC50/LD50 determination using R

Problem:
There are R codes for calculating LC50 or LD50 using 4 to 6 data points. However, none uses control mortality correction (e.g. Abbott correction). There is an ecotoxicology package which does this, but it seems hugely complex to use.

Solution:
In the code below, we put the last line as the control mortality (dose=0), which is used to calculate the correction
Adjusted mortality = (mortality – control.mortality)/ (1 – mortality)

Then the last line of data was dropped by using subset.

Data structure:

dose 	dead 	alive
1000.00   15     0
500.00   13     2
250.00   12     5
125.00    8     8
62.50    6     9
31.25    3    12
0.00    2    13

New data after mortality correction:

the adjusted data would have this
     dose      dead     alive
1 1000.00 15.000000  0.000000
2  500.00 12.692308  2.307692
3  250.00 11.230769  5.769231
4  125.00  6.769231  9.230769
5   62.50  4.615385 10.384615
6   31.25  1.153846 13.846154

Code:

data = read.csv("LD50-2.csv")
c=nrow(data)  
N=data$dead+data$alive

mort=data$dead/N
ad.mort=(mort-mort[c])/(1-mort[c])
data$dead=N*ad.mort
data$alive=N-data$dead

data=subset(data, dose!=0)

y = cbind(data$alive,data$dead)
model.results = glm(y~log(data$dose+1),binomial)
summary(model.results)
library(MASS)

ld50=exp(dose.p(model.results,p=0.5))
ld50

The calculated LD50 here is 141.7363 +- 0.1774705 (SE).
The unit will be the same as the dose in the first file. i.e. if it is mg, then it is 141.7 mg here.

One thought on “Adjusting control mortality in LC50/LD50 determination using R

  1. Stephen Young

    try this script with the same data:
    cf<-read.csv("abbotseg.csv")
    attach(cf)
    str(cf)
    resp<-cbind(dead,alive)
    ####searches for best x offest value
    ac<-vector()
    for (i in 1:1000)
    {
    m1<-glm(resp~log(dose+.2*i),binomial)
    ac<-c(ac,AIC(m1))
    }
    xoff<-.2*which.min(ac)
    #########fits binomial glm, including dose=0
    m2<-glm(resp~log(dose+xoff),binomial)
    ####calculates ld50
    cc<-coef(m2)
    ld50<-exp(-cc[1]/cc[2])-xoff
    ######plots data and model curve
    ppn<-dead/(dead+alive)
    plot(ppn~log(dose+xoff),xaxt="n",xlab="Dose",ylab="Proportion dead")
    xx<-c(0,20,50,100,200,500,1000)
    axis(1,at=log(xx+xoff),labels=xx)
    xp<-seq(4.6,7.1,.01)
    pr<-xp*cc[2]+cc[1]
    py<-exp(pr)/(1+exp(pr))
    lines(xp,py)
    ###adds ld50 to graph
    points(log(ld50+xoff),.5,col="red",cex=2)

    It lets the glm do its stuff with raw proportions – no warnings!, and leaves the "control mortality" 0 dose in the model, by choosing an optimumx offest.

Leave a Reply

Your email address will not be published. Required fields are marked *