Temporal dynamics of whole body residues of the neonicotinoid insecticide imidacloprid in live or dead honeybees

In cases of acute intoxication, honeybees often lay in front of their hives for several days, exposed to sunlight and weather, before a beekeeper can take a sample. Beekeepers send samples to analytical laboratories, but sometimes no residues can be detected. Temperature and sun light could influence the decrease of pesticides in bee samples and thereby residues left for analysis. Moreover, samples are usually sent via normal postal services without cooling. We investigated the temporal dynamics of whole-body residues of imidacloprid in live or dead honeybees following a single-meal dietary exposure of 41 ng/bee under various environmental conditions, such as freezing, exposure to UV light or transfer of individuals through the mail system. Immobile, “dead” looking honeybees recovered from paralysis after 48 hours. The decrease of residues in living but paralysed bees was stopped by freezing (= killing). UV light significantly reduced residues, but the mode of transport did not affect residue levels. Group feeding increased the variance of residues, which is relevant for acute oral toxicity tests. In conclusion, elapsed time after poisoning is key for detection of neonicotinoids. Freezing before mailing significantly reduced the decrease of imidacloprid residues and may increase the accuracy of laboratory analysis for pesticides.

The ("contributed") R-packages that will be used in the following are loaded: > library( Hmisc) > library( car) > library( lattice) > library( multcomp) Storing Rs current options for their recovery after generating this report: > startoptions <-options()

Import and conversion
After passing an initial, minimal format check these 3 CSV-files are imported into a 'data frame' in R, Version 3.3.2 [1, R Development Core Team]. These data frames are combined into a list with named components.

Imidachloprid
Question: Does Imidachloprid degrade in dead bees at room temperature?
Asked a bit more precisely: Does the average concentration of Imidachloprid (in dead bees at room temperature) decrease along time?
To answer this question we perform a (log-)linear regression of Imidachloprid on time t (with data at the times t ∈ {1, 24, 48}).
Points to consider: Repeated measurements on the same hive usually induce dependencies (and hence correlation) within each group of measurements from the same hive. However, thorough model diagnostics (in particular, inspection of the distribution of the residuals) revealed no evidence for a correlation structure here. In addition, the small number of only three groups (hives) did not warrant the use of a mixed-effects regression model (which could typically be employed for the analysis of data with a hierarchical grouping structure like the present).
Graphical Exploratory Data Analysis (EDA) Fig. 1 presents a first exploratory graph of the Imidachloprid data. It displays -without consideration of a potential influence of time -the distribution of the Imidachloprid values separately for each hive using "strip plots" (also called one-dimensional scatter plots) on different scales: the left panel uses the original measurement scale, the right one the decimal logarithmic scale, i.e., log 10 -scale (for reasons given below). Fig. 2 shows -color-coded for all hives overlaid -the distribution of Imidachloprid values versus time on the original scale (left) and on the decimal logarithmic scale (right). The strip plots on the original scale appear to present a slight heteroscedasticity (i. e., inhomogeneity of variances) of Imidachloprid along time: the variability (spread) of the data seems to decline with falling average Imidachloprid values (i.e., with advancing time), across hives. This is (maybe a bit over-)compensated by the (typically variance-stabilizing) logarithmic transformation. (For details regarding the log-transformation see the following remarks.) Remarks: Any logarithmic transformation -due to its non-linearity -"stretches" the lower end of the measurement scale relatively stronger than the upper end. Consequently, "cluttered" data points at the lower end of the scale become separated and wide-spread data points at the upper end get pushed together. This has typically two statistically useful consequences: it may symmetrize an asymmetric data distribution and it may homogenize variances. In addition, it ensures that the log-linearly modelled response values yield positive values on the original measurement scale (after re-transforming the model with the antilogarithm, i.e., exponentiation with the base 10).

