##chargement des packages----
library(questionr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tableone)
library(labelled)
library(gtsummary)
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(readxl)
library(effects)
## Le chargement a nécessité le package : carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(survival)
library(survminer)
## Le chargement a nécessité le package : ggpubr
## 
## Attachement du package : 'survminer'
## 
## L'objet suivant est masqué depuis 'package:survival':
## 
##     myeloma
library(ggplot2)
library(dplyr)
library(knitr)
library(cowplot)
## 
## Attachement du package : 'cowplot'
## 
## L'objet suivant est masqué depuis 'package:ggpubr':
## 
##     get_legend
## 
## L'objet suivant est masqué depuis 'package:lubridate':
## 
##     stamp
##chargement des données 

install.packages("readxl")
## Warning: le package 'readxl' est en cours d'utilisation et ne sera pas installé
library("readxl")
library(devtools)
## Le chargement a nécessité le package : usethis
library(openxlsx)
library(dplyr)

data <- read_excel("E:/Internes/INTERNES NOV 2024 - MAI 2025/Camille/sfpo/prediction_cyp3a_inh_m/analyses stats/20250317/datacamilleef.xlsx")
## New names:
## • `biblio` -> `biblio...32`
## • `biblio` -> `biblio...33`
## • `` -> `...56`
## • `ratio_AUC` -> `ratio_AUC...60`
## • `ratio_AUC` -> `ratio_AUC...62`
## • `ratio_AUC` -> `ratio_AUC...64`
## • `ratio_AUC` -> `ratio_AUC...66`
install.packages("devtools")
## Warning: le package 'devtools' est en cours d'utilisation et ne sera pas
## installé
data_im<-filter(data, c(eligible=="ok"))
databuilding<-filter(data_im, c(population=="building"))
validation<-filter(data_im, c(population=="validation"))

##recodage des variables et bases de données le cas échéant----
databuilding$seuilinibition_1_5<-ifelse(databuilding$ratio_auc_inh_m>=1.5, 1, 0)
databuilding$seuilinibition_2<-ifelse(databuilding$ratio_auc_inh_m>=2, 1, 0)
databuilding$seuilpredictif<-ifelse(databuilding$ratio_auc_inh_p>=3.25, 1, 0)

validation$seuilinibition_1_5<-ifelse(validation$ratio_auc_inh_m>=1.5, 1, 0)
validation$seuilinibition_2<-ifelse(validation$ratio_auc_inh_m>=2, 1, 0)
validation$seuilpredictif<-ifelse(validation$ratio_auc_inh_p>=3.25, 1, 0)

data_im$seuilinibition_2<-ifelse(data_im$ratio_auc_inh_m>=2, 1, 0)

##renommer des variables pour présentation dans les tableaux de résultats
library(labelled)
var_label(data_im$ratio_auc_inh_p) <- "Ratio AUC avec inhibiteur 3A puissants"
var_label(data_im$ratio_auc_inh_m) <- "Ratio AUC avec inhibiteur 3A modérés"
var_label(data_im$biodisp) <- "Facteur de biodisponibilité absolue"
var_label(data_im$substrat_transporteur) <- "Médicaments substrats de transporteurs (P-gp ou BCRP)"
var_label(data_im$substrat_autre_cyp) <- "Médicaments substrats d'autres isoenzyme que cyp3A"

var_label(databuilding$ratio_auc_inh_p) <- "Ratio AUC avec inhibiteur 3A puissants"
var_label(databuilding$ratio_auc_inh_m) <- "Ratio AUC avec inhibiteur 3A modérés"
var_label(databuilding$biodisp) <- "Facteur de biodisponibilité absolue"
var_label(databuilding$substrat_transporteur) <- "Médicaments substrats de transporteurs (P-gp ou BCRP)"
var_label(databuilding$substrat_autre_cyp) <- "Médicaments substrats d'autres isoenzyme que cyp3A"

var_label(validation$ratio_auc_inh_p) <- "Ratio AUC avec inhibiteur 3A puissants"
var_label(validation$ratio_auc_inh_m) <- "Ratio AUC avec inhibiteur 3A modérés"
var_label(validation$biodisp) <- "Facteur de biodisponibilité absolue"
var_label(validation$substrat_transporteur) <- "Médicaments substrats de transporteurs (P-gp ou BCRP)"
var_label(validation$substrat_autre_cyp) <- "Médicaments substrats d'autres isoenzyme que cyp3A"

