LC_probit {ecotox}

Lethal Concentration Probit

Description

Calculates lethal concentration (LC) and its fiducial confidence limits (CL) using a probit analysis according to Finney 1971, Wheeler et al. 2006, and Robertson et al. 2007.

Usage

LC_probit(formula, data, p = NULL, weights = NULL,
          subset = NULL, log_base = NULL, log_x = TRUE,
          het_sig = NULL, conf_level = NULL, conf_type = NULL,
          long_output = TRUE)

Arguments

formula

an object of class formula or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under Details.

data

an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which LC_probit is called.

p

Lethal Concentration (LC) values for given p, example will return a LC50 value if p equals 50. If more than one LC value wanted specify by creating a vector. LC values can be calculated down to the 1e-16 of a percentage (e.g. LC99.99). However, the tibble produced can round to nearest whole number.

weights

vector of 'prior weights' to be used in the fitting process. Only needs to be supplied if you are taking the response / total for your response variable within the formula call of LC_probit. Otherwise if you use cbind(response, non-response) method you do not need to supply weights. If you do the model will be incorrect. If you don't supply weights there is a warning that will help you to make sure you are using one method or the other. See the following StackExchanges post about differences cbind() function in R for a logistic regression and input format for binomial glm.

subset

allows for the data to be subseted if desired. Default set to NULL.

log_base

default is 10 and will be used to calculate results using the anti of log10() given that the x variable has been log10 transformed. If FALSE results will not be back transformed.

log_x

default is TRUE and will calculate results using the antilog of determined by log_base given that the x variable has been log() transformed. If FALSE results will not be back transformed.

het_sig

significance level from person's chi square goodness-of-fit test (pgof) that is used to decide if a heterogeneity factor is used. NULL is set to 0.15.

conf_level

adjust confidence level as necessary or NULL set at 0.95.

conf_type

default is "fl" which will calculate fudicial confidence limits per Finney 1971. If set to "dm" the delta method will be used instead.

long_output

default is TRUE which will return a tibble with all 17 variables. If FALSE the tibble returned will consist of the p level, n, the predicted LC for given p level, lower and upper confidence limits.

Value

Returns a tibble with predicted LC for given p level, lower CL (LCL), upper CL (UCL), Pearson's chi square goodness-of-fit test (pgof), slope, intercept, slope and intercept p values and standard error, and LC variance.

References

Finney, D.J., 1971. Probit Analysis, Cambridge University Press, Cambridge, England, ISBN: 052108041X

Wheeler, M.W., Park, R.M., and Bailey, A.J., 2006. Comparing median lethal concentration values using confidence interval overlap or ratio tests, Environ. Toxic. Chem. 25(5), 1441-1444.10.1897/05-320R.1

Robertson, J.L., Savin, N.E., Russell, R.M. and Preisler, H.K., 2007. Bioassays with arthropods. CRC press. ISBN: 9780849323317

Examples

head(lamprey_tox)

# within the dataframe used, control dose, unless produced a value
# during experimentation, are removed from the dataframe,
# as glm cannot handle values of infinite. Other statistical programs
# make note of the control dose but do not include within analysis

# calculate LC50 and LC99

m <- LC_probit((response / total) ~ log10(dose), p = c(50, 99),
               weights = total,
               data = lamprey_tox[lamprey_tox$nominal_dose != 0, ],
               subset = c(month == "May"))
# OR

m1 <- LC_probit(cbind(response, survive) ~ log10(dose), p = c(50, 99),
               data = lamprey_tox[lamprey_tox$nominal_dose != 0, ],
               subset = c(month == "May"))

# view calculated LC50 and LC99 for seasonal toxicity of a pisicide,
# to lamprey in 2011

m

# these are the same

m1

# dose-response curve can be plotted using 'ggplot2'

# library(ggplot2)

# lc_may <- subset(lamprey_tox, month %in% c("May"))

# p1 <- ggplot(data = lc_may[lc_may$nominal_dose != 0, ],
#              aes(x = log10(dose), y = (response / total))) +
#   geom_point() +
#   geom_smooth(method = "glm",
#               method.args = list(family = binomial(link = "probit")),
#               aes(weight = total), colour = "#FF0000", se = TRUE)

# p1

# calculate LC50s and LC99s for multiple toxicity tests, June, August, and September

j <- LC_probit((response / total) ~ log10(dose), p = c(50, 99),
        weights = total,
        data = lamprey_tox[lamprey_tox$nominal_dose != 0, ],
        subset = c(month == "June"))

a <- LC_probit((response / total) ~ log10(dose), p = c(50, 99),
        weights = total,
        data = lamprey_tox[lamprey_tox$nominal_dose != 0, ],
        subset = c(month == "August"))

s <- LC_probit((response / total) ~ log10(dose), p = c(50, 99),
        weights = total,
        data = lamprey_tox[lamprey_tox$nominal_dose != 0, ],
        subset = c(month == "September"))

# group results together in a dataframe to plot with 'ggplot2'

results <- rbind(m[, c(1, 3:8, 11)], j[,c(1, 3:8, 11)],
                 a[, c(1, 3:8, 11)], s[, c(1, 3:8, 11)])
results$month <- factor(c(rep("May", 2), rep("June", 2),
                          rep("August", 2), rep("September", 2)),
                        levels = c("May", "June", "August", "September"))

# p2 <- ggplot(data = results, aes(x = month, y = dose,
#                              group = factor(p), fill = factor(p))) +
#   geom_col(position = position_dodge(width = 0.9), colour = "#000000") +
#   geom_errorbar(aes(ymin = LCL, ymax = UCL),
#                 size = 0.4, width = 0.06,
#                 position = position_dodge(width = 0.9))

# p2