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.

Author: Zachary Huang

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

  1. 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 *