Estimation de modèle d’absorption-élimination à un compartiment par maximum de vraisemblance

Le code peut être trouvé dans ce repo GitHub.

Paquêts

library(docstring)
library(tidyverse)

Lecture des données

fichier <- here("data", "PK.tsv")
data <- read_delim(fichier, delim = " ")
sujet11 <- data %>% filter(sujet == 11) %>% select(temps, concentration, dose)

head(sujet11)

Définition de fonctions

fonc.regr.dose <- function(D){
  #' Fixer la dose initiale 
  #' 
  #' @description Créer une fonction pour un modèle 
  #' d'absorption-élimination à 1 compartiment, en 
  #' fixant la dose initiale.
  #' 
  #' \deqn{Y_{i} = \frac{D}{V}\frac{k_a}{k_a - k_e}
  #' (e^{-k_e t_i} - e^{-k_a t_i}) + \varepsilon _i}
  #' Ce modèle est déterminé par quatre paramètres
  #'   D, V, ka, ke
  #' et la covariable explicative, le temps t.
  #' 
  #' Étant donné que la dose initiale D est connue,
  #' les paramètres sont réduits à :
  #'   V, ka, ke
  #'
  #' @param D réel. La dose initiale.
  #'
  #' @usage fonc.regr.dose(dose)
  #'
  #' @return Une fonction `concentration(t, V, ka, ke)`
  #'
  #' @note Le premier argument de la fonction "concentration"
  #' est la covariable `t` 
  #'
  
  concentration <- function(t, V, ka, ke){
    #stopifnot("ka == ke leads to zero-division"= ka != ke)
    (D/V) * (ka/(ka - ke)) * (exp(-ke * t) - exp(-ka * t))
  }
   
  return(concentration)
}

closest.to.mean <- function(df, group_var) {
  group_var <- enquo(group_var)

  df %>%
    filter(
      abs(!!group_var - mean(!!group_var)) == 
      min(abs(!!group_var - mean(!!group_var)))
    )
}

closest.to.median <- function(df, group_var) {
  group_var <- enquo(group_var)

  df %>%
    filter(
      abs(!!group_var - median(!!group_var)) == 
      min(abs(!!group_var - median(!!group_var)))
    )
}

less.more <-  function(a, b){
  c(a - b, a + b)
}

more.less <-  function(a, b){
  c(a + b, a - b)
}

Estimation des valeurs initiales

dose.by.subject <- data %>% group_by(sujet) %>% select(dose) %>%
  summarise( dose = unique(dose)) 

## Adding missing grouping variables: `sujet`

f.sujet11 <- fonc.regr.dose(unique(sujet11$dose))

v_i <- sujet11 %>%
   mutate(V = dose / concentration) %>% select(V) %>% min()

t.half_life <- sujet11 %>% closest.to.median(concentration) %>% 
  select(temps) %>% as.numeric()
ke_i <- log(2.0) / t.half_life

t.optim <- sujet11 %>% filter(concentration == max(concentration)) %>% 
  select(temps) %>% as.numeric()
dydt_0 <- (sujet11$concentration[1] - 0) / (sujet11$temps[1] - 0)
ka_i <- (dydt_0 * v_i ) / unique(sujet11$dose)

theta0 <- list(V=v_i, ka=ka_i, ke=ke_i)

Estimation du modèle

y.sujet11 <- nls(
  concentration ~ f.sujet11(temps, V, ka, ke), 
  data = sujet11, start = theta0,
  trace = TRUE, model = TRUE
)

## 34.35709 :  7.02857143 0.43428571 0.01925409
## 12.6584 :  7.38257696 0.73086183 0.01863309
## 9.639168 :  6.9315292 0.7554925 0.0203235
## 9.601533 :  6.90821105 0.73924867 0.02064228
## 9.600813 :  6.90675854 0.74071249 0.02067514
## 9.600797 :  6.90564574 0.74028554 0.02068265
## 9.600797 :  6.90563630 0.74032490 0.02068326
## 9.600797 :  6.90561174 0.74031476 0.02068342

s.y.sujet11 <- summary(y.sujet11)

plt1.sj11 <- sujet11 %>%
  mutate(
    `Concentration du médicament` = concentration,
    Temps = temps
  ) %>%
  ggplot(mapping = aes(x = Temps, y = `Concentration du médicament`)) +
    geom_point() +
    geom_function(aes(colour = "nls"), 
                  fun = f.sujet11, args = as.list(coef(y.sujet11))) +
    geom_function(aes(colour = "Valeurs Initiales"), 
                  fun = f.sujet11, args = theta0) +
  ggtitle(
 label = "Estimation de la concentration du médicament dans le plasma sanguin"
  )

print(plt1.sj11)

sujet11.var <- vcov(y.sujet11)

IC_V  <- less.more(coef(y.sujet11)[[1]], qnorm(0.975) * sqrt(sujet11.var[1,1])) 
IC_ka <- less.more(coef(y.sujet11)[[2]], qnorm(0.975) * sqrt(sujet11.var[2,2])) 
IC_ke <- less.more(coef(y.sujet11)[[3]], qnorm(0.975) * sqrt(sujet11.var[3,3])) 

y.coeffs <- tibble(as.data.frame(s.y.sujet11$coef)) %>%
  mutate(
    `Std. Estimate` = Estimate / `Std. Error`,
    `p(N)` = 2*(1 - pnorm(`Std. Estimate`)),
    param = rownames(as.data.frame(s.y.sujet11$coef))
  ) %>%
    select(param, everything())

y.coeffs