Copyright (C) 2018 Crista Moreno

Bootstrap Residuals 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/>.
source("dose_response_models.R")
source("extra_functions.R")
source("log_likelihood.R")

Data

data <- read.csv("./data/cryptosporidium_iowa_data.csv", header = TRUE)
d <- data$dose
n <- data$subjects
y <- data$positives

Residuals

\[\pi_i^{0} = \dfrac{y_i}{n_i}\]

\[\epsilon_i = \dfrac{\pi_i^{0} - \pi_i}{\sqrt{\pi_i(1 - \pi_i)/n_i}} \]

model_estimates <- function(N, d, model, theta_0) {

  pi <- numeric(N)
  
  for (i in 1:N) {
    pi_temp <- model(d[i], theta_0)
    pi[i] <- clamp(pi_temp)
  }
  return(pi)
}
epsilon_residuals <- function(N, n, pi_0, pi) {
  
  epsilon <- numeric(N)
  
  for (i in 1:N) {
    epsilon[i] <- (pi_0[i] - pi[i]) / sqrt(pi[i] * (1 - pi[i]) / n[i])
  }
  return(epsilon)
}
N <- length(d)
pi_0 <- y / n
theta_0 <- 0.00419109099824652
model <- exponential

pi <- model_estimates(N, d, model, theta_0)
epsilon <- epsilon_residuals(N, n, pi_0, pi)

\[\pi_i^{m} = \pi_i + \epsilon_q \sqrt{\pi(1 - \pi_i)/n_i}\]

\[y_i^{m} = \text{Bin}(\pi_i^{m}, n_i)\]

Bootstrap Replicates

bootstrap_replicates <- function(M, N, model, theta1, epsilon, pi, n) {
  
  theta_M <- numeric(M)
  
  for(m in 1:M) { 
    
    pi_m <- numeric(N)
    yi_m <- numeric(N)
    
    for(i in 1:N) {
      
      pi_m_temp <- pi[i] + sample(epsilon, 1)*sqrt(pi[i] * (1 - pi[i]) / n[i])
      pi_m[i] <- clamp(pi_m_temp)
      yi_m[i] <- rbinom(1, n[i], pi_m[i])
      
    }
    
    mle_exp_M = nlminb(theta1, function(theta)
      log_likelihood(theta, d, n, yi_m, N, model), lower = 0)
    
    theta_M[m] <- mle_exp_M$par
  }
  return(theta_M)
}
M <- 100
theta_1 <- 0.0000001
model <- exponential
theta_M <- bootstrap_replicates(M, N, model, theta_1, epsilon, pi, n)

Plot Bootstrap Models

colors <- c("chartreuse4","palevioletred1", "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)

# bootstrap exponential models
for(m in 1:M) {
  curve(exponential(x,theta_M[m]), 0.000001, max(d), add=TRUE, col=colors[2], lwd = 0.8)
}
points(d, y / n, pch = 19, col = colors[1], lwd=4)

# grid
grid(lwd = 1.2)

# legend
legend("bottomright", legend_names, fill = colors, title="Models & Data")