#' --- #' title: 'Chronic and cumulative toxicity of boscalid' #' author: ['Gilles San Martin', 'Louis Hautier', 'Noa Simon'] #' date: "`r Sys.setlocale('LC_ALL', 'en_GB.UTF-8'); format(Sys.time(), '%d %B %Y')`" #' #' output: #' #pdf_document: #' bookdown::pdf_document2: #' toc: true #' toc_depth: 5 #' fig_caption: true #' number_sections: true #' latex_engine: pdflatex #' # latex_engine : xelatex #' highlight: default # default # espresso # zenburn #' geometry: margin=0.7in #' documentclass: article # report #' fontsize: 11pt #' fontfamily: lmodern # mathpple # times # utopia # lmodern # bookman # mathpazo # for pdflatex engine #' header-includes: #' - \usepackage{float} # probably unnecessary with fig.pos = "H" #' - \usepackage{fancyvrb} #' - \usepackage{graphicx} # needed for fig.align = 'center' option #' - \usepackage{pdflscape} # to have some pages in landscape orientation #' --- # /* # ---------------- This part is just for report production ------------------- # */ #+ echo = FALSE, warning = FALSE, error = FALSE, message = FALSE # source("/home/gilles/stats/makereport.R") # makereport("Chronic_toxicity.R", output = c("latex"), odtdpi = 300, wd= "/home/gilles/stats/Rprojects/20170727_chronic_toxicity") library(knitr) opts_chunk$set(echo=FALSE, results= "markup", cache = TRUE, dev = c("cairo_pdf"), fig.width = 10/2.54, fig.height= 10/2.54, fig.align='center', dpi = 500, dev.args=list(bg='white'), fig.cap = ' ', fig.pos = "H", warning = FALSE, error = FALSE, message = FALSE, tidy = FALSE, tidy.opts=list(width.cutoff=95)) library(pander) panderOptions("table.style" , "rmarkdown") panderOptions("table.split.table" , Inf) panderOptions("table.alignment.rownames" , "right") options(width = 95) #' #' #' ********************************************************************************* #' #' \newpage #' # /* # ---------------- Load packages and scripts ------------------- # */ #+ echo = FALSE source("/home/gilles/stats/mytoolbox.R") setwd("/home/gilles/stats/Rprojects/20170727_chronic_toxicity") # load ggplot, change the default theme and change the locale (language = English) library(ggplot2) #+ results = "hide", echo = FALSE Sys.setlocale("LC_ALL", 'en_GB.UTF-8') mytheme <- theme_bw(10) + theme(axis.text.x=element_text(size=8), legend.key = element_rect(color = NA)) theme_set(mytheme) library(drc) library(reshape) library(lme4) library(car) library(multcomp) #' #' \newpage #' # /* # ---------------- Introduction - context ------------------- # */ #' # Introduction - Context #' #' ## Context #' #' The full datasets, R code, and raw results output related to this document are #' available in a public figshare repository : https://figshare.com/s/865e87feaad34c095bbd #' #' This report presents the detailed analyses performed for the following paper : #' Chronic and cumulative toxicity for honey bees of the fungicide boscalid revealed #' with a time-to-death approach but not with the standard 10 days test #' by Noa Simon-Delso, Gilles San Martin, Etienne Bruneau, Louis Hautier #' #' The context of this study is twofold : #' #' 1. the evaluation of the methodologies used to assess the (chronic) toxicity of #' pesticides on adult bees #' 1. a specific fungicide (boscalid) whose toxicity for honey bees has started to attract #' attention in recent years #' #' **Methods to estimate toxicity of pesticides on bees** #' #' * It is important to find international agreements on the methodologies in order to #' make the ecotoxicological studies used for risk assessment comparable between labs and #' between products. #' * However existing methods might need evaluation, validation and improvement. #' * It might be argued that it is better (at least from a regulatory point of view) #' to have a standardized, stable, but maybe imperfect method than to gradually improve #' existing methods with a continuously changing protocol that makes the studies not #' comparable. #' However, even in that case, it might be useful to understand the limitations of #' existing protocols. #' #' Many questions can arise, in particular concerning oral chronic toxicity tests : #' #' * Is the standard duration of 10 days for chronic tests on adults enough to detect #' chronic troxicity ? #' * How to detect cumulative toxicity or "time reinforced toxicity" caused by bioaccumulation or other mechanisms? #' * Instead of a fixed time evaluation (eg mortality after 10 days) how can the #' kinetics of the toxicity improve our evaluation of the environmental risk of a pesticide ? #' * How does the syrup consumption change over time and between treatments and what is the impact #' on the ecotoxicological results ? #' * What is the impact of evaporation on the evaluation of the dose consumed by the bees ? #' #' **Boscalid** #' #' * Boscalid is a persistent, systemic and widely used fungicide. #' * It is one of the most frequently found pesticides within beehives in Belgium. #' * It has also been found in pollen pellets collected by bees very late in the season #' (October) at periods of the year where the pesticide is unlikely to be used in the field. #' #' --> Hence, it seems that bees are frequently exposed and during extended periods to #' this pesticide. The question of its chronic toxicity is therefore legitimate. #' #' * A former observative field study has shown a positive relationship between unexplained #' winter colony health problems and the fungicides load (mainly boscalid) found within the beehive. #' * It is however only a correlation and it remained unclear whether boscalid itself has #' a direct, causal effect in the observed relationship. #' #' A few lab studies have been recently performed that allow to explore a potential #' causal relationship between the presence of boscalid and the health problem observed #' in the field : #' #' * Only the acute toxicity of boscalid on bees has been evaluated before commercialization #' * A recent study has shown no (or limited) chronic toxicity of pure boscalid on honey bees larvae #' * A recent study has shown a synergist effect of boscalid on a neo-nicotinoid. Previous #' monitoring studies done in the region to know the level of contamination of beehives #' did not detect many residues of neonicotinoids, but the sensitivity of the analyses #' used was not good enough. Therefore, bearing in mind the amounts of neonicotinoids sold #' in Belgium (official data, not shown), a potential exposure cannot be excluded. #' * Several studies have shown no or limited synergist effect of boscalid on acaricides #' on honey bees or other arthropods #' #' ## Aims of this project #' #' 1. Test the chronic oral toxicity of the fungicide boscalid (as a formulated product) #' on adult bees during a longer period than the classical 10 days duration of standard chronic tests #' 1. Test potential bioacculation of the product using different methodological approaches #' proposed in the literature #' 1. Compare the feeding rate over time, between treatments and evaluate the levels of #' syrup evaporation and their consequences on the results #' #' #' ## Methods and data analysis #' #' We used Cantus as commercial formulation of boscalid (50% w/w). #' Hence, the toxicity we observe here could be due to the co-formulants, #' not boscalid itself... It was not possible to use pure boscalid in this study because #' it was impossible to dilute it efficiently. In a former study on larvae, we used pure #' boscalid diluted in royal jelly that showed no or limited toxicity on larvae. #' #' In the whole study, we make the distinction between #' #' * `Concentration` which is the concentration of the #' product within the syrup provided to the bees. The units are `mg a.i./l syrup`. #' * `Dose` which is the quantity of product really consumed by the bees through their food. #' It can be expressed as `average mg a.i./(bee*day)` for the daily average dose #' or `total mg a.i./bee` for the total cumulative dose. #' #' We computed 4 types of dose or time to effect (mortality here) statistics. #' #' * `LTx (days)`: Lethal Time based on the model : Mortality vs Time #' * `LCx (mg a.i./l syrup)` : Lethal Concentration based on the model : Mortality vs Concentration #' * `LDDx (average mg a.i./(bee*day))`: Lethal Dietary Dose based on the model : Mortality #' vs Average Dose consumed per bee and per day #' * `LCDx (total mg a.i./bee)`: Lethal Cumulative Dose based on a model of Mortality vs #' Cumulative (total) Dose consumed per bee since D0 #' #' x represent the intended mortality rate (10%, 20%,...). #' NB : for the sake of simplicity, we will some time use "EDx" (Effect Dose) as a generic describer of #' LCx, LDDx, LCDx and even LTx, in line with the terminology used in the `drc` package. #' #' We chose a 3 parameters Weibull 2 model (_sensu_ Ritz et al. 2010) to model the #' uncorrected mortality vs time of concentration/doses. Three other types of models #' have been tested for all relationships : #' logistic, log-logistic and Weibull 1 (_sensu_ Ritz et al. 2010), all with 3 parameters. #' The Weibull 2 and logistic models performed most of the time better (better fit, lower AIC) #' than the two other types of models. #' The Weibull 2 model was more stable than the logistic models for which the standard errors of the EDx #' were often impossible to compute. Both models provided very similar estimates of EDx in any cases. #' See section \@ref(ModelChoice) for more details on the comparison of the 4 types of models #' and on the structure of these models. The `raw_output` directory provides the detailed results #' for the 4 types of models at each time, each concentration and for ED10 up to ED90 (by step of 10%). #' #' The LCx, LDDx and LCDx estimates are corrected for the mortalities in the control because we use #' a 3 parameters model that estimates a parameter for the lower asymptote of the sigmoïd #' dose response curve (estimating the mortality when the concentration/dose is 0) #' and we compute the LCx, LDDx and LCDx relative to this asymptote. #' #' For LTx we also use a 3 parameters model, however the lower asymptote in that case estimates #' the mortality at D0 which is always 0. Hence, the LTx values as we calculate them here #' are not corrected for mortality in the control. #' #' Applying correctly a mortality correction might be difficult because the total number #' of bees are not exactly equal at the start of the experiment. Abbott's formula is #' applicable in such unbalanced cases. However binomial models require that the #' response variable (the % of dead bees - corrected or not) is a ratio of two integers #' which will most of the time not be the case after Abbott's correction when the sample #' sizes are unequal. This is the reason why we prefer to use 3 parameters models on #' uncorrected mortalities rather than the more classical 2 parameters models on corrected mortalities. #' NB : a 3 parameters model on uncorrected mortalities gives very similar results #' (in terms of Effect Dose estimates) than a 2 parameters model on corrected mortalities #' (see section \@ref(cormortality) for more details and an empiric demonstration). #' #' #' #' #' \newpage #' # /* # ---------------- Raw data and computations ------------------- # */ #' #' #' \newpage #' #' #' # Raw data and computations #' #' The data from the last day (D33) have been removed because there is no data about #' syrup consumption on this last day #' #' Raw data : #' #' * Treat : Treatment : CANTUS (= commercial name of the fungicide containing boscalid) - #' 2CNA (not used here) - TS (Toxic Standard) - CONTROL #' * Rep : Replicate #' * Conc (mg a.i./l syrup) : concentration of the a.i. in the syrup (sugar solution) #' The syrup is a 50% w/v sugar solution (e.g 1 kg sugar in 1 l water). #' This solution has a density of 1.23 g/ml for : 50% sugar with a purity of 100% at ~25°C, #' [source](http://www.sugartech.co.za/density/index.php). #' This Conc value is used to compute the LCx (Lethal Concentration x%). #' * Eff (bees) : Total number of Bees on D0 (Effectif in French). Generally Eff=10. #' In few cases eff = 9 because some bees died btween the moment they were collected for #' the test and the moment to start the test (and there were no other available bees of the same age) #' * Day (days): Experiment day starting at Day 0 #' * Alive (bees): number of living bees on day d #' * Dead (bees): number of dead bees since day D0 #' * Abnormal (bees) : number of bees showing an abnormal behavior #' * Weightt0 (g) : measure on day d of the weight of the new syrup seringue provided #' * Weightt1 (g) : measure on day d+1 of the Weight of the old syrup seringue provided on day d #' #' The rest of the columns are calculated on the basis of these first columns #' #' * Time (days) : the day number starting at 0 (numeric version of Day) #' * MortRate : mortality rate = Dead/Eff #' * AliveRate : "living rate" = Alive/Eff (not used) #' * Conso (g/cage): syrup consumption per cage on a specific day = Weightt0-Weightt1 #' (and = 0 when the data are not available, because all bees were already dead) #' * ConsoBee : syrup consumption by the living bees on a specific day = Conso/Alive (g/living bee) #' * Dose = (mg a.i./cage) Dose of ai ingested per cage by the living bees on a specific day (mg ai/cage) #' = `Conc*Conso/(1.23*1000)` (taking into account the syrup density) #' * DoseBee (mg ai/living bee) = Dose of ai ingested per living bees on a specific day #' = `Conc*ConsoBee/(1.23*1000)` #' * CumDose (total mg ai/ cage) = total amount of a.i. consumed by the bees taht are dead on day d since D0 = #' `cumsum(Dose) - Dose` for each combination of Treatment, concentration and replicate #' The dose consumed on day d is consumed by the living bees but not by the dead bees noted on that day. #' This is the reason why we remove the dose from dau d in the calculations. Otherwise said #' for day d, we compute the cumulative dose up to day -1 and the dose on D0 is 0 #' * CumDoseBee (total mg ai/ dead bee) : the same per bee = `cumsum(DoseBee) - DoseBee`. #' these values are used for the computations of LCDx (Lethal Cumulated Doses x%) #' * MeanDose (average mg ai / (cage * day) ): Average Dose in a whole cage consumed by #' the dead bees on day d since D0 (ie this is a kind of cumulated average). #' Basically this is = `CumDose/Day` for each #' replicate and each treatment. However the first value is not divided by 0 but by 1 to #' obtain an average dose = 0 at D0. And once the mortality reaches 100% the mean dose consumed #' remains constant. #' * MeanDoseBee (average mg ai / (bee * day) ) : The same but per bee i.e. average dose #' consumed per bees that are dead on day d since D0. Basically this is = CumDoseBee/Day with the #' same restrictions as MeanDose. #' MeanDoseBee is used to compute the LDDx (Lethal Dietary Dose x%) #' d <- read.csv2("data/RAW_Chronic_toxicity-Noa.csv") d <- d[d$Day != "D33",] d$Time <- as.numeric(substr(d$Day, 2,3)) d$MortRate <- d$Dead/d$Eff d$AliveRate <- d$Alive/d$Eff # syrup consumption per cage (g) d$Conso <- ifelse(d$Alive == 0, 0, d$Weightt0-d$Weightt1) # syrup consumption per bee (g/bee) d$ConsoBee <- ifelse(d$Alive == 0, 0, d$Conso/d$Alive) # Dose of ai ingested per cage (mg ai/cage) on a specific day d$Dose <- d$Conc*d$Conso/(1.23*1000) # Dose of ai ingested per living bee (mg ai/ alive bee) on a specific day d$DoseBee <- d$Conc*d$ConsoBee/(1.23*1000) # cumulative consumption of a.i. (mg ai/ dead bee) # = total amount of a.i. consumed by the dead bees d$CumDose <- ave(d$Dose, d$Treat,d$Conc,d$Rep, FUN = cumsum) - d$Dose # per cage d$CumDoseBee <- ave(d$DoseBee, d$Treat,d$Conc,d$Rep, FUN = cumsum) - d$DoseBee # per bee # Average consumption of a.i. up to day d # # Compute the denominator to divise the CumDose and CumDoseBee # Basically it is the day number but the first value is set to 1 (divided by 1) to # avoid a division per 0 and once the mortality reaches 100% the denominator is fixed # to obtain a constant mean dose denominator <- ave(d$Alive, d$Treat,d$Conc,d$Rep, # apply the following function on d$Alive separately for each combination of Treat, Conc, Rep FUN = function(x) { c(1, # 1 for the first value to obtain a 0 as result (0/1) seq_along(x[x>0]), # Sequence from 1 to the number of days with living bees rep(sum(x>0), sum(x==0)))[-length(x)] # repeat the day when you reached 100% mortality (and remove the last value) }) d$MeanDose <- d$CumDose/denominator d$MeanDoseBee <- d$CumDoseBee/denominator # comparison with the excell file discussed with Louis & Noa # to check if we aggree on the calculations # d[d$Treat == "TS" & d$Rep == 2,] # d[d$Treat == "CONTROL" & d$Rep == 2,] # d[d$Treat == "CANTUS" & d$Rep == 2 & d$Conc == 18000,] #' \newpage #' summary(d) store <- d # store the main result # write raw results on the disc if(all(!grepl("^raw_results$", list.dirs()))) {dir.create("raw_results")} tmp <- d tmp <- tmp[tmp$Treat != "2CNA",] tmp$Treat <- factor(tmp$Treat, levels = c("CONTROL", "CANTUS", "TS")) tmp <- tmp[order(tmp$Treat, tmp$Conc, tmp$Rep, tmp$Time),] write.csv2(tmp, "raw_results/raw_data_with_computed_doses.csv", row.names = FALSE) # Main data taking into account the evaporation # NB in the end the object `d`` contains the data without evarporation correction # and the object `devap` contains the same dataset bu after correction for evaporation # using the average evaporation across all days and all cages (after removing one outlier) # # Compute average evaporation after removing one outlier evap <- read.csv2("data/Evaporation.csv") evap$Evaporation <- evap$weightt0 - evap$weightt1 avg_evap <- mean(evap[evap$Evaporation <0.20, "Evaporation"]) d <- read.csv2("data/RAW_Chronic_toxicity-Noa.csv") d <- d[d$Day != "D33",] d$Time <- as.numeric(substr(d$Day, 2,3)) d$MortRate <- d$Dead/d$Eff d$AliveRate <- d$Alive/d$Eff # syrup consumption per cage (g) d$Conso <- ifelse(d$Alive == 0, 0, d$Weightt0-d$Weightt1 - avg_evap) # For cantus/boscalid, 2 values are negative, they are set to 0 d[d$Conso < 0,"Conso"] <- 0 # syrup consumption per bee (g/bee) d$ConsoBee <- ifelse(d$Alive == 0, 0, d$Conso/d$Alive) # Dose of ai ingested per cage (mg ai/cage) on a specific day d$Dose <- d$Conc*d$Conso/(1.23*1000) # Dose of ai ingested per living bee (mg ai/ alive bee) on a specific day d$DoseBee <- d$Conc*d$ConsoBee/(1.23*1000) # cumulative consumption of a.i. (mg ai/ dead bee) # = total amount of a.i. consumed by the dead bees d$CumDose <- ave(d$Dose, d$Treat,d$Conc,d$Rep, FUN = cumsum) - d$Dose # per cage d$CumDoseBee <- ave(d$DoseBee, d$Treat,d$Conc,d$Rep, FUN = cumsum) - d$DoseBee # per bee # Average consumption of a.i. up to day d # # Compute the denominator to divise the CumDose and CumDoseBee # Basically it is the day number but the first value is set to 1 (divided by 1) to # avoid a division per 0 and once the mortality reaches 100% the denominator is fixed # to obtain a constant mean dose denominator <- ave(d$Alive, d$Treat,d$Conc,d$Rep, # apply the following function on d$Alive separately for each combination of Treat, Conc, Rep FUN = function(x) { c(1, # 1 for the first value to obtain a 0 as result (0/1) seq_along(x[x>0]), # Sequence from 1 to the number of days with living bees rep(sum(x>0), sum(x==0)))[-length(x)] # repeat the day when you reached 100% mortality (and remove the last value) }) d$MeanDose <- d$CumDose/denominator d$MeanDoseBee <- d$CumDoseBee/denominator devap <- d # change the name d <- store # set back the main dataset (without correction to the original name) # # plot(devap$MeanDoseBee ~ d$MeanDoseBee) # abline(0,1) # Function to compare different models (AIC, goodness of fit test) mselect2 <- function(l, dec = 3) { l <- l[!is.na(l)] res <- lapply(X = l, FUN = function(x){ try(cbind(modelFit(x)[2,3:5], AIC = AIC(x), ED50 = ED(x, 50, interval = "delta", display = FALSE, od = FALSE)), TRUE) }) res <- res[sapply(res, class) == "data.frame"] res <- do.call(rbind, res) res <- res[order(res[,"AIC"]),-6] colnames(res) <- c("df", "Chisq", "p", "AIC", "ED50", "LowerCI", "UpperCI") round(res,dec) } # /* # ---------------- Compute Lethal concentrations ------------------- # */ # In this section we compute 4 type of dose response curves for each day between # D8 and D25 and for each type of effect : LC, LDD, LCD and LT. # Days for which we want to estimate models and LCx # NB at lower and higher values of days the models are unstable or do not # converge because you have too few dead or to much dead... timevalues <- 8:25 LC_AIC <- vector("list", length = length(timevalues)) LC <- vector("list", length = length(timevalues)) for(i in 1:length(timevalues)) { time <- timevalues[i] PctMort <- seq(10,90,10) dataset <- d[d$Treat %in% c("CONTROL", "CANTUS") & d$Time == time,] # 4 type of models mLC_L <- drm(MortRate ~ Conc, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, NA, 1, NA, 1), names = c("b", "c", "d", "e", "f")) , type = "binomial") mLC_LL <- drm(MortRate ~ Conc, weights = Eff, data = dataset, fct = LL.3u(), type = "binomial") mLC_W1 <- drm(MortRate ~ Conc, weights = Eff, data = dataset, fct = W1.3u(), type = "binomial") mLC_W2 <- drm(MortRate ~ Conc, weights = Eff, data = dataset, fct = W2.3u(), type = "binomial") mLC <- list(mLC_L, mLC_LL, mLC_W1, mLC_W2) names(mLC) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLC) res <- data.frame( Time = time, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LC_AIC[[i]] <- res # Compute Effect doses (here LC) res <- lapply(mLC, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Time = time, PctMort = PctMort, Model = rep(names(mLC), each = length(PctMort)), res ) LC[[i]] <- res } LC_AIC <- do.call(rbind, LC_AIC) LC <- do.call(rbind, LC) # /* # ---------------- Compute Lethal Dietary Dose ------------------- # */ # days for which we want to estimate models and LCx # NB lower and higher values of days the models are unstable or do not # converge because you have too few dead or to much dead... timevalues <- 8:25 LDD_AIC <- vector("list", length = length(timevalues)) LDD <- vector("list", length = length(timevalues)) for(i in 1:length(timevalues)) { time <- timevalues[i] PctMort <- seq(10,90,10) dataset <- d[d$Treat %in% c("CONTROL", "CANTUS") & d$Time == time,] # 4 type of models mLC_L <- drm(MortRate ~ MeanDoseBee, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, NA, 1, NA, 1), names = c("b", "c", "d", "e", "f")) , type = "binomial") mLC_LL <- drm(MortRate ~ MeanDoseBee, weights = Eff, data = dataset, fct = LL.3u(), type = "binomial") mLC_W1 <- drm(MortRate ~ MeanDoseBee, weights = Eff, data = dataset, fct = W1.3u(), type = "binomial") mLC_W2 <- drm(MortRate ~ MeanDoseBee, weights = Eff, data = dataset, fct = W2.3u(), type = "binomial") mLC <- list(mLC_L, mLC_LL, mLC_W1, mLC_W2) names(mLC) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLC) res <- data.frame( Time = time, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LDD_AIC[[i]] <- res # Compute Effect doses (here LDD) res <- lapply(mLC, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Time = time, PctMort = PctMort, Model = rep(names(mLC), each = length(PctMort)), res ) LDD[[i]] <- res } LDD_AIC <- do.call(rbind, LDD_AIC) LDD <- do.call(rbind, LDD) # /* # ---------------- Compute Lethal cumulative Dose ------------------- # */ # days for which we want to estimate models and LCx # NB lower and higher values of days the models are unstable or do not # converge because you have too few dead or to much dead... timevalues <- 8:25 LCD_AIC <- vector("list", length = length(timevalues)) LCD <- vector("list", length = length(timevalues)) for(i in 1:length(timevalues)) { time <- timevalues[i] PctMort <- seq(10,90,10) dataset <- d[d$Treat %in% c("CONTROL", "CANTUS") & d$Time == time,] # 4 type of models mLC_L <- drm(MortRate ~ CumDoseBee, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, NA, 1, NA, 1), names = c("b", "c", "d", "e", "f")) , type = "binomial") mLC_LL <- drm(MortRate ~ CumDoseBee, weights = Eff, data = dataset, fct = LL.3u(), type = "binomial") mLC_W1 <- drm(MortRate ~ CumDoseBee, weights = Eff, data = dataset, fct = W1.3u(), type = "binomial") mLC_W2 <- drm(MortRate ~ CumDoseBee, weights = Eff, data = dataset, fct = W2.3u(), type = "binomial") mLC <- list(mLC_L, mLC_LL, mLC_W1, mLC_W2) names(mLC) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLC) res <- data.frame( Time = time, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LCD_AIC[[i]] <- res # Compute Effect doses (here LC) res <- lapply(mLC, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Time = time, PctMort = PctMort, Model = rep(names(mLC), each = length(PctMort)), res ) LCD[[i]] <- res } LCD_AIC <- do.call(rbind, LCD_AIC) LCD <- do.call(rbind, LCD) # /* # ---------------- Compute Lethal Time ------------------- # */ # concentration values concvalues <- c(0, 1125, 2250, 4500, 9000, 18000) LT_AIC <- vector("list", length = length(concvalues)) LT <- vector("list", length = length(concvalues)) for(i in 1:length(concvalues)) { conc <- concvalues[i] PctMort <- seq(10,90,10) dataset <- d[d$Treat %in% c("CONTROL", "CANTUS") & d$Conc == conc, ] # 4 type of models mL <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, NA, 1, NA, 1), names = c("b", "c", "d", "e", "f")) , type = "binomial") mLL <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = LL.3u(), type = "binomial") mW1 <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = W1.3u(), type = "binomial") mW2 <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = W2.3u(), type = "binomial") mLT <- list(mL, mLL, mW1, mW2) names(mLT) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLT) res <- data.frame( Conc = conc, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LT_AIC[[i]] <- res # Compute Effect doses (here LT) res <- lapply(mLT, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Conc = conc, PctMort = PctMort, Model = rep(names(mLT), each = length(PctMort)), res ) LT[[i]] <- res } LT_AIC <- do.call(rbind, LT_AIC) LT <- do.call(rbind, LT) # Lethal Time for the Toxic Standard # conc <- 1.5 PctMort <- seq(10,90,10) dataset <- d[d$Treat %in% c("TS") , ] # 4 type of models mL <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, NA, 1, NA, 1), names = c("b", "c", "d", "e", "f")) , type = "binomial") mLL <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = LL.3u(), type = "binomial") mW1 <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = W1.3u(), type = "binomial") mW2 <- drm(MortRate ~ Time, weights = Eff, data = dataset, fct = W2.3u(), type = "binomial") mLT <- list(mL, mLL, mW1, mW2) names(mLT) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLT) res <- data.frame( Conc = conc, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LT_TS_AIC <- res # Compute Effect doses (here LT) res <- lapply(mLT, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Conc = conc, PctMort = PctMort, Model = rep(names(mLT), each = length(PctMort)), res ) LT_TS <- res # Write all raw results on the disc # if(all(!grepl("raw_results$", list.dirs()))) {dir.create("raw_results")} write.csv2(LT, "raw_results/LT.csv", row.names = FALSE) write.csv2(LT_AIC, "raw_results/LT_AIC.csv", row.names = FALSE) write.csv2(LT_TS, "raw_results/LT_TS.csv", row.names = FALSE) write.csv2(LT_TS_AIC, "raw_results/LT_TS_AIC.csv", row.names = FALSE) write.csv2(LC, "raw_results/LC.csv", row.names = FALSE) write.csv2(LC_AIC, "raw_results/LC_AIC.csv", row.names = FALSE) write.csv2(LDD, "raw_results/LDD.csv", row.names = FALSE) write.csv2(LDD_AIC, "raw_results/LDD_AIC.csv", row.names = FALSE) write.csv2(LCD, "raw_results/LCD.csv", row.names = FALSE) write.csv2(LCD_AIC, "raw_results/LCD_AIC.csv", row.names = FALSE) #' #' \newpage #' # /* # ---------------- Syrup consumption ------------------- # */ #' #' # Syrup consumption and evaporation #' #' #' ## Evaporation data {#evaporation} #' #' We collected partial data about the difference of syrup weight that might be due #' to simple evaporation by adding syringes with syrup into cages without bees or with 10 #' dead bees during some days of the experiment (when cages were free because all bees were dead). #' #' One rather extreme data seems to be a transcription or measurement error. #' If we don't take this outlier into consideration, it appears that the weight loss #' due to evaporation is between 0.04 and 0.1 g per cage and per day. There is clear variation #' from day to day so it appears that a daily measurement would be better to correct for evaporation #' for each day. The presence of dead bees does not seem to affect the values #' however the presence of living bees might change the evaporation. #' #' We compared the results with and without correction for evaporation (using the average #' evaporation value after excluding the outlier). The impact of this correction of the daily #' consumption is discussed in section \@ref(consumption). For the rest of the analyses #' (ea LCx, LDDx, LCDx,... and cumulative toxicity) there were little differences when we applied #' the correction or not (results not shown here). In the rest of the document we show the #' results without evaporation consumption because such kind of ecotoxicilogical studies #' have historically not corrected for evaporation. #' A more detailed study about the daily evaporation with different number of living be should #' be performed separately. #' #' Note : As the consumption is generally measured as a consumption per bee, the measurement error #' caused by evaporation will have a higher impact when the mortality is high (ie when #' there are few living bees). #' evap <- read.csv2("data/Evaporation.csv") evap$Evaporation <- evap$weightt0 - evap$weightt1 #+ fig.width = 8/2.54, fig.height = 8/2.54 # dev.new(width = 8/2.54, height = 8/2.54) ggplot(evap, aes(y = Evaporation, x = Day, shape = Empty_cage)) + geom_point(size = 3) + ylim(0, 0.15) + ylab("Evaporation (g/cage)") + scale_shape_manual(name = "Cage", values = c(1,19)) + scale_x_continuous(breaks = seq(6, 18, 2)) + theme(legend.position = "top") #' Graph including the extreme value #+ fig.width = 8/2.54, fig.height = 8/2.54 # dev.new(width = 8/2.54, height = 8/2.54) ggplot(evap, aes(y = Evaporation, x = Day, shape = Empty_cage)) + geom_point(size = 3) + ylab("Evaporation (g/cage)") + scale_shape_manual(name = "Cage", values = c(1,19)) + scale_x_continuous(breaks = seq(6, 18, 2)) + theme(legend.position = "top") #' \newpage #' #' ## Consumption of syrup between doses and over time {#consumption} #' #' ### Kinetics of the consumption without evaporation correction #' #' In the following graph, we look at the consumption of syrup (not the dose of a.i.) consumed #' by the bees (without evaporation correction). #' Each gray line represent the data of one cage (replicate) for each concentration #' (each facet of the graph) of boscalid/cantus. The red line is a loess smoother #' (locally weighted polynomial regression) hat shows the general trend. #' The data displayed are the mean syrup consumption per bee and per day without evaporation #' correction. #' #' Main results : #' #' * There is a very large variation of consumption from day to day #' * At most concentrations (but not in the control), there is a strong increase of syrup #' consumption per bee at the end of the experiment, above 0.08 g/bee #' #+ fig.width = 16/2.54, fig.height = 10/2.54 # dev.new(width = 16/2.54, height = 10/2.54) tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"),] tmp <- tmp[tmp$Alive > 0,] tmp$Conc <- factor(tmp$Conc) tmp$Rep <- factor(paste0(tmp$Conc, "_", tmp$Rep)) ggplot(tmp, aes(y = ConsoBee, x = Time) )+ geom_line(aes(group = Rep), color = "gray50") + stat_smooth(aes(color = Conc), se = FALSE)+ # geom_point() + facet_wrap(~ Conc) + ylab("Syrup consumed per bee and per day\n(g/(bee*day))")+ xlab("Time (days) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid concentration\n(mg a.i./l syrup)", color="Boscalid concentration\n(mg a.i./l syrup)") + theme(legend.position = "top") + guides(color = guide_legend(byrow = TRUE)) #' #' \newpage #' #' Here we look at the same response (syrup consumption per bee) but against the number of #' living bees. The very high levels of syrup consumption (above 0.08 g/bee) occur only when #' there are only 1 or 2 living bee(s) left (and this case never occurs in the control). #' So this increase of consumption at the end of the test could be due to the stress of #' the only bee left. An important contributing factor could also be the evaporation that #' has been measured so be between 0.04 and 0.1 g per cage (in empty cages - #' see section \@ref(consumption)). #' When there are 10 bees in the cage, the maximum measurement error (over estimation) of #' syrup consumption **per bee** due to evaporation would then be #' between 0.004 and 0.01 g per bee but when there are only one or 2 bee(s) left this can have #' a stronger effect on the measurement. #' #+ fig.width = 16/2.54, fig.height = 10/2.54 # dev.new(width = 16/2.54, height = 10/2.54) # ggplot(tmp, aes(y = ConsoBee, x = Alive) )+ # geom_line() + geom_point(shape = 1) + facet_wrap(~ Conc) + stat_smooth(aes(color = Conc), se = FALSE)+ scale_x_continuous(breaks = seq(1, 10, 2))+ ylab("syrup consumed per bee and per day\n(g/(bee*day))")+ xlab("Number of bees alive") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid concentration\n(mg a.i./l syrup)", color="Boscalid concentration\n(mg a.i./l syrup)") + theme(legend.position = "top") + guides(color = guide_legend(byrow = TRUE)) #' #' ### Kinetics of the consumption with evaporation correction #' #' If we use the data corrected for evaporation, most of the peak consumption at the end #' of the test disappear. However there are still a few cases of peaking consumption at the #' doses 4500 mg/l and 18000 mg/l associated with low number of bees. #' So it seems that these observed peaks might be due to a combination of measurement error #' due to evaporation (and more important when there are few bees) and the mere fact that the #' bees are in low number. #' #' #+ fig.width = 16/2.54, fig.height = 10/2.54 # dev.new(width = 16/2.54, height = 10/2.54) tmp <- devap[devap$Treat %in% c("CANTUS", "CONTROL"),] tmp <- tmp[tmp$Alive > 0,] tmp$Conc <- factor(tmp$Conc) tmp$Rep <- factor(paste0(tmp$Conc, "_", tmp$Rep)) ggplot(tmp, aes(y = ConsoBee, x = Time) )+ geom_line(aes(group = Rep), color = "gray50") + stat_smooth(aes(color = Conc), se = FALSE)+ # geom_point() + facet_wrap(~ Conc) + ylab("Syrup consumed per bee and per day\n after a global correction for evaporation \n(g/(bee*day))")+ xlab("Time (days) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid concentration\n(mg a.i./l syrup)", color="Boscalid concentration\n(mg a.i./l syrup)") + theme(legend.position = "top") + guides(color = guide_legend(byrow = TRUE)) #+ fig.width = 16/2.54, fig.height = 9/2.54 # dev.new(width = 16/2.54, height = 9/2.54) # ggplot(tmp, aes(y = ConsoBee, x = Alive) )+ # geom_line() + geom_point(shape = 1) + facet_wrap(~ Conc) + stat_smooth(aes(color = Conc), se = FALSE)+ scale_x_continuous(breaks = seq(1, 10, 2))+ ylab("Syrup consumed per bee and per day\n after a global correction for evaporation \n(g/(bee*day))")+ xlab("Number of living bees") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid concentration\n(mg a.i./l syrup)", color="Boscalid concentration\n(mg a.i./l syrup)") + theme(legend.position = "top") + guides(color = guide_legend(byrow = TRUE)) #' #' \newpage #' #' ### Kinetics of the consumption up to 50% of mortality #' #' The question of the difference of consumption over time and between doses is mainly important #' up to 50% mortality (eg for application of the cumulative toxicity test proposed by EFSA, #' for ED50 calculations ...). #' We can then look at the consumption data after removing the data corresponding to more #' than 50% of mortality. #' NB : the results hereafter are almost identical with or without correction for evaporation. #' We show the data without evaporation correction. #' #' There is here a clear pattern of increased consumption at the beginning of the test followed #' by a decrease. The colored lines are loess smoothers (locally weighted polynomial regression). #' tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"),] tmp <- tmp[tmp$Alive > 0,] tmp$Conc <- factor(tmp$Conc) tmp$Rep <- factor(paste0(tmp$Conc, "_", tmp$Rep)) #+ fig.width = 16/2.54, fig.height = 10/2.54 # dev.new(width = 16/2.54, height = 10/2.54) ggplot(tmp[tmp$Alive >= 5,], aes(y = ConsoBee, x = Time) )+ geom_line(aes(group = Rep), color = "gray50") + stat_smooth(aes(color = Conc), se = FALSE)+ # geom_point() + facet_wrap(~ Conc) + ylab("Syrup consumed per bee and per day\n(g/(bee*day)")+ xlab("Time (days) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid concentration\n(mg a.i./l syrup)", color="Boscalid concentration\n(mg a.i./l syrup)") + theme(legend.position = "top") + guides(color = guide_legend(byrow = TRUE)) #' \newpage #' #' We can group the loess curves on the same graph. The kinetics of the consumption is clearly #' different with a peak at different moments and with a trend toward a higher consumption at lower doses #' #' #+ fig.width = 8.5/2.54, fig.height = 8/2.54 # dev.new(width = 8.5/2.54, height = 8/2.54) # tmp2 <- tmp # tmp2$Time <- scale(tmp2$Time, scale = FALSE) # ggplot(tmp[tmp$Alive >= 5,], aes(y = ConsoBee, x = Time) )+ # stat_smooth(aes(color = Conc), se = FALSE)+ # ylab("Syrup consumed per bee and per day\n(g/bee)")+ # xlab("Time (days)") + # scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + # labs(shape="Boscalid concentration (mg a.i./l syrup)", # color="Boscalid concentration (mg a.i./l syrup)") + # theme(legend.position = "top") + # guides(color = guide_legend(byrow = TRUE, title.position = "top")) #+ fig.width = 11/2.54, fig.height = 7/2.54 # dev.new(width = 11/2.54, height = 7/2.54) ggplot(tmp[tmp$Alive >= 5,], aes(y = ConsoBee, x = Time) )+ stat_smooth(aes(color = Conc), se = FALSE)+ ylab("Syrup consumed per bee and per day\n(g/(bee*day)")+ xlab("Time (days)") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration\n(mg a.i./l syrup)", color="Boscalid \nconcentration\n(mg a.i./l syrup)") #' This visual observation is confirmed by statistical testing. #' We used a gaussian mixed model with the consumption per bee per day as response and as fixed #' explanatory variable the Concentration, a second order polynomial of the time (to take #' into account the bell shaped relationship) and their first level interactions. The #' cage was used as random effect. The Time was centered before the analysis (mean day = 12.5) to #' limit problems of correlations between coefficients. We used only the consumption data #' when there were at least 5 living bees in the cage (i.e. 50% mortality). #+ echo = TRUE library(lme4) tmp2 <- tmp tmp2$Time <- scale(tmp2$Time, scale = FALSE) m <- lmer(ConsoBee ~ Conc + Time + I(Time^2) + Time:Conc + I(Time^2):Conc + (1 |Rep), data = tmp2[tmp2$Alive >=5,]) #' \newpage #' #' The model quality is correct as revealed by residual plots, eg a second order polynomial #' is sufficient to model the bell shaped relationship. #' A model without the second order polynomial shows clear problem of non linearity #' between the consumption and time (as can be expected based the above plots). #+ fig.width = 18/2.54, fig.height = 4.5/2.54 # dev.new(width = 18/2.54, height = 4.5/2.54) diagplot(m, layout = c(1,4), pt.cex = 0.5) diagplot2(m, layout = c(1,4), pt.cex = 0.5) #' **Likelihood ratio test (Type II)** #' The time x concentration interactions are significant meaning that the kinetics of the #' consumption (consumption vs time) is different between the concentrations. #' The fact that the main Conc (concentration) effect is highly significant means #' that even without considering these interactions #' (type II test) the consumption at the average day (~12.5 days) is different between #' the treatments (concentrations). #' #' NB a model without centering the time still shows clear differences of consumption at day 0 #' (ie the intercept for non centered data) at least between the control and the 2 highest #' concentrations. #' res <- Anova.lmer(m) res$"" <- signif(as.numeric(as.character(res[,3]))) pander(res) #' \newpage summary(m, corr = F) #' However the day to day variability is quite high as can be seen here when we compare #' the consumption just during the first 4 days of the experiment. #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) ggplot(tmp[tmp$Time <= 3,], aes(y = ConsoBee, x = Conc, color = Conc) )+ geom_point() + facet_wrap(~ Day, nrow = 1) + ylab("Syrup consumed per bee and per day\n(g/bee)")+ xlab("Boscalid concentration (mg a.i./l syrup)") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration\n(mg a.i./l syrup)", color="Boscalid \nconcentration\n(mg a.i./l syrup)") + theme(legend.position = "top", axis.text.x=element_text(angle=45, hjust=1))+ guides(color = guide_legend(byrow = TRUE)) #' #' ### Average consumption (without taking the time into account) #' #' In most bee tests we compare simply the average consumption between the treatment #' (without taking into account the time / kinetics). #' #' NB on the following graph the black dots represent the mean and the bars their standard #' deviation. #+ fig.width = 7.5/2.54, fig.height = 7/2.54 # dev.new(width = 7.5/2.54, height = 7/2.54) ggplot(tmp, aes(y = ConsoBee, x = Conc) )+ geom_point(aes(color = Conc), position = position_jitter(height = 0, width = 0.2), alpha = 0.2) + stat_summary(aes(group = Conc), fun.data = "mean_sdl", fun.args = list(mult = 1), geom = "pointrange", shape = 20, size = 0.5) + ylab("Syrup consumed per bee and per day\n(g/bee)")+ xlab("Boscalid concentration (mg a.i./l syrup) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration\n(mg a.i./l syrup)", color="Boscalid \nconcentration\n(mg a.i./l syrup)") + theme(legend.position = "top")+ guides(color = FALSE) # guides(color = guide_legend(byrow = TRUE)) #' We can test if there is a significant difference between the treatments with a gaussian #' mixed model with the concentration as fixed effect and the cage as random effect. #' The overall difference test is not significant : m <- lmer(ConsoBee ~ Conc + (1|Rep), data = tmp) pander(Anova(m, test = "F")) #' The absence of significant difference might be due to the somewhat extreme consumption #' values observed when the mortality is high. #' #' Here we use only the consumption data up to 10 days as in a classical chronic test. #' There is still no difference (NB choosing other time point does not change the results). #' This is due to the high variability including the variability due to time that is not #' taken into account here. #+ fig.width = 7.5/2.54, fig.height = 7/2.54 # dev.new(width = 7.5/2.54, height = 7/2.54) ggplot(tmp[tmp$Time <= 10,], aes(y = ConsoBee, x = Conc) )+ geom_point(aes(color = Conc), position = position_jitter(height = 0, width = 0.2), alpha = 0.2) + stat_summary(aes(group = Conc), fun.data = "mean_sdl", fun.args = list(mult = 1), geom = "pointrange", shape = 20, size = 0.5) + ylab("Syrup consumed per bee and per day\n up to day 10 (g/bee)")+ xlab("Boscalid concentration (mg a.i./l syrup) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration\n(mg a.i./l syrup)", color="Boscalid \nconcentration\n(mg a.i./l syrup)") + theme(legend.position = "top")+ guides(color = FALSE) # guides(color = guide_legend(byrow = TRUE)) m <- lmer(ConsoBee ~ Conc + (1|Rep), data = tmp[tmp$Time <= 10,]) pander(Anova(m, test = "F")) #' #' If you keep only the data up to 50% of mortality, then you can see a significant difference #' between the concentrations. A post-hoc all pairwise comparison test (`multcomp` package #' - single step method for p-values adjustment) shows that the consumption is lower in #' all doses relative to the control excepted at 4500 mg/l. #' The significant differences are displayed on the graph with letters (concentrations #' sharing the same letter are not significantly different). #' m <- lmer(ConsoBee ~ Conc + (1|Rep), data = tmp[tmp$Alive >= 5,]) pander(Anova(m, test = "F")) mc <- glht(m, linfct = mcp(Conc = "Tukey")) summary(mc) compact_letters <- cld(mc)$mcletters$Letters #+ fig.width = 7.5/2.54, fig.height = 7/2.54 # dev.new(width = 7.5/2.54, height = 7/2.54) y <- max(tmp[tmp$Alive >= 5,"ConsoBee"]) # position of the compact letters y <- y + 0.1*y ggplot(tmp[tmp$Alive >= 5,], aes(y = ConsoBee, x = Conc) )+ geom_point(aes(color = Conc), position = position_jitter(height = 0, width = 0.2), alpha = 0.2) + stat_summary(aes(group = Conc), fun.data = "mean_sdl", fun.args = list(mult = 1), geom = "pointrange", shape = 20, size = 0.5) + annotate("text", y = y, x = as.factor(names(compact_letters)), label = compact_letters) + ylab("Syrup consumed per bee and per day\nup to 50% of mortality (g/bee)")+ xlab("Boscalid concentration (mg a.i./l syrup) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration\n(mg a.i./l syrup)", color="Boscalid \nconcentration\n(mg a.i./l syrup)") + theme(legend.position = "top")+ guides(color = FALSE) # guides(color = guide_legend(byrow = TRUE)) #' \newpage #' #' ## Conclusions #' #' Note : the raw consumption data are provided in the file #' `raw_results/raw_data_with_computed_doses.csv` #' #' **Evaporation** #' #' * The evaporation is quite variable from day to day and may also be influenced by the presence #' living bees (here we measured the evaporation in empty cages or cages with dead bees for #' only a subset of days). #' * The effect of evaporation is mainly important in the evaluation of the daily consumption #' when there are only 1 or 2 living bees left. However the impact of the evaporation #' (as measured here, ie not very precisely) on the overall results (kinetics of the toxicity and #' kinetics of the consumption during most of the test) seem very minor. #' * Taking into account the evaporation might however improve the estimation of the dose to effect #' statistics and reduce the day to day measured variability in consumption. #' #' **Consumption** #' #' * The daily consumption seems to increase strongly when there are only 1 or 2 living bee left. #' This might be partially due to the measurement error induced by evaporation but even after #' correcting for evaporation we still observe some large peaks of consumption when there are only #' a few bees left #' * there is a huge day to day variation in consumption even in the control but some general #' patterns are nevertheless visible. #' * Even when there are more than 5 bees per cage, the consumption is not stable over time. #' In the control, the consumption slowly increases to reach a maximum between 15 and 20 #' days and then decreases slowly. #' * A similar pattern is observed at the different concentrations but with a different kinetics : #' the higher the concentration the sooner the maximum and the steeper the increase and then decrease. #' * The bees tend to have a lower overall consumption in the higher concentrations #' at a given point in time #' * These differences of consumption between treatments are only visible if you take #' into account the time and/or if you do not use the consumption data after 50% of mortality. #' Otherwise the variability of the consumption between days masks the differences. #' #' #' #' #' #' \newpage #' # /* # ---------------- Descriptive ------------------- # */ #' #' # Kinetics of the toxicity : LTx, LCx, LDDx, LCDx at different time or concentration #' #' ## Raw mortality rates #' #' Mortality rate over time for the different doses of boscalid/cantus, the control and toxic standard. #' The first vertical dashed line at day 10 shows when the test should have been stopped according #' to the standard protocol. #' The second vertical dashed line shows the time (day 20) where the mortality rate in the control reached #' 15% i.e. the validity criterion of the 10 day test. #' The test was stopped at 33 days when the mortality in the control reached 50%. #' #' The mortality at the highest concentration is higher than 50% only at Day 8. Before that date #' it is complicated to estimate any LC50 or LDD50. On day 10 the mortality has reached #' almost 100% at the highest concentration while all the other doses are below 50% and most #' of hem close to 0%. #' On Day 20 when the mortality in the control has reached 15%, the mortality in the 3 highest doses #' is 100% or close to 100%. And on day 31, the mortality of all concentrations has reached #' 100% while the control is still just below 50%. #' tmp <- d[d$Treat != "2CNA", ] tmp <- aggregate(tmp[c("Eff", "Dead")], tmp[c("Treat", "Conc", "Time")], sum) tmp$Treatment <- factor(paste(tmp$Treat, tmp$Conc), levels = c("CONTROL 0", "CANTUS 1125", "CANTUS 2250", "CANTUS 4500", "CANTUS 9000", "CANTUS 18000", "TS 1.5")) levels(tmp$Treatment) <- c("Control", "Boscalid 1125", "Boscalid 2250", "Boscalid 4500", "Boscalid 9000", "Boscalid 18000", "Toxic Standard") tmp$MortRate <- tmp$Dead/tmp$Eff #+ fig.width = 16/2.54, fig.height = 8/2.54 # dev.new(width = 16/2.54, height = 8/2.54) ggplot(tmp, aes(y = MortRate, x = Time, group = Treatment, color = Treatment, shape = Treatment)) + geom_line() + geom_point() + ylab("Mortality rate") + xlab("Time (days)") + geom_vline(xintercept = c(10, 20), linetype = 2, color = "gray50")+ scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(7)) + labs(color = "", shape = "") #' \newpage #' #' ## LTx vs Concentration #' #+ fig.width = 9/2.54, fig.height = 8/2.54 # dev.new(width = 9/2.54, height = 8/2.54) tmp <- LT[LT$Model == "Weibull2" & LT$PctMort %in% c(10, 50, 90) ,] tmp$PctMort <- factor(paste0("LT", tmp$PctMort)) ggplot(tmp, aes(y = Estimate, x = Conc, group = PctMort, color = PctMort, shape = PctMort)) + geom_line(aes(linetype = PctMort)) + geom_point() + geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) + ylab("LTx and Confidence Interval \n(days)")+ xlab("Boscalid concentration (mg a.i. / l syrup)") + scale_color_manual(values = c("dodgerblue", "gold", "orangered")) + scale_linetype_manual(values = c(2, 1, 3))+ labs(shape="", color="", linetype = "") + theme(legend.position = "top") #' #' ## LCx vs time #' #+ fig.width = 9/2.54, fig.height = 8/2.54 # dev.new(width = 9/2.54, height = 8/2.54) tmp <- LC[LC$Model == "Weibull2" & LC$PctMort %in% c(10, 50, 90) ,] tmp$PctMort <- factor(paste0("LC", tmp$PctMort)) ggplot(tmp, aes(y = Estimate, x = Time, group = PctMort, color = PctMort, shape = PctMort)) + geom_line(aes(linetype = PctMort)) + geom_point() + geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) + ylab("LCx and Confidence Interval \n(mg a.i./l syrup)")+ xlab("Time (days) ") + scale_color_manual(values = c("dodgerblue", "gold", "orangered")) + scale_linetype_manual(values = c(2, 1, 3))+ labs(shape="", color="", linetype = "") + theme(legend.position = "top") #' #' ## LDDx vs Time #' #+ fig.width = 9/2.54, fig.height = 8/2.54 # dev.new(width = 9/2.54, height = 8/2.54) tmp <- LDD[LDD$Model == "Weibull2" & LDD$PctMort %in% c(10, 50, 90) ,] tmp$PctMort <- factor(paste0("LDD", tmp$PctMort)) ggplot(tmp, aes(y = Estimate, x = Time, group = PctMort, color = PctMort, shape = PctMort)) + geom_line(aes(linetype = PctMort)) + geom_point() + geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) + ylab("LDDx and Confidence Interval \n(average mg a.i./(bee*day))")+ xlab("Time (days) ") + scale_color_manual(values = c("dodgerblue", "gold", "orangered")) + scale_linetype_manual(values = c(2, 1, 3))+ labs(shape="", color="", linetype = "") + theme(legend.position = "top") #' #' ## LCDx vs Time #' #+ fig.width = 9/2.54, fig.height = 8/2.54 # dev.new(width = 9/2.54, height = 8/2.54) tmp <- LCD[LCD$Model == "Weibull2" & LCD$PctMort %in% c(10, 50, 90) ,] tmp$PctMort <- factor(paste0("LCD", tmp$PctMort)) ggplot(tmp, aes(y = Estimate, x = Time, group = PctMort, color = PctMort, shape = PctMort)) + geom_line(aes(linetype = PctMort)) + geom_point() + geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) + ylab("LCDx and Confidence Interval \n(total mg a.i./bee)")+ xlab("Time (days) ") + scale_color_manual(values = c("dodgerblue", "gold", "orangered")) + scale_linetype_manual(values = c(2, 1, 3))+ labs(shape="", color="", linetype = "") + theme(legend.position = "top") #' #' \newpage #' #' #' ## Conclusion #' #' The Lethal Times are shorter when the concentration is higher and the Lethal Concentrations #' and Lethal Daily Doses are decreasing when the exposure time is increasing. #' This is not very surprising. The question here is to know if this decrease is as expected #' according to Haber's rule (ie when there is no cumulative toxicity). This question #' is explored in section \@ref(loglog). However it seems to be not very pratical that these #' standard toxicity estimates depend so much on the exposure time... #' #' The Lethal Cumulated Dose show however a different and interesting pattern : #' The LCDx is more or less stable up to ~ day 17 and then suddenly starts to decrease. #' This is a first clue that there might be some cumulative toxicity or change in the capacities #' of the bees to detoxify the product. #' The fact that the LCDx is stable during the first part of the tes seems however to #' be an interesting feature #' because it means that you can have an estimate of toxicity that is relatively independent #' of the test duration providing that there is no cumulative toxicity #' (you need however to maintain the test long enough to reach 50% #' of mortality in some of the doses). #' #' #' #' \newpage #' # /* # ---------------- cumulative toxicity ------------------- # */ #' #' # Testing for cumulative toxicity (= "time reinforced toxicity") #' #' We have explored three main ways to check for cumulative toxicity properties of boscalid #' as described in the literature : #' #' 1. Log-Log relationships between concentration and time. If there is no cumulative toxicity, the slope #' of this linear relationship is expected to be = -1. There are 3 main possibilities : #' #' + 1a log(Concentration) vs log(LTx) #' + 1b log(LCx) vs log(Time) #' + 1c log(LDx) vs log(Time) #' #' 2. Comparing cumulative doses between concentrations (EFSA protocol - EFSA 2013). If there is no #' cumulative toxicity, you expect to observe the same level of cumulative dose to reach a #' given level of mortality whatever the concentration. #' 3. Estimation of $\alpha$ and $\beta$ : two exponents describing the kinetics between the #' concentration and the time (method proposed by Miller 2000). Their ratio is expected #' to be = 1 when there is no cumulative toxicity #' #' #' #' #' ## Log-log relationship between concentration and time {#loglog} #' #' According to Haber's rule, if there is no cumulative toxicity, when the concentration/dose #' is divided by 2 the time of exposure to reach the same level of mortality should be doubled. #' #' There are two main traditional ways to see if the concentration ~ time relationship #' follows Haber's rule : #' #' 1. fit a model of log(Concentration) vs log(LTx) #' 2. fit a model of log(LCx) vs log(Time) or log(LDDx) vs log(Time) #' #' If the toxicity follows Haber's rule, the slope of these models should be ~= -1. #' If there is cumulative toxicity, the slope should be lower than -1. #' #' NB : strictly speaking it would be more correct to fit a model of #' log(LTx) vs log(Concentration) because the LTx is the value estimated with a certain #' error and this error should be on the y axis. However then the interpretation #' of the coefficient is reversed : the slope should be higher than -1 when there is #' cumulative toxicity. This would make the comparison and interpretation less straightforward and this #' is why we keep the time and LT as X axis for all approaches. #' #' We also estimated a log(LCDx) vs log(Time) model for comparison however this should not #' be used in the classical way to test if the Haber's rule holds because LCD already #' includes time accumulation (it is a cumulative dose). #' With this model, you expect a slope = 0 under Haber's rule. #' On the contrary, if there is cumulative toxicity, you expect a slope significantly lower than #' 0 (and not -1). #' #' We have estimated the slope for each type of model and each level of expected mortality #' (by steps of 10%). The following graph shows the slope estimates (black dots) and their #' confidence interval. The vertical dashed line shows the value expected under Haber's rule. #' (-1). #' The confidence intervals are larger (lower precision) for the log(Concentration) #' vs log(LTx) models because they are based on a lower number of points : 5 points (one for #' each dose). The control dose has been excluded. #' The slope is always significantly different from -1 excepted for LT90 #' (confidence interval : -3.378 -0.995) however the upper bound is very close to -1 and #' the CI is the largest of all estimates. #' #' As expected the slopes of the log(LCDx) vs log(Time) models are not significantly #' less than -1 but they are significantly lower than 0. #' # Estimate of the log-log slope and the intercept under Haber'srule for each mortality rate # and based on LC, LCD and LT estimates ConcvsLT <- matrix(ncol = 7, nrow = length(PctMort)) colnames(ConcvsLT) <- c("PctMort", "Intercept", "Slope", "CIlower", "CIupper","Haber_int", "Haber_slope") for (i in 1:length(PctMort)) { tmp <- LT[LT$Model == "Weibull2" & LC$PctMort == PctMort[i], ] haber <- lm(log(Conc) ~ 1 + offset(-1* log(Estimate)), data = tmp[-1,]) m <- lm(log(Conc) ~ log(Estimate), data = tmp[-1,]) ConcvsLT[i,] <- c(PctMort[i], coef(m), confint(m)[2,], coef(haber), -1) } LCvsTime <- matrix(ncol = 7, nrow = length(PctMort)) colnames(LCvsTime) <- c("PctMort", "Intercept", "Slope", "CIlower", "CIupper","Haber_int", "Haber_slope") for (i in 1:length(PctMort)) { tmp <- LC[LC$Model == "Weibull2" & LC$PctMort == PctMort[i], ] haber <- lm(log(Estimate) ~ 1 + offset(-1* log(Time)), data = tmp) m <- lm(log(Estimate) ~ log(Time), data = tmp) LCvsTime[i,] <- c(PctMort[i], coef(m), confint(m)[2,], coef(haber), -1) } LDDvsTime <- matrix(ncol = 7, nrow = length(PctMort)) colnames(LDDvsTime) <- c("PctMort", "Intercept", "Slope", "CIlower", "CIupper","Haber_int", "Haber_slope") for (i in 1:length(PctMort)) { tmp <- LDD[LDD$Model == "Weibull2" & LDD$PctMort == PctMort[i], ] haber <- lm(log(Estimate) ~ 1 + offset(-1* log(Time)), data = tmp) m <- lm(log(Estimate) ~ log(Time), data = tmp) LDDvsTime[i,] <- c(PctMort[i], coef(m), confint(m)[2,], coef(haber), -1) } LCDvsTime <- matrix(ncol = 7, nrow = length(PctMort)) colnames(LCDvsTime) <- c("PctMort", "Intercept", "Slope", "CIlower", "CIupper","Haber_int", "Haber_slope") for (i in 1:length(PctMort)) { tmp <- LCD[LCD$Model == "Weibull2" & LCD$PctMort == PctMort[i], ] haber <- lm(log(Estimate) ~ 1 + offset(-1* log(Time)), data = tmp) m <- lm(log(Estimate) ~ log(Time), data = tmp) LCDvsTime[i,] <- c(PctMort[i], coef(m), confint(m)[2,], coef(haber), -1) } slope <- rbind(ConcvsLT, LCvsTime, LDDvsTime, LCDvsTime) slope <- data.frame( Method = rep(c("Conc vs LTx", "LCx vs Time", "LDDx vs Time", "LCDx vs Time"), each = length(PctMort)), EDx = paste0(rep(c("LT", "LC", "LDD", "LCD"), each = length(PctMort)), slope[,"PctMort"]), slope) slope$Method <- factor(slope$Method, levels = c("Conc vs LTx", "LCx vs Time", "LDDx vs Time", "LCDx vs Time")) #+ fig.width = 16/2.54, fig.height = 8/2.54 # dev.new(width = 16/2.54, height = 8/2.54) ggplot(slope, aes(y = EDx, x = Slope)) + geom_point() + geom_errorbarh(aes(xmin = CIlower, xmax = CIupper), height = 0.2)+ facet_wrap(~Method, scale = "free_y", nrow = 1) + geom_vline(xintercept = -1, linetype = "22", color = "gray60") + ylab("") #' In conclusion : there is little doubt based on these results that the concentration vs time #' relationship does not follow Haber's rule and this implies that there is some level of #' cumulative toxicity. #' #' \newpage #' #' However, even if all the slopes are significantly smaller than -1, #' these slopes are just a rough summary of the results and again an interesting pattern #' appears when we plot for example the log(LDDx) vs log (Time) and compare to the line #' expected under Haber's rule. #' NB : the Haber's rule line is a linear regression with a fixed -1 slope #' (the model estimates only the intercept). #' #' #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) tmp <- LDD[LDD$Model == "Weibull2" & LDD$PctMort == c(10,50,90), ] tmp$LDD <- factor(paste0("LDD", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LDD, scales = "free") + stat_smooth(aes(linetype = "LDDx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LDD50) (average mg a.i./(bee*day))") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top", legend.margin = margin(0, 0, 0, 0))+ guides(linetype=guide_legend(ncol=1)) #' The log-log LDDx vs Time relationship seems to follow closely Haber's rule up to day #' 17-18 with a slope of ~ -1 (i.e. when the Dietary Dose is divided by 2, the time to reach the same #' mortality is doubled). At day 17-18 the slope abruptly decreases and the relationship #' clearly deviates from Haber's rule. This pattern is more marked for lower levels of mortality #' (LDD10, LDD20, etc) and tend to disappear at higher levels of mortality (LDD90). #' For the LDD90 the points are almost perfectly aligned on a straight line without #' inflection point. However even for the LDD90 where the observed regression line seems #' to be quite close to the theoretic Haber's slope, the slope is significantly lower #' than the expected -1 (estimate = -1.359 with a 95% confidence interval of #' [-1.492,-1.226]). #' #' This pattern is not visible on the LTx vs Concentration Log-Log graphs and is less marked on #' the LCx vs Time graphs. However the LTx vs Concenration regression are based on only 5 #' points (one for each concentration) and it might be therefore difficult to visualise #' the sudden change visible on the LDDx plot. #' #' All the graphs are provided below #' #' \newpage #' #' ### log(Concentraton) vs log(LTx) #' #+ fig.width = 14/2.54, fig.height = 12/2.54 # dev.new(width = 14/2.54, height = 12/2.54) tmp <- LT[LT$Model == "Weibull2" & LT$Conc > 0, ] tmp$LT <- factor(paste0("LT", tmp$PctMort)) ggplot(tmp, aes(x = log(Estimate), y = log(Conc))) + geom_point(shape = 1) + facet_wrap(~ LT) + stat_smooth(aes(linetype = "Concentr. vs LTx (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(Concentration) (mg a.i./l)") + xlab("log(LTx) (days)") + scale_linetype_manual(name = "", values = c("solid","22")) + theme(legend.position = "top") #' #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) tmp <- LT[LT$Model == "Weibull2" & LT$Conc > 0 & LT$PctMort == c(10,50,90), ] tmp$LT <- factor(paste0("LT", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Conc))) + geom_point(shape = 1) + facet_wrap(~ LT, scales = "free") + stat_smooth(aes(linetype = "Concentration vs LTx (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(Concentration) (mg a.i./l)") + xlab("log(LTx) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top", legend.margin = margin(0, 0, 0, 0))+ guides(linetype=guide_legend(ncol=1)) #' #' \newpage #' #' ### log(LCx) vs log(Time) #' #+ fig.width = 14/2.54, fig.height = 12/2.54 # dev.new(width = 14/2.54, height = 12/2.54) tmp <- LC[LC$Model == "Weibull2", ] tmp$LC <- factor(paste0("LC", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LC) + stat_smooth(aes(linetype = "LCx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LCx) (mg a.i./l)") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top") #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) tmp <- LC[LC$Model == "Weibull2" & LC$PctMort == c(10,50,90), ] tmp$LC <- factor(paste0("LC", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LC, scales = "free") + stat_smooth(aes(linetype = "LCx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LCx) (mg a.i./l)") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top", legend.margin = margin(0, 0, 0, 0))+ guides(linetype=guide_legend(ncol=1)) #' \newpage #' #' ### log(LDDx) vs log(Time) #' #' #+ fig.width = 14/2.54, fig.height = 12/2.54 # dev.new(width = 14/2.54, height = 12/2.54) tmp <- LDD[LDD$Model == "Weibull2", ] tmp$LDD <- factor(paste0("LDD", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LDD) + stat_smooth(aes(linetype = "LDDx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LDDx) (average mg a.i./(bee*day))") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top") #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) tmp <- LDD[LDD$Model == "Weibull2" & LDD$PctMort == c(10,50,90), ] tmp$LDD <- factor(paste0("LDD", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LDD, scales = "free") + stat_smooth(aes(linetype = "LDDx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LDDx) (average mg a.i./(bee*day))") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top", legend.margin = margin(0, 0, 0, 0))+ guides(linetype=guide_legend(ncol=1)) #' #' \newpage #' #' ### log(LCDx) vs log(Time) #' #' #+ fig.width = 14/2.54, fig.height = 12/2.54 # dev.new(width = 14/2.54, height = 12/2.54) tmp <- LCD[LCD$Model == "Weibull2", ] tmp$LCD <- factor(paste0("LCD", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LCD) + stat_smooth(aes(linetype = "LCDx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LCDx) (mg a.i./bee)") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top") #' Only with LCD10, LCD50 and LCD90 with free y scales #' #+ fig.width = 18/2.54, fig.height = 8/2.54 # dev.new(width = 18/2.54, height = 8/2.54) tmp <- LCD[LCD$Model == "Weibull2" & LCD$PctMort == c(10,50,90), ] tmp$LCD <- factor(paste0("LCD", tmp$PctMort)) ggplot(tmp, aes(y = log(Estimate), x = log(Time))) + geom_point(shape = 1) + facet_wrap(~ LCD, scales = "free") + stat_smooth(aes(linetype = "LCDx vs Time (estimated b)"), show.legend = TRUE, method = "lm", formula = y ~ x, se = FALSE, color = "gray30", size = 0.75) + stat_smooth(aes(linetype = "Haber's rule (b = -1)"), show.legend = TRUE, method = "lm", formula = y ~ 1 + offset(-1*x), se = FALSE, color = "gray60", size = 0.75)+ ylab("log(LCDx) (total mg a.i./bee)") + xlab("log(Time) (days)") + scale_linetype_manual(name = "", values = c("22","solid")) + theme(legend.position = "top", legend.margin = margin(0, 0, 0, 0))+ guides(linetype=guide_legend(ncol=1)) #' \newpage #' Raw results of the model (slope and confidence intervals) : pander(slope[, c(1,2,5:7)]) #' \newpage #' #' ## Comparing cumulative dose between concentrations (EFSA protocol) #' #' To test for cumulative toxicity the EFSA (2013) propose to : #' #' * determine the LC50 at 48h #' * launch a test with bees fed at a high dose = LC50 and others fed at a low dose = 0.25*LC50 #' * measure the total active substance consumed when 50% mortality occured for both the #' high dose and the low dose #' * compare these values with a t test (with a power of 80% to detect a difference of 35%) #' * You can express the difference as a % of the high dose consumption : #' 100* (high dose consumption - low dose consumption)/high dose consumption #' * if there is cumulative toxicity, you expect that the low dose total consumption to reach #' 50% of mortality will be lower than the total consumption of the higher dose #' #' Here the LC50 at 48h is impossible to estimate. However we can apply the same idea by #' comparing the cumulated dose consumed when the mortality reaches 50%. #' #' We compare between the concentration the cumulative dose consumed per bee once the cage #' reached 50% mortality. We use a simple analysis of variance folowed by all pairwise #' post-hoc comparisons (similar to a Tukey test here using `multcomp` package single-step #' p-value correction method). #' Here the LC50 at 10 days is ~ 10000 mg ai/l (close to our 9000 dose). However #' if Haber's rule holds and that you want to test a quarter of this dose you expect to have #' to wait ~ 40 days to reach the same mortality... #' #' If there is no cumulative toxicity, no differences in total doses between #' the treatments (concentrations) are expected. #' #' ### With raw mortalities #' #' The global differences are highly significant. If you compare the 3 highest concentrations #' they are all significantly different from concentrations four times lower. #' NB : a student test would also be significant as the student t test is more powerfull than #' this corrected p-value test. #' So applying the EFSA protocol could have worked here but the toxicity was undetectable #' at 2 days. #' #' NB : these are non corrected mortalities. This could impact the results because #' the honey bees will finally die even with very low doses of the product. #' #' Another potential problem is that some times (quite often in fact) you do not observe #' an exact mortality of 50%. For example you have 4 dead bees one day and then you have #' 6 or 7 dead bees the next day. You could expect that the cumulated dose from the cages #' with more than 5 dead bees might be higher but it does not seem to be the case here #' (see on the graph below) #' #' \newpage #' #' Anova table : #' tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] tmp$Mortality <- factor(ifelse(tmp$MortRate >=0.5, "More than 50%", "Less than 50%")) tmp <- tmp[tmp$MortRate >=0.5,] tmp <- split(tmp, f = list(tmp$Treat, tmp$Rep, tmp$Conc)) for (i in 1: length(tmp)) { tmp[[i]] <- tmp[[i]][1,] } tmp <- do.call(rbind, tmp) tmp <- tmp[!is.na(tmp$Treat),] tmp$Conc_f <- factor(tmp$Conc) m <- lm(CumDoseBee ~ Conc_f, data = tmp[tmp$Treat != "CONTROL",]) pander(anova(m)) # summary(m) library(multcomp) mc <- glht(m, linfct = mcp(Conc_f = "Tukey")) # summary(mc) compact_letters <- cld(mc)$mcletters$Letters #' Graph : concentrations with different letters are significantly different #' (post-hoc Tukkey like test). Some horizontal noise has been added to the points position #' to avoid overplotting. #' #+ fig.width = 10/2.54, fig.height = 6/2.54 # dev.new(width = 10/2.54, height = 6/2.54) set.seed(123) y <- max(tmp[tmp$Treat != "CONTROL","CumDoseBee"]) y <- y + 0.1*y # ggplot(tmp[tmp$Treat != "CONTROL",], aes(y = CumDoseBee, x = Conc, fill = factor(Dead))) + # geom_point(shape = 21, size = 3, position = position_jitter(width = 600, height = 0)) + # scale_fill_manual(values = c("gray95", "gray50", "gray10"), # name = "Number of \ndead bees") + # annotate("text", y = y, x = as.numeric(names(compact_letters)), # label = compact_letters)+ # xlab("Concentration (mg a.i./l)") + # ylab("Cumulative dose consumed \nper bee when mortality \nreaches 50% (total mg a.i./bee)") # tmp$Conc <- as.factor(tmp$Conc) set.seed(1) ggplot(tmp[tmp$Treat != "CONTROL",], aes(y = CumDoseBee, x = Conc, fill = factor(Dead))) + geom_point(shape = 21, size = 3, position = position_jitter(width = 0.15, height = 0)) + scale_fill_manual(values = c("gray95", "gray50", "gray10"), name = "Number of \ndead bees") + annotate("text", y = y, x = 1:5, label = compact_letters)+ xlab("Concentration (mg a.i./l)") + ylab("Cumulative dose consumed \nper bee when mortality \nreaches 50% (total mg a.i./bee)") #' Raw data for the record : #' row.names(tmp) <- NULL pander(tmp[,c(1:5, 6, 19)]) #' #' \newpage #' #' ### With corrected mortalities #' #' Applying a mortality correction adequately is not straightforward. #' Ideally we should group the data from the 3 replicates and have identical initial number #' of bees. However as here we use only the mortality rate to select when we should stop looking #' at cumulated consumption, one of the forms of Abbot's correction seems appropriate : #' `CorMortRate = (MortRate - MortRateControl)/(1-MortRateControl)` where `MortRate` is the #' mortalitay rate on a given day in a given replicate and `MortRateControl` is the #' global mortality rate of the 3 control replicates pooled on the same day. #' #' The results are quite similar excepted that the difference of total consumption #' between concentration 18000 and 4500 is no more significant. #' Note however that this difference would be borderline significant with a single #' student t test (without correction of the p-values - as recommended by the EFSA protocol) #' with p-value = 0.057. #' #' Anova table : #' tmp <- d[d$Treat %in% c("CONTROL"),] MortRateControl <- aggregate(tmp["MortRate"], tmp["Day"], mean ) colnames(MortRateControl)[2] <- "MortRateControl" tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] tmp <- merge(tmp, MortRateControl, all.x = TRUE, sort = FALSE) tmp$CorMortRate <- (tmp$MortRate - tmp$MortRateControl)/(1-tmp$MortRateControl) tmp <- tmp[tmp$CorMortRate >=0.5,] tmp <- split(tmp, f = list(tmp$Treat, tmp$Rep, tmp$Conc)) for (i in 1: length(tmp)) { tmp[[i]] <- tmp[[i]][1,] } tmp <- do.call(rbind, tmp) tmp <- tmp[!is.na(tmp$Treat),] tmp$Conc_f <- factor(tmp$Conc) m <- lm(CumDoseBee ~ Conc_f, data = tmp[tmp$Treat != "CONTROL",]) pander(anova(m)) # summary(m) library(multcomp) mc <- glht(m, linfct = mcp(Conc_f = "Tukey")) # summary(mc, test = adjusted("none")) # without p-value correction # summary(mc) compact_letters <- cld(mc)$mcletters$Letters #' Graph : concentrations with different letters are significantly different #' (post-hoc Tukey like test). Some horizontal noise has been added to the points position #' to avoid overplotting. #' #+ fig.width = 10/2.54, fig.height = 6/2.54 # dev.new(width = 10/2.54, height = 6/2.54) set.seed(123) y <- max(tmp[tmp$Treat != "CONTROL","CumDoseBee"]) y <- y + 0.1*y # ggplot(tmp[tmp$Treat != "CONTROL",], aes(y = CumDoseBee, x = Conc, fill = CorMortRate)) + # geom_point(shape = 21, size = 3, position = position_jitter(width = 600, height = 0)) + # scale_fill_gradient(low = "gray90", high = "gray10", # name = "Corrected\nMortality") + # annotate("text", y = y, x = as.numeric(names(compact_letters)), # label = compact_letters)+ # xlab("Concentration (mg a.i./l)") + # ylab("Cumulative dose consumed \nper bee when corrected mortality \nreaches 50% (total mg a.i./bee)") # tmp$Conc <- as.factor(tmp$Conc) set.seed(1) ggplot(tmp[tmp$Treat != "CONTROL",], aes(y = CumDoseBee, x = Conc, fill = CorMortRate)) + geom_point(shape = 21, size = 3, position = position_jitter(width = 0.15, height = 0)) + scale_fill_gradient(low = "gray90", high = "gray10", name = "Corrected\nMortality") + annotate("text", y = y, x = 1:5, label = compact_letters)+ xlab("Concentration (mg a.i./l)") + ylab("Cumulative dose consumed \nper bee when corrected mortality \nreaches 50% (total mg a.i./bee)") #' #' \newpage #' #' Raw data for the record : #' row.names(tmp) <- NULL pander(tmp[,c(2:4, 1, 5, 7, 19, 22, 23)]) #' Mean cumulative dose for each treatment : tmp <- tmp[,c(2:4, 1, 5, 7, 19, 22, 23)] pander(aggregate(tmp["CumDoseBee"], tmp["Conc"], mean)) #' #' #' #' Another way to look at these data is to plot the cumulative doses consumed by the bees #' and visualize when you reach e.g. 50% mortality. The black dots show the LT50 for each concentration. #' If there were no cumulative toxicity, one would expect that the LT50 will be reached #' at similar cumulative doses (or at similar time as the control). #' It is not the case here : we observe that the cumulative dose leading to 50% of mortality #' at the lowest concentration (1125 mg a.i./ml) is less than half the #' cumulative dose killing 50% of the bees at 4500, 9000 and 18000 (mg a.i/ml). #' #' #+ fig.width = 9/2.54, fig.height = 8/2.54 # dev.new(width = 9/2.54, height = 8/2.54) tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"),] tmp <- aggregate(tmp["CumDoseBee"], tmp[c("Conc", "Time")], mean) tmp$Conc <- factor(tmp$Conc) tmpLT <- LT[LT$PctMort == 50 & LT$Model == "Weibull2",] tmpLT$Time <- round(tmpLT$Estimate, 0) tmp <- merge(tmp, tmpLT, all.x = TRUE) ggplot(tmp, aes(y = CumDoseBee, x = Time, group = Conc, color = Conc, shape = Conc)) + geom_line() + geom_point() + geom_point(aes(x = Estimate), shape = 20, color = "black", size = 3) + geom_errorbarh(aes(xmin = Lower, xmax = Upper), color = "black") + ylab("Cumulative dose consumed per bee\naverage of the 3 cages \n(total mg a.i./bee)")+ xlab("Time (days) ") + scale_color_manual(values = colorRampPalette(c("dodgerblue", "gold", "orangered"))(6)) + labs(shape="Boscalid \nconcentration", color="Boscalid \nconcentration") + theme(legend.position = "top") + guides(shape = guide_legend(byrow = TRUE), color = guide_legend(byrow = TRUE)) #' \newpage #' #' ## Estimation of $\alpha$ and $\beta$ #' #' #' #' The idea here is to use a more general form of the C vs t relationship : #' #' $$C^\alpha t^\beta=k$$ #' #' That relationship can be simplified into #' #' $$C = k^{1/\alpha} t^{-\beta/\alpha} = k' t^\gamma$$ #' #' where $\gamma$ is equivalent to the $b$ of the classical equation and should be = 1 #' when Haber's rule holds. #' According to Miller et al. (2000) we could estimate $\alpha$ and $\beta$ with a #' probit model : $Y = m + \alpha log(C) + \beta log(t)$ #' where $Y$ are the observed mortalities. #' #' We used a binomial GLM with a probit link #' (estimated by quasilikelihood to take into account overdispersion). Because of the presence #' of 0 values we used a log(x+1) transformation for both Concentration and Time. #' #' We tested other link functions (logit, cloglog which is equivalent to a weibull, cauchit) #' but the probit was the one that provided the best fit (lowest QAICc). #' #+ echo = TRUE m <- glm(MortRate ~ log(Conc+1) + log(Time+1), weights = Eff, family = quasibinomial(link = "probit"), data = d[d$Treat %in% c("CONTROL", "CANTUS"),]) summary(m) # confint(m) # diagplot(m) # diagplot2(m) # coef(m)[3] / coef(m)[2] # b = beta / alpha #' #' The regression coefficient for the concentration ($\alpha$) is `r round(coef(m)[2],3)` #' and the regression coefficient for time ($\beta$) is `r round(coef(m)[3],3)` and $\gamma = b = \beta/\alpha =$ #' `r round(coef(m)[3] / coef(m)[2],3)`. This would indicate a high level of cumulative toxicity, #' however this value is very different from the value calculated #' with the more classical approach above (~ 1.5 - 2). #' #' However this model does not fit the real data very well. #' It is quite clear that this model tends to underestimate toxicity at higher concentrations #' and to overestimate toxicity at low concentrations. The black line should be #' less steep. #' #' In the graph below the two predictive variables, Concentration and time, are represented on #' y and x axes and the color represent the predicted values for the response (mortality). #' The line represent the 50% mortality predicted values. The points #' (for the 3 replicates) show wen we observed 50% mortality in each cage. #' #' # prediction of mortalities according to this model X <- expand.grid(seq(min(d$Conc), max(d$Conc), by = 100), seq(min(d$Time), max(d$Time), by = 0.1)) pred <- family(m)$linkinv(as.matrix(cbind(1, log(X+1))) %*% coef(m)) pred <- data.frame(Conc = X[,1], Time = X[,2], MortRate = pred) mypal <- colorRampPalette(colors = adjustcolor(c("steelblue", "grey95", "orangered"), 0.5)) tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] tmp$Mortality <- factor(ifelse(tmp$MortRate >=0.5, "More than 50%", "Less than 50%")) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(pred, aes(y=Time, x=Conc, fill=MortRate)) + geom_raster() + geom_point(data = tmp, aes(color = Mortality, group = Rep), shape = 20, position = position_dodge(width = 700))+ geom_line(data = pred[round(pred$MortRate,3) == 0.5, ]) + # scale_fill_gradientn(colours = mypal(100)) + scale_fill_gradient2(low = "steelblue", mid = "grey95", high = "red", midpoint = 0.5, name = "Predicted \nmortality rate") + scale_color_manual(name = "Observed mortality", values = c("gray90", "gray40")) + ylim(0,35) + ylab("Time (days)") + xlab("Concentration (mg a.i./l)") + coord_flip() #' We can apply the same method on the average daily dose instead of the concentration. #' The fit is better (lower QAICc than the previous model) #+ echo = TRUE m <- glm(MortRate ~ log(MeanDoseBee+1) + log(Time+1), weights = Eff, family = quasibinomial(link = "probit"), data = d[d$Treat %in% c("CONTROL", "CANTUS"),]) summary(m) #' #' * Concentration coefficient : $\alpha$ = `r round(coef(m)[2],3)` #' * Time coefficient : $\beta$ = `r round(coef(m)[3],3)` #' * $\gamma = b = \beta/\alpha =$ `r round(coef(m)[3] / coef(m)[2],3)` #' #' This is now much lower than the b value estimated before but also much lower than b=1 #' expected under Haber's rule... #' #' Compared to the real data this model seems to do a better job. #' NB in the following graph the points might be slightly moved vertically to limit overplotting. #' # prediction of mortalities according to this model tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] X <- expand.grid(seq(min(tmp$MeanDoseBee), max(tmp$MeanDoseBee), by = 0.01), seq(min(tmp$Time), max(tmp$Time), by = 0.1)) pred <- family(m)$linkinv(as.matrix(cbind(1, log(X+1))) %*% coef(m)) pred <- data.frame(MeanDoseBee = X[,1], Time = X[,2], MortRate = pred) mypal <- colorRampPalette(colors = adjustcolor(c("steelblue", "grey95", "orangered"), 0.5)) tmp$Mortality <- factor(ifelse(tmp$MortRate >=0.5, "More than 50%", "Less than 50%")) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(pred, aes(y=Time, x=MeanDoseBee, fill=MortRate)) + geom_raster() + geom_point(data = tmp, aes(color = Mortality, group = Rep), shape = 20, position = position_dodge(width = 0.01))+ geom_line(data = pred[round(pred$MortRate,2) == 0.5, ]) + # scale_fill_gradientn(colours = mypal(100)) + scale_fill_gradient2(low = "steelblue", mid = "grey95", high = "red", midpoint = 0.5, name = "Predicted \nmortality rate") + scale_color_manual(name = "Observed mortality", values = c("gray90", "gray40")) + ylim(0,35) + ylab("Time (days)") + xlab("Average daily dose (mg a.i./ (bee*day))") + coord_flip() #' Interestingly if instead of using the observed mortalities, we use the LCx data #' with the intended mortality rate x as the response and Time and the LC estimate #' instead of the concentration, the beta calculated is more similar to what could be expected. #' #+ echo = TRUE m <- glm(PctMort/100 ~ log(Estimate+1) + log(Time+1), family = quasibinomial(link = "probit"), data = LC[LC$Model == "Weibull2",]) summary(m) #' #' * Concentration coefficient : $\alpha$ = `r round(coef(m)[2],3)` #' * Time coefficient : $\beta$ = `r round(coef(m)[3],3)` #' * $\gamma = b = \beta/\alpha =$ `r round(coef(m)[3] / coef(m)[2],3)` #' # prediction of mortalities according to this model LCtemp <- LC[LC$Model == "Weibull2",] X <- expand.grid(seq(min(LCtemp$Estimate), max(LCtemp$Estimate), by = 100), seq(0, 32, by = 0.1)) pred <- family(m)$linkinv(as.matrix(cbind(1, log(X+1))) %*% coef(m)) pred <- data.frame(Conc = X[,1], Time = X[,2], MortRate = pred) mypal <- colorRampPalette(colors = adjustcolor(c("steelblue", "grey95", "orangered"), 0.5)) tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] tmp$Mortality <- factor(ifelse(tmp$MortRate >=0.5, "More than 50%", "Less than 50%")) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(pred, aes(y=Time, x=Conc, fill=MortRate)) + geom_raster() + geom_point(data = tmp, aes(color = Mortality, group = Rep), shape = 20, position = position_dodge(width = 500))+ geom_line(data = pred[round(pred$MortRate,2) == 0.5, ]) + # scale_fill_gradientn(colours = mypal(100)) + scale_fill_gradient2(low = "steelblue", mid = "grey95", high = "red", midpoint = 0.5, name = "Predicted \nmortality rate") + scale_color_manual(name = "Observed mortality", values = c("gray90", "gray40")) + ylim(0,32) + ylab("Time (days)") + xlab("Concentration (mg a.i./l)") + coord_flip() #' #' **Conclusion** #' The results obtained with this method are puzzling. Using the concentration or the #' dose gives very different results and both results seem to be unlikely. #' #' There is also something strange from the theoretical point of view. #' If $$C^\alpha t^\beta=k$$ then $$\alpha log(C) + \beta log(t) = log(k)$$ (where k #' is a given level of mortality) but here we model instead #' $$\alpha log(C) + \beta log(t) = k$$. This would be OK if used a log link but here #' we use a probit (after the recommendation of Miller et al. 2000) #' #' #' \newpage #' # /* # ---------------- Models comparisons ------------------- # */ #' #' # Comparison of the 4 types of models {#ModelChoice} #' #' NB : the aim of this rather long section is to explain why we choose a Weibull 2 model #' in all the analyses shown above. #' In most of the studies the type of model is chosen _a priori_ and without justification. #' #' We compare 4 types of dose response curves (logistic, loglogistic, #' weibull 1 and weibull 2) for each day between #' D8 and D25 and for each type of effect : LC, LDD, LCD. #' A similar approach is used to compute the LT (using the data from all days). #' The aim is to choose the model with the best fit (most of the time) and best #' statistical properties. #' #' We compare the quality of the fit of each model by extracting the AIC value and computing #' the difference between the best model (with lowest AIC) and the other AIC values for a #' given time (for LC, LDD, LCD) or a given concentration (for LT). We have also computed #' a goodness of fit test (`modelFit` function from `drc` package) for each model (available in the #' raw outputs). #' We compare also graphically the differences of LCx, LDDx, LCDx and LTx computed with #' each model. #' #' Here are the formulas of the 4 models (see Ritz et al. 2010 for more details): #' #' 3 parameters Logistic Model (with d = 1) : #' $$f(x) = c + \frac{d-c}{(1+\exp(b(x - e)))}$$ #' #' 3 parameters Log-Logistic Model (with d = 1) #' $$f(x) = c + \frac{d-c}{(1+\exp(b(\log(x)-\log(e))))}$$ #' #' 3 parameters Weibull 1 Model (with d = 1) #' $$f(x) = c + (d-c) \exp(-\exp(b(\log(x)-\log(e))))$$ #' #' 3 parameters Weibull 2 Model (with d = 1) #' $$f(x) = c + (d-c) (1 - \exp(-\exp(b(\log(x)-\log(e)))))$$ #' #' ## Main results #' #' Here are the main conclusion for the data analysis below. #' #' The logistic model and Weibull2 model are almost always the best models (the ones with #' the lowest AIC). The logistic has some times a much lower AIC than the Weibull2 model #' (particularly after Day 20) however logistic model is often unable to estimate the standard #' errors of the parameters while the Weibull2 models has almost always standard errors. #' In addition LCx, LDDx, LCDx and LTx estimates from both models are often very close and #' the standard errors are quite similar. #' #' So we decided to the the Weibull2 model as a good compromise for the this study. #' #' \newpage #' # /* # ---------------- Model comparison : Lethal concentrations ------------------- # */ #' #' ## Lethal concentrations (LCx) at each time #' #' First lines of the AIC and goodness of fit results (note that ED50 stands for LC50 here): # summary(LC_AIC) pander(LC_AIC[1:12,]) #' First lines of the LCx results # summary(LC) pander(LC[1:15,]) #' #' \newpage #' #' Comparison of the delta AIC values (on a log10(x+1) scale) for the 4 types of models at #' each day after treatment (between D8 and D25). #' The horizontal dotted red line is the classical threshold of difference of AIC = 2. #' When a dot is plotted, it means that the LC50 was estimated at that day but that no #' standard error (and hence no confidence interval) could be estimated by the model). #' #' Interpretation : The logistic model is almost always the best model or very close to it however #' from day 8 to day 15 the standard error of the LC50 could not be estimated by the model. #' The Weibull2 model has generally an AIC very close to the logistic model excepted at #' days >15 where it regularly peaks far away from the best model. tmp <- LC_AIC tmp$no_SE <- is.nan(tmp$LowerCI) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(tmp, aes(y = log10(DeltaAIC+1), x = Time, color = Model)) + geom_line() + geom_point(aes(alpha = no_SE))+ geom_hline(yintercept = log10(2+1), color = "orangered", linetype = 2 ) + xlab("Time (number of days since D0)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' If you compare the LCx estimates for each rate (LC10, LC20,...) for the 4 types of models #' you can see that the estimates of the Weibull2 model and Logistic models are very close #' to each other. For the LC50 all four models provide similar estimates. #' The Weibull1 and LogLogistic models can deviate quite strongly from the others particularly #' for the estimates of LC10, LC20, LC80 and LC90. #' # #+ fig.width = 18/2.54, fig.height = 15/2.54 # dev.new(width = 18/2.54, height = 15/2.54) ggplot(LC, aes(y = Estimate, x = Time, color = Model)) + geom_line() + facet_wrap(~PctMort, scales = "free") + ylab("LCx estimate") + xlab("Time (Days)") + theme(legend.position = "top") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' Here we compare the standard errors of 3 of the models types for the LC50 estimates. #' You can see that even after day 15 where the Weibull2 models had some times #' AIC clearly higher that the Logistic model, both the estimates and their confidence #' intervals are very similar (and confidence intervals are not available for logistic models #' for days 8-15). #' #' The Weibull 2 models seem to be a good compromise that works well in almost situations. #' The confidence interval is not available on day 14 and on day 13 it abnormally small. #' However the estimates at these dates are quite similar to the estimates of other models and within the #' range of the confidence interval of the next best model (Weibull1) #' #+ fig.width = 12/2.54, fig.height = 8/2.54 # dev.new(width = 12/2.54, height = 8/2.54) ggplot(tmp[tmp$Model %in% c("Logistic", "Weibull2", "Weibull1"),], aes(y = ED50, x = Time, color = Model)) + geom_line() + geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1, position = position_dodge(0.4)) + ylab("LCx estimate") + xlab("Time (Days)") + coord_cartesian(ylim=c(0, 17000)) + scale_color_manual(name = "Model type", values = c("dodgerblue", "orangered", "gold")) #' **Conclusion** : the Weibull2 model seem to be a good compromise that can be used to #' estimates LC50 and other LCx at all days between D8 and D25. #' #' \newpage #' # /* # ---------------- Model comparison : Lethal Dietary Dose ------------------- # */ #' #' ## Lethal dietary dose (LDDx) at each time #' #' Similar approach but instead of the nominal concentration of each treatment, #' we use the effective average dose of a.i. consumed by the bees since D0. #' #' First lines of the AIC and goodness of fit results (note that ED50 stands for LDD50 here): # summary(LDD_AIC) pander(LDD_AIC[1:12,]) #' First lines of the LDDx results # summary(LDD) pander(LDD[1:15,]) #' #' \newpage #' #' #' Comparison of the delta AIC values (on a log10(x+1) scale) for the 4 types of models at #' each day after treatment (between D8 and D25). #' The horizontal dotted red line is the classical threshold of difference of AIC = 2. #' When a (colored) dot is plotted, it means that the LDD50 was estimated at that day but that no standard error #' (and hence no confidence interval) could be estimated. #' #' Interpretation : The logistic model is almost always the best model or very close to it however #' for several days, the standard error of the LCD50 could not be estimated by the model. #' The Weibull2 model has generally an AIC very close to the logistic model excepted at #' days >20 where it far away from the best model. tmp <- LDD_AIC tmp$no_SE <- is.nan(tmp$LowerCI) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(tmp, aes(y = log10(DeltaAIC+1), x = Time, color = Model)) + geom_line() + geom_point(aes(alpha = no_SE))+ geom_hline(yintercept = log10(2+1), color = "orangered", linetype = 2 ) + xlab("Time (number of days since D0)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' The different models show very similar patterns excepted the model Weibull1 that deviates #' some times from the others. For the LDD10 estimates there is also more variability between #' the models with model Weibull1 being the closest to # #+ fig.width = 18/2.54, fig.height = 15/2.54 # dev.new(width = 18/2.54, height = 15/2.54) ggplot(LDD, aes(y = Estimate, x = Time, color = Model)) + geom_line() + facet_wrap(~PctMort, scales = "free") + ylab("LDDx estimate") + xlab("Time (Days)") + theme(legend.position = "top") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' Comparison of the confidence intervals #' #+ fig.width = 12/2.54, fig.height = 8/2.54 # dev.new(width = 12/2.54, height = 8/2.54) ggplot(tmp[tmp$Model %in% c("Logistic", "Weibull2", "Weibull1"),], aes(y = ED50, x = Time, color = Model)) + geom_line() + geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1, position = position_dodge(0.4)) + ylab("LDDx estimate") + xlab("Time (Days)") + # coord_cartesian(ylim=c(0, 17000)) + scale_color_manual(name = "Model type", values = c("dodgerblue", "orangered", "gold")) #' **Conclusion** : the Weibull2 model seems to be a good compromise that can be used to #' estimate LDD50 and other LDDx at all days between D8 and D25. The logistic model has #' often a better AIC than the Weibull2 model but at several time points, the logistic #' model is unable to estimate the standard error of the LDD50. In addition the #' estimates of both models are very close to each other so using one or the other model #' should not change the final results. #' #' \newpage #' # /* # ---------------- Model comparison : Lethal cumulative Dose ------------------- # */ #' #' ## Lethal cumulative dose (LCDx) at each time #' #' Similar approach but instead of the nominal concentration we use the effective #' total dose of a.i. consumed by the bees. #' #' First lines of the AIC and goodness of fit results (note that ED50 stands for LCD50 here): # summary(LCD_AIC) pander(LCD_AIC[1:12,]) #' First lines of the LCDx results # summary(LCD) pander(LCD[1:15,]) #' #' \newpage #' #' #' Comparison of the delta AIC values (on a log10(x+1) scale) for the 4 types of models at #' each day after treatment (between D8 and D25). #' The horizontal dotted red line is the classical threshold of difference of AIC = 2. #' When a dot is plotted, it means that the LCD50 was estimated at that day but that no standard error #' (and hence no confidence interval) could be estimated by the model). #' #' Interpretation : The logistic model is almost always the best model or very close to it however #' for several days, the standard error of the LCD50 could not be estimated by the model. #' The Weibull2 model has generally an AIC very close to the logistic model excepted at #' days >20 where it far away from the best model. tmp <- LCD_AIC tmp$no_SE <- is.nan(tmp$LowerCI) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(tmp, aes(y = log10(DeltaAIC+1), x = Time, color = Model)) + geom_line() + geom_point(aes(alpha = no_SE))+ geom_hline(yintercept = log10(2+1), color = "orangered", linetype = 2 ) + xlab("Time (number of days since D0)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' The different models show similar patterns. There is a very interseting plateau of the #' LCDx followed by a drop. Without cumulative toxicity, you would expect to have #' a LCDx independent of time (flat line). # #+ fig.width = 18/2.54, fig.height = 15/2.54 # dev.new(width = 18/2.54, height = 15/2.54) ggplot(LCD, aes(y = Estimate, x = Time, color = Model)) + geom_line() + facet_wrap(~PctMort, scales = "free") + ylab("LCDx estimate") + xlab("Time (Days)") + theme(legend.position = "top") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' Comparison of the confidence intervals #' #+ fig.width = 12/2.54, fig.height = 8/2.54 # dev.new(width = 12/2.54, height = 8/2.54) ggplot(tmp[tmp$Model %in% c("Logistic", "Weibull2", "Weibull1"),], aes(y = ED50, x = Time, color = Model)) + geom_line() + ylab("LCDx estimate") + xlab("Time (Days)") + geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1, position = position_dodge(0.4)) + # coord_cartesian(ylim=c(0, 17000)) + scale_color_manual(name = "Model type", values = c("dodgerblue", "orangered", "gold")) #' **Conclusion** : the Weibull2 model seems to be a good compromise that can be used to #' estimate LCD50 and other LCDx at all days between D8 and D25. #' #' \newpage #' # /* # ---------------- Model comparison : Lethal Time ------------------- # */ #' #' ## Lethal time (LTx) at each concentration #' #' Same approach as for the Lethal Concentration but we fit 4 types of models #' of the mortality vs time for each concentration #' First lines of the AIC and goodness of fit results (note that ED50 stands for LT50 here): # summary(LT_AIC) pander(LT_AIC[1:12,]) #' First lines of the LTx results # summary(LT) pander(LT[1:15,]) #' #' \newpage #' #' Comparison of the delta AIC values for the 4 types of models at #' each concentration. #' The horizontal dotted red line is the classical threshold of difference of AIC = 2. #' When a dot is plotted, it means that the LT50 was estimated at that day but that no standard error #' (and hence no confidence interval) could be estimated by the model). #' #' Interpretation : The weibull2 model is almost always the best model. tmp <- LT_AIC tmp$no_SE <- is.nan(tmp$LowerCI) #+ fig.width = 14/2.54, fig.height = 8/2.54 # dev.new(width = 14/2.54, height = 8/2.54) ggplot(tmp, aes(y = DeltaAIC, x = Conc, color = Model)) + geom_line() + geom_point(aes(alpha = no_SE))+ geom_hline(yintercept = log10(2+1), color = "orangered", linetype = 2 ) + xlab("Boscalid concentration (mg a.i. / l syrup)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' The LT50 estimates are very similar for all 4 models. There are some differences #' for other mortality thresholds particularly for LT10. #' #+ fig.width = 18/2.54, fig.height = 15/2.54 # dev.new(width = 18/2.54, height = 15/2.54) ggplot(LT, aes(y = Estimate, x = Conc, color = Model)) + geom_line() + facet_wrap(~PctMort, scales = "free") + ylab("LTx estimate (Days)") + xlab("Boscalid concentration (mg a.i. / l syrup)") + theme(legend.position = "top") + scale_color_manual(name = "Model type", values = c("dodgerblue", "darkblue", "orangered", "gold")) #' Here we compare the standard errors of 3 of the models types for the LT50 estimates. #' #+ fig.width = 12/2.54, fig.height = 8/2.54 # dev.new(width = 12/2.54, height = 8/2.54) ggplot(tmp[tmp$Model %in% c("Logistic", "Weibull2", "Weibull1"),], aes(y = ED50, x = Conc, color = Model)) + geom_line() + geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), width = 0.1, position = position_dodge(500)) + ylab("LT50 and Confidence Interval (Days)")+ xlab("Boscalid concentration (mg a.i. / l syrup)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "orangered", "gold")) #' The same with LT10 #' #+ fig.width = 12/2.54, fig.height = 8/2.54 # dev.new(width = 12/2.54, height = 8/2.54) ggplot(LT[LT$Model %in% c("Logistic", "Weibull2", "Weibull1") & LT$PctMort == 10 ,], aes(y = Estimate, x = Conc, color = Model)) + geom_line() + geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1, position = position_dodge(500)) + ylab("LT10 and Confidence Interval (Days)")+ xlab("Boscalid concentration (mg a.i. / l syrup)") + scale_color_manual(name = "Model type", values = c("dodgerblue", "orangered", "gold")) #' **Conclusion** : the Weibull2 model seems to be a good compromise : it is always the best #' model with only one exception and the estimates are close to the other models in most circumstances #' with narrow Confidence Intervals. # /* # ---------------- Comparison models 2-3 parameters ------------------- # */ #' \newpage #' #' # Comparison of models with 2 or 3 parameters using corrected mortality or not {#cormortality} #' #' In the previous analyses we used 3 parameters models and uncorrected mortalities. #' The results of such models are identical or very similar to estimates of more #' classical 2 parameters models based on corrected mortalities. #' #' We compare here the results based on these two approaches. #' #' We have used the corrected mortalities (Abbott's formula on the % of dead) and #' performed 2 parameters models. #' This approach is not ideal for at least two reasons : #' #' 1. the corrected mortality is no more a ratio of to integers (ie because of the #' unbalance in the design) #' 2. the total number of individuals used in the model (as weights) #' remains constant while it should decrease when the correction is higher. #' However this should have mainly an impact on the standard errors and not on the #' estimate itself which is our main interest here. #' #' It could be possible to correct directly the number of dead bees and the total and then #' compute the corrected ratio. This means however that we have to group the data #' from the 3 replicates, slightly change de way the consumption is calculated and the correction #' will be unfair for the doses with a different initial number of bees. #' #' The comparison of the 2 approaches shows that the estimates are very close to #' each other particularly for the Weibull2 model used in the previous analyses. #' NB1 : the black line shows a line with slope = 1 and intercept =0 #' (ie the line showing the perfect match). #' NB2 : the few points that deviates in the LogLogistic models and Weibull1 models #' are due to 2 parameters models with a clear lack of fit and very large standard errors. #' #' tmp <- d[d$Treat %in% c("CONTROL"),] MortRateControl <- aggregate(tmp["MortRate"], tmp["Day"], mean ) colnames(MortRateControl)[2] <- "MortRateControl" tmp <- d[d$Treat %in% c("CANTUS", "CONTROL"), ] tmp <- merge(tmp, MortRateControl, all.x = TRUE, sort = FALSE) tmp$CorMortRate <- (tmp$MortRate - tmp$MortRateControl)/(1-tmp$MortRateControl) dcor <- tmp timevalues <- 8:25 LCcor_AIC <- vector("list", length = length(timevalues)) LCcor <- vector("list", length = length(timevalues)) for(i in 1:length(timevalues)) { time <- timevalues[i] PctMort <- seq(10,90,10) dataset <- dcor[dcor$Treat %in% c("CONTROL", "CANTUS") & dcor$Time == time,] # 4 type of models with 2 parameters, based on the mortality rate mLC_L <- drm(CorMortRate ~ Conc, weights = Eff, data = dataset, fct = logistic(fixed = c(NA, 0, 1, NA, 1), names = c("b", "c", "d", "e", "f")), type = "binomial") mLC_LL <- drm(CorMortRate ~ Conc, weights = Eff, data = dataset, fct = LL.2(), type = "binomial") mLC_W1 <- drm(CorMortRate ~ Conc, weights = Eff, data = dataset, fct = W1.2(), type = "binomial") mLC_W2 <- drm(CorMortRate ~ Conc, weights = Eff, data = dataset, fct = W2.2(), type = "binomial") mLC <- list(mLC_L, mLC_LL, mLC_W1, mLC_W2) names(mLC) <- c("Logistic", "LogLogistic", "Weibull1", "Weibull2") # Compare AIC and goodness of fit statistics res <- mselect2(mLC) res <- data.frame( Time = time, Model = row.names(res), res[,1:4], DeltaAIC = res$AIC - res$AIC[1], res[,5:7]) row.names(res) <- NULL LCcor_AIC[[i]] <- res # Compute Effect doses (here LC) res <- lapply(mLC, ED, respLev = PctMort, interval = "delta", display = FALSE) res <- do.call(rbind, res) row.names(res) <- NULL res <- data.frame( Time = time, PctMort = PctMort, Model = rep(names(mLC), each = length(PctMort)), res ) LCcor[[i]] <- res } LCcor_AIC <- do.call(rbind, LCcor_AIC) LCcor <- do.call(rbind, LCcor) tmp <- rbind(data.frame(LC, Mortality = "Corrected"), data.frame(LCcor, Mortality = "Uncorrected")) tmp <- data.frame(LC, corEstimate = LCcor[,4]) #+ fig.width = 12/2.54, fig.height = 10/2.54 # dev.new(width = 12/2.54, height = 10/2.54) ggplot(tmp, aes(y = Estimate, x = corEstimate)) + geom_abline(slope = 1, intercept = 0) + geom_point( shape = 1, color = "gray50", size = 2, alpha = 0.5) + facet_wrap(~Model, scales = "free") + ylab("LCx estimates based on 3 parameters \nmodels and uncorrected mortality") + xlab("LCx estimates based on 2 parameters \nmodels and corrected mortality") # /* # ---------------- Session Info & Citations ------------------- # */ #' #' \newpage #' #' # References #' #' Bretz, F., Hothorn, T., Westfall, P., Westfall, P.H., 2010. Multiple Comparisons Using R. #' CRC Press. #' #' European Food Safety Authority, 2013. Guidance on the risk assessment of plant protection #' products on bees (Apis mellifera, Bombus spp. and solitary bees). EFSA Journal 2013;11(7):3295, 266 pp. doi:10.2903/j.efsa.2013.3295. #' #' Miller, F.J., Schlosser, P.M., Janszen, D.B., 2000. Haber’s rule: a special case in a #' family of curves relating concentration and duration of exposure to a fixed level of #' response for a given endpoint. Toxicology 149, 21–34. #' #' Ritz, C., 2010. Toward a unified approach to dose–response modeling in ecotoxicology. #' Environmental Toxicology and Chemistry 29, 220–229. #' #' #' **Automatic citation of R and all packages used :** #' #+ results = 'asis', echo = FALSE cat("**R** \n") print(citation(), style = "text") cat("\n ") pkgs <- sort(names(sessionInfo()$otherPkgs)) for(i in 1:length(pkgs)) { cat(paste0("**", pkgs[i], "** \n")) print(citation(pkgs[i]), style = "text") cat("\n ") } #' \pagebreak #' #' # Session Information #' opts_chunk$set(results= "markup") sessionInfo()