Dear all,
I am trying to apply the logistic regression to determine the limit of
detection (LOD) of a molecular biology assay, the polymerase chain reaction
(PCR). The aim of the procedure is to identify the value (variable
"dilution") that determine a 95% probability of success, that is
"positive"/"total"=0.95. The procedure I have implemented seemed to work
looking at the figure obtained from the sample set 1; however the figure
obtained from the sample set 2 shows that interpolation is not correct.
Admittedly, these are not the best sample sets to determine the LOD -there
are too many 100% successes - but the regression should still work.
Does anybody have some tip to solve this incongruence?
Many thanks,
Luigi Marongiu, MSc
----------------------------------------------------------------------------
----------------
### SAMPLE SET No. 1
### define data
labs<-c("dilution", "total", "positive")
p.1<-c(10, 20, 19)
p.2<-c(100, 20, 20)
p.3<-c(1000, 20, 20)
p.4<-c(10000, 20, 20)
p.5<-c(100000, 20, 20)
p.6<-c(1000000, 20, 20)
### create data frame from matrix
my.mat<-matrix(c(p.1, p.2, p.3, p.4, p.5, p.6), nrow=6, byrow=TRUE)
dimnames(my.mat)<-list(c(1:6),labs)
my.data<-as.data.frame(my.mat)
attach(my.data)
my.data
### define frequency data
Y<-cbind(positive, total-positive)
### fit model LOGIT
model<-glm(Y ~ dilution, family=binomial(link=logit), data=my.data)
### plot data and regression line
plot(dilution, positive/total, pch=16, ylim=c(0,1), log = "x",
xaxt="n", xlim=c(1, 1000000),
main="Positive amplification by PCR", ylab="Fraction amplified",
xlab=expression(paste("Standard dilutions (c/",mu,"l)")))
### add x axis (William Dunlap, personal communication)
xlim <- par("usr")[1:2]
lab.x <- seq(ceiling(xlim[1]), floor(xlim[2]))
axis(side=1, at=10^lab.x, lab=as.expression(lapply(lab.x,
function(x)bquote(10^.(x)))))
lines(dilution, fitted(model), lty=1)
### set significance level and logit
# define significance
level<-0.95
# set logit value
logit.LOD<-log(level/(1-level))
### set function to determine LOD (Vito Ricci, personal communication)
LOD<-function(mod) as.vector((logit.LOD-coef(mod)[1])/coef(mod)[2])
### computation of the LOD value
x.LOD<-LOD(model)
### compute S.E.
# set variables
co<-model$coef
model.sum<-summary(model)
SE.co<-model.sum$coef[,2]
COV.co<-model.sum$cov.scaled[1,2]
# compute SE
SE.LOD<-abs(x.LOD) * sqrt( (SE.co[1]/co[1])^2 +
(SE.co[2]/co[2])^2 -
2*(COV.co/(co[1]*co[2])) )
log.SE.LOD<-log10(SE.LOD)
# define boundaries
bottom<-x.LOD-(1.96*log.SE.LOD)
ceiling<-x.LOD+(1.96*log.SE.LOD)
### REPORT
# mean value
x.LOD
# bottom
bottom
# ceiling
ceiling
### plot LOD
lines(c(0.1, x.LOD), c(0.95, 0.95), lty=2)
lines(c(x.LOD, x.LOD), c(0.95, 0), lty=2)
text(x.LOD, 0.02, labels = round(x.LOD, digits = 2), pos=4, offset =
0.3, cex = 0.7)
text(1, 0.93, label="95%", pos=4, offset = 0.3, cex = 0.7)
lines(c(bottom, ceiling), c(0,0), lty=1, lend="butt", lwd=0.7)
points(x.LOD, 0, pch=21, cex=0.8, bg="white")
detach(my.data)
############################################################################
###############################################################
############################################################################
###############################################################
### SAMPLE SET No. 2
### define data
labs<-c("dilution", "total", "positive")
p.1<-c(10, 10, 4)
p.2<-c(100, 10, 10)
p.3<-c(1000, 10, 10)
p.4<-c(10000, 10, 10)
p.5<-c(100000, 10, 10)
p.6<-c(1000000, 10, 10)
### create data frame from matrix
my.mat<-matrix(c(p.1, p.2, p.3, p.4, p.5, p.6), nrow=6, byrow=TRUE)
dimnames(my.mat)<-list(c(1:6),labs)
my.data<-as.data.frame(my.mat)
attach(my.data)
my.data
### define frequency data
Y<-cbind(positive, total-positive)
### fit model LOGIT
model<-glm(Y ~ dilution, family=binomial(link=logit), data=my.data)
### plot data and regression line
plot(dilution, positive/total, pch=16, ylim=c(0,1), log = "x",
xaxt="n", xlim=c(1, 1000000),
main="Positive amplification by PCR", ylab="Fraction amplified",
xlab=expression(paste("Standard dilutions (c/",mu,"l)")))
### add x axis (William Dunlap, personal communication)
xlim <- par("usr")[1:2]
lab.x <- seq(ceiling(xlim[1]), floor(xlim[2]))
axis(side=1, at=10^lab.x, lab=as.expression(lapply(lab.x,
function(x)bquote(10^.(x)))))
lines(dilution, fitted(model), lty=1)
### set significance level and logit
# define significance
level<-0.95
# set logit value
logit.LOD<-log(level/(1-level))
### set function to determine LOD (Vito Ricci, personal communication)
LOD<-function(mod) as.vector((logit.LOD-coef(mod)[1])/coef(mod)[2])
### computation of the LOD value
x.LOD<-LOD(model)
### compute S.E.
# set variables
co<-model$coef
model.sum<-summary(model)
SE.co<-model.sum$coef[,2]
COV.co<-model.sum$cov.scaled[1,2]
# compute SE
SE.LOD<-abs(x.LOD) * sqrt( (SE.co[1]/co[1])^2 +
(SE.co[2]/co[2])^2 -
2*(COV.co/(co[1]*co[2])) )
log.SE.LOD<-log10(SE.LOD)
# define boundaries
bottom<-x.LOD-(1.96*log.SE.LOD)
ceiling<-x.LOD+(1.96*log.SE.LOD)
### REPORT
# mean value
x.LOD
# bottom
bottom
# ceiling
ceiling
### plot LOD
lines(c(0.1, x.LOD), c(0.95, 0.95), lty=2)
lines(c(x.LOD, x.LOD), c(0.95, 0), lty=2)
text(x.LOD, 0.02, labels = round(x.LOD, digits = 2), pos=4, offset =
0.3, cex = 0.7)
text(1, 0.93, label="95%", pos=4, offset = 0.3, cex = 0.7)
lines(c(bottom, ceiling), c(0,0), lty=1, lend="butt", lwd=0.7)
points(x.LOD, 0, pch=21, cex=0.8, bg="white")
detach(my.data)
############################################################################
###############################################################
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