Model fitting
Tab. 1 contains a linear regression model of Imidachloprid on time with (arbitrary) regression lines along time that vary with hive in their vertical position and in their slope ("interaction"). This allows to assess if a significant time effect is present (across hives) and if there is a difference between hives with respect to the hive-specific Imidachloprid-level or with respect to the hive-specific Imidachloprid-trend along time. (Such models are also known as "analyis of covariance" (or ANCOVA) models.) Table 1: Regression of Imidachloprid on time and hive with interaction.
The R-output in short: the Estimate-column of the Coefficients-block in tab. 1 contains the estimated values for the regression coefficients β 0 , the α h s, β 1 , and the γ h s in its rows (starting with (Intercept)).
The last column Pr(>|t|) presents the p-values of the significance tests for the respective coefficients. The Multiple R-squared value is a goodness-of-fit measure and quantifies the proportion of the observed total variability of the response variable that can be "explained" by the regression model.
Tab. 2 contains almost the same model as tab. 1. The only difference is that the response variable is log 10transformed (so that this could be called a log-linear regression model). Fig. 3 visualizes in its left part the hive-specific estimated regression lines of the fitted model on the log 10scale in separate so-called panels, one per hive, augmented by the raw data (as they are already presented in fig. 2 on the right). In its right part fig. 3 shows the respective estimated regression functions re-transformed (using the antilogarithm, i.e., exponentiation with the base 10) onto the original scale (and also augmented by the raw data as they are seen in fig. 2 on the left).  Short formal representation of the model in tab. 2: For hive h ∈ {1, 2, 3} in bee-group i ∈ {1, . . . , 5} the observed residue level at time t j with j ∈ {1, 2, 3} is denoted Y hi (t j ) and modelled as: with both α 1 and γ 1 set to 0 (zero). This model is described in short together with the corresponding R-output after equation (1).
The right part of fig. 3 shows the re-transformed estimated hive-specific regression functionŝ with the estimated regression coefficientsβ 0 ,α h ,β 1 , andγ h , whereα 1 = 0 andγ 1 = 0.   (1) and (2): in its left column for the model with Imidachloprid on its original scale, and in its right column for the log 10 -transformed response (see column titles).
Summary: None of the two models indicates serious violations of the typical model assumptions of homoscedasticity and normality of errors, as can be seen in the plots of the residuals vs. the fitted values (top row) and in the normal q-q plots for the (studentized) residuals (bottom row). Also, there are no unduly influential observations (outliers) according to Cook's distances (middle row).
Comparing the multiple R 2 -values (see Multiple R-squared in tables 1 and 2), which measure the goodnessof-fit of the models to the data, the non-transformed model appears to fit slightly better. However, its technical "advantage" to produce positive response values leads us to prefer the log 10 -transformed model.   table (ANOVA table) of the so-called sequential tests (also known as "type-I tests") for the model terms of the regression model of the log-transformed response Imidachloprid on time and hive with interaction. The ANOVA table has one row per model term and the respective term is given at the beginning of the row. It contains in its last column (Pr(>F)) the p-value of the test of the hypothesis of no influence of the respective model term, given that the model already contains the terms "above" the row's term. This is the reason why the tests are called sequential. (Warning: If the design is unbalanced -which is not the case here -the results depend on the order of the terms in the model formula.) Table 3: Sequential ("type-I tests") ANOVA table for the log-linear regression model of Imidachloprid on time and hive with interaction.

