Uncategorized (Page 19)

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<-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è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<-b1exp(-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) irisRead More →

i<-0 while (i<10){ print(i) i<-i+1 } print(« on sort de la boucle »)Read More →

for ( i in 1:10) { print(i) }Read 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 →