In alex-conanec/optisure:

knitr::opts_chunk$set(echo = TRUE)

Simulation

library(tidyverse)
library(MOOVaR)
library(latex2exp)

path = "~/Seafile/PLMbox/These/Redaction/manuscript_these/figure_manuscript/"
set.seed(1234)

n = 300
d1 = 11
d2 = 2
n_lev = 2:4
p = 3

if (d2 >= 2){
  lev = sample(n_lev, d2, replace = TRUE)
}else{
  lev = c(2, 5)
}


quali_design = purrr::cross(
  lapply(lev, seq, from = 1)) %>%
  lapply(function(x){
    df = as.data.frame(x)
    colnames(df) = paste0("X", 1:NCOL(df))
    df
  }) %>% bind_rows() %>% mutate_if(is.integer, as.factor)

beta = matrix(runif(n = p*(d1 + sum(lev)), -1, 1), ncol = p)


f = function(X, beta){
  is_fac = sapply(X, is.factor)
  if (any(is_fac)){
    X_dummy = lapply(which(is_fac), function(j){
      as.data.frame(model.matrix(~. + 0, data = X[,j, drop = FALSE]))
    }) %>% bind_cols()

    X = purrr::quietly(bind_cols)(X[!is_fac], X_dummy)$result
  }

  as.matrix(X) %*% beta

}
formals(f)$beta = beta


X_quanti = matrix(runif(n = n * d1), ncol = d1)
colnames(X_quanti) = paste0("X", 1:NCOL(X_quanti))
idx = sample(1:NROW(quali_design), n, replace = TRUE)
X_quali = quali_design[idx,]
colnames(X_quali) = paste0("X", (NCOL(X_quanti)+1):(NCOL(X_quanti) + NCOL(X_quali)))
X = cbind(X_quanti, X_quali)
Y = f(X)

#3 scenario
rr = c(-0.5, 0, 0.5)
names(rr) = rr
tau_list = c(0.05, 0.15, 0.3, 0.5)
names(tau_list) = tau_list
X = X %>% mutate(id=1:NROW(.)) %>% select(id, everything())
res_all = lapply(rr, function(r){
  C = matrix(r, 3, 3)
  diag(C) = 1
  epsilon = mvtnorm::rmvnorm(n=n, mean=rep(0, 3), sigma = C, method="chol")
  cor(epsilon)

  Y = Y + epsilon

  res_list = lapply(tau_list, function(tau){
    lapply(c(FF=FALSE, TT=TRUE), function(globale_tau){
      MOOVaR(X, Y, tau = rep(tau, NCOL(Y)),
             globale_tau = rep(globale_tau, NCOL(Y)),
             X_space_csrt = FALSE,
             optim_method = "nsga",
             TT = 20, sens = rep("max", p))
    })
  })

  #calcul d'un HI mais avec le meme espace de MC
  N_MC=30^p

  Y_all_front = lapply(tau_list, function(tau){
    lapply(c(FF="FF", TT="TT"), function(globale_tau){
      res_list[[tau]][[globale_tau]]$Y
    }) %>% bind_rows()
  }) %>% bind_rows()

  Y_MC = sapply(1:p, function(j){
    runif(N_MC,
          min = min(Y_all_front[,j]),
          max = max(Y_all_front[,j]))
  })

  HI_esti = lapply(tau_list, function(tau){
    tau = as.character(tau)
    lapply(c(FF="FF", TT="TT"), function(globale_tau){
      HI(Y_front = res_list[[tau]][[globale_tau]]$Y,
         Y_MC = Y_MC, A_MC=1, N_MC = N_MC)
    })
  })

  list(res_list = res_list,
       HI_esti = HI_esti)

})
# saveRDS(res_all, "simu/res_tau_ajustement.RDS")
# res_all = readRDS("simu/res_tau_ajustement.RDS")

