R/ecotoxicology.R

In ecotoxicology: Methods for Ecotoxicology

#' @name Dunnett.t.Statistic
#' @title Critical Values of Dunnett's t Statistic
#' @description  Critical Values of Dunnett's t Statistic, Two-Tailed Comparisons
#' @docType data
#' @usage Dunnett.t.Statistic
#' @details Critical Values of Dunnett's t Statistic - data columns
#' \itemize{
#'   \item df. Degress of freedom.
#'   \item alpha. significance level.
#'   \item 2. k=2, Number of Treatment Means,  Including Control.
#'   \item 3. k=3, Number of Treatment Means,  Including Control.
#'   \item 4. k=4, Number of Treatment Means,  Including Control.
#'   \item 5. k=5, Number of Treatment Means,  Including Control.
#'   \item 6. k=6, Number of Treatment Means,  Including Control.
#'   \item 7. k=7, Number of Treatment Means,  Including Control.
#'   \item 8. k=8, Number of Treatment Means,  Including Control.
#'   \item 9. k=9, Number of Treatment Means,  Including Control.
#'   \item 10. k=10, Number of Treatment Means,  Including Control.
#' }
#' @references C. W. Dunnett, 1964. 
#' New tables for multiple comparisons with a control. 
#' Biometrics 20. 482–491.
#' @author Jose Gama
#' @keywords data
NULL

#' @name AphisRumicisDerrisMalaccensis
#' @title data on the toxicity to Aphis rumicis of an ether extract of Derris malaccensis
#' @description data on the toxicity to Aphis rumicis of an ether extract of Derris malaccensis
#' @docType data
#' @usage AphisRumicisDerrisMalaccensis
#' @details 
#' \itemize{
#'   \item concentration. concentration
#'   \item n. number of insects
#'   \item r. number of observed affected
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve. pp 238.
#' Cambridge University Press
#' 
#' Martin, J. T ., 1940
#' The problem of the evaluation of rotenone-containing plants. V. The relative 
#' toxicities of different species of derris. Ann. Appl. Biol. 27, 274-94.
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table1Finney1964
#' @title Transformation of Percentages to Probits, table I of Finney, 1964
#' @description Transformation of Percentages to Probits, table I of Finney, 1964
#' @docType data
#' @usage Table1Finney1964
#' @details Transformation of Percentages to Probits - data columns
#' \itemize{
#'   \item Percentage. Percentage.
#'   \item Col0.0. Column for 0.0
#'   \item Col0.1. Column for 0.1
#'   \item Col0.2. Column for 0.2
#'   \item Col0.3. Column for 0.3
#'   \item Col0.4. Column for 0.4
#'   \item Col0.5. Column for 0.5
#'   \item Col0.6. Column for 0.6
#'   \item Col0.7. Column for 0.7
#'   \item Col0.8. Column for 0.8
#'   \item Col0.9. Column for 0.9
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table2Finney1964
#' @title The Weighting Coefficient and Q/Z, table II of Finney, 1964
#' @description The Weighting Coefficient and Q/Z, table II of Finney, 1964
#' @docType data
#' @usage Table2Finney1964
#' @details The Weighting Coefficient and Q/Z - data columns
#' \itemize{
#'   \item Y. expected probit
#'   \item Q/Z. 
#'   \item C=0. 0% of natural mortality
#'   \item C=1. 1% of natural mortality
#' ...
#'   \item C=89. 89% of natural mortality
#'   \item C=90. 90% of natural mortality
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table3Finney1964
#' @title Maximum and Minimum working probits and Range, table III of Finney, 1964
#' @description Maximum and Minimum working probits and Range, table III of Finney, 1964
#' @docType data
#' @usage Table3Finney1964
#' @details Maximum and Minimum working probits and Range - data columns
#' \itemize{
#'   \item Ymin. Minimum working probit - expected
#'   \item Y0. Minimum working probit - Y0 = Y-P/Z
#'   \item Yrange. Range 1/Z
#'   \item Y100. Maximum working probit - Y100 = Y+Q/Z
#'   \item Ymax. Maximum working probit - expected
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table4Finney1964
#' @title Working probits, table IV of Finney, 1964
#' @description Working probits, table IV of Finney, 1964
#' @docType data
#' @usage Table4Finney1964
#' @details Working probits - data columns
#' \itemize{
#'   \item Kill%. % kill
#'   \item Col2 Expected probit 2.0
#'   \item Col2.1 Expected probit 2.1
#' ...
#'   \item Col7.8 Expected probit 7.8
#'   \item Col7.9 Expected probit 7.9
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table5Finney1964
#' @title The Probability, P, the Ordinate, Z, and Z^2, table V of Finney, 1964
#' @description Probability, P, the Ordinate, Z, and Z^2, table V of Finney, 1964
#' @docType data
#' @usage Table5Finney1964
#' @details The Probability, P, the Ordinate, Z, and Z^2 - data columns
#' \itemize{
#'   \item Y. Expected probit
#'   \item P. Probability P of expected probit
#'   \item Z. Ordinate to the normal distribution corresponding to the probability P
#'   \item Z^2. Z^2
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table8Finney1964
#' @title The Weighting Coefficient in Wadley's Problem, table VIII of Finney, 1964
#' @description The Weighting Coefficient in Wadley's Problem, table VIII of Finney, 1964
#' @docType data
#' @usage Table8Finney1964
#' @details The Weighting Coefficient in Wadley's Problem - data columns
#' \itemize{
#'   \item Y. Expected probit
#'   \item w. Weighting Coefficient
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name Table9Finney1964
#' @title Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling, table IX of Finney, 1964
#' @description Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling, table IX of Finney, 1964
#' @docType data
#' @usage Table9Finney1964
#' @details Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling - data columns
#' \itemize{
#'   \item Y. Expected probit
#'   \item MinWorkProbit. Minimum working probit
#'   \item Range. Range 1/Z
#'   \item WeightingCoef. Weighting Coefficient
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @keywords data
NULL

#' @name SheepsheadMinnow40SK
#' @title Mortality data from a fathead minnow larval survival and growth test (40 organisms per concentration)
#' @description Mortality data from a fathead minnow larval survival and growth test (40 organisms per concentration)
#' @docType data
#' @usage SheepsheadMinnow40SK
#' @details Mortality data from a fathead minnow larval survival and growth test - data columns
#' \itemize{
#'   \item Concentration. Concentration.
#'   \item Mortality. Mortality
#' }
#' @references USEPA, 2002
#' Short-term Methods for Estimating the Chronic Toxicity of Effluents and Receiving Waters to Freshwater Organisms.
#' 4th Edition,USEPA,Office of Water,October 2002,EPA 821-R-02-013
#' TABLE J1. pp 312
#' @author Jose Gama
#' @keywords data
NULL


#' Calculate working probit
#' @description Returns the working probit
#' @param Y numeric, expected probit
#' @param p numeric, kill percentage
#' @return the working probit
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' # Example from page 50 of Finney 1964:
#' # kill p = 72.3%, expected probit Y = 6.2
#' Y <- 6.2
#' p <- 72.3/100
#' # working probit = 5.366
#' ProbitWorkingP(Y,p)
ProbitWorkingP<-function(Y,p) {# calculate working probit
P <- pnorm(Y-5)
Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)
Y+(p-P)/Z
}

