Initial set-up

# Load packages
library(tidyverse)
# library(biodosetools)
devtools::load_all()

# Additional options
theme_set(theme_light())

Dose estimation

fit_coeffs <- data.frame(
  estimate = c(0.0096, 0.0172, 0.0110),
  std.error = c(0.0012, 0.0032, 0.0010)
) %>%
  `rownames<-`(c("coeff_C", "coeff_alpha", "coeff_beta")) %>% 
  as.matrix()

fit_var_cov_mat <- matrix(
  0,
  nrow = nrow(fit_coeffs),
  ncol = nrow(fit_coeffs)
) %>%
  `colnames<-`(rownames(fit_coeffs)) %>%
  `rownames<-`(rownames(fit_coeffs))

for (x_var in rownames(fit_var_cov_mat)) {
  fit_var_cov_mat[[x_var, x_var]] <- fit_coeffs[x_var, "std.error"] * fit_coeffs[x_var, "std.error"]
}
case_data <- data.frame(
  C0 = 997, C1 = 3, C2 = 0, C3 = 0, C4 = 0, C5 = 0
) %>%
  dplyr::rename_with(
    .fn = toupper,
    .cols = dplyr::everything()
  ) %>%
  dplyr::mutate(
    dplyr::across(
      .cols = dplyr::starts_with("C"),
      .fns = as.integer
    )
  ) %>%
  dplyr::select(
    dplyr::starts_with("C")
  ) %>% 
  calculate_aberr_table(
    type = "case",
    assessment_u = 1
  ) %>%
  dplyr::rename(y = mean, y_err = std_err)

case_data
results_whole_merkle <- estimate_whole_body_merkle(
  case_data, 
  fit_coeffs,
  fit_var_cov_mat,
  protracted_g_value = 1,
  aberr_module = "dicentrics"
)
results_whole_delta <- estimate_whole_body_delta(
  case_data, 
  fit_coeffs,
  fit_var_cov_mat,
  conf_int = 0.95, 
  protracted_g_value = 1, 
  cov = TRUE, 
  aberr_module = "dicentrics"
)
case_data
fit_coeffs
fit_var_cov_mat
results_whole_merkle$est_doses
results_whole_delta$est_doses
plot_estimated_dose_curve(
  est_doses = list(whole = results_whole_delta),
  fit_coeffs,
  fit_var_cov_mat,
  protracted_g_value = 1,
  conf_int_curve = 0.95,
  aberr_name = "Dicentrics"
)

Check curve

# Generalised fit coefficients and variance-covariance matrix
general_fit_coeffs <- biodosetools:::generalise_fit_coeffs(fit_coeffs[, "estimate"])
general_fit_var_cov_mat <- biodosetools:::generalise_fit_var_cov_mat(fit_var_cov_mat)

coeff_C <- general_fit_coeffs[[1]]
coeff_alpha <- general_fit_coeffs[[2]]
coeff_beta <- general_fit_coeffs[[3]]
protracted_g_value <- 1

lambda_est <- case_data[["y"]]

lambda_est_corr <- correct_yield(lambda_est, "estimate", general_fit_coeffs, general_fit_var_cov_mat, conf_int = 0)
# lambda_est <- biodosetools:::correct_yield(lambda_est, "estimate", general_fit_coeffs, general_fit_var_cov_mat, conf_int = 0.95)
lambda_est <- 0.5
# Update coeff_beta to correct for protracted exposures
coeff_beta <- coeff_beta * protracted_g_value

# Auxiliary variable
z <- coeff_alpha^2 + 4 * coeff_beta * (lambda_est - coeff_C)

# Get estimate for dose
dose_est <- (-coeff_alpha + sqrt(z)) / (2 * coeff_beta)
dose_est
biodosetools:::project_yield(
  yield = 0,
  type = "estimate",
  general_fit_coeffs = general_fit_coeffs,
  general_fit_var_cov_mat = NULL,
  protracted_g_value = protracted_g_value,
  conf_int = 0
)
    dose_est <- project_yield(
      yield = lambda_est,
      type = "estimate",
      general_fit_coeffs = general_fit_coeffs,
      general_fit_var_cov_mat = NULL,
      protracted_g_value = protracted_g_value,
      conf_int = 0
    )


biodosimetry-uab/biodosetools documentation built on Jan. 26, 2024, 5:36 p.m.