##tableau descriptif population globales construction modèle ----
tbl_summary(
  databuilding, include = c("ratio_auc_inh_p", "seuilpredictif", "ratio_auc_inh_m", "biodisp", "substrat_transporteur","substrat_autre_cyp"),
  type = c(biodisp) ~ "continuous2",
  statistic = all_continuous() ~ "{median} [{min} - {max}]",
  digits=all_categorical()~ c(0,1)
)
Characteristic N = 171
Ratio AUC avec inhibiteur 3A puissants 2.7 [1.1 - 24.0]
seuilpredictif 8 (47.1%)
Ratio AUC avec inhibiteur 3A modérés 1.83 [1.03 - 5.00]
Facteur de biodisponibilité absolue
    Median [Min - Max] 0.12 [0.04 - 0.46]
    Unknown 10
Médicaments substrats de transporteurs (P-gp ou BCRP) 12 (70.6%)
Médicaments substrats d'autres isoenzyme que cyp3A 4 (23.5%)
1 Median [Min - Max]; n (%) 
##tableau descriptif population construction modele par ratio auc modéré ----
tbl_summary(
  databuilding, include = c("ratio_auc_inh_p", "seuilpredictif", "ratio_auc_inh_m", "biodisp", "substrat_transporteur","substrat_autre_cyp"),
  statistic = all_continuous() ~ "{median} [{min} - {max}]",
  type = c(biodisp) ~ "continuous2",
  by="seuilinibition_2",
  digits=all_categorical()~ c(0,1)
)
Characteristic 0
N = 91		 1
N = 81
Ratio AUC avec inhibiteur 3A puissants 1.9 [1.1 - 2.7] 6.6 [3.5 - 24.0]
seuilpredictif 0 (0.0%) 8 (100.0%)
Ratio AUC avec inhibiteur 3A modérés 1.50 [1.03 - 1.83] 3.99 [2.00 - 5.00]
Facteur de biodisponibilité absolue

    Median [Min - Max] 0.12 [0.05 - 0.46] 0.22 [0.04 - 0.46]
    Unknown 6 4
Médicaments substrats de transporteurs (P-gp ou BCRP) 6 (66.7%) 6 (75.0%)
Médicaments substrats d'autres isoenzyme que cyp3A 3 (33.3%) 1 (12.5%)
1 Median [Min - Max]; n (%) 
##analyses univariée prédiction ratio modéré >2 par ratio puissant en valeur continue
mod<-glm(  seuilinibition_2 ~ ratio_auc_inh_p, databuilding, family="binomial")
## Warning: glm.fit : l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
summary(mod)
## 
## Call:
## glm(formula = seuilinibition_2 ~ ratio_auc_inh_p, family = "binomial", 
##     data = databuilding)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)      -151.80   87094.55  -0.002    0.999
## ratio_auc_inh_p    49.23   28272.78   0.002    0.999
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2.3508e+01  on 16  degrees of freedom
## Residual deviance: 5.3772e-09  on 15  degrees of freedom
## AIC: 4
## 
## Number of Fisher Scoring iterations: 25
exp(coefficients(mod))
##     (Intercept) ratio_auc_inh_p 
##    1.190414e-66    2.400585e+21
exp(confint(mod, level=0.95))
## Waiting for profiling to be done...
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit : l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit : l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
##                 2.5 % 97.5 %
## (Intercept)         0    Inf
## ratio_auc_inh_p     0    Inf
mod%>%tbl_regression(intercept = TRUE, exponentiate = TRUE)
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit : l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit : l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Characteristic OR 95% CI p-value
(Intercept) 0.00 0.00, Inf >0.9
Ratio AUC avec inhibiteur 3A puissants 2,400,584,838,719,154,946,088 0.00, Inf >0.9
Abbreviations: CI = Confidence Interval, OR = Odds Ratio 
##estimation de la valeur seuil de Ratio inh puissant

###courbe ROC pour ratio auc p et seuil ratio modéré à 2
library(ROCR)