Analysis of Variance
Summary: There is no significant difference (p = 0.157) between hives in the overall log 10 (Imidachloprid)levels (averaged across time). After adjusting for a main effect of hive there is a significant time main effect (p = 1.14 × 10 −13 ). After adjusting for the main effects of hive and time there is no significant interaction effect between hive and time (p = 0.937).
Consequently, we fit the simpler log-linear regression model of Imidachloprid on only time without any hiveeffect (see tab. 4), and compare it with the previous, more complicated model to find out if all hive-effects are statistically negligeable simultaneously (see tab. 5):  Remark: The not shown qualitative diagnostic plots for the model neither indicate violations of the typical model assumptions of homoscedasticity and normality of errors, neither appear overly influential observations (outliers). The slightly smaller multiple R 2 -value is nearly fully compensated by the simplicity of this model. Summary: There is no significant difference (p = 0.417) between the two models. This means that the terms in which the two models differ provide no significant contribution to the distributional "behaviour" of Imidachloprid, i.e., no hive-effect is statistically significant. So, we end up with the following model: Short formal representation of the model in tab. 4: For hive h ∈ {1, 2, 3} in bee-group i ∈ {1, . . . , 5} the observed residue level at time t j with j ∈ {1, 2, 3} is denoted Y hi (t j ) and modelled as: Note: the regression function now does (of course) not depend on hive! So, the re-transformed estimated regression function reads: with the estimated regression coefficientsβ 0 andβ 1 , whose values are provided in the rows (Intercept) and time, respectively, of the Estimate-column of the Coefficients-table in the corresponding R-output in tab. 4. In addition, tab. 6 presents the estimated regression coefficients together with their lower and upper confidence limits for a confidence level of 95 %. Remark: From the re-transformed fitted model in (5) we deduce that for any two time points which are ∆ time units apart the ratio q of their pertaining Imidachloprid concentrations is estimated by

Some additional considerations and visualizations
So, for a given relative change q in the Imidachloprid concentration the required time to reach it can be estimated by∆ For example, the time required to reach a reduction to 50 % of any original Imidachloprid concentration, i.e., the half time, is estimated to be log 10   > Response <-c( "5-hydroxyimidacloprid" = "Hyd")
Graphical EDA Figures 6 and 7 show exploratory displays analogous to figures 1 and 2. For technical details confer with the beginning of section 2.1.1 on page 8, and regarding the R code that creates the following respective figures or regression models see the corresponding parts in section 2.1.1.
Remark: Fig. 7 indicates that a transformation appears to be unnecessary here, and, in particular, a log-transformation even counterproductive.  Tables 7 and 8 contain the linear and the log-linear regression model, respectively, of 5-hydroxyimidacloprid on time and hive with interaction, details for which are given on pages 10 and 11, and whose corresponding short formal representations are identical to those given in (1) and (2), respectively. Comparison of the goodness-of-fit measure multiple R 2 for the models (Multiple R-squared in tables 7 and 8) suggests to prefer the not-transformed model of tab. 7 (regardless of the technical advantage of the log-linear model to produce positive response values). In addition, the qualitative diagnostic plots in fig. 8 for each of the models (left: for the untransformed model; right: for the model with log 10 -transformed response) give also reason to prefer the untransformed model because its normal q-q plot for the (studentized) residuals (bottom row) looks better than the one for the log-linear model, while the plots of the residuals vs. the fitted values (top row) and the Cook-distances (bottom row) do not present noteable differences between the two models.
Summary: The linear regression model of (untransformed) 5-hydroxyimidacloprid on time and hive with interaction appears to describe the trend in 5-hydroxyimidacloprid along time well enough over the timerange observed (so that a log-linear transformation seems not necessary here).
Summary: There is a significant difference (p = 1.58×10 −5 ) between hives in the overall 5-hydroxyimidaclopridlevels (averaged across time). After adjusting for a main effect of hive there is a significant time main effect (p = 3.29 × 10 −7 ). After adjusting for the main effects of hive and time there is a significant interaction effect between hive and time (p = 9.64 × 10 −10 ). (The latter significance, in turn, means that each main effect hypothesis has to be interpreted on average across the levels/values of the respective other factor/variable.) Tab. 10 presents the estimated regression coefficients of the model in tab. 7 together with their lower and upper confidence limits for a confidence level of 95 %.  > Response <-c( "Oelefin" = "Oel")

