IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Voir le flux RSS

Professeur Leger (l. g...)

un peu de modélisation pour enrichir la discussion

Noter ce billet
par , 29/04/2020 à 16h25 (243 Affichages)
Je ne suis pas épidémiologiste, et loin de moi l'idée que le modèle proposé permette d'effectuer des prévisions.

Voici pour résumer le code R nécessaire pour afficher la jolie courbe en pièce jointe.

D'abord, on établit les données cumulées pour le nombre de cas, en France (faire copier/coller du code source)

# gestion du coronavirus avec R
cat("coronavirus:\n");
cat("edit(file=\"coronavirus.R\");\n");
# modification de l'option max.print de façon à permettre l'affichage de la totalité du fichier csv
options(max.print="200000");
donnees_corona<-read.csv("coronavirus_data.csv");

ls(donnees_corona);
typeof(donnees_corona$dateRep);
str(donnees_corona);

donnees_corona<-read.csv("coronavirus_data.csv",colClasses=c("dateRep"="character"));
donnees_corona$dateRep2<-as.Date(donnees_corona$dateRep,"%d/%m/%y");
# affichage des dates pour la France
donnees_corona$dateRep2[which(donnees_corona$countriesAndTerritories=="France")];

donnees_france<-data.frame(seq(donnees_corona$dateRep[which(donnees_corona$countriesAndTerritories=="France")]));

l0<-length(donnees_corona$dateRep[which(donnees_corona$countriesAndTerritories=="France")])
donnees_france<-data.frame(seq(l0))
donnees_france$dateRep<-donnees_corona$dateRep[which(donnees_corona$countriesAndTerritories=="France")];
donnees_france$cases<-donnees_corona$cases[which(donnees_corona$countriesAndTerritories=="France")];
donnees_france$death<-donnees_corona$deaths[which(donnees_corona$countriesAndTerritories=="France")];

edit(donnees_france);

donnees_france$num<-seq(length(donnees_france$death),by=-1);

edit(donnees_france);




plot.new();
plot(donnees_france$num,donnees_france$cases,col="red",main="nb de cas");
lines(donnees_france$num,donnees_france$death,col="black",main="nb de deces");


# representation graphique de donnees cumuleees

gl_reverse <-function(X) {
# X colonne,
j=length(X)
gl_r<-c(0)
for(i in 1:j)
{gl_r[j-i+1]<-X[i]}
gl_r
};

gl_cumul <-function(X) {
# X colonne,
j=length(X)
gl_r<-c(0)
gl_r[1]=0
for(i in 2:j)
{gl_r[i]<-gl_r[i-1]+X[i]}
gl_r
};

cumul_cas_france<-gl_cumul(gl_reverse(donnees_france$cas))
cumul_cas_france

plot(cumul_cas_france,main="coronavirus:nombre de cas cumules")

plot(cumul_cas_france,
main="coronavirus: nombre de cas cumules",
xlab="nb de jours",
ylab="nb de cas"
);


#----------------------------------------------------------------------------------------------------------------------------------------------------------------------
# calcul de la fonction logistique
#----------------------------------------------------------------------------------------------------------------------------------------------------------------------

# solution
solution_logistic<-function(t,a,b,c,d){
sol<-a+b/(1+exp(-(t-c)/d))
sol
}

# on reprend l'exemple pour l'Espagne :
# logis.m1 <- nls(Casos ~ logis(dia, a, b, c,d), data = agregados_por_fecha, start = list(a = 0, b = 180000, c = 40, d=5))

# on générera une séquence T croissante
n_size=length(cumul_cas_france)
T=seq(n_size)
logis.m1 <- nls(cumul_cas_france ~ solution_logistic(T, a,b, c,d), start = list(a = 0, b = 180000, c = 40, d=5))
# erreur le pas 0.000488281 est devenu inférieur à 'minFactor' de 0.000976562
# on fait varier a,b,c ... pour reduire le pb
logis.m1 <- nls(cumul_cas_france ~ solution_logistic(T, a,b, c,d), start = list(a = 10, b = 380000, c = 80, d=5))
summary(logis.m1)
coef(logis.m1)
solution_logistic (T,-340.241381,123116.928901,95.784888,6.658315)

curve(solution_logistic (x,-340.241381,123116.928901,95.784888,6.658315),from=1,to=210,
add=FALSE)
lines(cumul_cas_france,col="blue")

#-----------------------------------------------------------------------------------------------------------------------------------------------------
# exporter le graphique obtenu avec R
#-----------------------------------------------------------------------------------------------------------------------------------------------------

# la documentation des fonctions de génération de fichiers graphiques est obtenue en tapant ces lignes
library(help ="grDevices")
help(png)

# jpeg(filename = "Rplot%03d.jpeg",
# width = 480, height = 480, units = "px", pointsize = 12,
# quality = 75,
# bg = "white", res = NA, ...,
# type = c("cairo", "Xlib", "quartz"), antialias)

# en appliquant la documentation on obtient le fichier image coronavirus-graphe.jpg ... disponible après dev.of()
jpeg(file = "coronavirus-graphe.jpg",quality=80,bg="grey")
curve(solution_logistic (x,-340.241381,123116.928901,95.784888,6.658315),from=1,to=210,
add=FALSE)
lines(cumul_cas_france,col="blue")
dev.off()
Miniatures attachées Images attachées  

Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Viadeo Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Twitter Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Google Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Facebook Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Digg Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Delicious Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog MySpace Envoyer le billet « un peu de modélisation pour enrichir la discussion » dans le blog Yahoo

Commentaires