(gH) -_- pot-pourri_r.txt ; TimeStamp (unix) : 24 Mai 2012 vers 15:04 ################################################################# # # # quelques "trucs" en R # # # ################################################################# # interaction utilisateur reponse <- readline(question) user.prompt # package Zelig par(ask=TRUE) # pour passer manuellement d'un graphique à l'autre # calcul et affichage dans la foulée : mettre des parenthèses ! x <- rnorm(10) # on ne voit rien ( x <- rnorm(10) ) # calcule et affiche # création automatique de noms de lignes rownames(x) <- paste("LIG",sprintf("%04d",1:(dim(x)[1])),sep="") # rownames est dans statgh.r, cela me parait plus simple que # row.names et cohérent avec colnames # affichage en colonnes des noms de colonnes : cbind( names(leDataFrame) ) # et pour les avoir par ordre alphabétique : # (mais je conseile plutot ma fonction lesColonnes !) cbind( sort(names(leDataFrame) )) # passage du format long au format court et reciproquement a2dm <- lit.dar("a2dm.dar") # lit.dar est dans statgh.r, c'est un read.table avec les "bonnes" options qt <- a2dm[,1] ql <- a2dm[,2] a2dmw <- as.matrix(unstack(a2dm)) a2dml <- stack( a2dmw) # conversion en format large (plus lisible) nbv <- mean(table(ql)) donw <- matrix(qt,nrow=nbv,byrow=FALSE) colnames(donw) <- levels(ql) row.names(donw) <- paste("suj",su[1:nbv]) ## changements des options de configuration options() # toutes les valeurs des options options()$width # juste celle de width # changement de la largeur d'affichage options(width=450) # changement du nombre de lignes d'affichage par commande options(max.print=5000) # changement des informations affichées pour connexion internet options(internet.info=0) # les informations sont renvoyées par warnings() # colorier un rectangle rect(0,0,3,100,col=rgb(red=180/256,green=210/256,blue=190/256)) # attendre quelques millisecondes sleep.milli(temps) # dans le package RandomFields # couleur par block de lignes coldc <- rep(c("red","green","blue","black","purple"),each=dim(donw)[1]) # un dotchart quand le beanplot produit une erreur via try bp <- try(beanplot(qt~ql,main=tit,col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6"),silent=TRUE) if (is.character(bp)) { bpp <- tolower(bp[1]) if (nchar(bpp)>0) { if (grep("error",bpp)) { coldc <- rep(c("red","green","blue","black","purple"),each=dim(donw)[1]) dotchart(donw,labels="",col=coldc) } ; # fin si } ; # fin si } ; # fin si # une boite à moustaches améliorée bam <- function(x,titre) { boxplot(x,main=titre,col="yellow") mx <- mean(x) ex <- sd(x) abline(h=median(x),col="green",lwd=2) abline(h=mx,col="red",lwd=2) abline(h=mx-ex,col="blue",lwd=1) abline(h=mx+ex,col="blue",lwd=1) } # fin de fonction bam # élimination des modalités non présentes pour une QL # après sélection de certaines valeurs newql <- factor(as.factor(ql)) # transformation de chaines paste(formatC(valc,format="f",width=3,dig=0)," %") paste(paste(chMod,"!"),collapse="") strtrim(x, width) chartr(old, new, x) casefold(x, upper = FALSE) nomficg <- paste("graph",i,".png",sep="") # toutes les qt d'un .dbf via statgh.r allqtdbf("chu825qt",chaineEnListe("G/l UI/l mmo/l mug/l % mg/dl an")) # connaissez-vous Sweave ? Sweave("demosweave.nw") # ceci est mon programme R # exemple de fichier LaTeX avec du R-Sweave : \begin{document} ~ \SweaveOpts{echo=true} \begin{center} {\LARGE PRODUCTION DYNAMIQUE DE \\[0.5cm] \textsc{DOCUMENTS,} \\[0.5cm] \textsc{STATISTIQUES} \\[0.5cm] \textsc{ET GRAPHIQUES} \\[3cm] }% \end{center} <>= dos <- dosInit # par exemple "iris" nlu <- nluInit # par exemple 50 usr <- usrInit# par exemple "GH" source("http://forge.info.univ-angers.fr/~gh/demosweave.r") nbl <- demoSw(dos,nlu,usr,1) @ \section{Quelques lignes de données du fichier : \Sexpr{dos } Le fichier elf contient \Sexpr{nbl} lignes. L'utilisateur \Sexpr{usr} veut traiter les \Sexpr{nlu} premières lignes de ce fichier. Cela produit donc \Sexpr{nlu*(nlu-1)/2} couples de lignes. <<>>= demoSw(dos,nlu,usr,2) @ \section{Quelques graphiques pour le fichier : iris } <>= demoSw(dos,nlu,usr,4) @ # fin de demo Sweave ################################################################# # # gestion de packages et data # ################################################################# install.packages("ape",dependencies=TRUE) # avec Sys.setlocale("LC_MESSAGES","C") ?? library(ape) data(package="ape") detach(package:ape) update.packages(ask=FALSE) # pour la version 12, avant install.packages : Sys.setlocale("LC_MESSAGES","C") ################################################################# # # courbes ROC et AUROC # ################################################################# library(limma,warn.conflicts=FALSE) auROC(qlbin,qt) library(ROCR) pred <- prediction(qt, qlbin) pred <- prediction(fr$score, fr$label) perf <- performance(pred, "tpr", "fpr") plot(perf) perf <- performance(pred, "auc") auRoc <- perf@y.values[[1]] detach(package:limma) ################################################################# # # corrélation intraclasse # ################################################################# library(psych) scores <- cbind(scorei1,scorei2) intraclass <- ICC(scores) cat("Coefficients de corrélation intraclasses (Shrout and Fleiss) \n") print(intraclass) ################################################################# # # corrélation/concordance de Lin # ################################################################# library(epiR) linc <- epi.ccc(scorei1,scorei2) print(names(linc)) cat("Concordance de Lin : ",linc$rho.c[,1]," intervalle de confiance [",linc$rho.c[,2],",",linc$rho.c[,3],"]\n") ################################################################# # # liste des symboles pour plot ou pchShow de statgh.r # ################################################################# example(pch) # 19 et 20 pour des ronds (gros et petits) ################################################################# # # production automatique d'un fichier de résultats (dans statgh.r) # ################################################################# source("test.r") # pour mettre au point sinksource("test.r") # produit test.sor, plus généralement un .sor de meme nom source("a.r",encoding="latin1") sink("a.sor",split=TRUE) ################################################################# ################################################################# # # exemples de session R pour un cours de L2 (initiation à R) # ################################################################# ################################################################# q() x <- 1 x+2 print(x) cat(" valeur = ", x , "\n") y <- 0:10 y*y y**4 length(y) sum(y) mean(y) var(y) sd(y) cv <- function(x) { sd(x)/mean(x) } cv(y) help.start() rbinom(10,5,0.5) sim <- function(n) { rbinom(n,5,0.5) } sim(100) table( sim(500) ) summary( y ) fivenum( y ) str( y ) ls() mode(x) dim(x) attributes(x) class(x) dimnames(x) is.na(x) paste("I",1:5,sep="") paste("I",1:5,sep="",collapse="+") source("k:/stat_ad/statgh.r") elfdata <- read.dbf("d:/elf.dbf")$dbf attach(elfdata) head(elfdata) tail(elfdata,n=5) str( elfdata ) ag <- AGE decritQT(" AGE ", ag, "ans ") decritQT(" AGE ", ag, "ans ", TRUE) decritQT(" AGE ", ag, "ans ", TRUE ," d:/elf_age.png") sx <- SEXE table(sx) table(sx)/sum(sx) round( table(sx)/sum(sx) ) decritQL(" SEXE ",sx, "homme femme") decritQL(" SEXE ",sx, "homme femme",TRUE) detach(elfdata) agehom <- ag[ sx==0 ] agehomjeune <- ag[ sx==0 & ag < 30 ] RSiteSearch("table") head(ag) head(ag,n=10) tail(ag) mean(ag) sort(ag) order(ag) stem(ag) plot(ag) plot(ag,main="Tracé de l'AGE (dossier ELF) ") plot(ag,pch=19,col="blue") boxplot(ag) boxplot(ag ~ sx) png("elf_sexe.png",width=800,height=600) coul <- c("blue","red") csex <- coul[1+sx] plot(1:length(ag),sort(ag),xlab="individu",ylab="age", col=csex,bg=csex,pch=21) text(-10,0.90,pos=4,"couleurs pour sx ") text(-05,0.85,pos=4,"homme=bleu femme=rouge") dev.off() nomficg <- paste("graph",i,".png",sep="") vinsdata <- read.dbf("d:/vins.dbf")$dbf attach(vinsdata) vinsdata <- vinsdata[ , -1 ] cbind( mean(vinsdata), sd(vinsdata) ) plot(RFA) par(new=TRUE) points(CANADA) cor(x,y) lm(y~x) anova( lm(y~x) ) summary(vins) cor( vins ) pairs( vins ) data() data(iris) demo("graphics") plot(ag) plot(ag,type="p") plot(ag,type="p",pch=19) plot(sort(ag),type="p",col="blue",pch=19,bg="orange", main="Age (dossier ELF)",ylim=c(0,100)) source("~/Bin/statgh.r") elf <- read.dar("elf.dar") attach(elf) sx <- SEXE ag <- AGE sy <- sx sy[sx==0] <- "H" sy[sx==1] <- "F" co <- sy co[sx==0] <- "blue" co[sx==1] <- "red" idr <- order(ag) plot(1:99,ag[idr],pch=sy[idr],col=co) plot(1:99,ag[idr],type="p",col="blue",pch=19, bg="orange", main="Age (dossier ELF)",ylim=c(0,100)) plot(1:99,ag[idr],type="p",col="blue",pch=sy, bg="orange", main="Age (dossier ELF)",ylim=c(0,100)) par(mfrow=c(2,1)) plot(...) numéro 1 plot(...) numéro 2 par(mfrow=c(1,1)) plot(...) par(mfrow=c(2,2)) plot(...) numéro 1 plot(...) numéro 2 plot(...) numéro 3 plot(...) numéro 4 boxplot(ag,col="yellow",main="Boite à moustaches d'AGE pour ELF") ma <- mean(ag) sa <- sd(ag) abline(h=median(ag),col="green",lwd=2) abline(h=ma,col="red",lwd=2) abline(h=ma-sa,col="blue",lwd=1.5) abline(h=ma+sa,col="blue",lwd=1.5) bam(ag,titre="Boite à moustaches d'AGEpour ELF") agh <- ag[sy=="H"] bam(agh,titre="Boite à moustaches d'AGE Homme pour ELF") agf <- ag[sy=="F"] bam(agf,titre="Boite à moustaches d'AGE Femme pour ELF") boxplot(ag~sx,main="Hommes/Femmes pour AGE dans ELF",col="yellow") rect(0,0,3,100,col=rgb(red=180/256,green=210/256,blue=190/256)) par(new=TRUE) boxplot(ag~sx,main="Hommes/Femmes pour AGE dans ELF", col="yellow",xaxt="n") axis(1,at=c(1,2),labels=c("Homme","Femme")) box() par(mfrow=c(2,1)) vins <- read.dbf("vins.dbf")$dbf head(vins) attach(vins) vins <- vins[,-1] boxplot(vins,col="yellow") boxplot(vins,col="yellow",ylim=c(0,30000)) par(mfrow=c(1,1)) hist(sx) hist(sx,freq=TRUE) hist(sx,freq=TRUE,col="blue") hist(sx,freq=TRUE,col="blue",xaxt="n",ylim=c(0,100)) axis(1,at=c(0.05,0.95),labels=c("Homme","Femme")) box() barplot(table(sx),ylim=c(0,100),col="blue") barplot(100*table(sx)/length(sx),ylim=c(0,100),col="blue") traceQL("SEXE, dossier ELF",sx,"Homme Femme") uk <- UK rf <- RFA rl <- lm(uk~rf) # régression linéaire cm <- anova(rl) # comparaison de moyennes names(rl) names(cm) ir <- order(rf) plot(rf[ir],uk[ir],pch=19,col="blue") abline(rl,col="red",lwd=2) plot(resid(rl)) abline(h=0,col="green",lwd=2) cor(rf,uk) rfo <- rf[ir] uko <- uk[ir] nld <- length(rfo)-2 rfo <- rfo[1:nld] uko <- uko[1:nld] rlo <- lm(uko~rfo) plot(rfo,uko,pch=19,col="blue") abline(rlo,col="red",lwd=2) res <- rl$residuals plot(res,pch=15,col="black") plot(res,pch=15,col="black") plot(uk,vaj,pch=15,col="black") abline(0,1,col="red",lwd=2) plot(rl) detach(elf) options(width=450) ronf <- read.dar("ronfle.dar") attach(ronf) head(ronf) ageronf <- AGE length(agelf) length(ageronf) ages <- c(agelf,ageronf) rep(1,10) dos <- c( rep(1,length(agelf)),rep(2,length(ageronf)) ) rect(0,0,3,100,col=rgb(red=180/256,green=210/256,blue=190/256)) par(new=TRUE) boxplot(age~dos,col="yellow",xaxt="n",main="Age ELF vs RONFLE") axis(1,at=c(1,2),labels=c("ELF","RONFLE")) lines(c(0.6,1.4),rep(mean(agelf),2),col="red",lwd=2) lines(c(0.6,1.4),rep(median(agelf),2),col="blue",lwd=2) lines(c(1.6,2.4),rep(mean(ageronf),2),col="red",lwd=2) lines(c(1.6,2.4),rep(median(ageronf),2),col="blue",lwd=2) cbind( length(agelf),mean(agelf),sd(agelf) ) cbind( length(ageronf),mean(ageronf),sd(ageronf) ) compMoyData("AGE Elf vs Ronfle",agelf,ageronf) rbind( cbind( length(agelf),mean(agelf),sd(ageronf) ), cbind( length(ageronf),mean(ageronf),sd(ageronf) ) ) -> mdr colnames(mdr) <- c("Taille","Moyenne","Ecart-type") row.names(mdr) <- c("Elf","Ronfle") mdr tapply(ag,sx,mean) elf <- read.dar("elf.dar") attach(elf) sx <- SEXE ag <- AGE sx <- as.factor(sx) levels(sx) <- c("homme","femme") boxplot(ag~sx) barplot(table(sx)) vins <- read.dbf("vins.dbf")$dbf vins <- vins[,-1] nomcol <- colnames(vins) unites <- rep("hl",length(nomcol)) allQT(vins,nomcol,unites) allqtdbf <- function(nomdbf,unites="???") { nomficdbf <- paste(nomdbf,".dbf",sep="") data <- read.dbf(nomficdbf)$dbf data <- data[,-1] nomcol <- colnames(data) if (length(unites)==1) { unites <- rep(unites,length(nomcol)) } # fin de si allQT(data,nomcol,unites) } # fin fonction allqtdbf unites <- rep("hl",8) allqtdbf("vins",unites) allqtdbf("vins","hl") allqtdbf("vins") allqtdbf("chu825qt") allqtdbf("chu825qt",chaineEnListe("G/l UI/l mmo/l mug/l % mg/dl an")) ronfle <- read.dbf("ronfleql.dbf")$dbf ronfcol <- 2:4 ronfqlm <- matrix(nrow=3,ncol=3) ronfqlm[1,1] <- "SEXE" ronfqlm[1,2] <- c(" Sexe de la personne ?") ronfqlm[1,3] <- lstMod(c("homme","femme")) ronfqlm[2,1] <- "RONFLE" ronfqlm[2,2] <- c(" Ronfleur ou pas ?") ronfqlm[2,3] <- lstMod(c("non","oui")) ronfqlm[3,1] <- "TABAC" ronfqlm[3,2] <- c(" Fumeur ou pas ?") ronfqlm[3,3] <- lstMod(c("oui","non")) print(ronfqlm) allQL(ronfle,ronfqlm,ronfcol) ronfqlm[1,3] <- lstMod("homme femme") allQL(ronfle,ronfqlm,ronfcol) rbinom(100,5,0.2) dbinom(1,5,0.2) dbinom(0:5,5,0.2) qbinom(0.5,5,0.2) si <- rbinom(100,5,0.5) table(si) barlplot( table(si) ) xi <- 0:5 ni <- table(si) pi <- ni/sum(ni) sum(pi*xi) sum(pi*xi**2) sqrt(sum(pi*xi**2)-sum(xi*pi)**2) t.test(x,y) compMoyData(" age H vs F ",ageh,agef) prop.test(c(ia,ib),c(na,nb)) compPourc(" Prop Hom/Fem ",10,80,30,70) rnorm(30) qnorm(0.975) qnorm(1-alpha/2) # fonctions sur chaines de caractères # avec le logiciel R modc <- unlist(strsplit(tabmoda,"!")) debutit <- substr(titrelong,1,10) sprintf("%-10s",intit) formatC(odl[nnt,1],format="d",width=3,dig=0) paste(formatC(valc,format="f",width=3,dig=0)," %") paste(paste(chMod,"!"),collapse="") paste("Voici les 10 premières lignes de données (il y en a ",nbl," en tout)",sep="") longueurChaine <- nchar( laChaine ) substr(x, start, stop) substring(text, first, last = 1000000) substr(x, start, stop) <- value substring(text, first, last = 1000000) <- value strtrim(x, width) deparse(args(lm)) strwidth(s, units = "user", cex = NULL) strheight(s, units = "user", cex = NULL) all.lett <- c(letters, LETTERS) chartr(old, new, x) tolower(x) toupper(x) casefold(x, upper = FALSE) # inclassables options(width=180) library(foreign) lbw <- read.dbf("lbw.dbf") # attention, homonyme de read.dbf ! sink("lbw.sor") lbw sink() library(gdata) dataCNV <- read.xls("s29.xls") update.packages() update.packages(ask=FALSE) mat[idl,,drop=FALSE] apply(mat,mARGIN=2,sum) # sommes en colonne # fin du pot-pourri