#' Calculate weighting coefficient from expected probit
#' @description Returns the weighting coefficient from expected probit
#' @param Y numeric, expected probit
#' @param C numeric, proportion of natural mortality
#' @return the weighting coefficient
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press. Formula 6.3.
#' @author Jose Gama
#' @examples
#' # Example from page 90 of Finney 1964:
#' # expected probit Y = 6.2, control mortality C = 59%
#' Y <- 6.2
#' C <- 0.59
#' # weighting coefficient = 0.141
#' Probitw(Y,C)
Probitw<-function(Y,C=0) {# weight w to be attached to the probit of P
# works with C (Z^2/(Q*(P+C/(1-C)))) and without it (default, w=Z^2/(PQ))
P <- pnorm(Y-5)
Q<-1-P
Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)
Z^2/(Q*(P+C/(1-C)))
}

#' Calculate the weighting coefficient
#' @description Returns the weighting coefficient
#' @param Z numeric, ordinate to the normal distribution corresponding to the probability P
#' @param Q numeric, 1-P
#' @param P numeric, Probability P of expected probit
#' @param C numeric, proportion of natural mortality
#' @return the weighting coefficient
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press. Formula 6.3.
#' @author Jose Gama
#' @examples
#' # Example from page 90 of Finney 1964:
#' # expected probit Y = 6.2, control mortality C = 59%
#' Y <- 6.2
#' C <- 0.59
#' P <- pnorm(Y-5)
#' Q <- 1-P
#' Z <- ProbitZ(Y)
#' # weighting coefficient = 0.141
#' ProbitWeightingCoef(Z,Q,P,C)
ProbitWeightingCoef<-function(Z,Q,P,C) Z^2/(Q*(P+C/(1-C))) # coefficient for the calculation of working probits

#' Calculate the ordinate to the normal distribution corresponding to the probability P
#' @description Returns the ordinate to the normal distribution corresponding to the probability P
#' @param Y numeric, expected probit
#' @return the ordinate to the normal distribution corresponding to the probability P
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press. Formula 3.5.
#' @author Jose Gama
#' @examples
#' # expected probit Y = 6.2
#' Y <- 6.2
#' ProbitZ(Y)
ProbitZ<-function(Y) 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)# calculate Z

#' Calculate the ordinate to the normal distribution corresponding to the probability P, exactly like Finney's
#' @description Returns the ordinate to the normal distribution corresponding to the probability P
#' with the exact same results as Finney's
#' @param Y numeric, expected probit
#' @return the ordinate to the normal distribution corresponding to the probability P
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press. Formula 3.5.
#' @author Jose Gama
#' @examples
#' # expected probit Y = 6.2
#' Y <- 6.2
#' ProbitZ4dec(Y)
ProbitZ4dec<-function(Y) round(1/round(sqrt(2*3.1416),4)*exp(-0.5*round((Y-5)^2,4)),4)# calculate Z exactly like Finney's

#' Estimate the column for Chi calculation
#' @description Estimates the column for Chi calculation
#' @param r numeric vector, number of observed affected
#' @param n numeric vector, number of insects
#' @param P numeric vector, Probability P of expected probit
#' @return numeric vector
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitChi<-function(r, n, P) (r-n*P)^2/(n*P*(1-P))

#' Variance of dosage
#' @description Variance of dosage
#' @param b numeric, rate of increase of probit value per unit increase in x 
#' @param m numeric, dosage
#' @param n numeric, number of insects
#' @param w numeric, weighting coefficients
#' @param x numeric, log concentration
#' @param xbar numeric, mean dosage
#' @return Variance of dosage
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitVarianceDosage<-function(b, m, n, w, x, xbar) 
{
Snw <- sum(n*w, na.rm =T)
(1/Snw+(m-xbar)^2/( sum(n*w*x^2, na.rm =T) - sum(n*w*x, na.rm =T)^2 /Snw  ))/b^2
}

#' Variance of rate of increase of probit value per unit increase in x 
#' @description Variance of rate of increase of probit value per unit increase in x 
#' @param n numeric, number of insects
#' @param w numeric, weighting coefficients
#' @param x numeric, log concentration
#' @param xbar numeric, mean dosage
#' @return Variance of rate of increase of probit value per unit increase in x 
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitVarianceRate<-function(n, w, x, xbar) 
{
Snw <- sum(n*w, na.rm =T)
1/( sum(n*w*x^2, na.rm =T) - sum(n*w*x, na.rm =T)^2 /Snw  )
}

#' Standard Error of rate of increase of probit value per unit increase in x 
#' @description Standard Error of rate of increase of probit value per unit increase in x 
#' @param n numeric, number of insects
#' @param w numeric, weighting coefficients
#' @param x numeric, log concentration
#' @param xbar numeric, mean dosage
#' @return Standard Error of rate of increase of probit value per unit increase in x 
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitStandardErrorRate<-function(n, w, x, xbar) sqrt(ProbitVarianceRate(n, w, x, xbar))

#' Standard Error of dosage
#' @description Standard Error of dosage
#' @param varianceDosage numeric, Variance of dosage
#' @return Standard Error of dosage
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitStandardErrorOfDosage<-function(varianceDosage) sqrt(varianceDosage)

#' Approximate Standard Error of dosage
#' @description Approximate Standard Error of dosage
#' @param b numeric, rate of increase of probit value per unit increase in x 
#' @param Snw numeric, sum of nw
#' @return Approximate Standard Error of dosage
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitApproxStandardErrorOfDosage<-function(b, Snw) sqrt(Snw)/b

#' Probit regression line
#' @description Probit regression line
#' @param x numeric, log concentration
#' @param n numeric, number of insects
#' @param r numeric, number of observed affected
#' @param adjAbbot logic, use Abbot adjustment
#' @param roundFinney logic, round as in Finney's book
#' @return Probit regression line a+bx
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitRegression<-function(x, n, r, adjAbbot=FALSE, roundFinney=FALSE)
{
p <- round(100*r/n)
if (adjAbbot) p <- round(AdjustAbbott(p,min(p),max(p)))
Y <- PercentageToProbit(p) # empirical probit
if (roundFinney) Y <- round(Y, 2)
posInf <- c(which(is.infinite(Y)), which(is.infinite(x)))
flagLast <- length(posInf)!=0 # keep track of -Inf at the end
if (flagLast) df.lm <- lm(Y[-posInf]~x[-posInf]) else df.lm <- lm(Y~x) # regression
list(a=as.numeric(coefficients(df.lm)[1]), b=as.numeric(coefficients(df.lm)[2]))
}

#' Probit value "g"
#' @description Probit value "g"
#' @param b numeric, rate of increase of probit value per unit increase in x
#' @param n numeric, number of insects
#' @param w numeric, weighting coefficients
#' @param x numeric, log concentration
#' @param xbar numeric, mean dosage
#' @param tPercent numeric, probability level
#' @return Probit value "g"
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitVALUEg<-function(b, n, w, x, xbar, tPercent) 
{
nw <- n*w
Vb <- 1/sum(nw*(x-xbar)^2, na.rm =T)
t5perct<-round(qt(1-tPercent/100/2, df=Inf, log.p=F),2)
t5perct^2*Vb/b^2
}

#' Probit Fiducial Limits
#' @description Probit Fiducial Limits
#' @param Vm numeric, variance of the logarithm
#' @param m numeric, logLD50
#' @param tPercent numeric, probability level
#' @param roundFinney logic, round as in Finney's book
#' @return Probit Fiducial Limits
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitFiducialLimits<-function(Vm, m, tPercent=5, roundFinney=FALSE)
{
SEm<-sqrt(Vm)
if (roundFinney) SEm<-round(SEm,3)
t5perct <- qt(1-tPercent/100/2, df=Inf, log.p=F)
t5perct*SEm*c(-1,1)+m
}

