---
title: "Motorette"
author: "La dream team MathApp"
format: html
editor: visual
---

On va analyser les données Motorette qui sont présentes dans le package "SMPracticals" (à installer donc).

```{r}
library(SMPracticals)
data(motorette)
head(motorette)
```

## Analyse descriptive

Je vous booste un peu la culture générale de stat. Ici on va voir Kaplan-Meier.

```{r}
library(survival)##le graal pour l'analyse de survie
library(survminer)##pour des graphiques plus jolis
fit <- survfit(Surv(y, cens) ~ x, data = motorette)
ggsurvplot(fit, conf.int = 0.95)
```

## Code Gibbs

```{r}
gibbs_motorette <- function(nchain, data, prop_sd, init){
  ## Gibbs pour exo 9
  
  chain <- matrix(NA, nchain + 1, 3)
  colnames(chain) <- c("beta0", "beta1", "shape")
  acc_rate <- rep(0, 3)
  chain[1,] <- init
  durations <- data$y
  covar <- cbind(1, log(data$x / 100))
  is_censored <- data$cens == 0
  
  for (iter in 1:nchain){
    ##Etat initial
    current <- chain[iter,]
    scale_current <- exp(covar %*% current[1:2]) 
    
    ## Mise a jour des beta
    for (j in 1:2){
      prop <- current[1:2]
      prop[j] <- rnorm(1, prop[j], prop_sd[j])
      
      scale_prop <- exp(covar %*% prop)
      
      top <- sum(dweibull(durations[!is_censored],
                          current["shape"],
                          scale = scale_prop[!is_censored],
                          log = TRUE)) +
        sum(pweibull(durations[is_censored],
                     current["shape"],
                     scale = scale_prop[is_censored],
                     lower.tail = FALSE, log.p = TRUE)) -
        0.5 * prop[j]^2 / 100
      
      bottom <- sum(dweibull(durations[!is_censored],
                          current["shape"],
                          scale = scale_current[!is_censored],
                          log = TRUE)) +
        sum(pweibull(durations[is_censored],
                     current["shape"],
                     scale = scale_current[is_censored],
                     lower.tail = FALSE, log.p = TRUE)) -
        0.5 * prop[j]^2 / 100
      
      acc_prob <- exp(top - bottom)
      
      if (is.na(acc_prob))
        browser()
      
      if (runif(1) < acc_prob){
        current[1:2] <- prop
        scale_current <- scale_prop
        acc_rate[j] <- acc_rate[j] + 1
      }
    }
    
    ## Mise à jour du paramètre de forme
    prop <- rlnorm(1, log(current["shape"]), prop_sd[3])
    
    top <- sum(dweibull(durations[!is_censored],
                        prop,
                          scale = scale_current[!is_censored],
                          log = TRUE)) +
        sum(pweibull(durations[is_censored],
                     prop,
                     scale = scale_current[is_censored],
                     lower.tail = FALSE, log.p = TRUE)) -
        prop / 2
      
      bottom <- sum(dweibull(durations[!is_censored],
                          current["shape"],
                          scale = scale_current[!is_censored],
                          log = TRUE)) +
        sum(pweibull(durations[is_censored],
                     current["shape"],
                     scale = scale_current[is_censored],
                     lower.tail = FALSE, log.p = TRUE)) -
        current["shape"] / 2
      
      acc_prob <- exp(top - bottom) * prop / current["shape"]
      
      if (is.na(acc_prob))
        browser()
      
      if (runif(1) < acc_prob){
        current["shape"] <- prop
        acc_rate[3] <- acc_rate[3] + 1
      }
      
      ## Mise à jour de la chaine
      chain[iter+1,] <- current
  }
  
  return(list(chain = chain, acc_rate = acc_rate / nchain))
  
}
```

Lançons notre application

```{r}
nchain <- 100000
prop_sd <- c(0.5, 0.6, 0.5)
init <- c(-log(4000), 0, 1)
ans <- gibbs_motorette(nchain, motorette, prop_sd, init)
```

Regardons la chaine

```{r}
library(coda)##utile pour analyse chaine issues MCMC
burnin <- 2000##burnin period
thin <- 5## élagage
chain <- ans$chain[seq(burnin, nchain, by = thin),]
plot(mcmc(chain[-(1:2000),]))
```

On va vérifier si le modèle est bon ou pas. Pour cela, on fait une estimation ponctuelle (Bayésienne) et l'on compare à Kaplan-Meier.

```{r}
post_mean <- colMeans(chain)
post_mean

## On peut se faire une fonction qui calcule la survie
## pour une température donnée
survie_pred <- function(duree){
  scale <- exp(post_mean["beta0"] + post_mean["beta1"] *
                 log(temp/100))
  pweibull(duree, post_mean["shape"], scale,
           lower.tail = FALSE)
}

cols <- RColorBrewer::brewer.pal(4, "Set2")##daltonien friendly
plot(fit, col = cols)##fit est estimation Kaplan-Meier

i <- 1
for (temp in unique(motorette$x)){
  curve(survie_pred, from = 0, to = 10000,
        col = cols[i], add = TRUE, lwd = 2, lty = 2)
  i <- i + 1
}
```

Ici notre modèle paramétrique a l'air plutôt raisonnable et je peux donc me permettre une folie : extrapoler pour la température 130 degrés.

```{r}
temp <- 130
plot(survie_pred, from = 0, to = 100000)
```

Si mon intérêt était sur du contrôle qualité via par exemple le quantile d'ordre 0.95 pour la température 130, je pourrais obtenir sa loi a posteriori

```{r}
shape <- chain[,"shape"]
scale <- exp(chain[,"beta0"] + chain[,"beta1"] * log(130/100))

pred <- qweibull(0.95, shape, scale)
hist(pred, freq= FALSE)
```