Oelefin
Question: Does Oelefin degrade in dead bees at room temperature?
Graphical EDA Figures 10 and 11 show exploratory displays analogous to figures 1 and 2. For technical details confer with the beginning of section 2.1.1 on page 8, and regarding the R code that creates the following respective figures or regression models see the corresponding parts in section 2.1.1.
Remark: Fig. 11 indicates that a transformation appears to be unnecessary here, and, in particular, a log-transformation even counterproductive.    Tables 11 and 12 contain the linear and the log-linear regression model, respectively, of Oelefin on time and hive with interaction, details for which are given on pages 10 and 11, and whose corresponding short formal representations are identical to those given in (1) and (2), respectively. Comparison of the goodness-of-fit measure multiple R 2 for the models (Multiple R-squared in tables 11 and 12) suggests a slightly worse fit of the log-transformed model of tab. 12. But, the technical advantage of the log-linear model to produce positive response values in combination with the qualitative diagnostic plots in fig. 12 for each of the models (left: for the untransformed model; right: for the model with log 10 -transformed response) give reason to prefer the transformed log-linear model: its diagnostic plots look a little better than the ones for the untransformed model.

Summary:
The log-linear regression model of Oelefin on time and hive with interaction appears to describe the hive-specific trends in Oelefin along time well enough over the time-range observed.      Summary: There is no significant difference (p = 0.268) between hives in the overall log 10 (Oelefin)-levels (averaged across time). After adjusting for a main effect of hive there is no significant time main effect (p = 0.229). After adjusting for the main effects of hive and time there is a significant interaction effect between hive and time (p = 8.7 × 10 −5 ). (The latter significance, in turn, means that each main effect hypothesis has to be interpreted on average across the levels/values of the respective other factor/variable.) Tab. 14 presents the estimated regression coefficients of the model in tab. 12 together with their lower and upper confidence limits for a confidence level of 95 %.  > lev.of.int <-c( "RT1", "RT24", "RT48", "RT1mail") # treatment levels of interest > reflev <-"RT1mail" # reference level for, e.g., multiple comparisons 1. Extracting the rows of Bees, whose corresponding elements in method are in the set {"RT1", "RT24", "RT48", "RT1mail"}, i.e., which belong to the treatments RT1, RT24, RT48 or RT1mail, and storing the "extract" in a new data frame named B0Mail: > B0Mail <-droplevels( subset( Bees, subset = method %in% lev.of.int)) 2. Changing method's level order, so that RT1mail is the first level and hence the reference level in later analyses (like multiple comparisons): 3. Generating variables containing some "values" which will be needed and repeatedly (!) used in the following paragraphs, so that flexible and passably efficiently modifiable code can be designed, which is also simply reusable for subsequent paragraphs: > Response <-c( Imidachloprid = "Imi"); Treatment <-"method"; Group <-"hive" > WorkData <-B0Mail

Imidachloprid
Question, stated a bit more precisely: Can the process of degradation of Imidachloprid be slowed down by short freezing in comparison to the "normal" process of degradation?
Even more precisely: Is the concentrarion of Imidachloprid in treatment-method RT1mail higher than in the other three?
To this end: Two-factorial ANOVA of Imidachloprid on the treatment (in the following called method) and on hive (with interaction) with subsequent "multiple comparisons with a control" (i.e., with a reference level) and calculation of simultaneous confidence intervals for the respective differences.
Points to consider: Repeated measurements on the same hive may induce dependencies within each group of measurements from the same hive, so analoguos remarks apply as in the "points to consider" on page 8.  The latter can also be interpreted as a so-called interaction plot to assess if an interaction may be present between hive and method. Indication for interaction would be given if the hive-specific "profiles" of line segments "along" method were not parallel. This is not strong enough the case here, esp. in view of the variability of the raw data around their means, and will be confirmed by a formal test yielding a nonsignificant interaction effect.   > # Orthogonal contrasts to increase numerical and statistical "stability" > oc <-options( contrasts = c( "contr.helmert", "contr.poly")) > # 'Building' the model formula using variables (for the sake of flexibility): > form <-formula( paste( Response, "~", Group, "*", Treatment)) > fit <-aov( form, data = WorkData) # Fitting the underlying linear model. > fit <-update( fit,~.) > anova( fit)  Table 16: Sequential ("type-I tests") two-factorial ANOVA of log 10 (Imidachloprid) on method and hive with interaction between method and hive.

