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 SIXlibrary(deSolve) # on utilise deSolve#on définit le système dans une fonction six icisix<-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èmeparms<-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 tempstimes<-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 lsodaout<-as.data.frame(lsoda(xstart, times, six,Read More →

On cherche souvent à modéliser un échantillon par une loi de probabilité. A partir d’un jeu de données, comment peut-on trouver les paramètres d’une loi préalablement fixée? Plusieurs méthodes peuvent être utilisées. On prend l’exemple ici du délai entre l’infection d’un individu et la détection de cet individu comme malade.On modélise ici ce délai (en jours) par une loi de Weibull (on peut aussi essayer les lois gamma et lognormale par exemple) La méthode la plus simple est d’utiliser la fonction fitdistr du package MASS.Cette fonction permet d’ajuster de nombreuses lois par maximum de vraisemblance. Regardons ce que ça donne pour une loi Weibull. #onRead More →

L’estimateur de Kaplan-Meier donne la fonction de survie non paramétrique.Pour l’obtenir sous R on peut utiliser le package survival. On se place ici dans un cas très simple où il n’y a ni censure ni troncature.Pour bien comprendre le code, je vous conseille vivement de regarder la documentation du package en question!! #survival analysisls()rm(list=ls())library(survival)#on crée un jeu de données correspondant à des durées (étudiées dans l’analyse de survie) z<-c(14,14,14,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,23)d<-data.frame(delay=z)#on crée une colonne status, ici tous les individus sont « morts » pendant l’expérience # mort au sens de l’analyse de survie d$status<-1 s<-survfit(Surv(d$delay,d$status)~1)plot(s,main= »survival function »)Read More →

Pour charger un fichier excel vous avez 3 façons de faire 1 Vous pouvez utiliser Rcmdr library(Rcmdr) # puis importer des données.. ça marche plutôt bien.. mais que sous windows En pratique Rcmdr utilise le package RODBC.Je le trouve peu intuitif en ligne de commande et vous propose d’utiliser le package xlsReadWrite 2 Le package xlsReadWrite library(xlsReadWrite)xls.getshlib() # indispensabledonnee<-read.xls(« data.xls ») Par contre cela ne marche jamais vraiment parfaitement bien lorqu’il y a des onglets ou pour des cas un peu spéciaux. Il existe une version pro, payante de ce package qui utilise des fichiers binaires propriétaires, mais il serait dommage d’utiliser cela. J’en arrive donc auRead More →

# on choisit le dossier dans lequel on veut sauver le fichiersetwd(« U:/simulations »)# on crée un jeu de données, #ici on simule une densité de probabilité d’une loi normale de moyenne 2 et d’écart type 0.5v<-seq(0,7,by=0.1)a<-dnorm(v,mean=2,sd=0.5)plot(v,a) # on vérifie visuellement#on crée un dataframe avant de le sauverdata<-data.frame(TL=v,F=a)# on utilise la fonction write.table, voir ?write.table# ne pas oublier le .csv à la fin du nom du fichier excel « data.csv »write.table(data, « data.csv », row.names=FALSE, sep= »t »,dec= », », na= » « )Read More →

Il existe plusieurs façons de faire un graphique avec deux ordonnées. En voici une qui utilise les outils graphiques de base # Données d’exemple (peu importe…)times<-seq(0,3000)p<- 0.002197451 * exp(- 0.0009076665 *times)b1<- 7.376812e-08b2<-0.2652811b3<- 1986.235s<-b1*exp(-0.5*(log(times/b3)/b2)^2)# On ouvre une nouvelle fenêtre plot.new()# On choisit les paramètres de la fenêtre, voir ?par, ici mar correspond aux marges par(mar=c(5,4,3,4))# On met le premier graphique en définissant les limites des axesplot.new()plot.window(xlim=c(0,3000),ylim=c(0,0.0022))lines(p~times,type=’l’,col=’burlywood1′,lwd=3)# on ajoute les axes et leurs légendesaxis(1)axis(2)title(xlab= »time »)title(ylab= »rp »)# On superpose le graphique avec un axe des ordonnées différent qui sera à droite du graphique (axis(4))plot.window(xlim=c(0,3000),ylim=c(0,8e-08))lines(s~times,type=’l’,col=’burlywood3′,lwd=3)axis(4)#titre du graphiquetitle(main= »force of infection »)#légende de l’ordonnée n°2mtext(« rs »,side=4,line=2.5)#on termine le graphiquebox()Read More →

Quand on travaille sur un grand jeu de données, on peut voir à quoi il ressemble en rentrant son nom directement dans la console : data(iris) iris Mais pour de gros objets ce n’est pas du tout pratique… et souvent votre ordinateur peut avoir du  mal à tout vous afficher. Utilisez plutôt : data(iris) View(iris) # pensez bien au V majuscule Cela ouvre une nouvelle fenêtre avec l’intégralité du jeu de données. Et vous pouvez facilement vous balader dedans. Par contre si seul l’aspect général du tableau vous intéresse, vous pouvez utiliser la fonction head : data(iris) head(iris) Elle vous affiche les 5 premières lignesRead More →

i<-0 while (i<10){ print(i) i<-i+1 } print(« on sort de la boucle ») While va réaliser ce qui est écrit entre les accolades {} tant que ce qui est dans les parenthèses () est vrai. Ce code peut donc se traduire de la façon suivante : i vaut 0 tant que (i est inférieur à 0){ afficher i augmenter i de 1 } Dans cet exemple au moment où i vaudra 10, on sortira de la boucle.Read More →

for ( i in 1:10) { print(i) } Cette commande peut se traduire par : Pour (i allant de 1 à 10) { affiche i} Il faut noter que les parenthèses () servent à définir la variable et les valeurs qu’elle va prendre successivement à chaque tour de boucle. Les accolades {} servent à délimiter les actions à effectuer pour chacune des valeurs prises par la variable. IMPORTANT : R n’aime pas vraiment les boucles for, il est beaucoup plus efficace d’utiliser apply. Tout particulièrement pour les très grandes et longues boucles, apply fait cela en une fraction de seconde… alors que for peut mettreRead More →

Ce message apparaît lorsque l’on essaye de lancer un script initialement écrit sous windows sur une station linux (et inversement). Il s’agit d’un souci d’encodage. il est possible de corriger cela grâce au script suivant (remplacer latin1 par l’encodage d’origine du fichier) options(encoding= »latin1″) Ou alors options(encoding= »utf-8″)Read More →