Styjun

Context

Often in dose-response models we regress some range of doses against a response variable, but we are really interested in identifying the dose required to elicit a particular response. Typically this is done with inverse regression techniques (i.e. after-fitting / reparameterisation). EDIT: To clarify - this is commonly done when you need to estimate the dose required to kill say 50%, or 99.99% for quarantine protocols. To derive these estimates people employ inverse regression techniques - the above link goes through this more carefully (see page 9).

Question

How can I carry out these inverse regression procedures using methods like robust linear models, quantile regression, or machine learning models (i.e. neural networks or support vector machines)? EDIT: To clarify, I want a programming solution to how I can estimate the dose required to elicit a response of 99.99 when the model I have fitted is one of the above mentioned. I have fitted example models below to these ends.

My data looks like this:

 df <- structure(list(Response = c(100, 91.1242603550296, 86.9822485207101, 
100, 0, 0, 90.5325443786982, 95.8579881656805, 88.7573964497041, 
96.4497041420118, 82.2485207100592, 99.4082840236686, 99.4082840236686, 
98.8165680473373, 91.7159763313609, 59.1715976331361, 44.9704142011834, 
0, 100, 95.2662721893491, 100, 82.8402366863905, 7.69230769230769, 
81.6568047337278, 62.7218934911243, 97.6331360946746, 73.9644970414201, 
8.87573964497041, 0, 98.8165680473373, 78.1065088757396, 98.2248520710059, 
52.6627218934911, 96.4497041420118, 52.0710059171598, 0, 62.043795620438, 
84.6715328467153, 97.8102189781022, 4.37956204379562, 89.051094890511, 
99.2700729927007, 99.2700729927007, 97.0802919708029, 81.7518248175183, 
80.2919708029197, 90.5109489051095, 99.2700729927007, 96.3503649635037, 
0, 0, 94.8905109489051, 79.5620437956204, 67.8832116788321, 73.7226277372263, 
100, 97.0802919708029, 93.4306569343066, 86.8613138686131, 33.5766423357664, 
32.1167883211679, 46.7153284671533, 98.5401459854015, 95.6204379562044, 
86.1313868613139, 14.5985401459854, 92.7007299270073, 86.1313868613139, 
0, 77.3722627737226, 89.051094890511, 80.2919708029197, 98.1818181818182, 
96.3636363636364, 30.9090909090909, 0, 60.9090909090909, 100, 
0, 83.6363636363636, 88.1818181818182, 97.2727272727273, 0, 0, 
99.0909090909091, 100, 100, 91.8181818181818, 88.1818181818182, 
46.3636363636364, 50.9090909090909, 99.0909090909091, 97.2727272727273, 
100, 0, 92.7272727272727, 60.9090909090909, 90.9090909090909, 
57.2727272727273, 76.3636363636364, 94.5454545454545, 50, 98.1818181818182, 
16.3636363636364, 87.2727272727273, 92.7272727272727, 87.2727272727273, 
88.1818181818182, 10.7438016528926, 91.7355371900827, 98.3471074380165, 
60.3305785123967, 95.8677685950413, 0, 63.6363636363636, 71.900826446281, 
0, 74.3801652892562, 76.8595041322314, 0, 61.9834710743802, 0, 
0, 0, 84.297520661157, 47.1074380165289, 69.4214876033058, 97.5206611570248, 
100, 61.1570247933884, 90.0826446280992, 78.5123966942149, 10.7438016528926, 
100, 98.3471074380165, 100, 98.3471074380165, 93.3884297520661, 
90.9090909090909, 57.8512396694215, 57.8512396694215, 92.5619834710744, 
77.6859504132231, 69.4214876033058), Covariate = c(20, 14, 14, 
20, 0, 0, 14, 14, 14, 16, 10, 20, 20, 20, 16, 10, 10, 0, 16, 
16, 16, 10, 0, 12, 10, 12, 12, 0, 0, 20, 12, 16, 10, 12, 12, 
0, 14, 14, 16, 0, 14, 20, 16, 20, 14, 12, 12, 20, 20, 0, 0, 14, 
12, 10, 10, 20, 16, 16, 14, 10, 10, 10, 20, 16, 10, 0, 12, 12, 
0, 12, 16, 14, 16, 14, 0, 0, 12, 20, 0, 12, 14, 14, 0, 0, 20, 
20, 20, 14, 14, 10, 10, 20, 16, 16, 0, 12, 10, 10, 10, 16, 16, 
12, 20, 10, 12, 12, 16, 14, 0, 16, 20, 12, 14, 10, 10, 0, 0, 
12, 12, 10, 10, 0, 0, 0, 14, 12, 12, 20, 20, 14, 14, 14, 12, 
20, 20, 20, 16, 16, 14, 10, 10, 16, 16, 16)), row.names = 433:576, class = "data.frame")