Analysis of Variance
and the log-transformed as with α 1 , β 1 , and γ 1 as well as α 1 , β 1 , and γ 1 set to 0 (zero). Fig. 16 displays three qualitative diagnostic plots for each of the models (8) and (9): in its left column for the model with Imidachloprid on its original scale, and in its right column for the log 10 -transformed response (see column titles).
Diagnostic summary: (Based on fig. 16.) None of the two models indicates serious violations of the model assumption of normality of errors, see the normal q-q plots for the (studentized) residuals (bottom row). However, the plot of residuals vs. fitted values for the untransformed model (top row, left) contains some hints against homoscedasticity which are not seen in the corresponding plot for the log-transformed model (top row, right). In neither model are there any markedly influential observations (outliers) according to Cook's distances (middle row), but for the log-transformed model they seem to look even a little bit more alike.
Comparing the multiple R 2 -values (goodness-of-fit measure) of 0.75 of the underlying (but not shown) untransformed linear model and 0.77 of the respective log-transformed model (also not shown), the latter appears to fit slightly better. Together with its technical "advantage" to produce positive response values this leads us to prefer the log 10 -transformed model. Statistical inference for the fitted models ANOVA summary: (Based on tab. 16.) There is no significant difference (p = 0.138) between hives in the overall log 10 (Imidachloprid)-levels (averaged across method). After adjusting for a main effect of hive there is a significant method main effect (p = 7.01 × 10 −15 ). After adjusting for the main effects of hive and method there is no significant interaction effect between hive and method (p = 0.879).

Additional considerations and visualizations: multiple comparisons
Tab. 17 contains the multiple tests for all pairwise comparisons with the reference level, i.e., control, RT1mail (MCC for short), based on the so-called Dunnett contrasts, and tab. 18 presents the pertaining simultaneous confidence intervals. Fig. 17 displays those simultaneous confidence intervals graphically.

Remark:
The following piece of R-code is for documentation purposes only and shows the technical preparations for the computations underlying tab. 17.

95% family−wise confidence level
Difference of log 10 −concentrations of Imidachloprid Further considerations and visualizations: MCC for the main effect of method Tab. 19 shows simultaneous confidence intervals for the interaction parameters of model tab. 16. According to the practical recommendation in [6, Hsu (1996)], p. 183, we consider the models with and without interaction term as practically equivalent since all simultaneous confidence intervals (and even the not shown unadjusted intervals!) of the interaction effects are not only close to, but in fact contain zero. Hence, we decide to compute multiple tests for all pairwise comparisons with the control and the pertaining simultaneous confidence intervals for the main effect of method ignoring the interaction terms.
Remark: The following piece of R-code is shown for documentation purposes only and contains the technical preparations to enable the computations of the simultaneous confidence intervals in tab. 19. Table 19: Imidachloprid (on log 10 -scale) by method and hive: simultaneous confidence intervals of the interaction parameters of model (9).

Simultaneous Confidence Intervals
Fit: aov(formula = log10(Imi)~hive + method + hive:method, data = WorkData) Tab. 20 contains for the two-way model without interaction the multiple tests for all pairwise comparisons with the control (here: RT1mail), based on the so-called Dunnett contrasts for the main effect of method, and tab. 21 presents the pertaining simultaneous confidence intervals. Fig. 18 displays those simultaneous confidence intervals graphically.     Hence, on a family-wise significance level of 95 %, RT1 is not significantly different from RT1mail, while RT24 and RT48 are both significantly different from RT1mail (with respect to their average log 10 (Imidachloprid)values).