#Calcul IC95
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attachement du package : 'pROC'
## 
## Les objets suivants sont masqués depuis 'package:stats':
## 
##     cov, smooth, var
roc1<-roc(seuilinibition_2~ratio_auc_inh_p, databuilding)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
coords(roc1, "best", ret=c("threshold", "specificity", "1-npv"))
##           threshold specificity 1-npv
## threshold     3.085           1     0
roc2 <- roc(databuilding$seuilinibition_2,
            databuilding$ratio_auc_inh_p, percent=TRUE,
            # arguments for auc
            partial.auc=c(100, 0), partial.auc.correct=TRUE,
            partial.auc.focus="sens",
            # arguments for ci
            ci=TRUE, boot.n=100, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
            print.auc=TRUE, show.thres=TRUE)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Warning in ci.auc.roc(roc, ...): ci.auc() of a ROC curve with AUC == 1 is
## always 1-1 and can be misleading.

##analyses univariée prédiction ratio modéré >2 par ratio puissant>=3.25 par chi2

chisq.test(databuilding$seuilinibition_2, databuilding$seuilpredictif, correct=FALSE)
## Warning in chisq.test(databuilding$seuilinibition_2,
## databuilding$seuilpredictif, : L’approximation du Chi-2 est peut-être
## incorrecte
## 
##  Pearson's Chi-squared test
## 
## data:  databuilding$seuilinibition_2 and databuilding$seuilpredictif
## X-squared = 17, df = 1, p-value = 3.738e-05
##Estimation VPP et VPN sur population building

tbl_summary(
  databuilding, include = c("seuilinibition_2"),
  by="seuilpredictif",
  digits=all_categorical()~ c(0,1)
)
Characteristic 0
N = 91		 1
N = 81
seuilinibition_2 0 (0.0%) 8 (100.0%)
1 n (%) 
tbl_summary(
  databuilding, include = c("seuilpredictif"),
  by="seuilinibition_2",
  digits=all_categorical()~ c(0,1)
)
Characteristic 0
N = 91		 1
N = 81
seuilpredictif 0 (0.0%) 8 (100.0%)
1 n (%) 
##Estimation VPP et VPN sur population de validation

tbl_summary(
  validation, include = c("seuilinibition_2"),
  by="seuilpredictif",
  digits=all_categorical()~ c(0,1)
)
Characteristic 0
N = 71		 1
N = 41
seuilinibition_2 0 (0.0%) 4 (100.0%)
1 n (%) 
tbl_summary(
  validation, include = c("seuilpredictif"),
  by="seuilinibition_2",
  digits=all_categorical()~ c(0,1)
)
Characteristic 0
N = 71		 1
N = 41
seuilpredictif 0 (0.0%) 4 (100.0%)
1 n (%) 
###graphiques de correlations

g1<-ggplot(data_im) +
  aes(x = ratio_auc_inh_p, y = ratio_auc_inh_m) +
  geom_smooth(method = "lm") +
  geom_point() +
  xlab("Ratio AUC avec inhibiteur 3A puissant") +
  ylab("Ratio AUC avec inhibiteur 3A modérés")

g1
## `geom_smooth()` using formula = 'y ~ x'

##analyses multivariée prédiction ratio modéré >2 autres facteurs sur population globale
data_im |>
  tbl_uvregression(
    y = seuilinibition_2,
    include = c(ratio_auc_inh_p, biodisp, substrat_transporteur, substrat_autre_cyp),
    method = glm,
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) |> 
  bold_labels()
## There was a warning constructing the model for variable "ratio_auc_inh_p". See
## message below.
## ! glm.fit : l'algorithme n'a pas convergé and glm.fit: des probabilités ont été
##   ajustées numériquement à 0 ou 1
## There was a warning running `tbl_regression()` for variable "ratio_auc_inh_p".
## See message below.
## ! glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1, glm.fit:
##   des probabilités ont été ajustées numériquement à 0 ou 1, glm.fit: des
##   probabilités ont été ajustées numériquement à 0 ou 1, glm.fit : l'algorithme
##   n'a pas convergé, glm.fit: des probabilités ont été ajustées numériquement à
##   0 ou 1, glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1,
##   and glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
Characteristic N OR 95% CI p-value
Ratio AUC avec inhibiteur 3A puissants 28 5,765,540,774,768,967,816,444,860,060 Inf, Inf >0.9
Facteur de biodisponibilité absolue 13 0.03 0.00, 2.22 0.2
Médicaments substrats de transporteurs (P-gp ou BCRP) 28 1.56 0.33, 7.88 0.6
Médicaments substrats d'autres isoenzyme que cyp3A 28 0.26 0.03, 1.40 0.14
Abbreviations: CI = Confidence Interval, OR = Odds Ratio