Lethal Concentration Probit — LC

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)
#> # A tibble: 6 × 7
#>   nominal_dose tank  month  dose response survive total
#>          <dbl> <chr> <chr> <dbl>    <dbl>   <dbl> <dbl>
#> 1          0   S     May    0.19        0      20    20
#> 2          0.7 A     May    0.79        0      20    20
#> 3          0.7 L     May    0.8         1      19    20
#> 4          0.7 M     May    0.76        0      20    20
#> 5          1.1 B     May    1.36       13       8    21
#> 6          1.1 K     May    1.22       10       9    19

# 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"))
#> Warning: Are you using cbind(response, non-response) method as your y variable, if so do not weight the model. If you are using (response / total) method, model needs the total of test organisms per dose to weight the model properly

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

m
#> # A tibble: 2 × 19
#>       p     n  dose   LCL   UCL    se chi_square    df pgof_sig     h slope
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>      <dbl> <int>    <dbl> <dbl> <dbl>
#> 1    50    18  1.25  1.18  1.31  1.02       14.1    16    0.590     1  10.3
#> 2    99    18  2.11  1.95  2.35  1.05       14.1    16    0.590     1  10.3
#> # ℹ 8 more variables: slope_se <dbl>, slope_sig <dbl>, intercept <dbl>,
#> #   intercept_se <dbl>, intercept_sig <dbl>, z <dbl>, var_m <dbl>,
#> #   covariance <dbl>

# these are the same

m1
#> # A tibble: 2 × 19
#>       p     n  dose   LCL   UCL    se chi_square    df pgof_sig     h slope
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>      <dbl> <int>    <dbl> <dbl> <dbl>
#> 1    50    18  1.25  1.18  1.31  1.02       14.1    16    0.590     1  10.3
#> 2    99    18  2.11  1.95  2.35  1.05       14.1    16    0.590     1  10.3
#> # ℹ 8 more variables: slope_se <dbl>, slope_sig <dbl>, intercept <dbl>,
#> #   intercept_se <dbl>, intercept_sig <dbl>, z <dbl>, var_m <dbl>,
#> #   covariance <dbl>

# 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