Imidachloprid
Question, stated a bit more precisely: Is the process of degradation of Imidachloprid speeding up by heavy UV exposure in comparison to the "normal" process of degradation? In this scenario an extreme sun exposure is simulated by means of an UV lamp.
Even more precisely: Is the average concentration of Imidachloprid in treatment RT24UV lower than in the other three and especially RT24?
To this end: Two-factorial ANOVA of Imidachloprid on the treatment (in the following called method) and on hive (with interaction) with subsequent "multiple comparisons with a control" and calculation of simultaneous confidence intervals for the respective differences.
Points to consider: Repeated measurements on the same hive may induce dependencies within each group of measurements from the same hive, so analoguos remarks apply as in the "points to consider" on page 8.  The latter can also be interpreted as a so-called interaction plot to assess if an interaction may be present between hive and method. Indication for interaction would be given if the hive-specific "profiles" of line segments "along" method were not parallel. This is not strong enough the case here, esp. in view of the variability of the raw data around their means, and will be confirmed by a formal test yielding a nonsignificant interaction effect.

Model fitting & model diagnostics
Tab. 22 presents the two-factorial analysis of variance (ANOVA) table of Imidachloprid on method and hive with interaction using sequential tests (aka "type-I tests") for the model terms. For details on structure and interpretation of the ANOVA table see the explanations regarding tab. 3. Table 22: Sequential ("type-I tests") two-factorial ANOVA of Imidachloprid on method and hive with interaction between method and hive. (update() is used here for the same reason as in tab. 1.) > # 'Building' the model formula using variables (for the sake of flexibility): > form <-formula( paste( Response, "~", Group, "*", Treatment)) > fit <-aov( form, data = WorkData) # Fitting the underlying linear model. > fit <-update( fit,~.) > anova( fit)  Table 23: Sequential ("type-I tests") two-factorial ANOVA of log 10 (Imidachloprid) on method and hive with interaction between method and hive.

Analysis of Variance
Short formal representation of the models in tab. 22 and 23 is completely analogue to the one in (8) and (9), respectively. Diagnostic summary: (Based on fig. 21.) None of the two models indicates serious violations of the model assumption of normality of errors, as can be seen in the normal q-q plots for the (studentized) residuals (bottom row). However, the plot of residuals vs. fitted values for the untransformed model (top row, left) contains hints against homoscedasticity which appear to be (only a little) weaker in the corresponding plot for the log-transformed model (top row, right). In neither model are there any markedly influential observations (outliers) according to Cook's distances (middle row), but for the log-transformed model they appear a bit more alike.
Comparing the multiple R 2 -values (goodness-of-fit measure) of 0.8 of the underlying (but not shown) untransformed linear model and 0.76 of the respective log-transformed model (also not shown), the first appears to fit better. In view of this and the diagnostic summary, we prefer the untransformed model (despite its technical "disadvantage" of possibly producing negative response values).  Statistical inference for the fitted models ANOVA summary: (Based on tab. 22.) There is no significant difference (p = 0.362) between hives in the overall Imidachloprid-levels (averaged across method). After adjusting for a main effect of hive there is a significant method main effect (p = 4.56 × 10 −16 ). After adjusting for the main effects of hive and method there is no significant interaction effect between hive and method (p = 8.4 × 10 −2 ).

Additional considerations and visualizations: multiple comparisons
Tab. 24 contains the multiple tests for all pairwise comparisons with the control, i.e., reference level RT24UV (MCC for short), based on the so-called Dunnett contrasts, and tab. 25 presents the pertaining simultaneous confidence intervals. Fig. 22 displays those simultaneous confidence intervals graphically.

Remark:
The following piece of R-code is for documentation purposes only and shows the technical preparations for the computations underlying tab. 24.