#' Probit estimation regression with Finney's method
#' @description Probit estimation regression with Finney's method
#' @param toxData numeric matrix, matrix with concentration, n ,r columns
#' @param tPercent numeric, probability level
#' @param showPlot logic, show regression line - plot
#' @param roundFinney logic, round as in Finney's book
#' @return Probit estimation regression
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitFinney<-function(toxData, tPercent=5, showPlot=FALSE, roundFinney=FALSE)
{
# 1 - sort descending by concentration

toxData <- toxData[order(toxData[,1]),]
concentration <- toxData[,1]

x<-log10(concentration)
n <- toxData[,2]
r <- toxData[,3]

# 5 - Calculate the percentage kill, p' = lOOr/n
pk <- round(100*r/n) # percentage kill
# 6 - Abbott adjustment
pAdj <- round(AdjustAbbott(pk,min(pk), max(pk)))
# 9 - empirical probit
empProbit <- PercentageToProbit(pAdj)
if (roundFinney) empProbit <- round(empProbit,2)
matR <- cbind(concentration,x,n,r,pk,pAdj,empProbit)

tmp <- ProbitRegression(x, n, r,TRUE,TRUE)
x0 <- tmp$a
m0 <- tmp$b

# 10 - Plot the empirical probits against x
if (showPlot)
{
plot(x, empProbit, pch=4)
abline(x0, m0) # regression line
}

# 11 - expected probits Y
Y <- x0+m0*x
if (roundFinney) Y <- round(Y,1)
matR<-cbind(matR, Y = Y)

# 11 - weighting coefficient for each Y
w <- Probitw(Y,0.02)
if (roundFinney) w <- round(w,5)
matR<-cbind(matR, w = w)
nw <- n*w
if (roundFinney) nw <- round(nw,1)
matR<-cbind(matR, nw=nw)

# 14 - working probit y
y <- ProbitWorkingP(Y,pAdj/100)
if (roundFinney) y <- round(y,2)
matR<-cbind(matR, y = y)
nwx <- nw*x
if (roundFinney) nwx <- round(nwx,3)
nwy <- nw*y
if (roundFinney) nwy <- round(nwy,3)
matR<-cbind(matR, nwx = nwx, nwy = round(nw*y,3))
matR[which(is.infinite(matR))]<-NA
matR
# 17, ..., 22 Snw, Snwxy
Snw<-sum(nw, na.rm =T)
Snwx<-sum(nwx, na.rm =T)
Snwy<-sum(nwy, na.rm =T)
xmean<-Snwx/Snw
ymean<-Snwy/Snw
Snwx2<-sum(nwx*x, na.rm =T)
Snwy2<-sum(nwy*y, na.rm =T)
Snwxy<-sum(nwx*y, na.rm =T)

S1<-Snwx^2/Snw
S2<-Snwx*Snwy/Snw
S3<-Snwy^2/Snw
Sxx<-Snwx2-S1
Sxy<-Snwxy-S2
Syy<-Snwy2-S3
b<-Sxy/Sxx
if (roundFinney) b<-round(b)

# 25 - equation for the probit regression line (expected numbers of insects killed at each concentration)
Px0<-b*-xmean+ymean
if (roundFinney) Px0<-round(Px0,2)
regPstr <- paste(Px0,'+',b,'*x',sep='')
#regP <- function(x) Px0+b*x

# 26 - chi^2
chi2 <- Syy-Sxy^2/Sxx

# 28 - standard error of b
sEb<-1/sqrt(Sxx)
cat('b=',b,'+/-',sEb,'\n')

# 30 - Find the logLD50
m <- (5-Px0)/b
if (roundFinney) m <- round(m, 3)
LD50 <- round(10^m)

# 31 - standard error of m
varianceDOSAGEm <- ProbitVarianceDosage(b, m, n, w, x, xmean) 
cat('Standard error of m = +/-',sqrt(varianceDOSAGEm),'\n')

# 32 - Calculate g for 5% probability level
ProbitVALUEg(b, n, w, x, xmean, tPercent)

# 33 - fiducial limits to m
fiduc <- ProbitFiducialLimits(varianceDOSAGEm, m, tPercent)

# 35 - Conclusions from the analysis
cat('maximum likelihood estimate of',LD50,'\n') 
cat('fiducial probability of',100-tPercent,'%\n') 
cat('The true value of the media.n lethal dose may be expected to lie between',round(10^fiduc),'\n') 

}


#' Spearman Karber estimation
#' @description Spearman Karber estimation
#' @param toxData numeric matrix, matrix with concentration, n ,r columns
#' @param N numeric, number of organisms
#' @param retData logic, return the results in a list
#' @param showOutput logic, show results in the console
#' @param showPlot logic, show regression line - plot
#' @return Spearman Karber estimation
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
SpearmanKarber<-function(toxData, N, retData=FALSE, showOutput=TRUE, showPlot=TRUE)
{
toxData <- as.matrix(toxData)
n <- dim(toxData)[1]
if (!IsMonotonicallyIncreasing(toxData[,2]/N)) 
{
toxData <- cbind(toxData,smoothed=MakeMonotonicallyIncreasing(cbind(rep(N,n),toxData[,2])))
AdjustMortality <- AdjustAbbott(toxData[,4])
} else AdjustMortality <- AdjustAbbott(toxData[,3])

toxData <- cbind(toxData,smoothedAdjusted=AdjustMortality)

# log10 LC50
m <- sum(diff(AdjustMortality)*(log10(toxData[-n,1])+log10(toxData[-1,1]))/2, na.rm =T)

# estimated variance of m
Vm <- sum(diff(AdjustMortality[c(-1,-n)]) * (1-diff(AdjustMortality[c(-1,-n)])) * ( log10(toxData[c(-(1:3)),1]) - log10(toxData[c(-1, -n,-(n-1)),1]) )^2 /(4*N-4))

# 95% confidence interval for m
confInt4m <- m+c(-2, 2)*sqrt(Vm)

# estimated LC50
LC50 <- 10^m
# estimated 95% confidence interval for the estimated LC50
confIntLC50 <- 10^confInt4m

if (showPlot)
{
plot(toxData[,c(1,3)], type='p',pch=4,ylim=c(0,1),ylab='Mortality Proportion',xlab='Effluent Concentration %')
par(new=TRUE)
plot(toxData[,c(1,4)], type='l',ylim=c(0,1),xlab='',ylab='')
par(new=TRUE)
plot(toxData[,c(1,5)], type='l',ylim=c(0,1),xlab='',ylab='', lty=2)
legend("topleft", c('Observed Mortality Proportion','Smoothed Mortality Proportion','Smoothed Adjusted Mortality Proportion'), lty = c(0, 1, 2), pch=c(4,NA,NA))
}

if (showOutput)
{
print(toxData)
cat('log10 LC50 =',m,'\n')
cat('estimated variance of m =',Vm,'\n')
cat('95% confidence interval for m =',paste(confInt4m,sep=','),'\n')
cat('estimated LC50 =',LC50,'\n')
cat('estimated 95% confidence interval for the estimated LC50 =',paste(confIntLC50,sep=','),'\n')
}
if (retData)
{
list(log10LC50 =m, varianceOfm=Vm, confidenceInterval95m =confInt4m, LC50=LC50, confidenceInterval95LC50=confIntLC50)
}
}

