featuring : ajout d'un script pour identifiabilité de alpha

refactor : changement du script identifiabilité (alpha et mu) pour mieux recup les résultats

src/identif_alpha.R

+source("generdata.R")
+library(doParallel)
+
+# Identifiabilité pratique
+
+
+
+control <- saemix::saemixControl(map = FALSE,fix.seed = T, fim = FALSE, ll.is = FALSE, nbiter.saemix = c(300, 50), nb.chains = 1, nbiter.burn = 5, nbiter.mcmc = c(2, 2, 2, 0), displayProgress = FALSE, print= FALSE, save = FALSE,  nbdisplay = 50)
+
+alpha_pop = exp(-6.34)
+mu_pop = exp(-26.8)
+sigma = 0.542
+omega_alpha = 1
+omega_mu = 0
+
+true_v = c(alpha_pop,mu_pop,sigma,omega_alpha,omega_mu)
+
+
+
+alpha <- c()
+mu <- c()
+omega_alpha <- c()
+omega_mu <- c()
+
+Ncpus <- parallel::detectCores() - 1
+cl <- parallel::makeCluster(Ncpus)
+doParallel::registerDoParallel(cl)
+
+res <- foreach::foreach(i=1:25 , .packages = c("tidyverse","mechanisticModel"))%dopar%{
+  
+  set.seed(i)
+  
+  control$seed <- i
+  
+  JD1 <- GDC(1000,22.38,1.33,exp(-6.34),exp(-26.8),0.542,1,0,20*365)
+  JD1 <- JD1 %>% filter(tps > 0)
+  
+  
+  m = mechanisticModel(
+    Surv(tps, event) ~ TTR(vdiag),
+    log(alpha) ~ 1, log(mu) ~ 1,
+    data = JD1,
+    alpha_0 = 0.001, mu_0 = 2e-12, sigma_0 = 0.5, omega_0 = c(3, 3),
+    beta_0 = list(alpha = numeric(), mu = numeric()),
+    options = control
+  ) 
+  m = fit(m, silent = F, N_sim_surv = 0)
+}
+
+
+parallel::stopCluster(cl)
+
+
+for(i in 1:25){
+  alpha <- append(alpha,unlist(res[[i]]$coefs[1])[1])
+  mu <- append(mu,unlist(res[[i]]$coefs[1])[2])
+  omega_alpha <- append(omega_alpha,unlist(res[[i]]$coefs[4])[1])
+  omega_mu <- append(omega_mu,unlist(res[[i]]$coefs[4])[4])
+}
+
+
+
+
+moy <- c(mean(alpha),mean(mu),mean(omega_alpha),mean(omega_mu)) 
+sd <- c(sd(alpha),sd(mu),sd(omega_alpha),sd(omega_mu)) 
+
+
+
+RSE_alpha <- 100*(sd[1]/moy[1])
+RSE_mu <- 100*(sd[2]/moy[2])
+RSE_omega_alpha <- 100*(sd[3]/moy[3])
+RSE_omega_mu <- 100*(sd[4]/moy[4])
+
+
+save(res,true_v,alpha,mu,omega_alpha,omega_mu,moy,sd,RSE_alpha,RSE_mu,RSE_omega_alpha,RSE_omega_mu,file = "identifi_alpha.Rdata")
\ No newline at end of file

src/identifiabilite.R

@@ -26,7 +26,7 @@ Ncpus <- parallel::detectCores() - 1
 cl <- parallel::makeCluster(Ncpus)
 doParallel::registerDoParallel(cl)
 
-res <- foreach::foreach( i = 1:25 , .packages = c("tidyverse","mechanisticModel"))%dopar%{
+res <- foreach::foreach(i=1:25 , .packages = c("tidyverse","mechanisticModel"))%dopar%{
   
   set.seed(i)
   
@@ -51,13 +51,11 @@ res <- foreach::foreach( i = 1:25 , .packages = c("tidyverse","mechanisticModel"
 parallel::stopCluster(cl)
 
 
-
-
 for(i in 1:25){
-  alpha <- append(alpha,unlist(res[i]$coefs[1])[1])
-  mu <- append(mu,unlist(res[i]$coefs[1])[2])
-  omega_alpha <- append(omega_alpha,unlist(res[i]$coefs[4])[1])
-  omega_mu <- append(omega_mu,unlist(res[i]$coefs[4])[4])
+  alpha <- append(alpha,unlist(res[[i]]$coefs[1])[1])
+  mu <- append(mu,unlist(res[[i]]$coefs[1])[2])
+  omega_alpha <- append(omega_alpha,unlist(res[[i]]$coefs[4])[1])
+  omega_mu <- append(omega_mu,unlist(res[[i]]$coefs[4])[4])
 }