--- title: "Simulation de donnees a baseline et de l'outcome" author: "Saint Pierre - Savy" date: "1/07/2022" output: html_document: default pdf_document: default word_document: default --- ```{r,include=FALSE} #knitr::opts_chunk$set(fig.width=10, fig.height=7) ``` ```{r,include=FALSE} library(rvinecopulib) # library(VineCopula) # library(ggraph) library(e1071) library(GGally) library(caret) library(MASS) library(tidyverse) library(corrr) library(lsr) library(cowplot) library(EnvStats) library(ggraph) #diabetes <- read.csv("C:/Users/nicolas.savy/Dropbox/PROJETS#RECHERCHE/7#BookChapter/Simulations/diabetes.csv") data <- read.csv("diabetes.csv") dim(data) # nettoyage des donn??es : les 0 sont des NA pour Glucose, BloodPressure,SkinThickness,Insulin,BMI diabetes<-data[data$Insulin!=0 & data$Glucose!=0 & data$BMI!=0,] summary(diabetes) dim(diabetes) # split Virtual Patient (VP) and virtual Outcome (VO) #separation var discrete et continue VP <- diabetes[,c(1:8)] VPDisc <- diabetes[,c(1,6)] VPDisc$Pregnancies = cut(VPDisc$Pregnancies,c(-0.1,0.9,100)) levels(VPDisc$Pregnancies)<-c('No','Yes') VPDisc$BMI = cut(VPDisc$BMI,c(0,25,35,200)) levels(VPDisc$BMI)<-c('<=25','25-35','>35') VPCont <- diabetes[,-c(1,6,9)] VPMod =as.data.frame(c(VPDisc,VPCont)) VO <- as.data.frame(as.factor(diabetes[,9])) names(VO) = c('Pronostic') Data_Originale<-cbind(VPMod,VO) ``` ## Description des variables explicatives et de l'outcome ```{r,echo=F,message = FALSE,fig.width=10, fig.height=10} # ###################### Estimation correlation network mixed data ##################### # #### Fonction trouvee sur internet # ############################################################################### # # # # cor2 # # -- Compute correlations of columns of a dataframe of mixed types # # # ############################################################################### # # # # author : Srikanth KS (talegari) # # license : GNU AGPLv3 (http://choosealicense.com/licenses/agpl-3.0/) # # # ############################################################################### # #' @title cor2 # #' # #' @description Compute correlations of columns of a dataframe of mixed types. # #' The dataframe is allowed to have columns of these four classes: integer, # #' numeric, factor and character. The character column is considered as # #' categorical variable. # #' # #' @details The correlation is computed as follows: \itemize{ # #' # #' \item integer/numeric pair: pearson correlation using `cor` function. The # #' valuelies between -1 and 1. # #' # #' \item integer/numeric - factor/categorical pair: correlation coefficient or # #' squared root of R^2 coefficient of linear regression of integer/numeric # #' variable over factor/categorical variable using `lm` function. The value # #' lies between 0 and 1. \item factor/categorical pair: cramersV value is # #' computed based on chisq test using `lsr::cramersV` function. The value lies # #' between 0 and 1. # #' } # #' # #' For a comprehensive implementation, use `polycor::hetcor` # #' # #' @param df input data frame # #' # #' @author Srikanth KS(gmail to sri dot teach) # #' # #' @keywords GNU AGPLv3 (http://choosealicense.com/licenses/agpl-3.0/) # #' # #' @examples # #' iris_cor <- cor2(iris) # #' corrplot::corrplot(iris_cor) # #' corrgram::corrgram(iris_cor) # #' # cor2 = function(df){ # # stopifnot(inherits(df, "data.frame")) # stopifnot(sapply(df, class) %in% c("integer" # , "numeric" # , "factor" # , "character")) # # cor_fun <- function(pos_1, pos_2){ # # # both are numeric # if(class(df[[pos_1]]) %in% c("integer", "numeric") && # class(df[[pos_2]]) %in% c("integer", "numeric")){ # r <- stats::cor(df[[pos_1]] # , df[[pos_2]] # , use = "pairwise.complete.obs" # ) # } # # # one is numeric and other is a factor/character # if(class(df[[pos_1]]) %in% c("integer", "numeric") && # class(df[[pos_2]]) %in% c("factor", "character")){ # r <- sqrt( # summary( # stats::lm(df[[pos_1]] ~ as.factor(df[[pos_2]])))[["r.squared"]]) # } # # if(class(df[[pos_2]]) %in% c("integer", "numeric") && # class(df[[pos_1]]) %in% c("factor", "character")){ # r <- sqrt( # summary( # stats::lm(df[[pos_2]] ~ as.factor(df[[pos_1]])))[["r.squared"]]) # } # # # both are factor/character # if(class(df[[pos_1]]) %in% c("factor", "character") && # class(df[[pos_2]]) %in% c("factor", "character")){ # r <- lsr::cramersV(df[[pos_1]], df[[pos_2]], simulate.p.value = TRUE) # } # # return(r) # } # # cor_fun <- Vectorize(cor_fun) # # # now compute corr matrix # corrmat <- outer(1:ncol(df) # , 1:ncol(df) # , function(x, y) cor_fun(x, y) # ) # # rownames(corrmat) <- colnames(df) # colnames(corrmat) <- colnames(df) # # return(corrmat) # } # # ### calcul des correlation grace ?? la fonction # res.cor <- cor2(VPMod) # res.cor # # ### representation du r??seau des correlations # network_plot(res.cor,min_cor = 0.2) # box and whisker plots for all attributes #ggpairs(VPMod) p<-ggpairs(Data_Originale, aes(colour=Pronostic), lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1))) for(i in 1:p$nrow) { for(j in 1:p$ncol){ p[i,j] <- p[i,j] + scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) + scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) } } #p ``` ## Taille de l'echantillon a simuler et nombre de simulations ```{r,echo=F} ## Tailel echantillon simule #################################################################### N<-500 cat("Taille de l'echantillon simule :", N,"\n") ## Nombre de simulations (nbre de p-value calculees) #################################################################### Nsimu<-100 cat("Nombre de simulations :", Nsimu,"\n") ``` ## Simulation des donnees a baseline avec la methode continue ```{r,echo=F} #### Fonction pour simuler avec la m??thode continue ######### Simulation_Continue<- function(DataBaseline) { ###################### LEARNING ############################### MC = NULL VC = NULL VOriC = DataBaseline VOriC[,1] = as.numeric(VOriC[,1]) VOriC[,2] = as.numeric(VOriC[,2]) MC = colMeans(VOriC) VC = var(VOriC) ##################### GENERATION of a SAMPLE of size N ################################# VGenC = mvrnorm(N,MC,VC) VGenC = as.data.frame(VGenC) # Discretisation de Pregnancies TAB1 = table(VPMod[,1]) P1 = prop.table(TAB1) CrV1 = MC[1]+sqrt(VC[1,1])*qnorm(P1[1],0,1) VGenC[,1]=as.factor(cut(VGenC[,1],c(min(VGenC[,1]),CrV1,max(VGenC[,1])),include.lowest=T)) # Discretisation de BMI TAB2 = table(VPMod[,2]) P2 = prop.table(TAB2) CrV21 = MC[2]+sqrt(VC[2,2])*qnorm(P2[1],0,1) CrV22 = MC[2]+sqrt(VC[2,2])*qnorm(P2[1]+P2[2],0,1) VGenC[,2]=as.factor(cut(VGenC[,2],c(min(VGenC[,2]),CrV21,CrV22,max(VGenC[,2])),include.lowest=T)) levels(VGenC$Pregnancies)<-c('No','Yes') levels(VGenC$BMI)<-c('<=25','25-35','>35') return(VGenC) } ``` ## Simulation des donnees a baseline avec la methode discrete ```{r,echo=F,message=F} #### Fonction pour simuler avec la m??thode discrete ######### Simulation_Discrete<- function(DataBaseline) { ###################### LEARNING ############################### VPDisc<-DataBaseline[,1:2] VPCont<-DataBaseline[,3:8] P = length(VPDisc$Pregnancies) L = NULL L[[1]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[1]) & (VPDisc$BMI == levels(VPDisc$BMI)[1])) L[[2]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[2]) & (VPDisc$BMI == levels(VPDisc$BMI)[1])) L[[3]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[1]) & (VPDisc$BMI == levels(VPDisc$BMI)[2])) L[[4]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[2]) & (VPDisc$BMI == levels(VPDisc$BMI)[2])) L[[5]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[1]) & (VPDisc$BMI == levels(VPDisc$BMI)[3])) L[[6]] = which((VPDisc$Pregnancies == levels(VPDisc$Pregnancies)[2]) & (VPDisc$BMI == levels(VPDisc$BMI)[3])) proba = NULL for (i in 1:6){ proba = c(proba,length(L[[i]])/P) } M = NULL V = NULL for (i in 1:6){ M[[i]] = colMeans(VPCont[L[[i]],]) V[[i]] = var(VPCont[L[[i]],]) } ##################### GENERATION of a SAMPLE of size N ################################# TEMP = sample(c("1","2","3","4","5","6"), N, replace=TRUE, prob = proba) VGen = matrix(0,ncol = 8, nrow = N) for (i in 1:N){ if ( TEMP[i] == "1" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[1] VGen[i,2] = levels(VPDisc$BMI)[1] VGen[i,3] = mvrnorm(1,M[[1]],V[[1]])[1] VGen[i,4] = mvrnorm(1,M[[1]],V[[1]])[2] VGen[i,5] = mvrnorm(1,M[[1]],V[[1]])[3] VGen[i,6] = mvrnorm(1,M[[1]],V[[1]])[4] VGen[i,7] = mvrnorm(1,M[[1]],V[[1]])[5] VGen[i,8] = mvrnorm(1,M[[1]],V[[1]])[6]} if ( TEMP[i] == "2" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[2] VGen[i,2] = levels(VPDisc$BMI)[1] VGen[i,3] = mvrnorm(1,M[[2]],V[[2]])[1] VGen[i,4] = mvrnorm(1,M[[2]],V[[2]])[2] VGen[i,5] = mvrnorm(1,M[[2]],V[[2]])[3] VGen[i,6] = mvrnorm(1,M[[2]],V[[2]])[4] VGen[i,7] = mvrnorm(1,M[[2]],V[[2]])[5] VGen[i,8] = mvrnorm(1,M[[2]],V[[2]])[6]} if ( TEMP[i] == "3" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[1] VGen[i,2] = levels(VPDisc$BMI)[2] VGen[i,3] = mvrnorm(1,M[[3]],V[[4]])[1] VGen[i,4] = mvrnorm(1,M[[3]],V[[4]])[2] VGen[i,5] = mvrnorm(1,M[[3]],V[[4]])[3] VGen[i,6] = mvrnorm(1,M[[3]],V[[4]])[4] VGen[i,7] = mvrnorm(1,M[[3]],V[[4]])[5] VGen[i,8] = mvrnorm(1,M[[3]],V[[4]])[6]} if ( TEMP[i] == "4" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[2] VGen[i,2] = levels(VPDisc$BMI)[2] VGen[i,3] = mvrnorm(1,M[[5]],V[[5]])[1] VGen[i,4] = mvrnorm(1,M[[5]],V[[5]])[2] VGen[i,5] = mvrnorm(1,M[[5]],V[[5]])[3] VGen[i,6] = mvrnorm(1,M[[5]],V[[5]])[4] VGen[i,7] = mvrnorm(1,M[[5]],V[[5]])[5] VGen[i,8] = mvrnorm(1,M[[5]],V[[5]])[6]} if ( TEMP[i] == "5" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[1] VGen[i,2] = levels(VPDisc$BMI)[3] VGen[i,3] = mvrnorm(1,M[[6]],V[[6]])[1] VGen[i,4] = mvrnorm(1,M[[6]],V[[6]])[2] VGen[i,5] = mvrnorm(1,M[[6]],V[[6]])[3] VGen[i,6] = mvrnorm(1,M[[6]],V[[6]])[4] VGen[i,7] = mvrnorm(1,M[[6]],V[[6]])[5] VGen[i,8] = mvrnorm(1,M[[6]],V[[6]])[6]} if ( TEMP[i] == "6" ) { VGen[i,1] = levels(VPDisc$Pregnancies)[2] VGen[i,2] = levels(VPDisc$BMI)[3] VGen[i,3] = mvrnorm(1,M[[1]],V[[1]])[1] VGen[i,4] = mvrnorm(1,M[[1]],V[[1]])[2] VGen[i,5] = mvrnorm(1,M[[1]],V[[1]])[3] VGen[i,6] = mvrnorm(1,M[[1]],V[[1]])[4] VGen[i,7] = mvrnorm(1,M[[1]],V[[1]])[5] VGen[i,8] = mvrnorm(1,M[[1]],V[[1]])[6]} } VGen = as.data.frame(VGen) for (i in 3:8){VGen[,i] = as.numeric(as.character(VGen[,i]))} colnames(VGen)<-colnames(DataBaseline) return(VGen) } ``` ## Simulation des donnees a baseline avec la methode Vine copula (structure data-driven) ```{r,echo=F,message=F} #### Fonction pour simuler avec la m??thode des r-vine copula ######### Simulation_Copule<- function(DataBaseline) { ###################### ESTIMATION ###################################### #Pseudo inverse variables continues U_cont<-pseudo_obs(DataBaseline[,3:8]) #Distribution des 2 variables discretes disc_1<-as.integer(DataBaseline[,1])-1 disc_2<-as.integer(DataBaseline[,2])-1 freq_disc1<-prop.table(table(DataBaseline[,1])) freq_disc2<-prop.table(table(DataBaseline[,2])) Freq_disc_t1 <- cbind(pdiscrete(disc_1+1,freq_disc1),pdiscrete(disc_2+1,freq_disc2)) Freq_disc_t0 <- cbind(pdiscrete(disc_1,freq_disc1),pdiscrete(disc_2,freq_disc2)) #densite :ddiscrete(x+1,freq) #distribution function : pdiscrete(x+1,freq) #quantile function : qdiscrete(u[, 1],freq) - 1 # Infos necessaire pour copule avec donnees mixte U_mixte<-cbind(Freq_disc_t1,U_cont,Freq_disc_t0) ## Model for mixed data fit_DataDriven <- vinecop(U_mixte, var_types = c("d","d", rep("c", 6))) #summary(fit_DataDriven) #plot(fit_DataDriven) #contour(fit_DataDriven) # definition de la distribution de la r-vine Fit_dist<-vinecop_dist(fit_DataDriven$pair_copulas,fit_DataDriven$structure,fit_DataDriven$var_types) ###################### SIMULATION ###################################### # Fonction pour passer des uniformes aux donnees avec l'??chelle d'orgine ### Fonction inverse de Pseudo Obs PseudoObsInverse<-function(Data,SimUniformData) { PsInverse<-list() for (j in 1:ncol(Data)){ ecdfj<-ecdf(Data[,j]) ECDFvar<-get("x",environment(ecdfj)) ECDFjump<-get("y",environment(ecdfj)) PsInverse[[j]]<-stepfun(ECDFjump[-length(ECDFjump)],ECDFvar) } SimData<-matrix(0,nrow(SimUniformData),ncol(SimUniformData)) for (j in 1:ncol(SimUniformData)){ SimData[,j]<-PsInverse[[j]](SimUniformData[,j]) } SimData<-as.data.frame(SimData) return(SimData) } ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine # taille de l'??chantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels Data_Simu<-PseudoObsInverse(DataBaseline,U_Simu) for (i in 1:2){Data_Simu[,i] = as.factor(Data_Simu[,i])} for (i in 3:8){Data_Simu[,i] = as.numeric(as.character(Data_Simu[,i]))} colnames(Data_Simu)<-colnames(DataBaseline) levels(Data_Simu[,1])<-c('No','Yes') levels(Data_Simu[,2])<-c('<=25','25-35','>35') return(Data_Simu) } # ## select by log-likelihood criterion from one-paramter families # fit <- vinecop(u, family_set = "onepar", selcrit = "bic") # summary(fit) # # ## Gaussian D-vine # fit <- vinecop(u, structure = dvine_structure(1:5), family = "gauss") # plot(fit) # contour(fit) # # ## Partial structure selection with only first tree specified # structure <- rvine_structure(order = 1:5, list(rep(5, 4))) # structure # fit <- vinecop(u, structure = structure, family = "gauss") # plot(fit) # # ## 1-truncated model with random structure # fit <- vinecop(u, structure = rvine_structure_sim(5), trunc_lvl = 1) # contour(fit) ``` ## Apprentissage via Random Forest pour predire l'outcome ```{r,echo=F,message=F} #################################################################### ## estimation du "modele" de prediction avec les donnees origine #################################################################### # Run algorithms using 10-fold cross validation control <- trainControl(method="cv", number=10) metric <- "Accuracy" # Random Forest set.seed(7) fit.rf <- train(Pronostic~., data=Data_Originale, method="rf", metric=metric, trControl=control) #fit.rf$finalModel$confusion #################################################################### #### Fonction pour simuler l'outcome ?? partir du "modele" et des donn??es simulees #### En tenant compte ou pas de l'erreur de prediction #################################################################### Simulation_DataOutcome<- function(Modele,Covariables) { # prediction sans erreur predictionsSE = NULL predictionsSE <- predict(Modele,Covariables) # prediction avec erreur n<-nrow(Covariables) predictionsAE = NULL for (i in 1:n) { RAN = runif(1) if (predictionsSE[i] == "1") { if (RAN < fit.rf$finalModel$confusion[1,3]) {predictionsAE[i] = "0" } else {predictionsAE[i] = "1" }} if (predictionsSE[i] == "0") { if (RAN < fit.rf$finalModel$confusion[2,3]) {predictionsAE[i] = "1" } else {predictionsAE[i] = "0" }} } predictionsAE = as.factor(predictionsAE) DataOutcome_Simu <- cbind(Covariables,predictionsSE,predictionsAE) return(DataOutcome_Simu) } #################################################################### #### Fonction pour calculer les p-values #################################################################### CalculPvalue<- function(DataOutcome) { # pvalues pour outcome sans erreur Pvalue_SE<-NULL Pvalue_SE[1]<-chisq.test(DataOutcome[,1], DataOutcome[,9])$p.value Pvalue_SE[2]<-chisq.test(DataOutcome[,2], DataOutcome[,9])$p.value for (i in 3:8) { Pvalue_SE[i]<-t.test(DataOutcome[,i] ~ DataOutcome[,9])$p.value } # pvalues pour outcome avec erreur Pvalue_AE<-NULL Pvalue_AE[1]<-chisq.test(DataOutcome[,1], DataOutcome[,10])$p.value Pvalue_AE[2]<-chisq.test(DataOutcome[,2], DataOutcome[,10])$p.value for (i in 3:8) { Pvalue_AE[i]<-t.test(DataOutcome[,i] ~ DataOutcome[,10])$p.value } return(list(p_SE=Pvalue_SE,p_AE=Pvalue_AE)) } ``` ## Representation d'une simulation avec la methode discrete ```{r,echo=F,message=F,fig.width=12, fig.height=12,eval=F} DataSimu_Discrete<-Simulation_Discrete(VPMod) Simu_Discrete<-Simulation_DataOutcome(fit.rf, DataSimu_Discrete) ggpairs(DataSimu_Discrete) p<-ggpairs(Simu_Discrete[-9], aes(colour=predictionsAE), lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1))) for(i in 1:p$nrow) { for(j in 1:p$ncol){ p[i,j] <- p[i,j] + scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) + scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) } } p ``` ## Representation d'une simulation avec la methode Continue ```{r,echo=F,message=F,fig.width=12, fig.height=12,eval=F} DataSimu_Continue<-Simulation_Continue(VPMod) Simu_Continue<-Simulation_DataOutcome(fit.rf, DataSimu_Continue) ggpairs(DataSimu_Continue) p<-ggpairs(Simu_Continue[-9], aes(colour=predictionsAE), lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1))) for(i in 1:p$nrow) { for(j in 1:p$ncol){ p[i,j] <- p[i,j] + scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) + scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) } } p ``` ## Representation d'une simulation avec la methode Copula ```{r,echo=F,message=F,fig.width=12, fig.height=12,eval=F} DataSimu_Copule<-Simulation_Copule(VPMod) Simu_Copule<-Simulation_DataOutcome(fit.rf, DataSimu_Copule) ggpairs(DataSimu_Copule) p<-ggpairs(Simu_Copule[-9], aes(colour=predictionsAE), lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1))) for(i in 1:p$nrow) { for(j in 1:p$ncol){ p[i,j] <- p[i,j] + scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) + scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) } } p ``` ## Representation de certaines correlations pour les 3 methode + originale ```{r,echo=F,message=F,fig.width=7, fig.height=6,dpi=600,eval=F} plot_grid( ggmatrix_gtable(ggpairs(VPMod[,c(3,5)],title ="Original data")), ggmatrix_gtable(ggpairs(DataSimu_Discrete[,c(3,5)],title ="Simulation : discrete approach")), ggmatrix_gtable(ggpairs(DataSimu_Continue[,c(3,5)],title ="Simulation : continuous approach")), ggmatrix_gtable(ggpairs(DataSimu_Copule[,c(3,5)],title ="Simulation : Copula approach")),nrow=2) plot_grid( ggmatrix_gtable(ggpairs(VPMod[,c(4,8)],title ="Original data")), ggmatrix_gtable(ggpairs(DataSimu_Discrete[,c(4,8)],title ="Simulation : discrete approach")), ggmatrix_gtable(ggpairs(DataSimu_Continue[,c(4,8)],title ="Simulation : continuous approach")), ggmatrix_gtable(ggpairs(DataSimu_Copule[,c(4,8)],title ="Simulation : Copula approach")),nrow=2) ``` ## Population d'origine ```{r,echo=F,message=F,fig.width=12, fig.height=12} ggpairs(Data_Originale[,c(1:9)]) ``` ## Individus virtuels issus de la population d'origine : ### simulations de 1000 individus suivant les estimations de la copule et des marginales (empirique) ```{r,echo=F,message=F,fig.width=12, fig.height=12} # simulation de donnees avec la m??thode copule N=500 N<-1000 # estimation de la copule sur les vraies donnees ###################### ESTIMATION ###################################### #Pseudo inverse variables continues U_cont<-pseudo_obs(Data_Originale[,3:8]) #Distribution des 2 variables discretes disc_1<-as.integer(Data_Originale[,1])-1 disc_2<-as.integer(Data_Originale[,2])-1 freq_disc1<-prop.table(table(Data_Originale[,1])) freq_disc2<-prop.table(table(Data_Originale[,2])) Freq_disc_t1 <- cbind(pdiscrete(disc_1+1,freq_disc1),pdiscrete(disc_2+1,freq_disc2)) Freq_disc_t0 <- cbind(pdiscrete(disc_1,freq_disc1),pdiscrete(disc_2,freq_disc2)) #densite :ddiscrete(x+1,freq) #distribution function : pdiscrete(x+1,freq) #quantile function : qdiscrete(u[, 1],freq) - 1 # Infos necessaire pour copule avec donnees mixte U_mixte<-cbind(Freq_disc_t1,U_cont,Freq_disc_t0) ## Model for mixed data fit_DataDriven <- vinecop(U_mixte, var_types = c("d","d", rep("c", 6))) # definition de la distribution de la r-vine Fit_dist<-vinecop_dist(fit_DataDriven$pair_copulas,fit_DataDriven$structure,fit_DataDriven$var_types) plot(fit_DataDriven,edge_labels="pair", tree=c(1,2)) plot(fit_DataDriven,edge_labels="family_tau", tree=1) # Simulation en fonction du modele de copule et des marginales ###################### SIMULATION ###################################### # Fonction pour passer des uniformes aux donnees avec l'echelle d'origine ### Fonction inverse de Pseudo Obs PseudoObsInverse<-function(Data,SimUniformData) { PsInverse<-list() for (j in 1:ncol(Data)){ ecdfj<-ecdf(Data[,j]) ECDFvar<-get("x",environment(ecdfj)) ECDFjump<-get("y",environment(ecdfj)) PsInverse[[j]]<-stepfun(ECDFjump[-length(ECDFjump)],ECDFvar) } SimData<-matrix(0,nrow(SimUniformData),ncol(SimUniformData)) for (j in 1:ncol(SimUniformData)){ SimData[,j]<-PsInverse[[j]](SimUniformData[,j]) } SimData<-as.data.frame(SimData) return(SimData) } ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar<-PseudoObsInverse(Data_Originale,U_Simu) for (i in 1:2){SimuVar[,i] = as.factor(SimuVar[,i])} for (i in 3:8){SimuVar[,i] = as.numeric(as.character(SimuVar[,i]))} colnames(SimuVar)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar[,1])<-c('No','Yes') levels(SimuVar[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu<-Simulation_DataOutcome(fit.rf, SimuVar) ############### representations graphiques ####################@ #ggpairs(SimuVar) #ggpairs(Simu) ``` ## Individus virtuels issus de la population d'origine : meme support (troncature) !! ### simulations d'individus suivant les estimations de la copule et des marginales (empirique, sauf une marginale modelisee par une loi log normale ```{r,echo=F,message=F,fig.width=12, fig.height=12} library(fitdistrplus) library(truncdist) library(truncnorm) ###################################### fit de la meilleure distribution pour Glucose ###################################### fit_Glucose_norm <- fitdist(Data_Originale$Glucose, "norm") fit_Glucose_lnorm <- fitdist(Data_Originale$Glucose, "lnorm") fit_Glucose_gamma <- fitdist(Data_Originale$Glucose, "gamma") fit_Glucose_weibull <- fitdist(Data_Originale$Glucose, "weibull") # par(mfrow=c(2,2)) # cdfcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # denscomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # qqcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # ppcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # gofstat(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), fitnames=c("normal", "lognormal", "gamma","Weibull")) ####################################### Simulation suivant le meilleur modele parametrique pour la distribution marginale ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique) GlucoseNew<-rlnorm(nrow(Data_Originale),meanlog=4.778,sdlog=0.249) # Simulation d'une distribution marginale via des loi parametrique tronquee (ici loi normale, valeur >140 ) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a=140,b=300 ,mean=122.63,sd=30.82) # attention ... si observations en dehors du support original (<56, >198), il faut en générer une autre !! #summary(Data_Originale$Glucose) #summary(GlucoseNew) for (i in 1:length(GlucoseNew)) { if (GlucoseNew[i]>198 | GlucoseNew[i]<56) { GlucoseNew[i]<-rlnorm(1,meanlog=4.778,sdlog=0.249) } } for (i in 1:length(GlucoseNew)) { if (GlucoseNew[i]>198 | GlucoseNew[i]<56) { GlucoseNew[i]<-rlnorm(1,meanlog=4.778,sdlog=0.249) } } #summary(GlucoseNew) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar1<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar1[,i] = as.factor(SimuVar1[,i])} for (i in 3:8){SimuVar1[,i] = as.numeric(as.character(SimuVar1[,i]))} colnames(SimuVar1)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar1[,1])<-c('No','Yes') levels(SimuVar1[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu1<-Simulation_DataOutcome(fit.rf, SimuVar1) ###################### Representations graphiques ###################################### #ggpairs(SimuVar1) ``` ## Individus virtuels issus de la population d'origine : meme support !! ### simulations d'individus suivant les estimations de la copule et des marginales (empirique, sauf une marginale modelisee par une loi normale tronquee ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### fit_Glucose_norm <- fitdist(Data_Originale$Glucose, "norm") fit_Glucose_lnorm <- fitdist(Data_Originale$Glucose, "lnorm") fit_Glucose_gamma <- fitdist(Data_Originale$Glucose, "gamma") fit_Glucose_weibull <- fitdist(Data_Originale$Glucose, "weibull") # par(mfrow=c(2,2)) # cdfcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # denscomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # qqcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # ppcomp(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), legendtext=c("normal", "lognormal", "gamma","Weibull")) # gofstat(list(fit_Glucose_norm, fit_Glucose_lnorm, fit_Glucose_gamma, fit_Glucose_weibull), fitnames=c("normal", "lognormal", "gamma","Weibull")) ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a=56,b=198 ,mean=mean(Data_Originale[,3]),sd=sd(Data_Originale[,3])) hist(GlucoseNew) # Simulation d'une distribution marginale via des loi parametrique tronquee (ici loi normale, valeur >140 ) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a=140,b=300 ,mean=122.63,sd=30.82) # attention ... si observations en dehors du support original (<56, >198), il faut en générer une autre !! #summary(Data_Originale$Glucose) #summary(GlucoseNew) # for (i in 1:length(GlucoseNew)) { # if (GlucoseNew[i]>198 | GlucoseNew[i]<56) { # GlucoseNew[i]<-rlnorm(1,meanlog=4.778,sdlog=0.249) # } # } # for (i in 1:length(GlucoseNew)) { # if (GlucoseNew[i]>198 | GlucoseNew[i]<56) { # GlucoseNew[i]<-rlnorm(1,meanlog=4.778,sdlog=0.249) # } # } #summary(GlucoseNew) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar2<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar2[,i] = as.factor(SimuVar2[,i])} for (i in 3:8){SimuVar2[,i] = as.numeric(as.character(SimuVar2[,i]))} colnames(SimuVar2)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar2[,1])<-c('No','Yes') levels(SimuVar2[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu2<-Simulation_DataOutcome(fit.rf, SimuVar2) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` ## Representation de certaines correlations pour donnees originales, simulees non param, simulee param 1 et simulee param 2 (dans les 2 dernier cas, la distrib est tronquee au support de la variable) ```{r,echo=F,message=F,fig.width=7, fig.height=6,dpi=600} colnames(Simu)[9]<-"Prediction" colnames(Simu1)[9]<-"Prediction" colnames(Simu2)[9]<-"Prediction" plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,6,9)],title ="Original data")), ggmatrix_gtable(ggpairs(Simu[,c(3,6,9)],title ="Simulation using empirical distribution for Glucose")), ggmatrix_gtable(ggpairs(Simu1[,c(3,6,9)],title ="Simulation using Lognormal distribution for Glucose")), ggmatrix_gtable(ggpairs(Simu2[,c(3,6,9)],title ="Simulation using Normal distribution for Glucose")), nrow=2) # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,6)],title ="Original data")), # ggmatrix_gtable(ggpairs(SimuVar[,c(3,6)],title ="Simulation using empirical distribution for Glucose")), # ggmatrix_gtable(ggpairs(SimuVar1[,c(3,6)],title ="Simulation using Lognormal distribution for Glucose")), # ggmatrix_gtable(ggpairs(SimuVar2[,c(3,6)],title ="Simulation using Normal distribution for Glucose")), # nrow=2) # # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,9)],title ="Original data")), # ggmatrix_gtable(ggpairs(Simu[,c(3,9)],title ="Standard simulation")), # ggmatrix_gtable(ggpairs(Simu1[,c(3,9)],title ="Simulation with lognormal distribution for Glucose")), # ggmatrix_gtable(ggpairs(Simu2[,c(3,9)],title ="Simulation with normal distribution for Glucose")),nrow=2) # # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,9)],title ="Original data")), # ggmatrix_gtable(ggpairs(Simu[,c(3,10)],title ="Standard simulation")), # ggmatrix_gtable(ggpairs(Simu1[,c(3,10)],title ="Simulation with lognormal distribution for Glucose")), # ggmatrix_gtable(ggpairs(Simu2[,c(3,10)],title ="Simulation with normal distribution for Glucose")),nrow=2) # c(4,8) # c(3,6) # c(3,8) ``` # Individus virtuels dans un intervalle de la population d'origine ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) ### filtering pour conserver les individus satisfaisant un critere d'inclusion ```{r,echo=F,message=F,fig.width=12, fig.height=12} # simulation de donnees avec la methode copule N=3000 NN<-6000 ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine # taille de l'echantillon definit avant U_Simu<-rvinecop(NN,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee Simu_afiltrer<-PseudoObsInverse(Data_Originale,U_Simu) for (i in 1:2){Simu_afiltrer[,i] = as.factor(Simu_afiltrer[,i])} for (i in 3:8){Simu_afiltrer[,i] = as.numeric(as.character(Simu_afiltrer[,i]))} colnames(Simu_afiltrer)<-colnames(Data_Originale[,c(1:8)]) levels(Simu_afiltrer[,1])<-c('No','Yes') levels(Simu_afiltrer[,2])<-c('<=25','25-35','>35') # Restriction de la population simulee par filtration : uniquement patient tel que glucose >= 140 # attention ... avec estimation empirique, toutes les simulation sont inférieures au max du support d'origine SimuVar3<-Simu_afiltrer[Simu_afiltrer$Glucose > 100 & Simu_afiltrer$Glucose <=150,] # on garde uniquement 500 individus au hasard SimuVar3<-SimuVar3[sample(1:nrow(SimuVar3), N), ] # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu3<-Simulation_DataOutcome(fit.rf, SimuVar3) ###################### Representations graphiques ###################################### #ggpairs(SimuVar3) ``` # Individus virtuels dans un intervalle de la population d'origine ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### ici une gaussienne ajuste sur les données et tronquee a 140 ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 101, b=150, mean=mean(Data_Originale[,3]),sd=sd(Data_Originale[,3])) #hist(GlucoseNew) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar4<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar4[,i] = as.factor(SimuVar4[,i])} for (i in 3:8){SimuVar4[,i] = as.numeric(as.character(SimuVar4[,i]))} colnames(SimuVar4)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar4[,1])<-c('No','Yes') levels(SimuVar4[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu4<-Simulation_DataOutcome(fit.rf, SimuVar4) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` # Individus virtuels dans un intervalle de la population d'origine ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### Une marginale avec une loi parametrique differente (par ex biblio) de la loi des données ... ici lognormale tronquee ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 140, mean=160,sd=50) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 101, b=150, mean=120,sd=30) GlucoseNew <- rlnormTrunc( nrow(Data_Originale), min= 101, max=150, meanlog=120,sdlog=30) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar5<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar5[,i] = as.factor(SimuVar5[,i])} for (i in 3:8){SimuVar5[,i] = as.numeric(as.character(SimuVar5[,i]))} colnames(SimuVar5)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar5[,1])<-c('No','Yes') levels(SimuVar5[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu5<-Simulation_DataOutcome(fit.rf, SimuVar5) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` ## Representation de certaines correlations pour donnees originales, simulees standard et filtering, simulee normal fit et filtering et simulee avec gaussienne tronquee ```{r,echo=F,message=F,fig.width=7, fig.height=6,dpi=600} colnames(Simu3)[9]<-"Prediction" colnames(Simu4)[9]<-"Prediction" colnames(Simu5)[9]<-"Prediction" plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[Data_Originale[,3]>100 & Data_Originale[,3]<=150,c(3,6,9)],title ="Filtering the original data")), ggmatrix_gtable(ggpairs(Simu3[,c(3,6,9)],title ="Filtering the simulated data (based on empirical distribution)")), ggmatrix_gtable(ggpairs(Simu4[,c(3,6,9)],title ="Filtering the simulated data (based on Normal distribution)")), ggmatrix_gtable(ggpairs(Simu5[,c(3,6,9)],title ="Filtering the simulated data (based on Lognormal distribution)")),nrow=2) # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,6)],title ="Original data")), # ggmatrix_gtable(ggpairs(SimuVar3[,c(3,6)],title ="Standard simulation and filtering")), # ggmatrix_gtable(ggpairs(SimuVar4[,c(3,6)],title ="Simulation with fitted normal distribution and filtering")), # ggmatrix_gtable(ggpairs(SimuVar5[,c(3,6)],title ="Simulation with a different normal distribution")),nrow=2) # # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,9)],title ="Original data")), # ggmatrix_gtable(ggpairs(Simu3[,c(3,9)],title ="Standard simulation and filtering")), # ggmatrix_gtable(ggpairs(Simu4[,c(3,9)],title ="Simulation with fitted normal distribution and filtering")), # ggmatrix_gtable(ggpairs(Simu5[,c(3,9)],title ="Simulation with a different normal distribution")),nrow=2) # # plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[,c(3,9)],title ="Original data")), # ggmatrix_gtable(ggpairs(Simu3[,c(3,10)],title ="Standard simulation and filtering")), # ggmatrix_gtable(ggpairs(Simu4[,c(3,10)],title ="Simulation with fitted normal distribution and filtering")), # ggmatrix_gtable(ggpairs(Simu5[,c(3,10)],title ="Simulation with a different normal distribution")),nrow=2) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>140 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) ### filtering pour conserver les individus satisfaisant un critere d'inclusion ```{r,echo=F,message=F,fig.width=12, fig.height=12} # simulation de donnees avec la methode copule N=3000 NN<-60000 ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine # taille de l'echantillon definit avant U_Simu<-rvinecop(NN,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee Simu_afiltrer<-PseudoObsInverse(Data_Originale,U_Simu) for (i in 1:2){Simu_afiltrer[,i] = as.factor(Simu_afiltrer[,i])} for (i in 3:8){Simu_afiltrer[,i] = as.numeric(as.character(Simu_afiltrer[,i]))} colnames(Simu_afiltrer)<-colnames(Data_Originale[,c(1:8)]) levels(Simu_afiltrer[,1])<-c('No','Yes') levels(Simu_afiltrer[,2])<-c('<=25','25-35','>35') # Restriction de la population simulee par filtration : uniquement patient tel que glucose >= 140 # attention ... avec estimation empirique, toutes les simulation sont inférieures au max du support d'origine SimuVar6<-Simu_afiltrer[Simu_afiltrer$Glucose > 140,] # on garde uniquement N individus au hasard SimuVar6<-SimuVar6[sample(1:nrow(SimuVar6), N), ] # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu6<-Simulation_DataOutcome(fit.rf, SimuVar6) ###################### Representations graphiques ###################################### #ggpairs(SimuVar3) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>140 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### ici une gaussienne ajuste sur les données et tronquee a 140 ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 140, mean=mean(Data_Originale[,3]),sd=sd(Data_Originale[,3])) #hist(GlucoseNew) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar7<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar7[,i] = as.factor(SimuVar7[,i])} for (i in 3:8){SimuVar7[,i] = as.numeric(as.character(SimuVar7[,i]))} colnames(SimuVar7)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar7[,1])<-c('No','Yes') levels(SimuVar7[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu7<-Simulation_DataOutcome(fit.rf, SimuVar7) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>140 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### Une marginale avec une loi parametrique differente (par ex biblio) de la loi des données ... ici lognormale tronquee ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 140, mean=160,sd=50) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 101, b=150, mean=120,sd=30) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 170, mean=175,sd=180) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar8<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar8[,i] = as.factor(SimuVar8[,i])} for (i in 3:8){SimuVar8[,i] = as.numeric(as.character(SimuVar8[,i]))} colnames(SimuVar8)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar8[,1])<-c('No','Yes') levels(SimuVar8[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu8<-Simulation_DataOutcome(fit.rf, SimuVar8) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` # simulation en dehors du support d'origine \>140 ## Representation de certaines correlations pour donnees originales, simulees standard et filtering, simulee normal fit et filtering et simulee avec lognormal tronquee ```{r,echo=F,message=F,fig.width=7, fig.height=6,dpi=600} colnames(Simu6)[9]<-"Prediction" colnames(Simu7)[9]<-"Prediction" colnames(Simu8)[9]<-"Prediction" plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[Data_Originale[,3]>140,c(3,6,9)],title ="Filtering the original data")), ggmatrix_gtable(ggpairs(Simu6[,c(3,6,9)],title ="Filtering the simulated data (based on empirical distribution)")), ggmatrix_gtable(ggpairs(Simu7[,c(3,6,9)],title ="Filtering the simulated data (based on best Normal distribution)")), ggmatrix_gtable(ggpairs(Simu8[,c(3,6,9)],title ="Filtering the simulated data (based on a given Normal distribution)")),nrow=2) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>170 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) ### filtering pour conserver les individus satisfaisant un critere d'inclusion ```{r,echo=F,message=F,fig.width=12, fig.height=12} # simulation de donnees avec la methode copule N=3000 NN<-60000 ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine # taille de l'echantillon definit avant U_Simu<-rvinecop(NN,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee Simu_afiltrer<-PseudoObsInverse(Data_Originale,U_Simu) for (i in 1:2){Simu_afiltrer[,i] = as.factor(Simu_afiltrer[,i])} for (i in 3:8){Simu_afiltrer[,i] = as.numeric(as.character(Simu_afiltrer[,i]))} colnames(Simu_afiltrer)<-colnames(Data_Originale[,c(1:8)]) levels(Simu_afiltrer[,1])<-c('No','Yes') levels(Simu_afiltrer[,2])<-c('<=25','25-35','>35') # Restriction de la population simulee par filtration : uniquement patient tel que glucose >= 140 # attention ... avec estimation empirique, toutes les simulation sont inférieures au max du support d'origine SimuVar6<-Simu_afiltrer[Simu_afiltrer$Glucose > 170,] # on garde uniquement N individus au hasard SimuVar6<-SimuVar6[sample(1:nrow(SimuVar6), N), ] # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu6<-Simulation_DataOutcome(fit.rf, SimuVar6) ###################### Representations graphiques ###################################### #ggpairs(SimuVar3) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>170 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### ici une gaussienne ajuste sur les données et tronquee a 140 ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 170, mean=mean(Data_Originale[,3]),sd=sd(Data_Originale[,3])) #hist(GlucoseNew) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar7<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar7[,i] = as.factor(SimuVar7[,i])} for (i in 3:8){SimuVar7[,i] = as.numeric(as.character(SimuVar7[,i]))} colnames(SimuVar7)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar7[,1])<-c('No','Yes') levels(SimuVar7[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu7<-Simulation_DataOutcome(fit.rf, SimuVar7) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` # Individus virtuels dans un intervalle en dehors du support de la population d'origine \>170 ### simulations d'individus suivant les estimations de la copule et des marginales (empirique) et une param ### Une marginale avec une loi parametrique differente (par ex biblio) de la loi des données ... ici lognormale tronquee ```{r,echo=F,message=F,fig.width=12, fig.height=12} ###################################### fit de la meilleure distribution pour Glucose ###################################### ###################### Preparation des donnees ###################################### # Modification des donnees univariee ... nouvelle variable glucose (simuler via modele parametrique normal tronquee) #GlucoseNew<-rweibull(nrow(Data_Originale),shape=4.225,scale=134.705) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 140, mean=160,sd=50) #GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 101, b=150, mean=120,sd=30) GlucoseNew <- rtruncnorm( nrow(Data_Originale), a= 170, mean=175,sd=180) # definition de la base de donnees Data_modif<-cbind(Data_Originale[,1:2],GlucoseNew,Data_Originale[,4:8]) # Simulation en fonction du modele de copule et des marginales (dont une modifiee) ###################### SIMULATION ###################################### ## Simulation d'un echantillon d'uniformes a partir de la distribution de la vine (estimee sur donnees originales) # taille de l'echantillon definit avant U_Simu<-rvinecop(N,Fit_dist) # Simulation de patients virtuels avec la copule sur les donnees d'origine et une marginale modifiee SimuVar8<-PseudoObsInverse(Data_modif,U_Simu) for (i in 1:2){SimuVar8[,i] = as.factor(SimuVar8[,i])} for (i in 3:8){SimuVar8[,i] = as.numeric(as.character(SimuVar8[,i]))} colnames(SimuVar8)<-colnames(Data_Originale[,c(1:8)]) levels(SimuVar8[,1])<-c('No','Yes') levels(SimuVar8[,2])<-c('<=25','25-35','>35') # prevision de l'outcome avec le modele random forest estimer (plus haut) ############### Prevision ####################@ Simu8<-Simulation_DataOutcome(fit.rf, SimuVar8) ###################### Representations graphiques ###################################### #ggpairs(SimuVar2) ``` # simulation en dehors du support d'origine \>170 ## Representation de certaines correlations pour donnees originales, simulees standard et filtering, simulee normal fit et filtering et simulee avec lognormal tronquee ```{r,echo=F,message=F,fig.width=7, fig.height=6,dpi=600} colnames(Simu6)[9]<-"Prediction" colnames(Simu7)[9]<-"Prediction" colnames(Simu8)[9]<-"Prediction" plot_grid( ggmatrix_gtable(ggpairs(Data_Originale[Data_Originale[,3]>170,c(3,6,9)],title ="Filtering the original data")), ggmatrix_gtable(ggpairs(Simu6[,c(3,6,9)],title ="Filtering the simulated data (based on empirical distribution)")), ggmatrix_gtable(ggpairs(Simu7[,c(3,6,9)],title ="Filtering the simulated data (based on best Normal distribution)")), ggmatrix_gtable(ggpairs(Simu8[,c(3,6,9)],title ="Filtering the simulated data (based on a given Normal distribution)")),nrow=2) ``` ## Calcul des p-values pour plusieurs simulations Le code R de cette section n'est pas evalue. ```{r,echo=F,message=F,eval = FALSE} #################################################################### ## P-value sur donnees originale (entre diabete et les autres variables) #################################################################### Pvalue_vraie<-NULL Pvalue_vraie[1]<-chisq.test(Data_Originale[,1], Data_Originale[,9])$p.value Pvalue_vraie[2]<-chisq.test(Data_Originale[,2], Data_Originale[,9])$p.value for (i in 3:8) { Pvalue_vraie[i]<-t.test(Data_Originale[,i] ~ Data_Originale[,9])$p.value } Pvalue_vraie #################################################################### #### Boucle Simulation methode discrete #################################################################### Mat_Pvalue_Disc_SE<-NULL Mat_Pvalue_Disc_AE<-NULL for (j in 1:Nsimu) { Simu_Cov<-NULL Simu_Data<-NULL Pvalue_SE<-NULL Pvalue_AE<-NULL Simu_Cov<-Simulation_Discrete(VPMod) Simu_Data<-Simulation_DataOutcome(fit.rf, Simu_Cov) Pvalue_SE<-CalculPvalue(Simu_Data)$p_SE Pvalue_AE<-CalculPvalue(Simu_Data)$p_AE Mat_Pvalue_Disc_SE<-rbind(Mat_Pvalue_Disc_SE,Pvalue_SE) Mat_Pvalue_Disc_AE<-rbind(Mat_Pvalue_Disc_AE,Pvalue_AE) } colnames(Mat_Pvalue_Disc_SE)<-colnames(VPMod) row.names(Mat_Pvalue_Disc_SE)<-c(1:nrow(Mat_Pvalue_Disc_SE)) Mat_Pvalue_Disc_SE<-as.data.frame(Mat_Pvalue_Disc_SE) colnames(Mat_Pvalue_Disc_AE)<-colnames(VPMod) row.names(Mat_Pvalue_Disc_AE)<-c(1:nrow(Mat_Pvalue_Disc_AE)) Mat_Pvalue_Disc_AE<-as.data.frame(Mat_Pvalue_Disc_AE) #################################################################### #### Boucle Simulation methode continue #################################################################### Mat_Pvalue_Cont_SE<-NULL Mat_Pvalue_Cont_AE<-NULL for (j in 1:Nsimu) { Simu_Cov<-NULL Simu_Data<-NULL Pvalue_SE<-NULL Pvalue_AE<-NULL Simu_Cov<-Simulation_Continue(VPMod) Simu_Data<-Simulation_DataOutcome(fit.rf, Simu_Cov) Pvalue_SE<-CalculPvalue(Simu_Data)$p_SE Pvalue_AE<-CalculPvalue(Simu_Data)$p_AE Mat_Pvalue_Cont_SE<-rbind(Mat_Pvalue_Cont_SE,Pvalue_SE) Mat_Pvalue_Cont_AE<-rbind(Mat_Pvalue_Cont_AE,Pvalue_AE) } colnames(Mat_Pvalue_Cont_SE)<-colnames(VPMod) row.names(Mat_Pvalue_Cont_SE)<-c(1:nrow(Mat_Pvalue_Cont_SE)) Mat_Pvalue_Cont_SE<-as.data.frame(Mat_Pvalue_Cont_SE) colnames(Mat_Pvalue_Cont_AE)<-colnames(VPMod) row.names(Mat_Pvalue_Cont_AE)<-c(1:nrow(Mat_Pvalue_Cont_AE)) Mat_Pvalue_Cont_AE<-as.data.frame(Mat_Pvalue_Cont_AE) #################################################################### #### Boucle Simulation methode copule #################################################################### Mat_Pvalue_Cop_SE<-NULL Mat_Pvalue_Cop_AE<-NULL for (j in 1:Nsimu) { Simu_Cov<-NULL Simu_Data<-NULL Pvalue_SE<-NULL Pvalue_AE<-NULL Simu_Cov<-Simulation_Copule(VPMod) Simu_Data<-Simulation_DataOutcome(fit.rf, Simu_Cov) Pvalue_SE<-CalculPvalue(Simu_Data)$p_SE Pvalue_AE<-CalculPvalue(Simu_Data)$p_AE Mat_Pvalue_Cop_SE<-rbind(Mat_Pvalue_Cop_SE,Pvalue_SE) Mat_Pvalue_Cop_AE<-rbind(Mat_Pvalue_Cop_AE,Pvalue_AE) #print(j) } colnames(Mat_Pvalue_Cop_SE)<-colnames(VPMod) row.names(Mat_Pvalue_Cop_SE)<-c(1:nrow(Mat_Pvalue_Cop_SE)) Mat_Pvalue_Cop_SE<-as.data.frame(Mat_Pvalue_Cop_SE) colnames(Mat_Pvalue_Cop_AE)<-colnames(VPMod) row.names(Mat_Pvalue_Cop_AE)<-c(1:nrow(Mat_Pvalue_Cop_AE)) Mat_Pvalue_Cop_AE<-as.data.frame(Mat_Pvalue_Cop_AE) ``` ## Boxplots divers ```{r,warning=F,echo=F,message=F,fig.width=8, fig.height=5,dpi=600,eval=F} # ## Boxplot P-value methode discrete # ############################################################ # Plot_Discrete<-function(Mat_SE,Mat_AE,IndexVariable) { # vectorize_Pvalue<-NULL # vectorize_Label<-NULL # vectorize_Pvalue<-c(Mat_SE[,IndexVariable],Mat_AE[,IndexVariable]) # vectorize_Label<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) # vectorize<-cbind(vectorize_Pvalue,vectorize_Label) # vectorize<-as.data.frame(vectorize) # vectorize[,2]<-as.factor(vectorize[,2]) # levels(vectorize[,2])<-c("Without error","With error") # # Plot_pvalue<-ggplot(vectorize) + # geom_jitter(aes(x=vectorize_Label, y=vectorize_Pvalue,color = vectorize_Label), width=0.2, alpha=0.5) + # geom_boxplot(aes(x=vectorize_Label, y=vectorize_Pvalue), width=0.1) + # scale_color_brewer(palette="Dark2") + theme(legend.position='none') + # stat_summary(geom = "point", mapping = aes(x=vectorize_Label, y=vectorize_Pvalue), fun = "mean", shape=17, size=2) + # labs(x ="Simulation Error", y = "P-value Cov. x Diabete") + # ggtitle("P-value for discrete method \n True p-value =",round(Pvalue_vraie[IndexVariable],8)) + # theme(plot.title = element_text(hjust = 0.5)) # } # # ## Boxplot P-value methode continue # ############################################################ # Plot_Continue<-function(Mat_SE,Mat_AE,IndexVariable) { # vectorize_Pvalue<-NULL # vectorize_Label<-NULL # vectorize_Pvalue<-c(Mat_SE[,IndexVariable],Mat_AE[,IndexVariable]) # vectorize_Label<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) # vectorize<-cbind(vectorize_Pvalue,vectorize_Label) # vectorize<-as.data.frame(vectorize) # vectorize[,2]<-as.factor(vectorize[,2]) # levels(vectorize[,2])<-c("Without error","With error") # # Plot_pvalue<-ggplot(vectorize) + # geom_jitter(aes(x=vectorize_Label, y=vectorize_Pvalue,color = vectorize_Label), width=0.2, alpha=0.5) + # geom_boxplot(aes(x=vectorize_Label, y=vectorize_Pvalue), width=0.1) + # scale_color_brewer(palette="Dark2") + theme(legend.position='none') + # stat_summary(geom = "point", mapping = aes(x=vectorize_Label, y=vectorize_Pvalue), fun = "mean", shape=17, size=2) + # labs(x ="Simulation Error", y = "P-value Cov. x Diabete") + # ggtitle("P-value for continuous method \n True p-value =",round(Pvalue_vraie[IndexVariable],8)) + # theme(plot.title = element_text(hjust = 0.5)) # } # # ## Boxplot P-value methode copule # ############################################################ # Plot_Copule<-function(Mat_SE,Mat_AE,IndexVariable) { # vectorize_Pvalue<-NULL # vectorize_Label<-NULL # vectorize_Pvalue<-c(Mat_SE[,IndexVariable],Mat_AE[,IndexVariable]) # vectorize_Label<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) # vectorize<-cbind(vectorize_Pvalue,vectorize_Label) # vectorize<-as.data.frame(vectorize) # vectorize[,2]<-as.factor(vectorize[,2]) # levels(vectorize[,2])<-c("Without error","With error") # # Plot_pvalue<-ggplot(vectorize) + # geom_jitter(aes(x=vectorize_Label, y=vectorize_Pvalue,color = vectorize_Label), width=0.2, alpha=0.5) + # geom_boxplot(aes(x=vectorize_Label, y=vectorize_Pvalue), width=0.1) + # scale_color_brewer(palette="Dark2") + theme(legend.position='none') + # stat_summary(geom = "point", mapping = aes(x=vectorize_Label, y=vectorize_Pvalue), fun = "mean", shape=17, size=2) + # labs(x ="Simulation Error", y = "P-value Cov. x Diabete") + # ggtitle("P-value for copula method \n True p-value = ",round(Pvalue_vraie[IndexVariable],8)) + # theme(plot.title = element_text(hjust = 0.5)) # } # # ## Boxplot P-value avec les 3 methodes # ############################################################ # Plot_3methode<-function(Mat_D,Mat_C,Mat_Cop,IndexVariable) { # vectorize_Pvalue<-NULL # vectorize_Label<-NULL # vectorize_Pvalue<-c(Mat_D[,IndexVariable],Mat_C[,IndexVariable],Mat_Cop[,IndexVariable]) # vectorize_Label<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu),rep(2,Nsimu))) # vectorize<-cbind(vectorize_Pvalue,vectorize_Label) # vectorize<-as.data.frame(vectorize) # vectorize[,2]<-as.factor(vectorize[,2]) # levels(vectorize[,2])<-c("Discrete","Continuous", "Copula") # # Plot_pvalue<-ggplot(vectorize) + # geom_jitter(aes(x=vectorize_Label, y=vectorize_Pvalue,color = vectorize_Label), width=0.2, alpha=0.5) + # geom_boxplot(aes(x=vectorize_Label, y=vectorize_Pvalue), width=0.1) + # scale_color_brewer(palette="Dark2") + theme(legend.position='none') + # stat_summary(geom = "point", mapping = aes(x=vectorize_Label, y=vectorize_Pvalue), fun = "mean", shape=17, size=2) + # labs(x ="P-value", y = "P-value Cov. x Diabete") + # ggtitle("P-value for the three methods \n True p-value = ",round(Pvalue_vraie[IndexVariable],8)) + # theme(plot.title = element_text(hjust = 0.5)) # } ############# Representation ############################ # P_Disc_1<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,1) # P_Disc_2<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,2) # P_Disc_3<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,3) # P_Disc_4<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,4) # P_Disc_5<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,5) # P_Disc_6<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,6) # P_Disc_7<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,7) # P_Disc_8<-Plot_Discrete(Mat_Pvalue_Disc_SE,Mat_Pvalue_Disc_AE,8) # # P_Cont_1<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,1) # P_Cont_2<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,2) # P_Cont_3<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,3) # P_Cont_4<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,4) # P_Cont_5<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,5) # P_Cont_6<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,6) # P_Cont_7<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,7) # P_Cont_8<-Plot_Continue(Mat_Pvalue_Cont_SE,Mat_Pvalue_Cont_AE,8) # # P_Cop_1<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,1) # P_Cop_2<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,2) # P_Cop_3<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,3) # P_Cop_4<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,4) # P_Cop_5<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,5) # P_Cop_6<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,6) # P_Cop_7<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,7) # P_Cop_8<-Plot_Copule(Mat_Pvalue_Cop_SE,Mat_Pvalue_Cop_AE,8) # # plot_grid(P_Disc_1,P_Disc_2,P_Disc_3,P_Disc_4,P_Disc_5,P_Disc_6,P_Disc_7,P_Disc_8,nrow=2) # plot_grid(P_Cont_1,P_Cont_2,P_Cont_3,P_Cont_4,P_Cont_5,P_Cont_6,P_Cont_7,P_Cont_8,nrow=2) # plot_grid(P_Cop_1,P_Cop_2,P_Cop_3,P_Cop_4,P_Cop_5,P_Cop_6,P_Cop_7,P_Cop_8,nrow=2) # 3 graphs erreur vs. sans erreur # Graph avec boxplot pour chaque methode (meme graph) pour chaque covariable # Tableau : vraie p-value et m??diane (+ moyenne ?) ##### grah simulation avec erreur et sans erreur ########### #### Fonction pour 'vectoriser' les p-values (on ajoute une variable pour la methode et une pour erreur ou pas) Vectorize<-function(IndexVariable) { vectorize_Pvalue_Disc<-NULL vectorize_Label_Disc<-NULL vectorize_Pvalue_Disc<-c(Mat_Pvalue_Disc_SE[,IndexVariable],Mat_Pvalue_Disc_AE[,IndexVariable]) vectorize_Label_Disc<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) vectorize_Disc<-cbind(vectorize_Pvalue_Disc,vectorize_Label_Disc,rep(1,length(vectorize_Label_Disc))) vectorize_Disc<-as.data.frame(vectorize_Disc) vectorize_Disc[,2]<-as.factor(vectorize_Disc[,2]) vectorize_Disc[,3]<-as.factor(vectorize_Disc[,3]) levels(vectorize_Disc[,2])<-c("Without error","With error") levels(vectorize_Disc[,3])<-c("Discrete method") colnames(vectorize_Disc)<-c("P_value","Error","Method") vectorize_Pvalue_Cont<-NULL vectorize_Label_Cont<-NULL vectorize_Pvalue_Cont<-c(Mat_Pvalue_Cont_SE[,IndexVariable],Mat_Pvalue_Cont_AE[,IndexVariable]) vectorize_Label_Cont<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) vectorize_Cont<-cbind(vectorize_Pvalue_Cont,vectorize_Label_Cont,rep(1,length(vectorize_Label_Cont))) vectorize_Cont<-as.data.frame(vectorize_Cont) vectorize_Cont[,2]<-as.factor(vectorize_Cont[,2]) vectorize_Cont[,3]<-as.factor(vectorize_Cont[,3]) levels(vectorize_Cont[,2])<-c("Without error","With error") levels(vectorize_Cont[,3])<-c("Continuous method") colnames(vectorize_Cont)<-c("P_value","Error","Method") vectorize_Pvalue_Cop<-NULL vectorize_Label_Cop<-NULL vectorize_Pvalue_Cop<-c(Mat_Pvalue_Cop_SE[,IndexVariable],Mat_Pvalue_Cop_AE[,IndexVariable]) vectorize_Label_Cop<-as.factor(c(rep(0,Nsimu),rep(1,Nsimu))) vectorize_Cop<-cbind(vectorize_Pvalue_Cop,vectorize_Label_Cop,rep(1,length(vectorize_Label_Cop))) vectorize_Cop<-as.data.frame(vectorize_Cop) vectorize_Cop[,2]<-as.factor(vectorize_Cop[,2]) vectorize_Cop[,3]<-as.factor(vectorize_Cop[,3]) levels(vectorize_Cop[,2])<-c("Without error","With error") levels(vectorize_Cop[,3])<-c("Copula method") colnames(vectorize_Cop)<-c("P_value","Error","Method") vectorize<-NULL vectorize<-rbind(vectorize_Disc,vectorize_Cont,vectorize_Cop) } ### Plot pour la variable 1 : pregnancie ggplot(Vectorize(1) ,aes(x=Error, y=P_value,fill=Method)) + geom_jitter(position=position_jitter(width=0.3), aes(colour=Method), alpha=0.9) + geom_boxplot(alpha = 0.5, show.legend = FALSE, aes(fill=Method)) + facet_grid(.~Method, labeller = as_labeller(Vectorize(1)$Method)) + theme(strip.text.x = element_text(size=9, color="black", face="bold")) + labs(y="P-value : Pronostic vs. Pregnancie", x="Error in simulation") + ggtitle(paste("Chi-2 test \n","Pronostic vs. Pregnancy \n","True p-value =",round(Pvalue_vraie[1],5))) + theme(plot.title = element_text(hjust = 0.5)) + scale_fill_manual(values=c("gainsboro", "gray47", "gray0")) + scale_color_manual(values=c("gainsboro", "gray47", "gray0")) # scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) + # scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) ### Plot pour la variable 4 : BloodPressure ggplot(Vectorize(4) ,aes(x=Error, y=P_value,fill=Method)) + geom_jitter(position=position_jitter(width=0.3), aes(colour=Method), alpha=0.9) + geom_boxplot(alpha = 0.5, show.legend = FALSE, aes(fill=Method)) + facet_grid(.~Method, labeller = as_labeller(Vectorize(4)$Method)) + theme(strip.text.x = element_text(size=9, color="black", face="bold")) + labs(y="P-value : Pronostic vs. BloodPressure", x="Error in simulation") + ggtitle(paste("Student test \n","Pronostic vs. BloodPressure \n","True p-value =",round(Pvalue_vraie[4],5))) + theme(plot.title = element_text(hjust = 0.5)) + scale_fill_manual(values=c("gainsboro", "gray47", "gray0")) + scale_color_manual(values=c("gainsboro", "gray47", "gray0")) ### Plot pour la variable 6 : Insulin ggplot(Vectorize(6) ,aes(x=Error, y=P_value,fill=Method)) + geom_jitter(position=position_jitter(width=0.3), aes(colour=Method), alpha=0.9) + geom_boxplot(alpha = 0.5, show.legend = FALSE, aes(fill=Method)) + facet_grid(.~Method, labeller = as_labeller(Vectorize(6)$Method)) + theme(strip.text.x = element_text(size=9, color="black", face="bold")) + labs(y="P-value : Pronostic vs. Insulin", x="Error in simulation") + ggtitle(paste("Student test \n","Pronostic vs. Insulin \n","True p-value =",round(Pvalue_vraie[6],5))) + theme(plot.title = element_text(hjust = 0.5)) + scale_fill_manual(values=c("gainsboro", "gray47", "gray0")) + scale_color_manual(values=c("gainsboro", "gray47", "gray0")) ### Plot pour la variable 7 : DiabetesPedigreeFunction ggplot(Vectorize(7) ,aes(x=Error, y=P_value,fill=Method)) + geom_jitter(position=position_jitter(width=0.3), aes(colour=Method), alpha=0.9) + geom_boxplot(alpha = 0.5, show.legend = FALSE, aes(fill=Method)) + facet_grid(.~Method, labeller = as_labeller(Vectorize(7)$Method)) + theme(strip.text.x = element_text(size=9, color="black", face="bold")) + labs(y="P-value : Pronostic vs. DiabetesPedigreeFunction", x="Error in simulation") + ggtitle(paste("Student test \n","Pronostic vs. DiabetesPedigreeFunction \n","True p-value =",round(Pvalue_vraie[7],5))) + theme(plot.title = element_text(hjust = 0.5)) + scale_fill_manual(values=c("gainsboro", "gray47", "gray0")) + scale_color_manual(values=c("gainsboro", "gray47", "gray0")) #### comparaison des 3 m??thodes (avec erreur) ########### # variable 4 : p-value = 0.0002 # Plot_3methode(Mat_Pvalue_Disc_AE,Mat_Pvalue_Cont_AE,Mat_Pvalue_Cop_AE,4) # cat("True p-value =",Pvalue_vraie[4],"\n") # cat("Median p-value (Discrete method) =",median(Mat_Pvalue_Disc_AE[,4]),"\n") # cat("Median p-value (Continuous method) =",median(Mat_Pvalue_Cont_AE[,4]),"\n") # cat("Median p-value (Copula method) =",median(Mat_Pvalue_Cop_AE[,4]),"\n") # # # variable 6 : p-value < 10^-7 # Plot_3methode(Mat_Pvalue_Disc_AE,Mat_Pvalue_Cont_AE,Mat_Pvalue_Cop_AE,6) # cat("True p-value =",Pvalue_vraie[6],"\n") # cat("Median p-value (Discrete method) =",median(Mat_Pvalue_Disc_AE[,6]),"\n") # cat("Median p-value (Continuous method) =",median(Mat_Pvalue_Cont_AE[,6]),"\n") # cat("Median p-value (Copula method) =",median(Mat_Pvalue_Cop_AE[,6]),"\n") ```