2
$\begingroup$

So from what I have been able to gather I don't think there is any way to compare the models in question but I figured I would put it to better minds than my own before I gave up. Here is the situation:

The data I have follows a sigmoidal dose-response curve something like this:

 df <- data.frame(Effect = c(0,0,0.025,0.05,0.1,0.15,0.2,0.3,0.4,0.55,0.7,0.8,0.9,0.95,0.97,0.98,0.99,1,1), Dose= seq(0:18))

 df$Effected <- df$Effect * 100
 df$Uneffected <- 100 - df$Effect

 plot(Effect ~ Dose, data = df)

From these data I want to determine the dose required to achieve a pre-determined effect. The effected and uneffected columns are the counts of individuals who have been effected by the dose. Classically, this has been done by fitting glms with binomial errors like so:

fit <- glm(cbind(Effected, Uneffected) ~ Dose, family = binomial(link = logit), data = df)

From here Fieller's formula is employed to back calculate the dosage required for a certain effected rate, along with confidence intervals. This formula though, as it currently stands, is incapable of dealing with quadratics meaning that we end up fitting straight lines to curved data. To make matters worse, these kinds of data are regularly over-dispersed relative to binomial dispersion - in some cases as much as 40 times the normal binomial dispersion so it is obviously an issue. There has been some consideration of using beta-binomial models to deal with this but that is something of a separate issue.

So as an alternative myself and some of my colleagues have been exploring simply swapping the terms around and modelling things like so:

fit2 <- lm(log(Dose) ~ Effect + I(Effect^2), data = df)  

This type of model gets around the use of fiellers formula as we can directly predict the variable of interest and as such can use polynomials. Is this appropriate though given that dose is technically independent and effect is dependent on the response?

The other issue we have come up against is how one can determine which approach is better in a particular case. Standard metrics aren't an option as the response changes. Could one compare assumption checks? If so how? We have also considered simulation experiments and seeing which type gets consistently closest to the true value. Are there any other alternatives?

In our case, the end goal is to find the best predicting model that will give the most accurate estimates, in particular for the upper end of the effect data (i.e. > 0.9). For a given response variable we have been comparing models using cross-validation to calculate prediction error metrics (I.e. MdAE). Is there any way of doing something similar when the response variable changes but the data included in the model altogether remains the same?

Lastly, solutions would ideally be possible in R as this is the system I am most proficient in; thank you in advance!

EDIT: Fixed mistake in the glm code

$\endgroup$
1
$\begingroup$

A couple of points:

  • There is a mismatch between the data you sent and your description, i.e., in the example, you have non-integer Effected, which you should not expect to have under the Binomial distribution. Do you perhaps have a continuous outcome in the $[0, 1]$ interval? If yes, then you could consider the Beta model or an extended Beta model to account for the zeroes and ones.
  • You could still obtain the dose for a target effect by calculating predictions for a fine grid of doses and then inverting the relationship. For an example, check the following:
df <- data.frame(Effect = c(0,0,0.025,0.05,0.1,0.15,0.2,0.3,0.4,0.55,0.7,0.8,
                            0.9,0.95,0.97,0.98,0.99,1,1), 
                 Dose = seq(0:18))

df$Effected <- round(df$Effect * 100)
df$Uneffected <- 100 - df$Effected 

############################################################

# Fit the model
fit <- glm(cbind(Effected, Uneffected) ~ poly(Dose, 2), family = binomial(), 
           data = df)

# Fine predictions for different doses
nDF <- with(df, data.frame(Dose = seq(min(Dose), max(Dose), length.out = 1000)))
nDF$preds <- predict(fit, newdata = nDF, type = "response")

# plot
plot(preds ~ Dose, data = nDF, type = "l")

# find dose for target effect
target_effect <- 0.835
dose_target_effect <- nDF$Dose[which.min(abs(nDF$preds - target_effect))]

# see it in the plot
lines(c(0, dose_target_effect), c(target_effect, target_effect), col = 2)
lines(c(dose_target_effect, dose_target_effect), c(0, target_effect), col = 2)

$\endgroup$
  • $\begingroup$ Thank you Dimitris! I made a mistake in my original code - the glm should have used the integer column Effected. It's good to see another strategy for predicting required dosage that circumvents Fieller's formula but I am still quite interested in making a direct comparison between models when the responses are swapped. Any ideas? $\endgroup$ – André.B Dec 4 '18 at 1:22
0
$\begingroup$

There are several issues for working with fit2 <- lm(log(Dose) ~ Effect + I(Effect^2), data = df) 1) The x-variable does not have homogeneous errors. How the errors might feed through to the y-scale would need to be sorted out. 2) Variability on the y-scale has an errors-in-x effect once transferred to the horizontal scale, and will strongly attenuate the slope. 3) One is using the wrong regression line, x on y rather than y on x. The line will on this account be flat relative to the way that the y on x line translates to the plot of x vs y. The errors-in-x effect will further flatten it. 4) On the y vs x plot, the effect will be a very steep line, much steeper than the y on x regression line. This will lead to unrealistically small LT99s etc.

On errors in x, see Maindonald and Braun: Data Analysis and Graphics Using R (CUP, 2010), pp.203-208. See also the help page for the errorsINx function in the DAAG package for R. Running the example code will show a simulation of errors-in-x effects.

One can adapt Fieller’s formula to work with quadratics — one needs to find roots of a nonlinear function, and it is not straightforward. One may get very wide CIs, because the way that the curve flexes is not very well defined. Often one can get a plausibly linear relationship by leaving off the first one or two points; working with that line then makes better sense. We really do not know how mortality extrapolates into high mortality regions, and that is perhaps what should worry us, more than some minor nonlinearity, if it really is minor.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.