#' Probit estimation similar to the EPA's Ecological Exposure Research Division (EERD) tool
#' @description Probit estimation similar to the EPA's Ecological Exposure Research Division (EERD) tool
#' @param toxData numeric matrix, matrix with concentration, n ,r columns
#' @param retData logic, return the results in a list
#' @param showOutput logic, show results in the console
#' @return Probit estimation regression
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
ProbitEPA<-function(toxData, retData=FALSE, showOutput=TRUE)
{
toxData <- toxData[order(toxData[,1]),]
concentration <- toxData[,1]

x<-log10(concentration)
# 		x <- AdjustAbbott(toxData[,3]/toxData[,2])
n <- toxData[,2]
r <- toxData[,3]
if (concentration[1]==0) NCTRL<-n[1] else NCTRL<-0
pobs<-r/n
MAXITER<-55
TOLERANCE<-1.0e-9
MINVALUE<-1e-9
MINP<-10.0*MINVALUE
INFINITESLOPE<-200
N<-length(n)
a<-5
b<-0
x50<-0
oldx50<-0
numiteration<-0
infSlope<-FALSE
bigC <-0
if (NCTRL != 0) bigC<-rep(pobs[1],N)
repeat {
oldx50<-x50
yp<-a+b*x
ypm5<-yp-5.0
z<-dnorm(ypm5)
py<-pnorm(ypm5)
xp <- (1-py)/z
if (any(abs(pobs[-which(is.na(py))]-py[-which(is.na(py))])<=MINP)) {w=0;ywp<-yp} else {w<-z^2/(py*(1-py));ywp<-yp+(pobs-py)/z}

y=ywp
Snw <-sum(n*w, na.rm =T)
Snwx<-sum(n*w*x, na.rm =T)
Snwy<-sum(n*w*y, na.rm =T)

Snwyy<-sum(n*w*y*y, na.rm =T)
Snwxy<-sum(n*w*x*y, na.rm =T)
Snwxx<-sum(n*w*x*x, na.rm =T)
Snwxp<-sum(n*w*xp, na.rm =T)
Snwxpp<-sum(n*w*xp^2, na.rm =T)
Snwxpy<-sum(n*w*xp*y, na.rm =T)
Snwxxp<-sum(n*w*xp*x, na.rm =T)

xbar<-Snwx/Snw
ybar<-Snwy/Snw
xpbar<-Snwxp/Snw
Syy=Snwyy-Snwy^2/Snw

Sxx=Snwxx-Snwx^2/Snw
Sxy=Snwxy-Snwx*Snwy/Snw
Sxpxp=Snwxpp-Snwxp*Snwxp/Snw
Sxxp=Snwxxp-Snwx*Snwxp/Snw
Sxpy=Snwxpy-Snwxp*Snwy/Snw
b=Sxy/Sxx
a=ybar-b*xbar

delC<-(1-bigC)*((Sxpy-Sxy/Sxx*Sxxp)/(Sxpxp - Sxxp^2/Sxx))
bigC<-bigC+delC

a<-ybar-b*xbar

varB<-1.0/Sxx
varYbar<-1.0/Snw
varA<-1.0/Snw+xbar^2/Sxx
chi2<-Syy-Sxy^2/Sxx
chi2<-(Syy - b * Sxy)
HET<-chi2/(N-2)
SEI = (1 / Snw) + (xbar * xbar * 1/(Sxx))
SEa<- sqrt(SEI*HET) # sqrt(1.0/Snw+xbar^2/Sxx)
SEb<-sqrt(1/(Sxx)*HET)

iSigma<-1.0/b
iMu<-(5-a)/b

x50<-(5.0-a)/b
if(is.infinite(x50)) x50=infSlope+1
if (a>INFINITESLOPE) infSlope<-TRUE
numiteration<-numiteration+1
if (numiteration>MAXITER)
{
numiteration<- -1
break
}
if (!((abs(oldx50-x50)>TOLERANCE) & !infSlope)) break
}
probLevel<-0.05
# t4P<-qt(1-probLevel/2, df=Inf, log.p=F) # distribution of t 5%=1.96
t4P<-c(12.706,4.303,3.182,2.776,2.571,2.447,2.365,2.306, 2.262,2.228,2.201,2.179,2.16,2.145,2.131,2.12,2.11,2.101)[N-2]

fLb<-b + c(-1,1)*sqrt(1/(Sxx) * HET)*t4P
fLy<-a + c(-1,1)*sqrt(SEI * HET)*t4P # Fiducial limits to Y
ECpoint<-c(1,5,10,15,50,85,90,95,99) # point
m<- (PercentageToProbit(ECpoint)-a)/b
EC<-10^m # concentration

g <- t4P^2/b^2*varB
fB<-(1-g)/Snw+(m-xbar)^2/varB
facSqrt<- sqrt(HET)*t4P/b/(1-g)
fB<-t4P/b/(1-g)*sqrt((1-g)/Snw+(m-xbar)^2/Sxx)
fA<- g * (m-xbar)/(1-g)

if(showOutput)
{
cat('						PROBIT ANALYSIS PROGRAM
					  USED FOR CALCULATING LC/EC VALUES')
cat('		Number	Number	Number
Number	Conc.	Resp.	Exposed')
print(cbind(1:N,concentration,n,r,pobs,pobs,py))
cat('Chi - Square for Heterogeneity = ',chi2,'\nMu = ',iMu,'\nSigma = ',iSigma)
cat('\n\nParameter       Estimate    Std. Err.         95% Confidence Limits
')
cat('Intercept',a,SEa,fLy,'\nSlope',b,SEb, fLb,'\n\n')
cat('					Estimated LC/EC Values and Confidence Limits\n')
cat('Exposure	95% Confidence Limits\n')
cat('Point	Conc.	Lower	Upper\n')
print(cbind('LC/EC\t',ECpoint, EC,' (',10^(m+fA-fB),', ', 10^(m+fA+fB),')'))
}

if(retData)
{
list(
tableData=cbind(concentration,n,r,pobs,pobs,py), chi2, mu=iMu, sigma=iSigma, intercept=c(a,SEa,fLy), slope=c(b,SEb, fLb), 
tableResult=cbind(point=ECpoint, concentration=EC,lower=10^(m+fA-fB), upper=10^(m+fA+fB))
)
}

}

#' Calculate LC for N between 0 (LC0) and 100 (LC100)
#' @description Returns the LC for n between 0 and 100
#' @param x numeric, log concentration
#' @param n numeric, number of insects
#' @param r numeric, number of observed affected
#' @param N numeric, Lethal Concentration "N"
#' @return the LC for n between 0 and 100
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
CalculateLCn<-function(x, n, r, N=50)
{
if (N<0) N <- 0
if (N>100) N <- 100
tmp <- ProbitRegression(x, n, r)
a <- tmp$a
b <- tmp$b
Y <- PercentageToProbit(N)
logLD <- (Y-a)/b
10^logLD
}







#' Calculate LC50 from a matrix with 3 columns: concentration, number of exposed subjects and number of deaths
#' @description Returns the LC50 from a matrix with 3 columns: concentration, number of exposed subjects and number of deaths
#' @param matrixConcExpoResp numeric vector
#' @return the LC50
#' @references 
#' Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' Trimmed spearman-karber method for estimating median
#' Lethal concentrations in toxicity bioassays.
#' Environ. Sci. Technol. 11(7): 714-719;
#' Correction 12(4):417 (1978).
#' @author Jose Gama
#' @examples
#' #Data from the example on page 5:
#' #Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' #Trimmed spearman-karber method for estimating median
#' #Lethal concentrations in toxicity bioassays.
#' #Environ. Sci. Technol. 11(7): 714-719;
#' #Correction 12(4):417 (1978).
#' concentration<-c(.5,1,2,4,8)
#' exposed<-c(10,10,10,10,10)
#' mortality<-c(0,2,4,9,10)
#' CalculateLC50(cbind(concentration, exposed, mortality))
CalculateLC50<-function(matrixConcExpoResp)
{# calculate LC50 from a matrix with 3 columns: concentration, number of exposed subjects and number of deaths
concentration<-matrixConcExpoResp[,1]
exposed<-matrixConcExpoResp[,2]
mortality<-matrixConcExpoResp[,3]
p<-mortality/exposed
x<-log(concentration)
if (!IsMonotonicallyIncreasing(p)) p<-MakeMonotonicallyIncreasing(matrixConcExpoResp)
m<-0
for (i in 1:(length(concentration)-1)) m<-m+(p[i + 1]-p[i])*(x[i] + x[i + 1])/2
exp(m)# LC50
}

#' Calculate corrected efficacy \% with Abbott's formula
#' @description Returns the corrected efficacy \% with Abbott's formula
#' @param smoothedObservedProportion numeric vector, treated population
#' @param ps0 numeric vector, control
#' @param p1 numeric vector, percentage 0 to1 or 0 to 100 (p1=1 or P1=100)
#' @return the corrected efficacy \%
#' @source
#' ehabsoft, last accessed 2015
#' \url{http://www.ehabsoft.com/ldpline/onlinecontrol.htm}
#' @references Puntener W., 1981 
#' Manual for field trials in plant protection second edition. Agricultural Division, Ciba-Geigy Limited.
#' @author Jose Gama
#' @examples
#' #same result as example on Short-term Methods for Estimating the Chronic Toxicity of 
#' #Effluents and Receiving Waters to Freshwater Organisms.TABLE J1. page 312
#' data(SheepsheadMinnow40SK)
#' IsMonotonicallyIncreasing(SheepsheadMinnow40SK[,2]/40)
#' mydata <- cbind(SheepsheadMinnow40SK,
#'   MakeMonotonicallyIncreasing(cbind(rep(40,6),SheepsheadMinnow40SK[,2])))
#' AdjustAbbott(mydata[,3])
AdjustAbbott<-function(smoothedObservedProportion, ps0=smoothedObservedProportion[1], p1=1)
{# ps0 is the control, p1 is the percentage 0 to1 or 0 to 100 (p1=1 or P1=100)
smoothedObservedProportion<-(smoothedObservedProportion-ps0)/(p1-ps0)
smoothedObservedProportion*p1
}

#' Calculate corrected efficacy \% with Henderson-Tilton's formula
#' @description Returns the corrected efficacy \% with Henderson-Tilton's formula
#' @param smoothedObservedProportion numeric vector, treated population
#' @param ps0 numeric vector, control
#' @param p1 numeric vector, percentage 0 to1 or 0 to 100 (p1=1 or P1=100)
#' @return the corrected efficacy \%
#' @source ehabsoft, last accessed 2015
#' \url{http://www.ehabsoft.com/ldpline/onlinecontrol.htm}
#' @references Puntener W., 1981 
#' Manual for field trials in plant protection second edition. Agricultural Division, Ciba-Geigy Limited.
#' @author Jose Gama
AdjustHendersonTilton<-function(smoothedObservedProportion, ps0=smoothedObservedProportion[1], p1=1)
{
ps0<-smoothedObservedProportion[1]
smoothedObservedProportion<-(smoothedObservedProportion-ps0)/(1-ps0)*100
smoothedObservedProportion
}

#' Calculate corrected efficacy \% with Sun-Shepard's formula
#' @description Returns the corrected efficacy \% with Sun-Shepard's formula
#' @param smoothedObservedProportion numeric vector, treated population
#' @param ps0 numeric vector, control
#' @param p1 numeric vector, percentage 0 to1 or 0 to 100 (p1=1 or P1=100)
#' @return the corrected efficacy \%
#' @source ehabsoft, last accessed 2015
#' \url{http://www.ehabsoft.com/ldpline/onlinecontrol.htm}
#' @references Puntener W., 1981 
#' Manual for field trials in plant protection second edition. Agricultural Division, Ciba-Geigy Limited.
#' @author Jose Gama
AdjustSunShepard<-function(smoothedObservedProportion, ps0=smoothedObservedProportion[1], p1=1)
{
ps0<-smoothedObservedProportion[1]
smoothedObservedProportion<-(smoothedObservedProportion-ps0)/(1-ps0)*100
smoothedObservedProportion
}

#' Calculate corrected efficacy \% with Schneider-Orelli's formula
#' @description Returns the corrected efficacy \% with Schneider-Orelli's formula
#' @param smoothedObservedProportion numeric vector, treated population
#' @param ps0 numeric vector, control
#' @param p1 numeric vector, percentage 0 to1 or 0 to 100 (p1=1 or P1=100)
#' @return the corrected efficacy \%
#' @source ehabsoft, last accessed 2015
#' \url{http://www.ehabsoft.com/ldpline/onlinecontrol.htm}
#' @references Puntener W., 1981 
#' Manual for field trials in plant protection second edition. Agricultural Division, Ciba-Geigy Limited.
#' @author Jose Gama
AdjustSchneiderOrelli<-function(smoothedObservedProportion, ps0=smoothedObservedProportion[1], p1=1)
{
ps0<-smoothedObservedProportion[1]
smoothedObservedProportion<-(smoothedObservedProportion-ps0)/(1-ps0)*100
smoothedObservedProportion
}

#' Determine if a series is monotonically decreasing
#' @description Returns TRUE if all proportions are in a monotonically decreasing sequence
#' @param p numeric vector
#' @return True is the series is monotonically decreasing
#' @references
#' Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' Trimmed spearman-karber method for estimating median
#' Lethal concentrations in toxicity bioassays.
#' Environ. Sci. Technol. 11(7): 714-719;
#' Correction 12(4):417 (1978).
#' @author Jose Gama
#' @examples
#' IsMonotonicallyDecreasing(1:10)
#' IsMonotonicallyDecreasing(6:2)
#' IsMonotonicallyDecreasing(c(1,3,2))
IsMonotonicallyDecreasing<-function(p)
{ # returns TRUE if all proportions are in a monotonically decreasing sequence
l1<-length(p)-1
l2<-l1+1
all(p[1:l1]>=p[2:l2])
}

#' Determine if a series is monotonically increasing
#' @description Returns TRUE if all proportions are in a monotonically increasing sequence
#' @param p numeric vector
#' @return True is the series is monotonically increasing
#' @references 
#' Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' Trimmed spearman-karber method for estimating median
#' Lethal concentrations in toxicity bioassays.
#' Environ. Sci. Technol. 11(7): 714-719;
#' Correction 12(4):417 (1978).
#' @author Jose Gama
#' @examples
#' #Data from the example on page 8:
#' #Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' #Trimmed spearman-karber method for estimating median
#' #Lethal concentrations in toxicity bioassays.
#' #Environ. Sci. Technol. 11(7): 714-719;
#' #Correction 12(4):417 (1978).
#' concentration<-c(1.1,2.3,4.5,8.8,17.1)
#' exposed<-c(10,10,9,10,10)
#' mortality<-c(1,5,4,2,7)
#' p<-mortality/exposed
#' x<-log(concentration)
#' IsMonotonicallyIncreasing(p)
IsMonotonicallyIncreasing<-function(p)
{
l1<-length(p)-1
l2<-l1+1
all(p[1:l1]<=p[2:l2])
}

#' Smoothed Mortality Proportion (monotonically increasing sequence)
#' @description Returns the Smoothed Mortality Proportion (monotonically increasing sequence)
#' @param matrixExpoResp numeric vector or matrix
#' @return The Smoothed Mortality Proportion (monotonically increasing sequence)
#' @references 
#' Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' Trimmed spearman-karber method for estimating median
#' Lethal concentrations in toxicity bioassays.
#' Environ. Sci. Technol. 11(7): 714-719;
#' Correction 12(4):417 (1978).
#' @author Jose Gama
MakeMonotonicallyIncreasing<-function(matrixExpoResp)
{# returns the Smoothed Mortality Proportion (monotonically increasing sequence)
ok<-FALSE
if(is.null(dim(matrixExpoResp))) if(is.numeric(matrixExpoResp)) { p<- matrixExpoResp;ok<-TRUE }
if(!is.null(dim(matrixExpoResp))) if(dim(matrixExpoResp)[2]==1) if(is.numeric(matrixExpoResp)) p<-c(matrixExpoResp)
if(!is.null(dim(matrixExpoResp))) if(is.numeric(matrixExpoResp)) {
ok<-TRUE
exposed<-matrixExpoResp[,1]
mortality<-matrixExpoResp[,2]
p<-mortality/exposed
}
if(!ok) stop('Input must be numeric, either vector or matrix')
while(TRUE){
if (IsMonotonicallyIncreasing(p)) break
l1<-length(p)-1
l2<-l1+1
w1<-which(p[1:l1]>p[2:l2])
p[c(w1,w1+1)]<-mean(p[c(w1,w1+1)])
}
p
}

#' Make monotonically decreasing sequence
#' @description Returns a monotonically decreasing sequence
#' @param matrixExpoResp numeric vector or matrix
#' @return monotonically decreasing sequence
#' @references 
#' Hamilton, m.a., R.c. Russo, and r.v. Thurston, 1977.
#' Trimmed spearman-karber method for estimating median
#' Lethal concentrations in toxicity bioassays.
#' Environ. Sci. Technol. 11(7): 714-719;
#' Correction 12(4):417 (1978).
#' @author Jose Gama
MakeMonotonicallyDecreasing<-function(matrixExpoResp)
{# returns a monotonically decreasing sequence
ok<-FALSE
if(is.null(dim(matrixExpoResp))) if(is.numeric(matrixExpoResp)) { p<- matrixExpoResp;ok<-TRUE }
if(!is.null(dim(matrixExpoResp))) if(dim(matrixExpoResp)[2]==1) if(is.numeric(matrixExpoResp)) p<-c(matrixExpoResp)
if(!is.null(dim(matrixExpoResp))) if(is.numeric(matrixExpoResp)) {
ok<-TRUE
exposed<-matrixExpoResp[,1]
mortality<-matrixExpoResp[,2]
p<-mortality/exposed
}
if(!ok) stop('Input must be numeric, either vector or matrix')
while(TRUE){
if (IsMonotonicallyDecreasing(p)) break
l1<-length(p)-1
l2<-l1+1
w1<-which(p[1:l1]<p[2:l2])
p[c(w1,w1+1)]<-mean(p[c(w1,w1+1)])
}
p
}

#' Convert percentages to Probit values
#' @description Converts percentages to Probit values
#' @param mypercentage numeric vector
#' @return Probit values
#' @references Statistical tests for significance, accessed October 2015
#' \url{http://archive.bio.ed.ac.uk/jdeacon/statistics/tress4.html}
#' @author Jose Gama
#' @examples
#' a<-c(.1,.5,1:10,50,96,97,98,99.5,99.99,99.999,99.9999)
#' b<-PercentageToProbit(a)
PercentageToProbit<-function(mypercentage) qnorm(mypercentage/100,5,1) # converts percentages to Probit values

#' Convert Probit values to percentages
#' @description Converts Probit values to percentages
#' @param myprobit numeric vector
#' @return percentages
#' @references Statistical tests for significance, accessed October 2015
#' \url{http://archive.bio.ed.ac.uk/jdeacon/statistics/tress4.html}
#' @author Jose Gama
#' @examples
#' a<-c(.1,.5,1:10,50,96,97,98,99.5,99.99,99.999,99.9999)
#' b<-PercentageToProbit(a)
#' d<-ProbitToPercentage(b)
ProbitToPercentage<-function(myprobit) pnorm(myprobit,5,1)*100 # converts Probit values to percentages

#' Convert percentages to Arcsin values
#' @description Converts percentages to Arcsin values
#' @param mypercentage numeric vector
#' @return Arcsin values
#' @references Statistical tests for significance, accessed October 2015
#' \url{http://archive.bio.ed.ac.uk/jdeacon/statistics/tress4.html}
#' @author Jose Gama
#' @examples
#' a<-c(.1,.5,1:10,50,96,97,98,99.5,99.99,99.999,99.9999)
#' b<-PercentageToProbit(a)
#' d<-ProbitToPercentage(b)
#' e<-PercentageToArcsin(d)
PercentageToArcsin<-function(mypercentage) asin(sqrt(mypercentage/100))*180/pi # converts percentages to Arcsin values

#' Convert Arcsin values to percentages
#' @description Converts Arcsin values to percentages
#' @param myarcsin numeric vector
#' @return percentages
#' @references Statistical tests for significance, accessed October 2015
#' \url{http://archive.bio.ed.ac.uk/jdeacon/statistics/tress4.html}
#' @author Jose Gama
#' @examples
#' a<-c(.1,.5,1:10,50,96,97,98,99.5,99.99,99.999,99.9999)
#' b<-PercentageToProbit(a)
#' d<-ProbitToPercentage(b)
#' e<-PercentageToArcsin(d)
#' f<-ArcsinToPercentage(e)
ArcsinToPercentage<-function(myarcsin) sin(myarcsin/180*pi)^2*100 # converts Arcsin values to percentages

#' Generate table I from Finney1964 "Transformation of percentages to probits"
#' @description Generates table I from Finney1964 "Transformation of percentages to probits"
#' @return table I from Finney1964 "Transformation of percentages to probits"
#' \itemize{
#'   \item Percentage. Percentage.
#'   \item Col0.0. Column for 0.0
#'   \item Col0.1. Column for 0.1
#'   \item Col0.2. Column for 0.2
#'   \item Col0.3. Column for 0.3
#'   \item Col0.4. Column for 0.4
#'   \item Col0.5. Column for 0.5
#'   \item Col0.6. Column for 0.6
#'   \item Col0.7. Column for 0.7
#'   \item Col0.8. Column for 0.8
#'   \item Col0.9. Column for 0.9
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableIFinney1964()
GenTableIFinney1964<-function(){
t1v <- c(seq(0,98,0.1)[-981],seq(98,100,0.01)[-201])
t1 <- matrix(t1v,98+20,10,byrow=T)
t2 <- t1[,1]
t1 <- round(PercentageToProbit(t1),4)
t1 <- cbind(t2, t1)
t1[1,2] <- NA
colnames(t1) <- c('Percentage','Col0.0','Col0.1','Col0.2','Col0.3','Col0.4','Col0.5','Col0.6','Col0.7','Col0.8','Col0.9')
t1
}

#' Generate table II from Finney1964 "The weighting coefficient and Q/Z"
#' @description Generates table II from Finney1964 "The weighting coefficient and Q/Z"
#' @return table II from Finney1964 "The weighting coefficient and Q/Z"
#' \itemize{
#'   \item Y. expected probit
#'   \item Q/Z. 
#'   \item C=0. 0% of natural mortality
#'   \item C=1. 1% of natural mortality
#' ...
#'   \item C=89. 89% of natural mortality
#'   \item C=90. 90% of natural mortality
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableIIFinney1964()
GenTableIIFinney1964<-function(){
Y <- seq(1.1,9.0,0.1);P <- pnorm(Y-5);Q<-1-P;Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2);n<-length(Y)
t1 <- matrix(NA,n,93)
t1[,1] <- Y
t1[,2] <- Q/Z
cnames <-c()
for (C in 0:90) { t1[,3+C] <- round(ProbitWeightingCoef(Z,Q,P,rep(C,n)/100),5); cnames <-c(cnames,paste('C=',C,sep='')) }
t1[which(t1==0)] <- NA
colnames(t1) <- c('Y','Q/Z',cnames)
t1[45:n,2]<-round(t1[45:n,2],4)
t1[16:44]<-round(t1[16:44,2],3)
t1[11:15,2]<-round(t1[11:15,2],2)
t1[6:10,2]<-round(t1[6:10,2],1)
t1[1:5,2]<-round(t1[1:5,2])
t1
}

#' Generate table III from Finney1964 "Maximum and minimum working probits and range"
#' @description Generates table III from Finney1964 "Maximum and minimum working probits and range"
#' @return table III from Finney1964 "Maximum and minimum working probits and range"
#' \itemize{
#'   \item Ymin. Minimum working probit - expected
#'   \item Y0. Minimum working probit - Y0 = Y-P/Z
#'   \item Yrange. Range 1/Z
#'   \item Y100. Maximum working probit - Y100 = Y+Q/Z
#'   \item Ymax. Maximum working probit - expected
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableIIIFinney1964()
GenTableIIIFinney1964<-function(){
Y <- seq(1.1,6.5,0.1);P <- pnorm(Y-5);Q<-1-P;Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2);
Y2 <- 10-Y;P2 <- pnorm(Y2-5);Q2<-1-P2;Z2 <- 1/sqrt(2*pi)*exp(-0.5*(Y2-5)^2);
y0 <- Y-P/Z
y100 <- Y2+Q2/Z2
#range <- as.numeric(formatC(1/Z, width = 5))
range <- 1/Z
t1<-cbind(Ymin=Y, y0=round(y0,4), range=range, y100=round(y100,4),Ymax=(Y2))
t1[1:5,'range']<-round(t1[1:5,'range'])
t1[6:10,'range']<-round(t1[6:10,'range'],1)
t1[11:15,'range']<-round(t1[11:15,'range'],2)
t1[16:19,'range']<-round(t1[16:19,'range'],3)
t1[20:55,'range']<-round(t1[20:55,'range'],4)
t1
}

#' Generate table IV from Finney1964 "Working probits"
#' @description Generates table IV from Finney1964 "Working probits"
#' @return table IV from Finney1964 "Working probits"
#' \itemize{
#'   \item Kill%. % kill
#'   \item Col2 Expected probit 2.0
#'   \item Col2.1 Expected probit 2.1
#' ...
#'   \item Col7.8 Expected probit 7.8
#'   \item Col7.9 Expected probit 7.9
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableIVFinney1964()
GenTableIVFinney1964<-function(){
Y <- seq(2.0,7.9,0.1);P <- pnorm(Y-5);Q<-1-P;Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2);n<-length(Y)
t1 <- matrix(0,101,(n+1))
cnames <- paste('Col',Y,sep='')
for (p in 1:101) t1[p,2:(n+1)]<- Y+((p-1)/100-P)/Z
t1[which(t1<0)]<-NA
t1[which(t1>=10)]<-NA
t1[,1]<-0:100#/100
t1[,2:61]<-round(t1[,2:61],3)
colnames(t1) <- c('Kill%', cnames)
t1
}

