Comment estimer les paramètres d'un modèle d'équations différentielles ordinaires par les moindres carrés avec R ?

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<-ap*exp(-bp*t)
    rs<-as*exp(-0.5*(log(t/gs)/bs)^2)
    
    dI<- (rp*X*S+rs*I*S)
    dS<- -(rp*X*S+rs*I*S)

    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…