De nombreux systèmes sont modélisés par des équations différentielles ordinaires. R peut permettre de résoudre certains de ces systèmes et aussi d’estimer leurs paramètres.
On prend ici l’exemple d’un modèle épidémiologique temporel SI.
#edo SIX
library(deSolve) # on utilise deSolve
#on définit le système dans une fonction six ici
six<-function(t,x,parms){
with( as.list(c(parms,x)),{
rp<-apexp(-bpt)
rs<-asexp(-0.5(log(t/gs)/bs)^2)
dI<- (rpXS+rsIS)
dS<- -(rpXS+rsIS)
res<-c(dI,dS)
list(res)})
}
#on définit les paramètres pour la simulation du système
parms<-c(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1396, X=1, N=1010)
#on crée un vecteur pour le temps
times<-seq(0:3000)
#valeurs initiales des variables (ici tous les individus sont sains au début)
y<- xstart <-(c(I = 0, S = 1010))
#on résout le système avec la fonction lsoda
out<-as.data.frame(lsoda(xstart, times, six, parms))
# toujours regarder le résultat pour détecter les incohérences
plot(out$S~out$time, type="l",col='blue')
plot(out$I~out$time, type="l",col='green',lwd=3)
#on crée un jeu de données en ajoutant un bruit blanc (gaussien)
noisy<-out$I+rnorm(nrow(out),sd=0.15*out$I)
plot(noisy)
#on ajuste le modèle sur ces données par les moindres carrés
# pour cela on utilise la fonction nls, il faut faire attention aux paramètres initiaux
fitnoisy<-nls(noisy~lsoda(xstart,times,six,c(ap=ap, bp=bp, as=as, bs=bs, gs=gs, X=1, N=1010))[,2]
,start=list(ap=0.002, bp=0.0084, as=5.9e-7, bs=0.25, gs=1000)
, control=list(minFactor=1e-20,maxiter=150))
summary(fitnoisy)
# un beau graphique pour visualiser tout ça!
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)
Vous pouvez définir différents systèmes d’équations différentielles en conservant la même structure.
Enfin pour estimer les paramètres il est possible d’utiliser d’autres fonctions et/ou d’autres méthodes (vraisemblance par exemple).
Pour que nls converge il faut donner des valeurs initiales le plus proche possibles des "vraies valeurs". Je vous conseille donc d’essayer de fitter votre système graphiquement. Si ça ne marche pas vous pouvez utiliser nls2 avec l’algorithme bruteforce…