lapply(rr, function(r){
  r = as.character(r)
  lapply(tau_list, function(tau){
    tau = as.character(tau)
    a = lapply(c("FF", "TT"), function(globale_tau){
      mu = apply(res_all[[r]][[1]][[tau]][[globale_tau]]$Y, 2, mean) %>% round(2)
      p_value = sapply(1:p, function(j){
        t.test(res_all[[r]][[1]][[tau]][["FF"]]$Y[,j],
               res_all[[r]][[1]][[tau]][["TT"]]$Y[,j])$p.value
      }) %>% signif(1)
      p_value = ifelse(p_value>0.05, p_value,
                       ifelse(p_value>0.01, "*",
                              ifelse(p_value>0.001, "**", "***")))


      if (globale_tau == "TT"){
        p_value = NA
      }

      cbind(
        `$1 - \\tau$` = 1 - as.numeric(tau),
        `Méthode` = globale_tau,
        `$\\Phi$` = res_all[[r]][[1]][[tau]][[globale_tau]]$globale_confiance %>% round(2),
        HI = as.data.frame(res_all[[r]][[2]][[tau]][[globale_tau]] %>% round(2) %>%
                             as.character()%>% t() ),
        as.data.frame(mu %>% round(2) %>%
                        as.character()%>% t() ),
        as.data.frame(res_all[[r]][[1]][[tau]][[globale_tau]]$tau_j %>%
                        round(2) %>%
                        as.character()%>% t() ),
        as.data.frame(p_value %>% t() ))

    }) %>% bind_rows()

    levels(a$Méthode) = c("Non ajustée", "Ajustée")
    colnames(a)[4] = "HI"
    colnames(a)[5:7] = paste0("$\\Bar{y}^{(", 1:3, ")}$")
    colnames(a)[8:10] = paste0("$\\tau + \\lambda^{(", 1:3,")}$")
    colnames(a)[11:13] = paste0("p_{value}")

    cbind(r = round(as.numeric(r), 1), a[,c(1:5, 11, 8, 6, 12, 9, 7, 13, 10)])

  }) %>% bind_rows()

}) %>% bind_rows() %>%
  xtable::xtable(type = "latex") %>%
  print(sanitize.text.function=identity, include.rownames=FALSE)

df = lapply(rr, function(r){
  r = as.character(r)
  lapply(tau_list, function(tau){
    tau = as.character(tau)
    lapply(c("FF", "TT"), function(globale_tau){
      data.frame(
        r = r,
        un_moins_tau = 1 - as.numeric(tau),
        Méthode = ifelse(globale_tau =="FF", "Non ajustée", "Ajustée"),
        HI = res_all[[r]][[2]][[tau]][[globale_tau]] %>% round(2),
        tau_moy = mean(res_all[[r]][[1]][[tau]][[globale_tau]]$tau_j),
        y_moy_agg = mean(apply(res_all[[r]][[1]][[tau]][[globale_tau]]$Y, 2, mean)),
        y_sd_agg = mean(apply(res_all[[r]][[1]][[tau]][[globale_tau]]$Y, 2, sd))
        )
    }) %>% bind_rows()

  }) %>% bind_rows()

}) %>% bind_rows()


ggplot(df) +
  geom_point(aes(x = un_moins_tau, y = HI, colour = r, shape = Méthode)) +
  geom_line(aes(x = un_moins_tau, y = HI, colour = r, linetype = Méthode)) +
  xlab(TeX("$1-\\tau$")) +
  theme_bw()

# ggsave(filename = paste0(path, "HI_simu_tau.png"))

ggplot(df %>% filter(Méthode == "Ajustée")) +
  geom_point(aes(x = un_moins_tau, y = tau_moy, colour = r), shape = 2) +
  geom_line(aes(x = un_moins_tau, y = tau_moy, colour = r), linetype = 2) +
  xlab(TeX("$1-\\tau$")) +
  ylab(TeX("$1/n \\sum_j \\tau + \\lambda^{(j)}$")) +
  theme_bw()
# ggsave(filename = paste0(path, "penalite_simu_tau.png"))

alex-conanec/optisure documentation built on Dec. 19, 2021, 12:27 a.m.