I conducted a toxicology study where I exposed fish to 6 concentrations of an effluent treatment for 10 days. At the end of the experiment, I measured 10 endpoints with the fish plasma to determine effects of treatment.
After collecting the data, I ran a Principal Component Analysis on all endpoints. PC1 here is a principal component that is heavily driven by 4 metabolic endpoints that impact growth, fat metabolism, and protein production. When compared to treatment, PC1 appears to show hormesis, a biphasic dose-response to an environmental agent or toxic substance characterized by a low dose stimulation or beneficial effect and a high dose inhibitory or toxic effect. All 4 of the metabolic endpoints I measured follow the same trend as PC1.
To evaluate hormesis, I am using a 5-parameter Brain-Cousens hormesis model* with the drc package in R. It fits the response $R$ as a function of concentration $x$ to the following nonlinear form:
$$R(x) = c + \frac{d-c+fx}{1+\exp[b(\log(x)-\log(e))]} .$$
Most of the parameters represent the concentration-response curve without hormesis: $c$ is the lower asymptote of response, $d$ is the upper asymptote, $b$ represents the steepness of the concentration-response curve, and $e$ is the concentration at the midpoint of the curve. $f$ is the size of the hormesis effect; if $f=0$ then there is no hormesis and there is a standard log-logistic concentration-response curve. Positive values of $f$ represent hormesis.
My model seems to be driven by one value in the highest treatment concentration** but I think the value is likely real, especially since this is likely below the EC50 (based on the way the model works). I hypothesize that the PC1 and 4 endpoint concentrations would continue to decrease (increase in the plot below) if a higher treatment concentration was applied.
The question: How should one handle outliers in this type of modified log-logistic model?
Analysis so far
The Brain-Cousens model in R would not work with my data unless I inverted the values to follow this shape. The "PC1" values are actually inverted here. In reality, PC1 and the 4 endpoints begin to decrease at the 20% treatment. So what I actually see in the raw data is improved (increased concentrations of) metabolic endpoints at the 0.4-5.3% treatments.
I will not be using this as a predictive model, just to describe how our exposure may be affecting the fish. I pasted a figure of the data from the model described. The points are the inverted PC1 values and the line is from the Brain-Cousens model. My code for the model with the below data is:
hm.m2 <- drm(PC1 ~ Treatment, data = data, fct = BC.5()) #BC.5 is Brain-Cousens
summary(hm.m2) #f parameter p-value is 0.0021
modelFit(hm.m2) #lack of fit is not significant
Without 20-3 data point, the $f$ parameter $p$-value in the summary(hm.m2) line is 0.77
Data:
| ID# | Treatment | PC1 |
|---|---|---|
| 0.1-1 | 0.1 | 0.576 |
| 0.1-2 | 0.1 | 0.705 |
| 0.1-3 | 0.1 | 0.518 |
| 0.1-4 | 0.1 | 1.939 |
| 0.4-1 | 0.4 | 0.588 |
| 0.4-2 | 0.4 | -0.926 |
| 0.4-3 | 0.4 | -0.639 |
| 0.4-4 | 0.4 | -3.321 |
| 0-1 | 0 | -1.312 |
| 0-2 | 0 | 0.949 |
| 0-3 | 0 | 0.569 |
| 0-4 | 0 | 0.615 |
| 1.4-1 | 1.4 | -1.679 |
| 1.4-2 | 1.4 | -1.557 |
| 1.4-3 | 1.4 | -1.273 |
| 1.4-4 | 1.4 | -0.994 |
| 20-1 | 20 | 0.626 |
| 20-2 | 20 | 0.142 |
| 20-3 | 20 | 7.039 |
| 20-4 | 20 | 1.598 |
| 5.3-1 | 5.3 | -2.498 |
| 5.3-2 | 5.3 | 0.704 |
| 5.3-3 | 5.3 | -1.480 |
| 5.3-4 | 5.3 | -0.889 |
*If you're interested in looking into the Brain Cousens and other Hormesis models, this website provides the very basic info about how to run the model in R.
There is not too much info about Brain Cousens models online but there is a good deal of info in the scientific literature.
This is the original article describing the model;
Ritz, et al. 2015 describe some dose response models in R; Knezevic, et al 2007. Weed Technology (21): 841 – 848 is another helpful resource for dose response models and is downloadable on google scholar. Lastly, this paper about hormetic responses seems to be freely downloadable
**To test for outliers, I used a rosnerTest() from the EnvStats package on the principal component values. The outlier is ID# 20-3 in the data. I tried to look at Cook's distance but it does not work with my hormesis model output.
Best Answer
The outlier problem is relatively easy to address. There is, however, an additional problem here that illustrates the risks of nonlinear curve fitting.
The
drcpackage provides many useful tools for fitting concentration-response curves. Some tools are designed to handle the apparent outliers that frequently are seen in such data.With your model for data that you display, you used the default least-squares fitting of mean values at each concentration. That necessarily puts a lot of weight on outliers. You can, however, choose from a variety of robust fitting methods that lessen the influence of apparent outliers on the model. For example, you could choose the
robust = "tukey"setting instead to have a model that applies a Tukey biweight to the residuals. Then you would get:Most of the coefficient values aren't much different from those in the model based on least squares.
The additional problem
All is not well, however. The Brain-Cousens model for hormesis only makes physical sense if both
b(related to curve steepness) andf(hormesis) are positive. In this model and in your model they have different signs. The inability to distinguish the value ofbfrom 0 should also be troubling. Ifb = 0in the Brain-Cousens model, then there is no log-logistic contribution to the concentration-response curve; the denominator of the second term of the equation has a value of 2 independent of concentrationx.This problem seems to come from the way that the
drm()function initializes its search for the coefficient values. This page has examples of what can go wrong if poor starting choices are made in fitting a nonlinear model. You might have found a particular fit of your data to the formula, but the result doesn't make sense in the context of positive hormesis. That's why you thought that you had to "invert the values" to fit the model. You were looking for a positive value offbut didn't recognize that you only got that by trading off for a negative value ofb. With nonlinear fits you need to take extra precautions to make sure that your model makes sense.If you think that a 5-parameter Brain-Cousens model fits your data, then you should be able to add a constant to all of your response values and still get an appropriate fit. In terms of the model parameters, that should just shift the values of
c(lower asymptote) andd(upper asymptote) in parallel. Let's do that while putting the PC1 values into their original direction so that this looks more like standard hormesis. The following makes all of the adjusted response values positive.With this adjustment to the response values, both
bandfare clearly positive andbis clearly distinguished from a value of 0. I suspect that the initialization withdrm()for the nonlinear fitting is designed to work well with the non-negative response values that are typical in this type of work and just didn't perform well with your combination of positive and negative values.Back to the outlier
Now try using the adjusted data to fit the model with least squares and plot the data and the predictions from both that model and the robust model:
The least-square and robust models don't disagree too badly at the lower concentrations, but they diverge dramatically at the concentration of 20 and at the higher concentrations where you expect that the trends would continue.
Other issues
Even though the robust fit seems to get around the apparent outlier and you have a "statistically significant" value of
f, I'd be hesitant to claim that your model actually documents hormesis.If you perform ANOVA on the response values at lower concentrations and then do a Dunnett test to examine differences from values at the baseline concentration of 0, you won't find any significant differences. Thus a claim of hormesis would depend very heavily on knowing that the Brain-Cousens model represented your data correctly.
Unfortunately, you don't have data at the higher concentrations that would be needed to document the reliability of that model near the lower asymptote
c. I'd particularly worry about this here, as your responses, principal-component values, don't necessarily have the asymptotic association with concentration you would expect for responses like mortality that the Brain-Cousens model was designed to handle.Furthermore, you barely have more concentrations (0 effluent, plus 5 concentrations) than you have parameters to fit (5). That can lead to overfitting and a model that doesn't describe anything beyond the data at hand, even though you have 4 replicates at each concentration. More data at higher concentrations would seem to be needed.