95% family−wise confidence level
Difference of Imidachloprid−concentrations Further considerations and visualizations: MCC for the main effect of method Tab. 26 shows simultaneous confidence intervals for the interaction parameters of model tab. 22. According to the practical recommendation in [6, Hsu (1996)], p. 183, we consider the models with and without interaction term as practically equivalent since all but two of the simultaneous confidence intervals of the interaction effects are not only close to, but in fact contain zero. The two that do not cover zero are anyhow very close to zero (and quite long). Hence, we decide to compute multiple tests for all pairwise comparisons with the control, i.e., the reference level, and the pertaining simultaneous confidence intervals for the main effect of method ignoring the interaction terms. Table 26: Imidachloprid by method and hive: simultaneous confidence intervals of the interaction parameters of model (8).

Simultaneous Confidence Intervals
Fit: aov(formula = Imi~hive + method + hive:method, data = WorkData) Tab. 27 contains for the two-way model without interaction the multiple tests for all pairwise comparisons with the control (here: the reference level RT24UV), based on the so-called Dunnett contrasts for the main effect of method, and tab. 28 presents the pertaining simultaneous confidence intervals. Fig. 23 displays those simultaneous confidence intervals graphically. > summary( mcc2 <-glht( fm, linfct = mcp( method = "Dunnett")))    Hence, on a family-wise significance level of 95 %, RT1 and RT24 are both significantly different from RT24UV, while RT48 is not significantly different from RT24UV (with respect to their average Imidachlopridvalues).

Imidachloprid
Question, stated a bit more precisely and with additional information: There are actually two questions here: Does the feeding method (single vs. group feeding) influence the measurement results of Imidachloprid with respect to their . . . Note 1: In the main experiment bees were fed individually to minimise "dilution"-effects by trophalaxis (feeding gustation drops from one bee to another). In the group feeding experiment 10 bees had access to the 10-fold volume of the same food source as in the other experiments. Every other factor was kept similar. Note 2: Variance is analysed first because the result typically influences the analysis of location.
The key questions even more precisely: a) Is the variance of Imidachloprid-concentrations in treatment RT24GF different from that in RT24 (taking hives into consideration!)?
b) Is the average concentration of Imidachloprid in treatment RT24GF different from that in RT1, RT24 or RT48 (taking hives into consideration!)?

Analysis of Variance
> anova( fitlog <-update( fit, log10( .)~.)) Short formal representation of the models in tab. 36 and 37 is completely analogue to the one in (8) and (9), respectively. Diagnostic summary: (Based on fig. 26.) Neither model indicates a violation of the model assumption of normality, as can be seen in the normal q-q plots for the (studentized) residuals (bottom row). Homoscedasticity of errors is doubtful for the model with untransformed values, though, readable from the plot of residuals vs. fitted values (top row) and from Bartlett's test in tab. 31. In neither model are there any markedly influential observations (outliers) according to Cook's distances (middle row), albeit they appear a little bit more alike in the log-transformed model.

Analysis of Variance
Comparing the multiple R 2 -values (goodness-of-fit measure) of 0.81 of the underlying (but not shown) untransformed linear model and 0.82 of the respective log-transformed model (also not shown), both appear to fit almost equally well. In view of this and of the diagnostic summary, we prefer the log-transformed model (further supported by its technical "advantage" of producing only positive response values). 5 10 15 20  ANOVA summary: (Based on tab. 37.) There is no significant difference (p = 0.616) between hives in the overall Imidachloprid-levels (averaged across method). After adjusting for a main effect of hive there is a significant method main effect (p = 2.51 × 10 −17 ). After adjusting for the main effects of hive and method there is no significant interaction effect between hive and method (p = 8.83 × 10 −2 ).  MCC summary: Fig. 27 displays the simultaneous confidence intervals for all pairwise comparisons with a control (here: the reference level RT24GF) graphically: Each horizontal line segment that does not intersect the dashed vertical line at zero indicates that the respective difference of log 10 -concentrations (indicated on the vertical axis on the left) is significantly different from zero. This means approximately, that the ratio of Imidachloprid's concentrations of the two pertaining treatments are significantly different from 1 (one) on the "original" scale. Hence, on a family-wise significance level of 95 %, RT48 and RT24GF are significantly different in all hives, while on the one hand RT24 and RT24GF are significantly different in hives 1 and 3, but not in hive 2, and on the other hand RT1 and RT24GF is only significantly different in hive 1, but not in hives 2 and 3 (with respect to their average log 10 (Imidachloprid)-values).