#' Generate table V from Finney1964 "The Probability, P, the Ordinate, Z, and Z^2"
#' @description Generates table V from Finney1964 "The Probability, P, the Ordinate, Z, and Z^2"
#' @return table V from Finney1964 "The Probability, P, the Ordinate, Z, and Z^2"
#' \itemize{
#'   \item Y. Expected probit
#'   \item P. Probability P of expected probit
#'   \item Z. Ordinate to the normal distribution corresponding to the probability P
#'   \item Z^2. Z^2
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableVFinney1964()
GenTableVFinney1964<-function(){
Y <- 50:90/10;P <- pnorm(Y-5);Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)
t1 <- cbind(Y=Y,P=P,Z=Z,"Z^2"=round(Z^2,5))
t1[38:41,2] <- round(t1[38:41,2],9)
t1[32:38,2] <- round(t1[32:38,2],8)
t1[25:31,2] <- round(t1[25:31,2],7)
t1[13:24,2] <- round(t1[13:24,2],6)
t1[1:13,2] <- round(t1[1:13,2],5)
t1[36:41,3] <- round(t1[36:41,3],8)
t1[29:35,3] <- round(t1[29:35,3],7)
t1[18:28,3] <- round(t1[18:28,3],6)
t1[1:17,3] <- round(t1[1:17,3],5)
t1
}