with my formula usually being something like:

Response ~ Covariate + I(Covariate^2)

Here is an example of the models I have fitted:

#Robust linear model
MASS::rlm(Response ~ Covariate + I(Covariate^2), data = df)

#Quantile regression
quantreg::rq(Response ~ Covariate + I(Covariate^2), data = df, tau = c(0.5, 0.95)) # In this case I want to predict the specified quantiles for the dose required to elicit a given response, although I realised this code doesn't do that...

#Machine learning algorithms were trained with caret
TRControl <- trainControl(method = "cv")

#Neural Network
caret::train(Response ~ Covariate, data = df, method = "neuralnet", trControl = TRControl)

#Support Vector Machine
caret::train(Response ~ Covariate, data = df, method = "polySVM", trControl = TRControl)

Further to my comments above, your data doesn't really resemble that of a typical dose-response measurement

library(ggplot2)
ggplot(df, aes(Covariate, log10(Response))) +
 geom_point()

Here I assume that Covariate is the dose/concentration.

Do the different measurements for every Covariate relate to different experiments/groups? Do you plan on fitting multiple dose response curves do different groups in order to compare them?

A possible analysis strategy

Here is something that might give you some ideas. I'm using drc here because it allows me to fit a "sensible" dose-response curve to your data. A sensible dose-response model has horizontal asymptotes for dose → 0 and dose → ∞.

1. 
   

   In this particular example we fit a four parameter Weibull function to your data.

   
   
   

   library(drc)
   model <- drm(Response ~ Covariate, data = df, fct = W2.4())
   

   
   


2. 
   

   Let's plot original data and model predictions (including confidence interval)

   
   
   

   library(tidyverse)
   df.pred <- data.frame(Covariate = 1.1 * seq(min(df$Covariate), max(df$Covariate), length.out = 20)) %>%
    bind_cols(as.data.frame(predict(model, data.frame(Covariate = Covariate), interval = "confidence"))) %>%
    rename(Response = Prediction)
   
   ggplot(df, aes(Covariate, Response)) +
    geom_point() +
    geom_line(data = df.pred, aes(Covariate, Response), color = "blue") +
    geom_ribbon(data = df.pred, aes(x = Covariate, ymin = Lower, ymax = Upper), fill = "blue", alpha = 0.2)
   

   
   
   

   

   
   


3. 
   

   We can now use uniroot to determine specific LDx values, which are defined as the dose required to reduce the maximum response by x / 100.

   
   
   

   getLDx <- function(model, x = 0.5) 
    maxResponse <- max(predict(model, data.frame(x = c(0, Inf))))
    uniroot(
    function(Covariate) predict(model, newdata = data.frame(Covariate = Covariate)) - x * maxResponse,
    interval = range(Covariate))$root
   
   

   
   
   

   This is basically an inversion of the model, so perhaps this is what the authors of the papers you link to in your original post refer to as "inverse regression techniques".

   
   


4. 
   

   Let's calculate the LD50 value (i.e. the dose required to reduce the response by 50%)

   
   
   

   getLDx(model, x = 0.5)
   #[1] 9.465188
   

   
   
   

   From an inspection of the plot you can see that this value indeed corresponds to the dose where the response is 50% of the maximum response value.