95% family−wise confidence level
Difference of log 10 −concentrations of Imidachloprid Further considerations and visualizations: MCC for the main effect of method Tab. 40 shows simultaneous confidence intervals for the interaction parameters of model tab. 37. According to the practical recommendation in [6, Hsu (1996)], p. 183, we consider the models with and without interaction term as practically equivalent since all but two of the simultaneous confidence intervals of the interaction effects are not only close to, but in fact contain zero. The two that do not cover zero are anyhow very close to zero (and quite long). Hence, we decide to compute multiple tests for all pairwise comparisons with the control, i.e., the reference level, and the pertaining simultaneous confidence intervals for the main effect of method ignoring the interaction terms. Table 40: Imidachloprid (on log 10 -scale!) by method and hive: simultaneous confidence intervals of the interaction parameters of model (9).

Simultaneous Confidence Intervals
Fit: aov(formula = log10(Imi)~hive + method + hive:method, data = WorkData) Tab. 41 contains for the two-way model without interaction the multiple tests for all pairwise comparisons with the control (here: the reference level RT24GF), based on the so-called Dunnett contrasts for the main effect of method, and tab. 42 presents the pertaining simultaneous confidence intervals. Fig. 28 displays those simultaneous confidence intervals graphically. Table 41: Imidachloprid (on log 10 -scale!) by method and hive: multiple tests for all pairwise comparisons with the control, i.e., the reference level RT24GF.

95% family−wise confidence level
Difference of log 10 −concentrations of Imidachloprid Figure 28: Imidachloprid (on log 10 -scale!) by method and hive: simultaneous confidence intervals for all pairwise comparisons with RT24GF ignoring interaction effects. (File name: MS-Q4 AOV Imi by method and hive SimCIplot2.pdf ) MCC Summary 2: Fig. 28 displays the simultaneous confidence intervals for all pairwise comparisons with a control (here: the reference level RT24GF) graphically: Each horizontal line segment that does not intersect the dashed vertical line at zero (which is invisible here because it is outside the presented part of the horizontal axis) indicates that the respective difference of log 10 -concentrations (indicated on the vertical axis on the left) is significantly different from zero. This means approximately, that the ratio of Imidachloprid's concentrations of the two pertaining treatments are significantly different from 1 (one) on the "original" scale.
Hence, on a family-wise significance level of 95 %, RT1, R24, and RT48 are all significantly different from RT24GF (with respect to their average log 10 (Imidachloprid)-values).

SFig1
Figure 29: Supplementary Fig.1: Dose-dependent effects on locomotory behaviour. Individual honeybees were fed with a single dose of sugar syrup (control) or a dosage of 3.7 and 41 ng/bee imidacloprid (20 bees per treatment group, single experiment). After one, 12, 24, 36, 48, 60, and 72 hours bees were examined for activity. The locomotory behaviour was categorised in (black) immobility, (shaded area) movements of body parts, including ventilation movements of the abdomen, and (white) coordinated activity (e.g. standing, hanging, or walking). After 12 hours, 68% of the bees of the 41 ng treatment group were immobile and 31% only showed movement of body parts. After 48 hours the 41 ng treatment group had the same number of completely immobile bees as the control but the percentage of bees with coordinated movements was lower in this group.