#' Generate table VI from Finney1964 "Distribution of chi^2"
#' @description Generates table VI from Finney1964 "Distribution of chi^2"
#' @return table VI from Finney1964 "Distribution of chi^2"
#' \itemize{
#'   \item Deg.freedom. Degrees of freedom
#'   \item 0.9. Probability 0.9
#'   \item 0.7. Probability 0.7
#'   \item 0.5. Probability 0.5
#'   \item 0.3. Probability 0.3
#'   \item 0.1. Probability 0.1
#'   \item 0.05. Probability 0.05
#'   \item 0.02. Probability 0.02
#'   \item 0.01. Probability 0.01
#'   \item 0.001. Probability 0.001
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableVIFinney1964()
GenTableVIFinney1964<-function(){
p<-c(.9,.7,.5,.3,.1,.05,.02,.01,.001)
t1 <- matrix(0,30,10)
t1[,1]<-1:30
for (n in 1:30) t1[n,2:10]<-qchisq(p, df=n, lower.tail =F)
l1 <- which(t1<1)
t1[l1] <- round(t1[l1],2)
t1[-l1] <- round(t1[-l1],1)
colnames(t1) <- c('Deg.freedom',p)
t1
}

#' Generate table VII from Finney1964 "Distribution of t"
#' @description Generates table VII from Finney1964 "Distribution of t"
#' @return table VII from Finney1964 "Distribution of t"
#' \itemize{
#'   \item Deg.freedom. Degrees of freedom
#'   \item 0.9. Probability 0.9
#'   \item 0.7. Probability 0.7
#'   \item 0.5. Probability 0.5
#'   \item 0.3. Probability 0.3
#'   \item 0.1. Probability 0.1
#'   \item 0.05. Probability 0.05
#'   \item 0.02. Probability 0.02
#'   \item 0.01. Probability 0.01
#'   \item 0.001. Probability 0.001
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableVIIFinney1964()
GenTableVIIFinney1964<-function(){
p<-c(.9,.7,.5,.3,.1,.05,.02,.01,.001)
m <- c(1:10,6:15*2,40,60,120, 1/0 )
t1 <- matrix(0,24,10)
t1[,1]<-m
for (n in 1:24) t1[n,2:10]<- qt(1-p/2, df=m[n], log.p=F)# 2 tails
colnames(t1) <- c('Deg.freedom',p)
t1[24,] <- round(t1[24,],3)
t1[1:23,] <- round(t1[1:23,],2)
t1[1,7:9] <- round(t1[1,7:9],1)
t1[2:4,10] <- round(t1[2:4,10],1)
t1[1,10] <- round(t1[1,10])
t1
}

