R/predict.LD50.R

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

#' predict.LD50 Estimate survival according to doses
#' @title Estimate survival according to doses
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return A data.frame with informations about survival
#' @param object A result file generated by LD50
#' @param doses A vector of temperatures
#' @param SE The standard error for doses, optional
#' @param range.CI The range of confidence interval for estimation, default=0.95
#' @param replicates Number of replicates to estimate CI
#' @param progressbar Logical. Does a progression bar must be shown
#' @param ... Not used
#' @description Estimate survival according to doses.\cr
#' The returned data.frame has the following components:\cr
#' doses, SE, survival, CI.minus.sexratio, CI.plus.sexratio, range.CI\cr
#' @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
#' }
#' @method predict LD50
#' @export


predict.LD50 <-
function(object, doses=NULL, SE=NULL, 
	range.CI=0.95, replicates=1000, progressbar=FALSE, ...) {
  
  # durations=NULL; SE=NULL; range.CI=0.95; replicates=1000; progressbar=FALSE
  

  if (is.null(doses)) doses <- object$doses
  
  if (!is.null(SE) & length(SE) != length(doses)) {
   stop("SE must be of the same length than doses")
  }
   
  x <- object
	equation <- x$equation
	par <- c(x$par, x$fixed.parameters)
	res <- x$SE
	

if (is.null(range.CI) | is.infinite(range.CI) | range.CI>1 | range.CI<0) {
  warning("The parameter range.CI is changed to 0.95")
  range.CI <- 0.95
}

    if (progressbar) pb<-txtProgressBar(min=1, max=length(doses), style=3)

    
    rep <- replicates-1
    df_par <- data.frame(P=unname(c(par["P"], rnorm(rep, par["P"], ifelse(is.na(res["P"]), 0, res["P"])))), 
                         S=unname(c(par["S"], rnorm(rep, par["S"], ifelse(is.na(res["S"]), 0, res["S"])))),
                         K=unname(ifelse(is.na(par["K"]), rep(NA, replicates), c(par["K"], rnorm(rep, par["K"], ifelse(is.na(res["K"]), 0, res["K"]))))),
                         K1=unname(ifelse(is.na(par["K1"]), rep(NA, replicates), c(par["K1"], rnorm(rep, par["K1"], ifelse(is.na(res["K1"]), 0, res["K1"]))))),
                         K2=unname(ifelse(is.na(par["K2"]), rep(NA, replicates), c(par["K2"], rnorm(rep, par["K2"], ifelse(is.na(res["K2"]), 0, res["K2"]))))))

    if (is.null(SE)) {
      SE <- rep(0, length(doses))
      replicates <- 1
    }
    
    df.out <- data.frame(doses=numeric(), 
                        SE=numeric(),
                        sexratio=numeric(), 
                        CI.minus.sexratio=numeric(), 
                        CI.plus.sexratio=numeric(), 
                        range.CI=numeric())
    
    for(i in 1:length(doses)) {
      if (progressbar) setTxtProgressBar(pb, i)
        doses_ec <- c(doses[i], rnorm(replicates-1, doses[i], SE[i]))
        out_tsd_plot <- apply(df_par, 1, function(par) {
        p <- .modelLD50(par, doses_ec, equation)
        return(list(sr=p))
      }
      )
      
      out_tsd3_plot <- unlist(out_tsd_plot)
      outquant <- quantile(out_tsd3_plot, probs = c((1-range.CI)/2, range.CI+((1-range.CI)/2)))
 #     outquant <- rbind(outquant, mean=out_tsd3_plot[1,], doses=doses_ec)
      df.out <- rbind(df.out, data.frame(doses=unname(doses[i]), 
                                SE=unname(SE[i]),
                                survival=unname(out_tsd3_plot[1]),
                                CI.minus.survival=unname(outquant[1]),
                                CI.plus.survival=unname(outquant[2]),
                                range.CI=range.CI))
      
    }

return(df.out)
}