20250317.R

##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 = 17¹
Ratio AUC avec inhibiteur 3A puissants	2.7 [1.1 - 24.0]
seuilpredictif	8 (47.1%)
Ratio AUC avec inhibiteur 3A modérés	2.00 [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%)
¹ 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 = 8¹	1 N = 9¹
Ratio AUC avec inhibiteur 3A puissants	2.0 [1.1 - 2.7]	6.5 [1.6 - 24.0]
seuilpredictif	0 (0.0%)	8 (88.9%)
Ratio AUC avec inhibiteur 3A modérés	1.55 [1.03 - 1.83]	3.80 [2.00 - 5.00]
Facteur de biodisponibilité absolue
Median [Min - Max]	0.26 [0.05 - 0.46]	0.12 [0.04 - 0.46]
Unknown	6	4
Médicaments substrats de transporteurs (P-gp ou BCRP)	6 (75.0%)	6 (66.7%)
Médicaments substrats d'autres isoenzyme que cyp3A	2 (25.0%)	2 (22.2%)
¹ 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: 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)       -5.552      2.858  -1.943   0.0521 .
## ratio_auc_inh_p    1.919      1.101   1.743   0.0813 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23.5081  on 16  degrees of freedom
## Residual deviance:  8.9845  on 15  degrees of freedom
## AIC: 12.984
## 
## Number of Fisher Scoring iterations: 8

exp(coefficients(mod))

##     (Intercept) ratio_auc_inh_p 
##     0.003878693     6.816051906

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: 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: 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: 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: 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: 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

##                        2.5 %      97.5 %
## (Intercept)     1.285809e-06   0.2140745
## ratio_auc_inh_p 1.661607e+00 170.0611030

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: 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: 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: 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: 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: 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

Characteristic	OR	95% CI	p-value
(Intercept)	0.00	0.00, 0.21	0.052
Ratio AUC avec inhibiteur 3A puissants	6.82	1.66, 170	0.081
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.1111111

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

##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 = 13.432, df = 1, p-value = 0.0002474

##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 = 9¹	1 N = 8¹
seuilinibition_2	1 (11.1%)	8 (100.0%)
¹ n (%)

tbl_summary(
  databuilding, include = c("seuilpredictif"),
  by="seuilinibition_2",
  digits=all_categorical()~ c(0,1)
)

Characteristic	0 N = 8¹	1 N = 9¹
seuilpredictif	0 (0.0%)	8 (88.9%)
¹ 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 = 7¹	1 N = 4¹
seuilinibition_2	0 (0.0%)	4 (100.0%)
¹ n (%)

tbl_summary(
  validation, include = c("seuilpredictif"),
  by="seuilinibition_2",
  digits=all_categorical()~ c(0,1)
)

Characteristic	0 N = 7¹	1 N = 4¹
seuilpredictif	0 (0.0%)	4 (100.0%)
¹ 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: 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: 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: 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: 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: 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: 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	9.86	2.52, 140	0.017
Facteur de biodisponibilité absolue	13	0.01	0.00, 0.74	0.10
Médicaments substrats de transporteurs (P-gp ou BCRP)	28	1.07	0.23, 5.04	>0.9
Médicaments substrats d'autres isoenzyme que cyp3A	28	0.45	0.08, 2.25	0.3
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

20250317.R

dujols_c

2025-03-20