#' Generate table VIII from Finney1964 "The Weighting Coefficient in Wadley's Problem"
#' @description Generates table VIII from Finney1964 "The Weighting Coefficient in Wadley's Problem"
#' @return table VIII from Finney1964 "The Weighting Coefficient in Wadley's Problem"
#' \itemize{
#'   \item Y. Expected probit
#'   \item w. Weighting Coefficient
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableVIIIFinney1964()
GenTableVIIIFinney1964<-function(){
Y <- 11:90/10;P <- pnorm(Y-5);Q<-1-P;Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)
#Y2 <- 51:90/10;P2 <- pnorm(Y2-5);Q2<-1-P2;Z2 <- 1/sqrt(2*pi)*exp(-0.5*(Y2-5)^2)
#cbind(Y=Y,WeightingCoef=Z^2/Q,Y=Y2,WeightingCoef=Z2^2/Q2)
t1 <- cbind(Y=Y,w=Z^2/Q)
t1[1:2,2] <- round(t1[1:2,2],8)
t1[3:5,2] <- round(t1[3:5,2],7)
t1[6:8,2] <- round(t1[6:8,2],6)
t1[9:80,2] <- round(t1[9:80,2],5)
t1
}

#' Generate table IX from Finney1964 "Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling"
#' @description Generates table IX from Finney1964 "Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling"
#' @return table IX from Finney1964 "Minimum Working Probit, Range, and Weighting Coefficient for Inverse Sampling"
#' \itemize{
#'   \item Y. Expected probit
#'   \item MinWorkProbit. Minimum working probit
#'   \item Range. Range 1/Z
#'   \item WeightingCoef. Weighting Coefficient
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' GenTableIXFinney1964()
GenTableIXFinney1964<-function(){
Y <- 50:90/10;P <- pnorm(Y-5);Q<-1-P;Z <- 1/sqrt(2*pi)*exp(-0.5*(Y-5)^2)
t1 <- cbind(Y=Y,MinWorkProbit=round(Y-Q/Z,4),Range=Q^2/Z,WeightingCoef=round(Z^2/P/Q^2,4))
t1[1:11,3] <- round(t1[1:11,3],5)
t1[12:20,3] <- round(t1[12:20,3],6)
t1[21:28,3] <- round(t1[21:28,3],7)
t1[29:34,3] <- round(t1[29:34,3],8)
t1[35:40,3] <- round(t1[35:40,3],9)
t1[41,3] <- round(t1[41,3],10)
t1
}

