Dose Response Modeling
Crista Moreno
October 14, 2018
Copyright (C) 2018 Crista Moreno
Dose Response Modeling is free software: you can redistribute
it and/or modify it under the terms of the GNU General Public
License as published by the Free Software Foundation, either
version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be
useful, but WITHOUT ANY WARRANTY; without even the implied
warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public
License along with this program. If not, see
<https://www.gnu.org/licenses/>.Dose Response Data
\[\textbf{x} = \{x_1, x_2, ..., x_N\}\]
where each data point \(x_i\) is a tripple \(x_i = (d_i, n_i, y_i)\).
\(d_i \text{ - dose of toxic agent}\)
\(n_i \text{ - number of subjects exposed to dose } d_i\)
\(y_i \text{ - number of subjects that test positive}\)
\(N \text{ - total number of dose levels}\)
Cryptosporidium (IOWA) Data Set
data <- read.csv("./data/cryptosporidium_iowa_data.csv", header = TRUE)
dose <- data$dose
subjects <- data$subjects
positives <- data$positivesModels \(F(d; \theta)\)
Beta-Poisson Model Approximation Furumoto and Mickey
\[F(d; \alpha, \beta) = \left\{ \begin{array}{lr} 1 - \left(1 + \dfrac{d}{\beta} \right)^{-\alpha}& : d \geq 0 \\ 0 & : d < 0 \end{array} \right.\]
where \(\theta = (\alpha, \beta)\).
beta_poisson_approximation <- function(dose, theta) {
d <- dose
alpha <- theta[1]
beta <- theta[2]
return (1 - (1 + (d / beta))^(-alpha))
}Exponential Model
\[F(d; \lambda) = \left\{ \begin{array}{lr} 1 - e^{- \lambda d}& : d \geq 0 \\ 0 & : d < 0 \end{array} \right.\]
where \(\theta = \lambda\).
exponential <- function(dose, theta) {
d <- dose
lambda <- theta
return(1 - exp(-lambda * d))
}Functions
clamp <- function(p) {
if (p > 1 - 10^(-15)) {
p = 1 - 10^(-15)
} else if (p <= 0){
p = 1e-15
} else {
# do nothing
}
return(p)
}Parameter Estimation by Minimizing the Negative Log Likelihood
For each dose level of an experiment indexed by \(i\) we define
\[\pi_i(\theta) := F(d_i; \theta)\]
to be the model prediction as a function of the parameter \(\theta\) and the dose \(d_{i}\).
Binomial Likelihood
\[P_i(\theta) := {n_i\choose y_i} \pi_i(\theta)^{y_i}(1 - \pi_i(\theta))^{n_i - y_i}\]
Likelihood Function
\[\mathcal{L}(\theta) := \Pi_i^N P_i(\theta) \]
The coefficient \({n_i \choose y_i}\) is dropped since it acts as a scalar.
Log Likelihood
\[\ln(\mathcal{L}(\theta)) = \sum_i^{N} \ln(P_i(\theta))\]
\[ = \sum_i^N y_i \ln(\pi_i(\theta)) + (n_i - y_i)\ln(1 - \pi_i(\theta))\]
To find the maximum of the log likelihood we take the negative of the result, and find the \(\theta\) that minimizes it.
log_likelihood <- function(theta, dose, subjects, positives, N, M) {
sum = 0
for (i in 1:N) {
y <- positives[i]
n <- subjects[i]
pi_theta <- M(dose[i], theta)
if(is.nan(pi_theta)) {
next
}
p <- clamp(pi_theta)
if (y > 0) {
term1 = y * log(p)
} else {
term1 = 0
}
if (n != y) {
term2 = (n - y) * log(1 - p)
} else {
term2 = 0
}
sum = sum + term1 + term2
}
return(-2*sum)
}Finds the parameters for which the negative log likelihood is a minimum.
Exponential Best Fit Parameter
Let \(\theta_{0} = 0.0001\) be the initial guess for \(\theta\).
# initial guess for theta
theta0 = 0.0001
d <- dose
n <- subjects
y <- positives
N <- length(d)
# Exponential Best Fit Parameter
mle_exp = nlminb(theta0, function(theta)
log_likelihood(theta, d, n, y, N, exponential), lower = 0)
# mle_exp
exponential_theta = mle_exp$par
print(paste("exponential parameter : ", exponential_theta))## [1] "exponential parameter : 0.00419109099824652"Graph of Exponential Model
colors <- c("magenta","orchid4", "deepskyblue3")
legend_names <- c("Data", "Exponential")plot(d, y / n, pch = 19, col = colors[1], lwd=4, ann=FALSE,
ylim=c(0,1), log="x",
xlim=c(10^floor(log(min(d), 10)), 10^ceiling(log(max(d), 10))),
xaxp=c(10^floor(log(min(d), 10)), ceiling(log(max(d), 10)), 1),
yaxp=c(0,1,10), cex.axis = 1)
title(xlab = "Dose", ylab = "Probability", main = "Dose Response Models", cex.lab = 1)
# exponential model
curve(exponential(x,exponential_theta), 0.000001, max(d), add=TRUE, col=colors[2], lwd = 2.2)
# grid
grid(lwd = 1.2)
# legend
legend("bottomright", legend_names, fill = colors, title="Models & Data")
Beta-Poisson Best Fit Parameter
\[P(d; \alpha, N_{50}) = \left\{ \begin{array}{lr} 1 - \left(1 + \dfrac{d}{N_{50}} ( 2^{1/\alpha} - 1) \right)^{-\alpha} & : d \geq 0 \\ 0 & : d < 0 \end{array} \right.\]
where \[N_{50} = \dfrac{\ln(1/2)}{- \lambda}\] (\(\lambda\) is the parameter estimate for the exponential model)
Let \(\theta_{0} = (\alpha_{0}, \beta_{0})\) be the initial guess for \(\theta\), where
\[\alpha_{0} = 0.5 \text{ and }\beta_{0} = \dfrac{N_{50}}{(2^{1/\alpha_{0}} - 1)}\]
N50 = log(0.5)/-exponential_theta
alpha0 = 0.5
beta0 = N50/(2^(1/alpha0) - 1)
#theta0 = c(log(alpha0), log(beta0))
theta0 = c(alpha0, beta0)
mle_beta_poisson = nlminb(theta0, function(theta)
log_likelihood(theta, d, n, y, N, beta_poisson_approximation))
# mle_beta_poisson
# beta_poisson_theta = exp(beta_poisson_fit$par)
beta_poisson_theta = mle_beta_poisson$par
print(paste("Beta-Poisson parameter: ", beta_poisson_theta))## [1] "Beta-Poisson parameter: 3.1601532399525"
## [2] "Beta-Poisson parameter: 604.292025645697"Graph of Models
legend_names <- c("Data", "Exponential", "Beta-Poisson_Approx")plot(d, y / n, pch = 19, col = colors[1], lwd=4, ann=FALSE, ylim=c(0,1), log="x",
xlim=c(10^floor(log(min(d), 10)), 10^ceiling(log(max(d), 10))),
xaxp=c(10^floor(log(min(d), 10)), ceiling(log(max(d), 10)), 1),
yaxp=c(0,1,10), cex.axis = 1)
title(xlab = "Dose", ylab = "Probability", main = "Dose Response Models", cex.lab = 1)
# exponential model
curve(pexp(x,exponential_theta), 0.000001, max(d), add=TRUE, col=colors[2], lwd = 2.2)
# beta-poisson approximate model
curve(beta_poisson_approximation(x, beta_poisson_theta),
0.000001, max(d), add=TRUE, col=colors[3], lwd = 2.2)
# grid
grid(lwd = 1.2)
# legend
legend("bottomright", legend_names, fill = colors, title="Models & Data")