R/plot.LD50.R

In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions

#' plot.LD50 plot result of IC50() that best describe IC50
#' @title Plot results of LD50() that best describe LD50
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Nothing
#' @param x A result file generated by IC50()
#' @param ... Parameters for plot()
#' @param las.x las parameter for x axis
#' @param las.y las parameter for y axis
#' @param at Position of ticks in x-axis
#' @param lab.PT Label to describe pivotal dose
#' @param lab.TRD Label to describe transitional range of dose
#' @param col.TRD The color of TRD
#' @param col.TRD.CI The color of CI of TRD based on range.CI
#' @param col.PT.CI The color of CI of PT based on range.CI
#' @param show.CI Do the CI for the curve should be shown
#' @description Plot the estimates that best describe lethality of doses.
#' @family LD50 functions
#' @examples
#' \dontrun{
#' data <- data.frame(Doses=c(80, 120, 150, 150, 180, 200),
#' Alive=c(10, 12, 8, 6, 2, 1),
#' Dead=c(0, 1, 5, 6, 9, 15))
#' LD50_logistic <- LD50(data, equation="logistic")
#' predict(LD50_logistic, doses=c(140, 170))
#' plot(LD50_logistic, xlim=c(0, 300))
#' }
#' @method plot LD50
#' @export



plot.LD50 <- function(x, ...,  
	las.x=1, las.y=1, lab.PT="LD50", at=NULL, 
	lab.TRD=paste0("Transitional range of doses l=",l*100,"%"), 
	col.TRD="gray", col.TRD.CI=rgb(0.8, 0.8, 0.8, 0.5), 
  col.PT.CI=rgb(0.8, 0.8, 0.8, 0.5), show.CI=TRUE) {


  range.CI.qnorm <- qnorm(1-((1-x$range.CI)/2))
  l <- x$l
  alive <- x$alive
  dead <- x$dead
  N <- x$N
  doses <- x$doses
  equation <- x$equation  
  par <- x$par
  res <- x$SE
  parP <- x$LD50
  resP <- x$SE_LD50
  
  result <- x
 
 
   L <- list(...)

		L1 <- modifyList(list(x=doses, y=alive/N, bty="n", type="n", xlab="Doses", 
		                      ylab="Alive proportion"), L)
  L1 <- modifyList(L1, list(ylim=c(0,1), xaxt="n", las=las.y))
  
  if (is.null(L$xlim)) {
    L1 <- modifyList(L1, list(xlim=c(floor(min(doses)), floor(1+max(doses)))))
  }
  
  a <- do.call(plot, L1) 
  
  x2 <- (par("usr")[1]+par("usr")[2]*26)/27
  x1 <- x2*26-par("usr")[2]/0.04
  
  axis(1, at=at, las=las.x)
  
  # je trace la TRD centree sur P
  
  
  polygon(c(result$TRD_limits[1], result$TRD_limits[1], result$TRD_limits[2], result$TRD_limits[2]), c(0,1,1,0), border=NA, col=col.TRD)  
  # limites de la limite basse de la TRD
  polygon(c(result$TRD_limits[1]+range.CI.qnorm*result$SE_TRD_limits[1], result$TRD_limits[1]+range.CI.qnorm*result$SE_TRD_limits[1], result$TRD_limits[1]-range.CI.qnorm*result$SE_TRD_limits[1], result$TRD_limits[1]-range.CI.qnorm*result$SE_TRD_limits[1]), c(0,1,1,0), border=NA, col=col.TRD.CI)
  # limites de la limite haute de la TRD
  polygon(c(result$TRD_limits[2]+range.CI.qnorm*result$SE_TRD_limits[2], result$TRD_limits[2]+range.CI.qnorm*result$SE_TRD_limits[2], result$TRD_limits[2]-range.CI.qnorm*result$SE_TRD_limits[2], result$TRD_limits[2]-range.CI.qnorm*result$SE_TRD_limits[2]), c(0,1,1,0), border=NA, col=col.TRD.CI)
  # limites de la PT
  polygon(c(parP-range.CI.qnorm*resP, parP-range.CI.qnorm*resP, parP+range.CI.qnorm*resP, parP+range.CI.qnorm*resP), c(0,1,1,0), border=NA, col=col.PT.CI)
  par(xpd=TRUE)
  segments(parP, 0, parP, 1.05, lty=4)
  segments(result$TRD_limits[1], 0, result$TRD_limits[1], 1.15, lty=3)
  segments(result$TRD_limits[2], 0, result$TRD_limits[2], 1.15, lty=3)
  text(x=parP, y=1.1, lab.PT)
  text(x=parP, y=1.2, lab.TRD)
  
  	  b <- getFromNamespace(".BinomialConfidence", ns="HelpersMG")(alive,N)
  	  L1 <- modifyList(list(x=doses, y=alive/N, bty="n", type="p", ylim=c(0,1), y.plus = b[,3], y.minus = b[,2]), L)
  L1 <- modifyList(L1, list(ylim=c(0,1), xlab="", ylab="", main="", axes=FALSE, xlim=c(x1, x2)))
  
  par(xpd=FALSE)

  par(new=TRUE)
  
  a <- do.call(plot_errbar, L1) 
  
  if (!is.null(result$CI)) {

x <- result$CI["doses", ]
p <- result$CI["mean", ]
   
  par(new=TRUE)
  L1 <- modifyList(list(x=x, y=p, bty="n"), L)
  L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", xlim=c(x1, x2)))

  a <- do.call(plot, L1)

if (show.CI) {
  
  pm <- result$CI["2.5%", ]
  pp <- result$CI["97.5%", ]
  
	par(new=TRUE)
    L1 <- modifyList(list(x=x, y=pm, bty="n"), L)
  L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", lty=2, xlim=c(x1, x2)))
  a <- do.call(plot, L1) 
  
  par(new=TRUE)
    L1 <- modifyList(list(x=x, y=pp, bty="n"), L)
  L1 <- modifyList(L1, list(ylim=c(0,1), axes=FALSE, xlab="", ylab="", type="l", main="", lty=2, xlim=c(x1, x2)))
  a <- do.call(plot, L1) 
}
}
}