#' Generate table 26 from Finney1964 "The Function for Planning Tests of Mixtures of Two Poisons"
#' @description Generates table 26 from Finney1964 "The Function for Planning Tests of Mixtures of Two Poisons"
#' @return table 26 from Finney1964 "The Function for Planning Tests of Mixtures of Two Poisons"
#' \itemize{
#'   \item rho. toxicity
#'   \item 0.1. distance 0.1 log rho in the left of the probit regression line
#' ...
#'   \item 0.9. distance 0.9 log rho in the left of the probit regression line
#' }
#' @references Finney D. J., 1964
#' Probit analysis: a statistical treatment of the sigmoid response curve.
#' Cambridge University Press
#' @author Jose Gama
#' @examples
#' TestMix2poisons()
TestMix2poisons<-function() #100*(p-p^th)/(p-1)
{
t1<-matrix(0,20,10)
p<-c(1.1,1.5,2:9,5*2:6,40,50,60,80,100)
th<-1:9/10
t1[,1]<-p
for (n in 1:20) t1[n,2:10]<-100*(p[n]-p[n]^th)/(p[n]-1)
t1 <- round(t1,1)
t1[1:5,3:10]<-round(t1[1:5,3:10])
t1[6:19,4:10]<-round(t1[6:19,4:10])
t1[11:15,5:10]<-round(t1[11:15,5:10])
t1[15:19,6:10]<-round(t1[15:19,6:10])
t1[20,7:10]<-round(t1[20,7:10])
colnames(t1)<-c('rho',1:9/10)
t1
}

#' Horsfall-Barratt Scale for Measuring Plant Disease
#' @description Horsfall-Barratt Scale for Measuring Plant Disease
#' @param percentAffected numeric vector
#' @return Horsfall-Barratt Scale for Measuring Plant Disease
#' @references Horsfall, J. G.; Barratt, R. W., 1945
#' An Improved Grading System for Measuring Plant Disease. Phytopathology.
#' @author Jose Gama
ScaleHorsfallBarratt<-function(percentAffected)
{#Horsfall-Barratt scale
if ((percentAffected>100)|(percentAffected<0)) stop('The input must be between 0 and 100')
if (percentAffected==0) return(1)
if (percentAffected<3) return(2)
if (percentAffected<6) return(3)
if (percentAffected<12) return(4)
if (percentAffected<25) return(5)
if (percentAffected<50) return(6)
if (percentAffected<75) return(7)
if (percentAffected<87) return(8)
if (percentAffected<94) return(9)
if (percentAffected<97) return(10)
if (percentAffected<100) return(11)
12
}

#' Archer Scale for assessment of leaf damage
#' @description Archer Scale for assessment of leaf damage
#' @param percentAffected numeric vector
#' @return Archer Scale for assessment of leaf damage
#' @references Archer, T.L., 1987
#' Techniques for screening maize for resistance to mites. pp.178-183. In:
#' Mihn, J.A., Wiseman, B.R. and Davis, F.M. (Eds.). Proceedings of the International
#' symposium on methodologies for developing host plant resistance to maize insects.
#' CIMMYT, Mexico.
#' @author Jose Gama
ScaleArcher<-function(percentAffected)
{# Assessment of leaf damage
# PROTOCOLS FOR THE BIOLOGICAL EVALUATION OF PESTICIDES.pdf
if ((percentAffected>100)|(percentAffected<1)) stop('The input must be between 1 and 100')
if (percentAffected<11) return(1)
if (percentAffected<21) return(2)
if (percentAffected<31) return(3)
if (percentAffected<41) return(4)
if (percentAffected<51) return(5)
if (percentAffected<61) return(6)
if (percentAffected<71) return(7)
if (percentAffected<81) return(8)
if (percentAffected<91) return(9)
10
}

#' Gauhl’s modification of Stover’s severity scoring system
#' @description Gauhl’s modification of Stover’s severity scoring system
#' @param percentShowingSymptoms numeric, proportion of the leaf area showing symptoms
#' @return Gauhl-Stover scale
#' @references Gauhl F., 1994
#' Epidemiology and ecology of black Sigatoka (Mycosphaerella fijiensis Morlet)
#' on plantain and banana (Musa spp.) in Costa Rica, Central America.
#' INIBAP, Montpellier, France. 120pp).
#' @author Jose Gama
ScaleGauhlStover<-function(percentShowingSymptoms)
{#Gauhl-Stover scale
percentShowingSymptoms <- round(percentShowingSymptoms)
if ((percentShowingSymptoms>100)|(percentShowingSymptoms<0)) stop('The input must be between 0 and 100')
if (percentShowingSymptoms==0) return(0)
if (percentShowingSymptoms < 1) return(1)
if (percentShowingSymptoms <= 5) return(2)
if (percentShowingSymptoms <= 15) return(3)
if (percentShowingSymptoms <= 33) return(4)
if (percentShowingSymptoms <= 50) return(5)
6
}

#' WAAPP Pest Count scoring system
#' @description WAAPP Pest Count scoring system
#' @param percentLeafDamage numeric, percentage of leaf damage
#' @return WAAPP Pest Count Score
#' @references Environmental Protection Agency Chemicals Control And Managemenet Centre (ACCRA),  2012
#' Protocols for the biological evaluation of pesticides on Selected crops grown 
#' in both the humid and sahel regions of West africa.
#' West Africa Agriculture Productivity Programme (WAAPP).
#' @author Jose Gama
WAAPPpestCount<-function(percentLeafDamage)
{
percentLeafDamage <- round(percentLeafDamage)
if ((percentLeafDamage>100)|(percentLeafDamage<0)) stop('The input must be between 0 and 100')
if (percentLeafDamage==0) return(0)
if (percentLeafDamage <= 10) return(1)
if (percentLeafDamage <= 25) return(2)
if (percentLeafDamage <= 50) return(3)
4
}