De nombreux syst\u00e8mes sont mod\u00e9lis\u00e9s par des \u00e9quations diff\u00e9rentielles ordinaires. R peut permettre de r\u00e9soudre certains de ces syst\u00e8mes et aussi d’estimer leurs param\u00e8tres.<\/p>\n
On prend ici l’exemple d’un mod\u00e8le \u00e9pid\u00e9miologique temporel SI.<\/p>\n
#edo SIX
library(deSolve) # on utilise deSolve
#on d\u00e9finit le syst\u00e8me dans une fonction six ici
six<-function(t,x,parms){
\u00a0\u00a0 \u00a0with( as.list(c(parms,x)),{
\u00a0\u00a0 \u00a0rp<-ap*exp(-bp*t)
\u00a0\u00a0 \u00a0rs<-as*exp(-0.5*(log(t\/gs)\/bs)^2)
\u00a0\u00a0 \u00a0
\u00a0\u00a0 \u00a0dI<- (rp*X*S+rs*I*S)
\u00a0\u00a0 \u00a0dS<- -(rp*X*S+rs*I*S)
\u00a0\u00a0 \u00a0res<-c(dI,dS)
\u00a0\u00a0 \u00a0list(res)})
}
#on d\u00e9finit les param\u00e8tres pour la simulation du syst\u00e8me
parms<-c(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1396, X=1, N=1010)
#on cr\u00e9e un vecteur pour le temps
times<-seq(0:3000)
#valeurs initiales des variables (ici tous les individus sont sains au d\u00e9but)
y<- xstart <-(c(I = 0, S = 1010))
#on r\u00e9sout le syst\u00e8me avec la fonction lsoda
out<-as.data.frame(lsoda(xstart, times, six, parms))
# toujours regarder le r\u00e9sultat pour d\u00e9tecter les incoh\u00e9rences
plot(out$S~out$time, type=\"l\",col='blue')
plot(out$I~out$time, type=\"l\",col='green',lwd=3)
#on cr\u00e9e un jeu de donn\u00e9es en ajoutant un bruit blanc (gaussien)
noisy<-out$I+rnorm(nrow(out),sd=0.15*out$I)
plot(noisy)
#on ajuste le mod\u00e8le sur ces donn\u00e9es par les moindres carr\u00e9s
# pour cela on utilise la fonction nls, il faut faire attention aux param\u00e8tres initiaux
fitnoisy<-nls(noisy~lsoda(xstart,times,six,c(ap=ap, bp=bp, as=as, bs=bs, gs=gs, X=1, N=1010))[,2]
\u00a0\u00a0 \u00a0 ,start=list(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1000)
\u00a0\u00a0 \u00a0, control=list(minFactor=1e-20,maxiter=150))
summary(fitnoisy)
# un beau graphique pour visualiser tout \u00e7a!
plot(noisy,pch=20,col='grey',main=\"SIX ode Model Fit\",xlab=\"time\",ylab=\"disease progress\")
lines(predict(fitnoisy,times),type='l',col='green',lwd=4.5)
<\/code><\/pre>\nVous pouvez d\u00e9finir diff\u00e9rents syst\u00e8mes d’\u00e9quations diff\u00e9rentielles en conservant la m\u00eame structure.
Enfin pour estimer les param\u00e8tres il est possible d’utiliser d’autres fonctions et\/ou d’autres m\u00e9thodes (vraisemblance par exemple).
Pour que nls converge il faut donner des valeurs initiales le plus proche possibles des \u00ab\u00a0vraies valeurs\u00a0\u00bb. Je vous conseille donc d’essayer de fitter votre syst\u00e8me graphiquement. Si \u00e7a ne marche pas vous pouvez utiliser nls2 avec l’algorithme bruteforce…<\/p>\n","protected":false},"excerpt":{"rendered":"
De nombreux syst\u00e8mes sont mod\u00e9lis\u00e9s par des \u00e9quations diff\u00e9rentielles ordinaires. R peut permettre de r\u00e9soudre certains de ces syst\u00e8mes et aussi d’estimer leurs param\u00e8tres. On prend ici l’exemple d’un mod\u00e8le \u00e9pid\u00e9miologique temporel SI. #edo SIXlibrary(deSolve) # on utilise deSolve#on d\u00e9finit le syst\u00e8me dans une fonction six icisix<-function(t,x,parms){\u00a0\u00a0 \u00a0with( as.list(c(parms,x)),{\u00a0\u00a0 \u00a0rp<-ap*exp(-bp*t)\u00a0\u00a0 \u00a0rs<-as*exp(-0.5*(log(t\/gs)\/bs)^2)\u00a0\u00a0 \u00a0\u00a0\u00a0 \u00a0dI<- (rp*X*S+rs*I*S)\u00a0\u00a0 \u00a0dS<- -(rp*X*S+rs*I*S)\u00a0\u00a0 \u00a0res<-c(dI,dS)\u00a0\u00a0 \u00a0list(res)})}#on d\u00e9finit les param\u00e8tres pour la simulation du syst\u00e8meparms<-c(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1396, X=1, N=1010)#on cr\u00e9e un vecteur pour le tempstimes<-seq(0:3000)#valeurs initiales des variables (ici tous les individus sont sains au d\u00e9but)y<- xstart <-(c(I = 0, S = 1010))#on r\u00e9sout le syst\u00e8me avec la fonction lsodaout<-as.data.frame(lsoda(xstart, times, six,Read More →<\/a><\/p>\n","protected":false},"author":6,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"content-type":"","rop_custom_images_group":[],"rop_custom_messages_group":[],"rop_publish_now":"initial","rop_publish_now_accounts":{"twitter_399453572_399453572":""},"rop_publish_now_history":[],"rop_publish_now_status":"pending","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[14],"tags":[],"class_list":{"0":"entry","1":"post","2":"publish","3":"author-melen","4":"post-564","6":"format-standard","7":"category-modelisation"},"acf":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p9O7Sx-96","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/564","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/users\/6"}],"replies":[{"embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/comments?post=564"}],"version-history":[{"count":2,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/564\/revisions"}],"predecessor-version":[{"id":4096,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/posts\/564\/revisions\/4096"}],"wp:attachment":[{"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/media?parent=564"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/categories?post=564"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/thinkr.fr\/abcdr\/wp-json\/wp\/v2\/tags?post=564"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}