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