# (gH) -_- statgh.r ; TimeStamp (unix) : 17 Décembre 2010 vers 15:19 # version_gH <- 3.93 ################################################################ # # +----------------------------------------------- (gH) --+ # | | # | statgh.r : quelques fonctions utiles | # | en Rstat pour les statistiques | # | descriptives, les lois et les tests. | # | | # +-------------------------------------------------------+ # ################################################################ # # derniers ajouts # --------------- # . décembre 2010 : couleurs et nombres puis pct dans traceQL # rchModeleLogistique, lesColonnes # auroc(ql,qt) # . novembre 2010 : ajout de couleurs_metaf(), modalites_metaf() et de gr() # reprise de allQT avec des valeurs NA (et donc de mdc) # ajout de pchShow(), trouvé dans example(pch) pour voir les # symboles utilisables dans plot # ajout de ucfirst (initiale en majuscule) # . octobre 2010 : ajoute de matId pour matrice identité, # matPen pour matrice de pénalité et discordance, # pairsi (la fonction pairs améliorée) # . septembre 2010 : ajout de encoding="latin1" pour source # . aout 2010 : bendist(data,col=TRUE), decritClass, lit.dac, verifSol, # analexies(txt), lit.texte(fn),versiongh() # . juin 2010 : sdmax(n,a,b), sdr(x,a,b), ifiab(x,a,b) # ajout du paramètre maxhist dans plotQTdet # bbpQT(titre,qt) # . mai 2010 : sigcodes, as.sigcode(pval) # # ################################################################ # # prévoir : nbdec pour decritQT # pctmax pour decritQL # enQL(v01,"homme femme") pour convertir en "vraie QL" au sens de R # ################################################################ # # # on charge ce fichier en écrivant # # source("statgh.r",encoding="latin1") # # avec éventuellement un chemin d'accès pour le fichier, # par exemple : # # - sous Dos on écrit / au lieu du \ habituel # source("D:/RstatDemo/statgh.r") # # - sous Linux # source("/home/gh/RstatDemo/statgh.r") # # ############################################################## # # Fonctions définies : # allQL analyse toutes les qualitatives (tris à plat et tris croisés) # allQLnonum analyse toutes les qualitatives (tris à plat) pour modalités non numériques # allQLrecap récapitulatif de tous les tris à plat # allQLtriap analyse toutes les qualitatives (tris à plat) # allQLtricr analyse toutes les qualitatives (tris croisés) # # allQT analyse toutes les quantitatives et analyse des corrélations linéaires # anaLin étude de la liaison linéaire entre 2 quantitatives # approxPoiss approxime suivant une loi de Poisson # boxplotQT trace la boite à moustaches avec sa moyenne # chi2 chi-deux d'adéquation # chi2Adeq calcule le chi-deux d'adéquation entre valeurs observées et théoriques # chi2Indep calcule le chi-deux d'indépendance pour 2 variables qualitatives # chi2IndepTable calcule le chi-deux d'indépendance pour un tric croisé # compMoy comparaison de moyennes # compMoyData comparaison de moyennes via les données (avec anova et test de Welch) # compMoyNoData comparaison de moyennes via les paramètres n,m,V et anova et test de Welch # compPourc comparaison de pourcentages # decritQT décrit une quantitative # decritQTparFacteur décrit une quantitative suivant les modalités d'une qualitative # ect calcule l'écart-type mathématique exact des données # histEffectifs affiche un histogramme via les effectifs # histProba affiche un histogramme via les probabilités # histQL affiche l'histogramme du tri à plat de la variable qualitative # identif identification probabiliste # identifgc identification probabiliste de groupes connus # moy calcule la moyenne des données # plotQT trace la courbe avec sa moyenne et sa médiane # lit.dar lit un fichier .dar (avec nom de lignes et de colonnes) # lit.dat lit un fichier .dat (avec nom de lignes seulement) # lit.dbf lit un fichier .dbf (Dbase) # lit.xls lit un fichier .xls (Excel) mais nécessite gdata # skku calcule skewness et kurtosis # triAplat effectue un tri à plat (valeurs numériques) # triAplatAvecOrdre effectue un tri à plat (valeurs numériques) et ordonne les résultats # triAplatNonum effectue un tri à plat (valeurs non numériques) # triCroise effectue un tri croisé sans marges # triCroiseAvecMarges effectue un tri croisé avec marges (en % du total) # trisCroises effectue les 4 tris croisés (effectifs, divisions ligne,colonne,total) # vvar calcule la variance mathématique exacte des données # simule # chaineEnListe # Syntaxe de ces fonctions : # allQL(tdata,tmoda,numCol) # allQLnonum(tdata,numCol,nomVars) # allQLrecap(tdata,tmoda,numCol) # allQLtriap(tdata,tmoda,numCol) # allQLtricr(tdata,tmoda,numCol) # allQT(dataMatrix,nomV,nomU) # anaLin(titreVar1,varQT1,unite1,titreVar2,varQT2,unite2) # approxPoiss(lambda,nbpoints,nbclasses) # boxplotQT(titreQt,vecteurQt) # chi2Adeq(vth,vobs,contr=FALSE) # chi2Indep(var1,var2,modal1,modal2) # chi2IndepTable(tcr,mo1="",mo2="") # compMoyData(titre,var1,var2) # compMoyNoData(titre,nb1,moy1,var1,nb2,moy2,var2) # compPourc(titre,ia,na,ib,nb) # decritQL(titreQL,nomVar,modalités,graphique,fichier_image) # decritQT(titreQT,nomVar,unité,graphique,fichier_image) # ect(nomVar) # histEffectifs(effectifs,valeurs,titre) ; # histProbas(probas,valeurs,titre) ; # histQL(titre,variableQL,modalités) ; # identif(nomv,vec,mat,affiche=1) # identifgc(fdac,fmdp,fngr) # matcor(dataMatrix,nomV,nomU) # moy(nomVar) # plotQT(nomVar,titreQT,unité="cm",tige=TRUE) # lit.dbf(nomFichierDbf) # skku(variableQT) # triAplat(titreQL,nomVar,labelMod) # triAplatAvecOrdre(titreQL,nomVar,labelMod) # triAplatNonum(titreQL,nomVar) # triCroise(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2) # triCroiseAvecMarges(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2) # trisCroises(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2) # vvar(nomVar) # Exemples d'utilisation de ces fonctions : # pbioCOLQL <- 2:3 # lesm <- lstMod(c("non réponse","oui","non")) # pbioQL <- matrix(nrow=length(pbioCOLQL),ncol=3) # pbioQL[1,1] <- c("Connaissez-vous les produits biologiques ?") # pbioQL[1,2] <- c("CONNAISS") # pbioQL[1,3] <- lesm # pbioQL[2,1] <- c("Y at-il une différence entre produit biologique et produit diététique ?") # pbioQL[2,2] <- c("DIFF") # pbioQL[2,3] <-lesm # # pbio <- lit.dbf("pbio.dbf") # pbioDATA <- pbio$dbf # # allQL(pbioDATA,pbioQL,pbioCOLQL) # allQLtriap(pbioDATA,pbioQL,pbioCOLQL) # allQLrecap(pbioDATA,pbioQL,pbioCOLQL) # allQLtricr(pbioDATA,pbioQL,pbioCOLQL) # allQT(lesQT,c("AGE","POIDS","TAILLE","ALCOOL"),c("ans","kg","cm","verres")) # triAplat("Sexe de la personne",sexeElf, c("homme","femme") ) # triAplatAvecOrdre("Sexe de la personne",sexeElf, c("homme","femme") ) # triAplatNonum("Métiers ",varM) # # triCroise( # "Connaissez-vous les produits bio ?", # connaitrePbio, # c("non réponse","oui","non"), # "Différenciez-vous bio et diététique ?", # diffPbioPdiet, # c("non réponse","oui","non") # ) ; # fin de fonction triCroise # chi2Adeq( c(20,20,20,140), c(24,11,20,145) ) # chi2Adeq( 150*rep(1/6,6), c(22,21,22,27,22,36), contr=TRUE ) # chi2Indep(SEXE,ETUD,c("Homme","Femme"), c("NR","Primaire","Bepc","Bac","Supérieur")) # chi2Indep( table(SEXE,ETUD) ) # chi2IndepTable(smo,row.names(smo),colnames(smo)) # # decritQT("nombre d'élèves intégrés en 1997",integr,"personne(s)") ; # decritQT("nombre d'élèves intégrés en 1997",integr,"personne(s)",graphique=TRUE,fichier_image="elv.ps") ; # allQT(ecoles,c("INT","NPP","NEC","EFT","BUD","DAS","NPM","NOI","MSE")) ; # anaLin("Age Hommes",ageh,"ans","Age femmes",agef,"ans") # # histProba( dbinom((0:5),5,0.5) ,c(0:5),"binomiale B(5,0.5)") # histEffectifs( " Filtrages ") , c( 70 , 60 , 20 , 20 , 8 , 70 ) , 0:5 ) # # approxPoiss(10,300,15) # # compPourc(" ORDINATEURS par Région",20,102,28,98) # compMoyData( "HOMMES vs FEMMES ",c(1:8), c(2:10) ) # compMoyNoData( "AGE FR/US ",97, 5.851, 3.978, 73, 4.877, 5.642 ) # # skku( ageElf ) # skku( c(-1,2) ) # plotQT(lng,"Longueur (en acides aminés) des chaines","aa",TRUE) # boxplotQT(lng,"Longueur (en acides aminés) des chaines") # identif("CCUG1",vec1,mat1) # identifgc("xmp1.dac","xmp1.mdp","xmp1.ngr") # identifgc("rch3.dac","rch3.mdp","xmp1.ngr") # Options modifiées : options(width=180) ; ################################################################ ################################################################ ## ## ## ghf FONCTIONS pour VARIABLES QUALITATIVES ## ## ## ################################################################ ################################################################ ################################################################ fql <- function() { ################################################################ cat("Pour les QL (variables qualitatives) vous pouvez utiliser :\n") cat(" decritQL() traceQL() plotQL() triAplats() histEffectifs() histProba() \n") cat(" anaTcr() traceTcr() triCroise() triCroiseAvecMarges() triCroises() \n") } ; # fin de fonction fql ################################################################ triAPlats <- function() { ################################################################ cat("Pour les tris à plat vous pouvez utiliser :\n") cat(" triAplat() triAplatAvecOrdre() triAplatNonum() \n") } ; # fin de fonction triAPlats ################################################################ decritQL <- function(titreQL,nomVar,labelMod="",graphique=FALSE,fichier_image="") { ################################################################ if (missing(titreQL)) { cat(" syntaxe : decritQL(titreQL,nomVar,labelMod [ ,graphique=FALSE,fichier_image=\"\" ] )\n") cat(" exemples d'utilisation : \n") cat(" decritQL(\"SEXE dossier ELF\",sx,\"homme femme\")\n") cat(" decritQL(\"SEXE dossier Cetop\",sexeCet,\"homme femme NR\",graphique=TRUE,fichier_image=\"CET_age.png\")\n") } else { if (missing(titreQL)) { labelMod <- as.character(sort(unique(nomVar))) } if (length(labelMod)==1) { labelMod <- chaineEnListe( labelMod ) } ##if (labelMod=="") { labelMod <- as.character(sort(unique(nomVar))) } triAplat(paste(titreQL,"(ordre des modalités)"),nomVar,labelMod) triAplatAvecOrdre(paste(titreQL,"(par fréquence décroissante)"),nomVar,labelMod) if (graphique) { traceQL(titreQL,nomVar,labelMod,grfile=fichier_image) } } # fin si cat("\n") } ; # fin de fonction decritQL ################################################################ allQL <- function(tdata,tmoda,numCol) { ################################################################ if (missing(tdata)) { cat(" syntaxe : allQL(tdata,tmoda,numCol) \n") } else { nbco <- length(numCol) tdataQL <- matrix(nrow=dim(tdata)[1],ncol=nbco) row.names(tdataQL) <- row.names(tdata) colnames(tdataQL) <- tmoda[,1] for (idc in (1:nbco)) { tdataQL[,idc] <- tdata[,numCol[idc]] } ; # fin pour print10(tdataQL) allQLtriap(tdata,tmoda,numCol) # contient le résumé : allQLrecap(tdata,tmoda,numCol) if (length(numCol)>1) { allQLtricr(tdata,tmoda,numCol) } ; # fin si } # fin si cat("\n") } ; # fin de fonction allQL ################################################################ allQLtriap <- function(tdata,tmoda,numCol) { ################################################################ allQLrecap(tdata,tmoda,numCol) cat("\n ANALYSE DE TOUTES LES VARIABLES QUALITATIVES \n") nc <- length(numCol) for (i in 1:nc) { numc <- numCol[i] varc <- tdata[,numc] titc <- paste(tmoda[i,1]," -- ", tmoda[i,2]) modc <- unlist(strsplit(tmoda[i,3],"!")) triAplat(titc,varc,modc) } ; # fin pour i } ; # fin de fonction allQLtriap ################################################################ allQLnonum <- function(tdata,numCol,nomVars) { ################################################################ nc <- length(numCol) for (i in 1:nc) { numc <- numCol[i] nomV <- nomVars[i] varc <- tdata[,numc] triAplatNonum(nomV,varc) } ; # fin pour i } ; # fin de fonction allQLnonum ################################################################ allQLrecap <- function(tdata,tmodal,numCol) { ################################################################ cat("\n TABLEAU RECAPITULATIF DES VARIABLES QUALITATIVES\n") nc <- length(numCol) # on construit le tableau des numéros de variable et mode tmodeQ <- matrix(nrow=nc,ncol=2) for (i in 1:nc) { numc <- numCol[i] varc <- tdata[,numc] titc <- tmodal[i,1] tap <- as.integer(table(varc)) tmodeQ[i,1] <- i tmodeQ[i,2] <- max(tap) } ; # fin pour i idc <- order(tmodeQ[,2],decreasing=TRUE) # affichage des correspondances court/long cat("\n Intitulé Question\n") ; cat( " -------- --------\n") ; for (i in 1:nc) { intit <- substr(paste(tmodal[i,1]," "),1,10) cat( " ",sprintf("%-10s",intit)," ",sprintf("%-70s",tmodal[i,2]),"\n") } ; # fin pour i # on peut alors construire le tableau récapitulatif # dans le bon ordre (fourni par idc) cat("\n Affichage par mode décroissant puis par effectifs décroissants\n") ; ncr <- 7 trecap <- matrix(nrow=nc,ncol=ncr) for (i in 1:nc) { j <- idc[i] numc <- numCol[j] varc <- tdata[,numc] tap <- as.integer(table(varc)) lemode <- round(100*tmodeQ[j,2]/sum(tap)) idx <- order(tap,decreasing=TRUE) modc <- unlist(strsplit(tmodal[j,3],"!")) trecap[i,1] <- paste(" ",tmodal[j,1]) trecap[i,2] <- paste(formatC(lemode,format="f",width=3,dig=0)," %") trecap[i,3] <- modc[ idx[1] ] valc <- round(100*tap[idx[2]]/sum(tap)) trecap[i,4] <- paste(formatC(valc,format="f",width=3,dig=0)," %") trecap[i,5] <- modc[ idx[2] ] if (length(idx)>2) { valc <- round(100*tap[idx[3]]/sum(tap)) trecap[i,6] <- paste(formatC(valc,format="f",width=3,dig=0)," %") trecap[i,7] <- modc[ idx[3] ] } else { trecap[i,6] <- "" trecap[i,7] <- "" } ; # fin si } ; # fin pour i colnames(trecap) <- rep(" ",ncr) rownames(trecap) <- rep(" ",nc) print(trecap,quote=FALSE) } ; # fin de fonction allQLrecap ################################################################ allQLtricr <- function(tdata,tmoda,numCol) { ################################################################ # # tous les tris croisés d'un coup ! # en hors-d'oeuvre : proposition d'un ordre de lecture # des tris croisés via les p-values du chi2 d'indépendance nc <- length(numCol) nbtcr <- nc*(nc-1)/2 if (nc>2) { cat("\n ORDRE CONSEILLE POUR LIRE LES ",nbtcr," TRIS CROISES \n\n") } odl <- matrix(nrow=nbtcr,ncol=8) lespv <- vector(length=nbtcr) alpha <- 5 nt <- 0 for (i in 1:(nc-1)) { numc1 <- numCol[i] nomc1 <- substr(paste(tmoda[i,1]," "),1,10) for (j in (i+1):nc) { nt <- nt + 1 odl[nt,1] <- numc1 odl[nt,2] <- nomc1 numc2 <- numCol[j] odl[nt,3] <- numc2 odl[nt,4] <- substr(paste(tmoda[j,1]," "),1,10) tchi <- chisq.test( table(tdata[,numc1],tdata[,numc2]), correct=TRUE) ddl <- as.integer(tchi$parameter) vc <- as.real(tchi$statistic) + 0.00 odl[nt,5] <- vc odl[nt,6] <- qchisq( 1 -alpha/100, ddl ) pv <- as.real(tchi$p.value) + 0.00 odl[nt,7] <- pv lespv[nt] <- pv odl[nt,8] <- ddl } ; # fin pour j } ; # fin pour i idc <- order(lespv) cat(" Variable 1 Variable 2 Chi2 Chi2Table ") cat(" p-value Signif. Ddl") cat("\n") for (nt in (1:nbtcr)) { nnt <- idc[ nt ] cat(" ",formatC(odl[nnt,1],format="d",width=3,dig=0)) cat(" ",odl[nnt,2]) cat(" ",formatC(odl[nnt,3],format="d",width=3,dig=0)) cat(" ",odl[nnt,4]) vc <- as.real(odl[nnt,5]) + 0.00 cat(" ",formatC(vc,format="f",width=8,dig=2)) vm <- as.real(odl[nnt,6]) + 0.00 cat(" ",formatC(vm,format="f",width=8,dig=2)) pv <- as.real(odl[nnt,7]) + 0.00 cat(" ",formatC(pv,format="f",width=15,dig=7)) cat(" ") if (pv<0.01) { cat(" ** ") } else { if (pv<0.05) { cat(" * ") } else { cat(" ") } } ; # fin de si cat(" ",formatC(odl[nnt,8],format="d",width=3,dig=0)) cat("\n") } ; # fin pour nt # puis tous les tris-croisés par ordre d'entrée for (i in 1:(nc-1)) { numc1 <- numCol[i] varc1 <- tdata[,numc1] titc1 <- tmoda[i,1] modc1 <- unlist(strsplit(tmoda[i,3],"!")) for (j in (i+1):nc) { numc2 <- numCol[j] varc2 <- tdata[,numc2] titc2 <- tmoda[j,1] modc2 <- unlist(strsplit(tmoda[j,3],"!")) triCroise(titc1,varc1,modc1,titc2,varc2,modc2) } ; # fin pour j } ; # fin pour i cat("\n") } ; # fin de fonction allQLtricr ################################################################ anaTcr <- function(titreVar1,varQL1,modaQL1,titreVar2,varQL2,modaQL2) { ################################################################ # trace le tri croisé de deux variables sous forme de 4 histogrammes # (les 2 histos de tris à plat et ceux croisés dans les deux sens) # puis analyse le tri croisé à l'aide d'un chi-deux # # exemple : # anaTcr("SEXE",sexe,"RONFLE",ronf,c("Homme","Femme"),c("ne ronfle pas","ronfle") ) if (missing(titreVar1)) { cat(" syntaxe : anaTcr(titreVar1,varQL1,modaQL1,titreVar2,varQL2,modaQL2) \n") } else { par(mfrow=c(2, 2)) traceQL(titreVar1,varQL1,modaQL1,couleur="blue") traceQL(titreVar2,varQL2,modaQL2,couleur="yellow") chi2Indep(varQL2,varQL1,modaQL2,modaQL1) traceTcr(titreVar1,varQL1,modaQL1,titreVar2,varQL2,modaQL2) traceTcr(titreVar2,varQL2,modaQL2,titreVar1,varQL1,modaQL1) par(mfrow=c(1, 1)) } # fin si cat("\n") }# fin de fonction anaTcr ################################################################ traceQL <- function(titreVar,varQL,modaQL,couleur="blue",grfile="",txtpct=FALSE,ymax=100) { ################################################################ # # histogramme du tri à plat en pourcentages if (missing(titreVar)) { cat(" syntaxe : traceQL(titreVar,varQL,modaQL,couleur=\"blue\",grfile=\"\") \n") } else { if (length(modaQL)==1) { modaQL <- chaineEnListe( modaQL ) } print(modaQL) if (!grfile=="") { png(grfile,width=1024,height=768) } eff <- table(varQL) pct <- 100.0*(eff/length(varQL)) pctf <- paste(sprintf("%3d",round(pct)),"%") modaQL2 <- modaQL for (idm in (1:length(modaQL))) { modaQL2[idm] <- paste("\n\n",modaQL[idm],pctf[idm],sep="\n") } # fin pour if (txtpct==FALSE) { barplot(pct,col=couleur,main=titreVar,names.arg=modaQL,space=0.8,ylim=c(0,ymax)) } else { barplot(pct,col=couleur,main=titreVar,names.arg=modaQL2,space=0.8,ylim=c(0,ymax)) pose <- barplot(pct,plot=FALSE) # la valeur 1.54 est expérimentale text(pose*1.54,pct+2,format(eff)) } # fin si if (!grfile=="") { dev.off() cat("\n vous pouvez utiliser ",grfile,"\n") } ; # fin de si } # fin si cat("\n") } # fin de fonction traceQL ################################################################ traceTcr <- function(titreVar1,varQL1,modaQL1,titreVar2,varQL2,modaQL2) { ################################################################ # # affiche et trace le tri croisé entre deux QL # dans le sens fourni. exemple d'utilisation : # # traceTcr("SEXE",sexe,"RONFLE",ronf,c("Homme","Femme"),c("ne ronfle pas","ronfle") ) # # on n'affiche aucun résultat statistique # la version "dans les deux sens" se nomme anaTcr # if (missing(titreVar1)) { cat(" syntaxe : traceTcr(titreVar1,varQL1,modaQL1,titreVar2,varQL2,modaQL2) \n") } else { tc <- table(varQL1,varQL2) chm <- paste(titreVar2," / ",titreVar1) nbr <- length(table(varQL1)) rgby <- c(rep("red",nbr),rep("green",nbr),rep("blue",nbr),rep("yellow",nbr),rep("brown",nbr)) modaQL3 <- modaQL2 nbm <- length(modaQL3) for (i in 1:nbm) { modaQL3[ i ] <- paste(modaQL3[ i ],"\n",paste(modaQL1,collapse=" ")) } ; # fin pour i print( modaQL3 ) barplot(tc,beside=TRUE,col=rgby,main=chm,names.arg=modaQL3) } # fin si cat("\n") } # fin de fonction traceTcr ################################################################ histEffectifs <- function(tit,e,v) { ################################################################ lv <- length(v) ; vmin <- min(v) vmax <- max(v) vv <- vmin:vmax vy <- rep(vv,e) ; pdv <- v brk <- c(vmin - 0.75 + 0.5*(0:(2*vmax+2))) xmin <- vmin - 1 ; xmax <- vmax + 1 ; emin <- min(e) ; emax <- max(e) ; ymin <- 0 ; nbvaly <- 4 ech <- 10**trunc(log10(emax)) ymax <- round(emax/ech) if (ech*ymax0) { minay <- 0 } maxay <- round(100*max(v))/100 ldvy <- c(0) hy <- 0 while (hymaxay){ maxay <- hy } nbax <- lv + 1 pdvemin <- min(pdv)-1 pdvemax <- max(pdv)+1 brk <- c(min(pdv) - 0.8 + 0.5*(0:(2*length(pdv)))) pdve <- c(pdvemin:pdvemax) lpdve <- c(" ",pdv," ") hist(vy,br=brk,col="red", freq = TRUE,xlim=c(xmin,xmax),ylim=100*c(minay,maxay),axes=F, main=c("Distribution de la Loi ",tit),xlab="valeurs",ylab="probabilites") axis(1,pdve,lpdve,col.axis="blue") axis(2,100*ldvy,ldvy,col.axis="blue") matH <- matrix(nrow=lv,ncol=2) ; matH[,1] <- pdv matH[,2] <- v colnames(matH) <- c("valeurs","probabilite") rownames(matH) <- rep(" ",lv) print(matH,quote=FALSE,right=TRUE) ; } ; # fin de fonction histProba ################################################################ histQL <- function(titre, vecteurQL,labelMod="") { ################################################################ # # un simple synonyme pour plotQL # plotQL(titre, vecteurQL,labelMod) } ; # fin de fonction histQL ################################################################ lstMod <- function(chMod) { ################################################################ # # fait une liste de toutes les modalités à partir # d'une chaine de caractères (les modalités sont séparées par le symbole !) # return( paste(paste(chMod,"!"),collapse="")) } ; # fin de fonction lstMod ################################################################ plotQL <- function(titre, vecteurQL,labelMod="") { ################################################################ # # un "hist" amélioré pour les QL # histEffectifs(titre,table(vecteurQL),c(min(vecteurQL),max(vecteurQL))) triAplat(titre,vecteurQL,labelMod) } ; # fin de fonction plotQL ################################################################ triAplat <- function(titreQL,nomVar,labelMod) { ################################################################ if (missing(titreQL)) { cat("\n") cat(" syntaxe : triAplat( titreQL , nomVar , labelMod ) \n") ; cat(" exemple : triAplat( \"SEXE\",sx,c(\"homme\",\"femme\") ) \n") ; } else { cat("\n TRI A PLAT DE : ",titreQL,"\n\n") ; if (length(labelMod)==1) { labelMod <- chaineEnListe( labelMod ) } as.vector(table(nomVar)) -> vetcp ; # effectifs length(vetcp) -> nbmod ; # nombre de valeurs sum(vetcp) -> efftot ; # effectif total 100*vetcp/efftot -> vftc ; # fréquences % round(vftc) -> vftcp ; # fréquences arrondies en % round(cumsum(vftc)) -> vdcum ; # cumul matrix(nrow=3,ncol=nbmod+1) -> mrtcp ; mrtcp[1,(1:nbmod)] <- vetcp ; mrtcp[1, nbmod+1] <- efftot ; mrtcp[2,(1:nbmod)] <- vftcp ; mrtcp[2, nbmod+1] <- sum(vftcp) ; mrtcp[3,(1:nbmod)] <- vdcum ; mrtcp[3, nbmod+1] <- sum(vftcp) if (length(labelMod)==1) { labelMod <- unlist( strsplit(labelMod," ") ) } colnames(mrtcp) <- c(labelMod," Total") ; rownames(mrtcp) <- c(" Effectif "," Frequence (en %)"," Cumul fréquences") ; print(mrtcp,quote=FALSE,right=TRUE) ; } ; # fin de si cat("\n") ; } # fin de fonction triAplat ################################################################ triAplatAvecOrdre <- function(titreQL,nomVar,labelMod) { ################################################################ cat("\n QUESTION : ",titreQL,"\n\n") ; as.vector(table(nomVar)) -> vetcp ; # effectifs length(vetcp) -> nbmod ; # nombre de valeurs sum(vetcp) -> efftot ; # effectif total 100*vetcp/efftot -> vftc ; # fréquences % round(vftc) -> vftcp ; # fréquences arrondies en % round(cumsum(vftc)) -> vdcum ; # cumul matrix(nrow=3,ncol=nbmod+1) -> mrtcp ; if (length(labelMod==1)) { labelMod <- unlist( strsplit(labelMod," ") ) } # on trie ... idx <- order( vetcp, decreasing=TRUE) mrtcp[1,(1:nbmod)] <- vetcp[idx] ; mrtcp[1, nbmod+1] <- efftot ; mrtcp[2,(1:nbmod)] <- vftcp[idx] ; round(cumsum(vftc[idx])) -> vdcum ; # cumul mrtcp[2, nbmod+1] <- sum(vftcp) ; mrtcp[3,(1:nbmod)] <- vdcum ; mrtcp[3, nbmod+1] <- sum(vftcp) colnames(mrtcp) <- c(labelMod[idx]," Total") ; rownames(mrtcp) <- c(" Effectif "," Frequence (en %)"," Cumul fréquences") ; print(mrtcp,quote=FALSE,right=TRUE) ; cat("\n") ; } # fin de fonction triAplatAvecOrdre ################################################################ triAplatNonum <- function(titreQL,nomVar) { ################################################################ cat("\n QUESTION : ",titreQL,"\n\n") ; as.vector(table(nomVar)) -> vetcp ; # effectifs length(vetcp) -> nbmod ; # nombre de valeurs sum(vetcp) -> efftot ; # effectif total round(100*vetcp/efftot) -> vftcp ; # fréquences en % matrix(nrow=2,ncol=nbmod+1) -> mrtcp ; mrtcp[1,(1:nbmod)] <- vetcp ; mrtcp[1, nbmod+1] <- efftot ; mrtcp[2,(1:nbmod)] <- vftcp ; mrtcp[2, nbmod+1] <- sum(vftcp) ; colnames(mrtcp) <- c(levels(nomVar)," Total") ; rownames(mrtcp) <- c(" Effectif "," Fréquence (en %)") ; print(mrtcp,quote=FALSE,right=TRUE) ; cat("\n") ; } # fin de fonction triAplatNoNum ################################################################ triCroiseSansMarges <- function(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2) { ################################################################ if (missing(titreQL1)) { cat("\n") cat(" syntaxe : triCroiseSansMarges( titreQL1 , nomVar1 , labelMod1,titreQL2,nomVar2,labelMod2) \n") ; cat(" exemple : triCroiseSansMarges( \"SEXE\",sx,c(\"homme\",\"femme\"), \"ETUD\",etud,metud ) \n") ; return(-1) } # fin si cat("\n TRI CROISE DES QUESTIONS : \n ",titreQL1," (en ligne) \n ",titreQL2," (en colonne)\n") ; cat("Effectifs\n") ; table(nomVar1,nomVar2) -> mtcp ; rownames(mtcp) <- labelMod1 ; colnames(mtcp) <- labelMod2 ; print(mtcp,quote=FALSE) ; cat("\n% du total\n") ; mtcp <- round(100*mtcp/sum(mtcp)) ; print(mtcp,quote=FALSE) ; } # fin de fonction triCroiseSansMarges ################################################################ triCroise <- function(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2,graphique=FALSE,grfile="") { ################################################################ if (missing(titreQL1)) { cat("\n") cat(" syntaxe : triCroise( titreQL1 , nomVar1 , labelMod1,titreQL2,nomVar2,labelMod2,graphique=FALSE,grfile) \n") ; cat(" exemple : triCroise( \"SEXE\",sx,c(\"homme\",\"femme\"), \"ETUD\",etud,metud ) \n") ; cat(" si graphique est TRUE, on effectue le test complet du CHI2\n") ; } else { cat("\n TRI CROISE DES QUESTIONS : \n ",titreQL1," (en ligne) \n ",titreQL2," (en colonne)\n") ; if (!grfile=="") { png(grfile,width=1024,height=768) ; graphique = TRUE } mtcp <- table(nomVar1,nomVar2) ; md <- dim(mtcp) ; gmd <- dim(mtcp) + 1 ; tmtcp <- matrix(nrow=md[1],ncol=md[2]) ; gmtcp <- matrix(nrow=gmd[1],ncol=gmd[2]) ; if (length(labelMod1)==1) { labelMod1 <- unlist( strsplit(labelMod1," ") ) } if (length(labelMod2)==1) { labelMod2 <- unlist( strsplit(labelMod2," ") ) } cat("Effectifs\n") ; tmtcp[(1:md[1]),(1:md[2])] <- mtcp[(1:md[1]),(1:md[2])] rownames(tmtcp) <- labelMod1 colnames(tmtcp) <- labelMod2 print(tmtcp,quote=FALSE) ; gmtcp[1:md[1],1:md[2]] <- mtcp ; gmtcp[gmd[1] ,1:md[2]] <- table(nomVar2) ; gmtcp[1:md[1], gmd[2]] <- table(nomVar1) ; tg <- sum(mtcp) ; gmtcp[gmd[1] , gmd[2]] <- tg ; rownames(gmtcp) <- c(labelMod1,"TOTAL") ; colnames(gmtcp) <- c(labelMod2,"TOTAL") ; cat("\n Valeurs en % du total\n") ; mtcp <- round(100*mtcp/tg) ; gmtcp[1:md[1],1:md[2]] <- mtcp ; gmtcp[gmd[1] ,1:md[2]] <- round(100*table(nomVar2)/tg) ; gmtcp[1:md[1], gmd[2]] <- round(100*table(nomVar1)/tg) ; gmtcp[gmd[1] , gmd[2]] <- 100 ; print(gmtcp,quote=FALSE) ; if (graphique) { anaTcr(titreQL1,nomVar1,labelMod1,titreQL2,nomVar2,labelMod2) } ; # fin de si if (!grfile=="") { dev.off() cat("\n vous pouvez utiliser ",grfile,"\n") } ; # fin de si } # fin de si } # fin de fonction triCroise ################################################################ ################################################################ ## ## ## ghf FONCTIONS pour VARIABLES QUANTITATIVES ## ## ## ################################################################ ################################################################ ################################################################ fqt <- function() { ################################################################ cat("Pour les QT (variables quantitatives) vous pouvez utiliser :\n") cat(" decritQT() boxplotQT() plotQT() allCalcQT() decritQTparFacteur() \n") cat(" anaLin() traceLin() mdc() allQT() allQTdf() \n") } ; # fin de fonction fqt ################################################################ allQT <- function(dataMatrix,nomV="",nomU="",graphique=FALSE,grfile="") { ################################################################ # # on passe en revue toutes les variables QT (quantitatives) # avec calculs de m,s,s/ma etc. puis de la matrice des ccl # (coefficients de corrélation linéiares) avec la liste # de ces ccl par ordre décroissant ; on donne enfin # les ccl pour la forte corrélation avec production d'un # graphique postscript pour cette plus forte corrélation # # les données sont une matrice numérique (ligne 1 numérique, # colonne 1 numérique car les identifiants de ligne et de # colonne doivent avoir été transmis à rownames et colnames # # nomV est le nom des variables # nomU est le nom des unités associées aux variable # # exemple d'utilisation : # # allQT(lesQT,c("AGE","POIDS","TAILLE","ALCOOL"),c("ans","kg","cm","verres")) # vins <- lit.dbf("vins.dbf")$dbf[,-1] # allQT(vins,colnames(vins),rep("hl",length(colnames(vins)))) if (missing(dataMatrix)) { cat(" syntaxe : allQT(dataMatrix,nomV,nomU) \n") } else { if (length(nomV)==1) { if (nomV=="") { nomV <- colnames(dataMatrix) } } # finsi if (length(nomV)==1) { nomV <- chaineEnListe( nomV ) } if (length(nomU)==1) { nomU <- chaineEnListe( nomU ) } dm <- dim(dataMatrix) ; nl <- dm[1] nc <- dm[2] resMat <- matrix(nrow=nc,ncol=9) # identifiants pour la matrice des résultats rownames(resMat) <- rep(" ",nc) ; colnames(resMat) <- c(" Num "," Nom "," Taille"," Moyenne"," Unite"," Ecart-type"," Coef. de var."," Minimum"," Maximum") tabCdv <- vector(length=nc) for (i in 1:nc) { vorg <- dataMatrix[,i] v <- vorg[!is.na(vorg)] resMat[i,1] <- i resMat[i,2] <- sprintf("%-8s",nomV[i]) resMat[i,3] <- length(v) mv <- mean(v) vv <- var(v) resMat[i,4] <- formatC(mv,"f",width=10,dig=3) resMat[i,5] <- nomU[i] resMat[i,6] <- formatC( sqrt(vv),"f",width=10,dig=3) lecdv <- abs(100*(sqrt(vv)/mv)) tabCdv[i] <- lecdv fcdv <- sprintf("%6.2f",lecdv) resMat[i,7] <- sprintf("%-8s %%",fcdv) resMat[i,8] <- min(v) resMat[i,9] <- max(v) # formatage éventuel resMat[i,8] <- formatC(round(1000*as.numeric(resMat[i,8]))/1000,"f",width=9,dig=3) resMat[i,9] <- formatC(round(1000*as.numeric(resMat[i,9]))/1000,"f",width=9,dig=3) } # fin pour i # tri par cdv décroissant tabCdv <- tabCdv * (-1) odx <- order(tabCdv) ; ordMat <- matrix(nrow=nc,ncol=9) ; rownames(ordMat) <- rep(" ",nc) ; colnames(ordMat) <- colnames(resMat) for (i in 1:nc) { ordMat[i,] <- resMat[ odx[i], ] } # fin pour i # tri par moyenne décroissante moy <- as.numeric(resMat[,4]) moy <- moy * (-1) odx <- order(moy) ; moyMat <- matrix(nrow=nc,ncol=9) ; rownames(moyMat) <- rep(" ",nc) ; colnames(moyMat) <- colnames(resMat) for (i in 1:nc) { moyMat[i,] <- resMat[ odx[i], ] } # fin pour i # affichages ##cat("\nDonnées\n") ; ##print(dataMatrix,rowlab=rep(" ",nl),quote=FALSE,right=TRUE) ; print10(dataMatrix) cat("\nPar cdv décroissant\n") ; print(ordMat,quote=FALSE,right=TRUE) ; cat("\nPar ordre d'entrée\n") ; print(resMat,quote=FALSE,right=TRUE) ; cat("\nPar moyenne décroissante\n") ; print(moyMat,quote=FALSE,right=TRUE) ; # matrice des corrélations mdc(dataMatrix,nomV,meilcor=TRUE) # graphique par paires éventuel if (graphique) { if (!grfile=="") { png(grfile,width=1024,height=768) ; graphique = TRUE } pairsi(dataMatrix) if (!grfile=="") { cat("\n vous pouvez utiliser ",grfile,"\n") dev.off() } ; # fin de si } # fin si sur graphique } # fin si cat("\n") } ; # fin de fonction allQT ################################################################ allQTdf <- function(dataFrame,nomU) { ################################################################ if (missing(nomU)) { nomUnite <- rep(" ",dim(dataFrame)[2]) ; } else { nomUnite <- nomU if (length(nomU)==1) { nomUnite <- chaineEnListe( nomUnite ) } } ; # fin si allQT(dataFrame,names(dataFrame),nomUnite) } ; # fin de fonction allQTdf ################################################################ allQT2 <- function(dataMatrix,nomV,nomU) { ################################################################ # # comme allQT mais avec cl, cov, comed et cdl # dm <- dim(dataMatrix) ; nl <- dm[1] nc <- dm[2] resMat <- matrix(nrow=nc,ncol=10) # identifiants pour la matrice des résultats rownames(resMat) <- rep(" ",nc) ; cn <- c(" Nom "," Num "," Taille"," Moyenne"," Unite") cn <- c(cn," Ecart-type"," Coef. de var."," Coef. de lat."," Minimum"," Maximum") colnames(resMat) <- cn tabCdv <- vector(length=nc) for (i in 1:nc) { v <- dataMatrix[,i] resMat[i,1] <- i resMat[i,2] <- sprintf("%-8s",nomV[i]) resMat[i,3] <- length(v) mv <- mean(v) vv <- vvar(v) resMat[i,4] <- formatC(mv,"f",width=9,dig=3) resMat[i,5] <- nomU[i] resMat[i,6] <- formatC( sqrt(vv),"f",width=9,dig=3) lecdv <- 100*(sqrt(vv)/mv) tabCdv[i] <- lecdv fcdv <- sprintf("%6.2f",lecdv) resMat[i,7] <- sprintf("%-8s %%",fcdv) fclv <- sprintf("%6.2f",cl(v)) resMat[i,8] <- sprintf("%-8s %%",fclv) resMat[i,9] <- min(v) resMat[i,10] <- max(v) # formatage éventuel resMat[i,9] <- formatC(round(1000*as.numeric(resMat[i,9]))/1000,"f",width=9,dig=3) resMat[i,10] <- formatC(round(1000*as.numeric(resMat[i,10]))/1000,"f",width=9,dig=3) } # fin pour i # tri par cdv décroissant tabCdv <- tabCdv * (-1) odx <- order(tabCdv) ; ordMat <- matrix(nrow=nc,ncol=10) ; rownames(ordMat) <- rep(" ",nc) ; colnames(ordMat) <- colnames(resMat) for (i in 1:nc) { ordMat[i,] <- resMat[ odx[i], ] } # fin pour i # tri par moyenne décroissante moy <- as.numeric(resMat[,4]) moy <- moy * (-1) odx <- order(moy) ; moyMat <- matrix(nrow=nc,ncol=10) ; rownames(moyMat) <- rep(" ",nc) ; colnames(moyMat) <- colnames(resMat) for (i in 1:nc) { moyMat[i,] <- resMat[ odx[i], ] } # fin pour i # affichages cat("\nDonnées\n") ; print(dataMatrix,rowlab=rep(" ",nl),quote=FALSE,right=TRUE) ; cat("\nPar cdv décroissant\n") ; print(ordMat,quote=FALSE,right=TRUE) ; cat("\nPar ordre d'entrée\n") ; print(resMat,quote=FALSE,right=TRUE) ; cat("\nPar moyenne décroissante\n") ; print(moyMat,quote=FALSE,right=TRUE) ; # matrice des corrélations mdc(dataMatrix,nomV,meilcor=TRUE) } ; # fin de fonction allQT2 ################################################################ panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...) { ################################################################ usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- cor(x, y,use="pairwise.complete.obs") txt <- format(c(r, 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") pval <- as.sigcode(cor.test(x, y,use="pairwise.complete.obs")$p.value) if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, paste(txt,pval,sep="") , cex = cex.cor * r / 1.4 ) } ; # fin de fonction panel.cor ################################################################ panel.corpvalue <- function(x, y, digits=2, prefix="", cex.cor, ...) { ################################################################ usr <- par("usr"); on.exit(par(usr)) par(usr = c(0, 1, 0, 1)) r <- cor(x, y,use="pairwise.complete.obs") txt <- format(c(r, 0.123456789), digits=digits)[1] txt <- paste(prefix, txt, sep="") pval <- as.sigcode(cor.test(x, y,use="pairwise.complete.obs")$p.value) if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt) text(0.5, 0.5, pval, cex = cex.cor * r ) } ; # fin de fonction panel.corpvalue ################################################################ panel.pointsreg <- function(x, y, digits=2, prefix="", cex.cor, ...) { ################################################################ points(x,y,pch=19) abline(lm(y~x),col="red") } ; # fin de fonction panel.pointsreg ################################################################ pairsi <- function(lesVar,opt=2,...) { ################################################################ if (opt==1) { pairs(lesVar,lower.panel=panel.pointsreg,...) } ; # fin si if (opt==2) { pairs(lesVar,upper.panel=panel.cor,lower.panel=panel.pointsreg,...) } ; # fin si ################################################################ } ; # fin de fonction pairsi ################################################################ anaLin <- function(titreVar1,varQT1,unite1="?",titreVar2,varQT2,unite2="?",graphique=FALSE,grfile="") { ################################################################ # affiche et trace la corrélation linéaire entre deux QT # dans les deux sens, avec les statistiques détaillées. exemple d'utilisation : # # analin("TAILLE",taille,"POIDS",poi,"cm","kg") # if (missing(titreVar1)) { cat(" syntaxe : analin(titreVar1,varQT1,unite1=\"?\",titreVar2,varQT2,unite2=\"?\",graphique=FALSE,grfile=\"\") \n") } else { if (!grfile=="") { png(grfile,width=1024,height=768) ; graphique = TRUE } if (graphique) { par(mfrow=c(2, 2)) plot(varQT1,main=titreVar1,col="blue",xlab="",ylab="") plot(varQT2,main=titreVar2,col="blue",xlab="",ylab="") #plotQT(titreVar1,varQT1,unite1,tigef=FALSE) #plotQT(titreVar2,varQT2,unite2,tigef=FALSE) } ; # fin de si cat("\n ANALYSE DE LA LIAISON LINEAIRE ENTRE ",titreVar1," ET ", titreVar2,"\n\n") cdc <- cor(varQT1,varQT2) cat("\n coefficient de corrélation : ",cdc," donc R2 = ",cdc*cdc,"\n\n") fcdc <- sprintf("%5.2f",cdc) cat(" p-value associée : ",cor.test(varQT1,varQT2)$p.value,"\n\n") traceLin(titreVar1,varQT1,titreVar2,varQT2,TRUE) traceLin(titreVar2,varQT2,titreVar1,varQT1,TRUE) if (!grfile=="") { cat("\n vous pouvez utiliser ",grfile,"\n") dev.off() } ; # fin de si if (graphique) { par(mfrow=c(1, 1)) } } # fin si cat("\n") } ; # fin de fonction anaLin ################################################################ boxplotQT <- function(titreQT, vecteurQT,ylims=0,moyenne=TRUE,lng="FR",unitQT="?") { ################################################################ # # un "boxplot" amélioré pour les QT # avec la moyenne en rouge if (unitQT=="?") { leTitreQT <- titreQT } else { leTitreQT <- paste(titreQT," (",unitQT,")",sep="") } # finsi if (length(ylims)==1) { boxplot(vecteurQT,main=leTitreQT,col="yellow",pch=21,bg="red") } else { boxplot(vecteurQT,main=leTitreQT,col="yellow",pch=21,bg="red",ylim=ylims) } ; # fin de si if (moyenne) { abline(c(mean(vecteurQT),0),col="red",pch=21,lw=2) if (lng=="FR") { legend(x="topright",c("moyenne","médiane"),col=c("red","black"),lty=c(1,1),bty="n") } else { legend(x="topright",c("mean","median"),col=c("red","black"),lty=c(1,1),bty="n") } # fin si } ; # fin de si } ; # fin de fonction boxplotQT ################################################################ beanplotQT <- function(titreQT, vecteurQT,ylims=0,moyenne=TRUE) { ################################################################ # # un "beanplot" amélioré pour les QT # avec la moyenne en rouge library(beanplot) if (length(ylims)==1) { beanplot(vecteurQT,main=titreQT,col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6") } else { beanplot(vecteurQT,main=titreQT,col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6",ylim=ylims) } ; # fin de si if (moyenne) { abline(c(mean(vecteurQT),0),col="red",pch=21,lw=2) } } ; # fin de fonction boxplotQT ################################################################ decritQT <- function(titreQT,nomVar,unite=" ",graphique=FALSE,fichier_image="",msg="oui",lng="FR",vmaxhist="") { ################################################################ if (missing(titreQT)) { cat(" syntaxe : decritQT(titreQT,nomVar [ ,unite=\"?\",graphique=FALSE,fichier_image=\"\" ] ) \n") cat(" exemples d'utilisation : \n") cat(" decritQT(\"AGE dans ELF\",ag,\"ans\")\n") cat(" decritQT(\"POIDS dans Ronfle\",poi,\"kg\",graphique=TRUE,fichier_image=\"RON_poi.png\")\n") } else { cat("\n") if (lng=="FR") { cat("DESCRIPTION STATISTIQUE DE LA VARIABLE ",titreQT,"\n") ; } else { cat("STATISTICAL DESCRIPTION OF THE VARIABLE ",titreQT,"\n") ; } # fin si mdsc <- matrix(nrow=11,ncol=2) colnames(mdsc) <- rep(" ",2) ; if (lng=="FR") { rownames(mdsc) <- c(" Taille "," Moyenne"," Ecart-type"," Coef. de variation", " 1er Quartile", " Mediane "," 3eme Quartile"," iqr absolu", " iqr relatif", " Minimum"," Maximum") } else { rownames(mdsc) <- c(" Size "," Mean"," Standard deviation"," Coef. of variation", " 1st Quartile", " Median "," 3d Quartile"," absolute iqr", " relative iqr", " Minimum"," Maximum") } # fin si taille <- length(nomVar) ; moyenne <- sum(nomVar)/taille ; ecartype <- ect(nomVar) cdv <- round( cv(nomVar) ) mdsc[1,1] <- formatC(taille,format="d",width=4) ; if (lng=="FR") { mdsc[1,2] <- " individus"; } else { mdsc[1,2] <- " individuals"; } # fin si mdsc[2,1] <- formatC(moyenne,format="f",width=9,dig=4) ; mdsc[2,2] <- unite ; mdsc[3,1] <- formatC(ecartype,format="f",width=9,dig=4) ; mdsc[3,2] <- unite ; mdsc[4,1] <- sprintf("%5.0f ",cdv) mdsc[4,2] <- "% " ; mdsc[5,1] <- formatC(quantile(nomVar,0.25),9,dig=4) ; mdsc[5,2] <- unite; mdsc[6,1] <- formatC(quantile(nomVar,0.50),9,dig=4) ; # ou median(nomVar) ; mdsc[6,2] <- unite; mdsc[7,1] <- formatC(quantile(nomVar,0.75),9,dig=4) ; mdsc[7,2] <- unite; mdsc[8,1] <- formatC(quantile(nomVar,0.75)-quantile(nomVar,0.25),9,dig=4) ; mdsc[8,2] <- unite; mdsc[9,1] <- sprintf("%5.0f ",round(100*(quantile(nomVar,0.75)-quantile(nomVar,0.25))/(quantile(nomVar,0.50)))) mdsc[9,2] <- "% "; mdsc[10,1] <- sprintf("%9.4f",min(nomVar)) ; mdsc[10,2] <- unite; mdsc[11,1] <- sprintf("%9.4f",max(nomVar)) ; mdsc[11,2] <- unite; print(mdsc,quote=FALSE,right=TRUE) ; if (graphique) { if (vmaxhist=="") { plotQTdet(titreQT,nomVar,unite,tigef=TRUE,grfile=fichier_image,msg=msg,lng=lng) } else { plotQTdet(titreQT,nomVar,unite,tigef=TRUE,grfile=fichier_image,msg=msg,lng=lng,maxhist=vmaxhist) } } # fin si sur graphique cat("\n") } # fin si sur missing(titreQT) } ; # fin de fonction decritQT ################################################################ tdn <- function(titreQT,nomVar,unite="?",graphique=FALSE,fichier_image="") { ################################################################ # tests de normalité cat("TESTS DE NORMALITE POUR LA VARIABLE y =",titreQT,"\n\n") # on teste respectivement y, log(y+1), sqrt(y), 1/y, y*y # par les test de ks et de shapiro y <- nomVar nbm <- 6 mdv <- matrix(nrow=length(y),ncol=nbm) mdv[,1] <- y mdv[,2] <- sqrt(y) mdv[,3] <- 1/y mdv[,4] <- y*y mdv[,5] <- log(y+1) mdv[,6] <- log(log(y+1)) noms <- vector(length=nbm) noms[1] <- "y" noms[2] <- "sqrt(y)" noms[3] <- "1/y" noms[4] <- "y^2" noms[5] <- "ln(y+1)" noms[6] <- "ln(ln(y+1))" mtn <- matrix(nrow=2*nbm,ncol=3) colnames(mtn) <- c("Variable","Test","p-value") row.names(mtn) <- rep("",dim(mtn)[1]) for (idm in (1:nbm)) { v <- mdv[,idm] jdm <- 2*idm-1 mtn[jdm,1] <- noms[idm] mtn[jdm+1,1] <- mtn[jdm,1] mtn[jdm,2] <- "Kolmogorov-Smirnov" mtn[jdm,3] <- sprintf("%16.10f",ks.test(v,pnorm,mean=mean(v),sd=sd(v))$p.value) mtn[jdm+1,2] <- "Shapiro" mtn[jdm+1,3] <- sprintf("%16.10f",shapiro.test(v)$p.value) } # fin pour idm # print(mtn,quote=FALSE,right=TRUE) ; idx <- rev(order(mtn[,3])) mtn <- mtn[idx,] print(mtn,quote=FALSE,right=TRUE) ; if (graphique) { if (!fichier_image=="") { png(fichier_image,width=1024,height=768) ; graphique = TRUE } par(mfrow=c(3,2)) for (idm in (1:nbm)) { vecteurQT <- mdv[,idm] nom <- sub("y",titreQT,noms[idm]) titreHDC <- paste("Histogramme des classes pour ",nom,sep="") hist(vecteurQT,col="light blue",probability = TRUE,main=titreHDC ,xlab="",ylab="") vnorm <- function(x) { dnorm(x,mean=mean(vecteurQT),sd=sd(vecteurQT)) } curve(vnorm,add=TRUE,col="red",lwd=2) } # fin pour idm par(mfrow=c(1,1)) if (!fichier_image=="") { dev.off() } } # fin si sur graphique } ; # fin de fonction tdn ################################################################ rch_VE <- function(titreQT,nomVar,unite="?") { ################################################################ # recherche de valeurs extremes cat("RECHERCHE DE VALEURS EXTREMES POUR LA VARIABLE y =",titreQT,"\n\n") lng <- length(nomVar) moy <- mean(nomVar) ect <- sd(nomVar) cdv <- 100*ect/moy vt <- sort(nomVar) cat(" Moyenne :",moy,"unite ; écart-type : ",ect," ; cdv : ",sprintf("%5.1f",cdv),"%\n\n") ; cat(" Les 5 premières valeurs arrondies : ",round(vt[1:5]),"\n") ; cat(" Les 5 dernières valeurs arrondies : ",round(rev(rev(vt)[1:5])),"\n") ; cat("\n") cat(" Les 10 premières valeurs arrondies : ",round(vt[1:10]),"\n") ; cat(" Les 10 dernières valeurs arrondies : ",round(rev(rev(vt)[1:10])),"\n") ; cat("\n") for (elg in c(2,3,5,7)) { vinf <- moy-elg*ect vsup <- moy+elg*ect nbinf <- sum(nomVarvsup) if (nbinf+nbsup==0) { cat(" Il n'y a aucune valeur à l'extérieur de l'intervalle [ m -",elg,"s , m + ",elg,"s ].\n") ; } else { if (nbinf==0) { cat(" Il y n'a aucune valeur inférieure à m -",elg,".\n") ; } else { cat(" Il y a",nbinf," valeurs, soit ",sprintf("%4.1f",100*nbinf/lng),"% inférieures à m -",elg,"s à savoir : ") ; cat(nomVar[nomVarvsup]," \n") ; } # fin si nbsup } # fin si nbinf+nbsup==0 cat("\n") } # fin pour cat("\n") m1 <- mean(nomVar,trim=0.1) m2 <- mean(nomVar,trim=0.2) cat(" Moyenne tronquée à 10 % : ",m1,"\n") ; cat(" Moyenne tronquée à 20 % : ",m2,"\n\n") ; v1 <- vt[vt>min(vt)] v2 <- vt[vt0) { cat("\n ",nbnas," valeurs manquantes supprimées sur ",length(nomVarQT),"\n") } ; # fin si nomVarQT <- nomVarQT[!nas] nomVarQL <- nomVarQL[!nas] if (unite==" ") { unite_ <- " " } else { unite_ <- paste(",unit=",unite,sep="") } # fin si cat("\n") if (lng=="FR") { cat("VARIABLE QT ",titreQT,unite_,"\n") ; cat("VARIABLE QL ",titreQL," labels : ",labels,"\n") ; } else { cat("QT VARIABLE: ",titreQT,unite_,"\n",sep="") ; cat("QL VARIABLE: ",titreQL,", labels =",labels,"\n") ; } ; # fin de si cat("\n") nbl <- length(nomVarQT) tdata <- matrix(nrow=nbl,ncol=2) tdata[,1] <- nomVarQT tdata[,2] <- nomVarQL if (length(labels)==1) { labels <- chaineEnListe( labels ) } ## print(labels) vglob <- list() vglob[[1]] <- tdata[,1] valQL <- sort(unique(nomVarQL)) nbm <- length(valQL) for (idm in (1:nbm)) { vglob[[1+idm]] <- nomVarQT[nomVarQL==valQL[idm]] } # fin pour idm vres <- matrix(nrow=1+nbm,ncol=11) vres[1,] <- allCalcQT("Global",vglob[[1]],unite) for (idm in (1:nbm)) { vres[1+idm,] <- allCalcQT(labels[idm],vglob[[1+idm]],unite) } # fin pour idm row.names(vres) <- c("Global",labels) if (lng=="FR") { colnames(vres) <- chaineEnListe("N Moy Unite Ect Cdv Q1 Med Q3 EIQ Min Max") } else { colnames(vres) <- chaineEnListe("N Mean Unit Std Cv Q1 Med Q3 IQR Min Max") } ; # fin de si print(vres,quote=FALSE,right=TRUE) ; if (!fichier_image=="") { png(fichier_image,width=1024,height=768) ; graphique = TRUE } if (graphique) { library(beanplot) titre <- paste(titreQT," vs ",titreQL) if (nbm>2) { par(mfrow=c(1,2)) } else { par(mfrow=c(1,3)) } # fin si coul <- rep("black",nbl) coul[nomVarQL==min(nomVarQL)] <- "blue" coul[nomVarQL==max(nomVarQL)] <- "red" titrex <- paste(min(nomVarQL),"en bleu, ", max(nomVarQL)," en rouge ") titrex <- paste(labels[1],"en bleu, ", labels[2]," en rouge ") titrey <- "" ql <- as.factor(nomVarQL) levels(ql) <- labels if (nbm<=2) { plot(1:nbl,nomVarQT,col=coul,pch=21,bg=coul,main=titre,xlab=titrex,ylab=titrey) } boxplot(nomVarQT~ql,main=titre,col="yellow",pch=21,bg="red",notch=FALSE) beanplot(nomVarQT~ql,main=titre,col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6") par(mfrow=c(1,1)) } # fin si if (!fichier_image=="") { cat("\n") if (lng=="FR") { cat(" vous pouvez utiliser ",fichier_image,"\n") } ; # fin de si dev.off() } ; # fin de si #cat("Résultats de l'anova : \n") print (anova(lm(nomVarQT~nomVarQL))) } # fin si } ; # fin de fonction decritQTparFacteur ################################################################ decritQTparFacteurTexte <- function(titreQT,nomVarQT,unite="?",titreQL,nomVarQL,lstMod,graphique=FALSE,fichier_image="",haricot=0) { ################################################################ if (missing(titreQT)) { cat(" syntaxe : decritQTparFacteurTexte(titreQT,nomVar,unite=\"?\",titreQL,nomVarQL,labels, [ graphique=FALSE,fichier_image=\"\" ] ) \n") cat(" exemples d'utilisation : \n") cat(" decritQTparFacteurTexte(\"AGE dans CETOP\",ag,\"ans\",\"SEXE\",sexT,\"hom fam\")\n") cat(" decritQTparFacteurTexte(\"POIDS dans SANTE\",poi,\"kg\",,\"SEXE\",sexeT,\"homme femme\"graphique=TRUE,fichier_image=\"RON_poi.png\")\n") } else { cat("\n") cat("VARIABLE QT ",titreQT," unité : ",unite,"\n") ; cat("VARIABLE QL ",titreQL," labels : ",lstMod,"\n") ; cat("\n") nbl <- length(nomVarQT) if (length(lstMod)==1) { lstMod <- chaineEnListe( lstMod ) } vglob <- list() vglob[[1]] <- nomVarQT nbm <- length(lstMod) for (idm in (1:nbm)) { vglob[[1+idm]] <- nomVarQT[nomVarQL==lstMod[idm]] } # fin pour idm vres <- matrix(nrow=1+nbm,ncol=11) vres[1,] <- allCalcQT("Global",vglob[[1]],unite) for (idm in (1:nbm)) { vres[1+idm,] <- allCalcQT(lstMod[idm],vglob[[1+idm]],unite) } # fin pour idm row.names(vres) <- c("Global",lstMod) colnames(vres) <- chaineEnListe("N Moy Unite Ect Cdv Q1 Med Q3 EIQ Min Max") print(vres,quote=FALSE,right=TRUE) ; if (graphique) { if (!fichier_image=="") { png(fichier_image,width=1024,height=768) } titre <- paste(titreQT," vs ",titreQL) par(mfrow=c(1,1)) if (haricot==0) { boxplot(nomVarQT~nomVarQL,main=titre,col="yellow",pch=21,bg="red",notch=TRUE) } else { beanplot(nomVarQT~nomVarQL,main=titre,col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6") } ; # fin si if (!fichier_image=="") { dev.off() } } # fin si print(anova(lm(nomVarQT ~ as.factor(nomVarQL)))) } # fin si } ; # fin de fonction decritQTparFacteurTexte ################################################################ compare2QT <- function(titreGen,modaQt1,varQt1,modaQt2,varQt2,unite,graph=FALSE,fichier_image="") { ################################################################ if (missing(titreGen)) { cat(" syntaxe : compare2QT \n\n") cat(" exemple : compare2QT(\"age vs sexe\",\"homme\",ah,\"femme\",af,\"ans\",TRUE) \n\n") } else { n1 <- length(varQt1) n2 <- length(varQt2) vaql <- c( rep(modaQt1,n1), rep(modaQt2,n2) ) moda <- c(modaQt1,modaQt2) vaqt <- c(varQt1,varQt2) print( shapiro.test(varQt1) ) print( shapiro.test(varQt2) ) print( var.test(varQt1,varQt2) ) print( t.test(varQt1,varQt2) ) print( wilcox.test(varQt1,varQt2) ) decritQTparFacteurTexte(titreGen,vaqt,unite,"Groupes",vaql,moda,graph,fichier_image) } # fin si } ; # fin de fonction compare2QT ################################################################ compare2QTappariees <- function(titreGen,modaQt1,varQt1,modaQt2,varQt2,unite,graph=FALSE,fichier_image="") { ################################################################ if (missing(titreGen)) { cat(" syntaxe : compare2QT \n\n") } else { n1 <- length(varQt1) n2 <- length(varQt2) vaql <- c( rep(modaQt1,n1), rep(modaQt2,n2) ) moda <- c(modaQt1,modaQt2) vaqt <- c(varQt1,varQt2) print( shapiro.test(varQt1) ) print( shapiro.test(varQt2) ) print( var.test(varQt1,varQt2) ) print( t.test(varQt1,varQt2,paired=TRUE) ) print( wilcox.test(varQt1,varQt2,paired=TRUE) ) decritQTparFacteurTexte(titreGen,vaqt,unite,"Groupes",vaql,moda,graph,fichier_image) } # fin si } ; # fin de fonction compare2QTappariees ################################################################ # le "vrai" écart-type, issu de la "vraie" variance, mathématique et "biaisée" ect <- function(vecteurQt) { ################################################################ return( sqrt( vvar(vecteurQt) ) ) } ; # fin de fonction ect ################################################################ # le coefficient de variation basé sur ect cv <- function(vecteurQt) { ################################################################ mdv <- mean(vecteurQt,na.rm=TRUE) if (abs(mdv)<10**(-9)) { fcv <- (-1) } else { fcv <- 100*ect(vecteurQt) / mdv } # fin de si return( fcv ) } ; # fin de fonction cv ################################################################ mdc <- function(dataMatrix,nomV="",meilcor=TRUE,method="pearson") { ################################################################ # matrice des corrélations # # nomV est le nom des variables # exemple d'utilisation : # # mdc(lesQT,c("AGE","POIDS","TAILLE","ALCOOL")) # dm <- dim(dataMatrix) ; nl <- dm[1] nc <- dm[2] resMat <- matrix(nrow=nc,ncol=9) if (length(nomV)==1) { if (nomV=="") { nomV <- colnames(dataMatrix) } } # finsi corMat <- matrix(nrow=nc,ncol=nc) ; pvalue <- matrix(nrow=nc,ncol=nc) ; rownames(corMat) <- as.vector(nomV) colnames(corMat) <- nomV meiCor <- (-3) iEtj <- "-99" for (i in 1:nc) { vi <- dataMatrix[,i] for (j in 1:nc) { if (i1) { if (method=="pearson") { if (meilcor) { # Meilleures corrélations xmeiCor <- iEtj[1] ymeiCor <- iEtj[2] nomX <- nomV[xmeiCor] nomY <- nomV[ymeiCor] lesX <- dataMatrix[,xmeiCor] lesY <- dataMatrix[,ymeiCor] #lesC <- lm( lesY ~ lesX,model=FALSE ) # à cause des arrondis, on recalcule les valeurs précises idX <- !is.na(lesX) idY <- !is.na(lesY) lesX <- lesX[ idX & idY ] lesY <- lesY[ idX & idY ] mX <- mean(lesX) eX <- sqrt( mean(lesX*lesX) - mX*mX) mY <- mean(lesY) eY <- sqrt( mean(lesY*lesY) - mY*mY) ro <- (mean(lesX*lesY) - mX*mY)/(eX*eY) vA <- ro*eY/eX vB <- mY - vA*mX wA <- ro*eX/eY wB <- mX - wA*mY if (vB<0) { sgn1 <- "" } else { sgn1 <- "+" } if (wB<0) { sgn2 <- "" } else { sgn2 <- "+" } pval <- formatC( cor.test(lesX,lesY,use="pairwise.complete.obs")$p.value, width=5,dig=3) cat("\nMeilleure corrélation ",meiCor," pour ",nomX," et ",nomY," p-value ",pval,"\n\n") ; fXY <- paste(nomY," = ",formatC(vA,"f",width=9,dig=3),"*",nomX," ",sgn1,formatC(vB,"f",width=9,dig=3)) ; fYX <- paste(nomX," = ",formatC(wA,"f",width=9,dig=3),"*",nomY," ",sgn2,formatC(wB,"f",width=9,dig=3)) ; cat("Formules ",fXY,"\n") ; cat(" et ",fYX,"\n") ; # tracé correspondant à la meilleure corrélation oX <- order(lesX) newX <- lesX[ oX] newY <- lesY[ oX] #par(lwd=2) #postscript("r.ps") ; #plot(newX,newY,"b",col="red",pch=21,bg="orange",fg="black",xlab=nomX,ylab=nomY) #title(paste(" Régression ",fXY," (rho=",meiCor,")")) #abline(c(vB,vA),col="blue") #dev.off() # affichage des corrélations par ordre décroissant # on en fait d'abord un vecteur de valeurs if (nc>2) { cat("\n Coefficients de corrélation par ordre décroissant \n\n") ; tvec <- nc*(nc-1)/2 vecMdc <- matrix(nrow=tvec,ncol=4) ; # initialisation for (idc in 1:tvec) { vecMdc[idc,1] <- 0 vecMdc[idc,2] <- 0 vecMdc[idc,3] <- 0 } # fin pour idc # remplissage idc <- 0 for (i in 2:nc) { for (j in 1:(i-1)) { idc <- idc + 1 vcor <- corMat[i,j] # pour debug : cat(i,j,idc,vcor,"\n") ; vecMdc[idc,1] <- vcor vecMdc[idc,2] <- i vecMdc[idc,3] <- j vecMdc[idc,4] <- pvalue[i,j] } # fin pour j } # fin pour i # tri de vecMdc if (tvec>1) { for (idc in 1:(tvec-1)) { for (jdc in (idc+1):tvec) { vecMdc[idc,1] -> vi vecMdc[jdc,1] -> vj if (vi vi2 ; vecMdc[idc,3] -> vi3 ; vecMdc[idc,4] -> vi4 vecMdc[jdc,2] -> vj2 ; vecMdc[jdc,3] -> vj3 ; vecMdc[jdc,4] -> vj4 vecMdc[idc,1] <- vj ; vecMdc[idc,2] <- vj2 ; vecMdc[idc,3] <- vj3 ; vecMdc[idc,4] <- vj4 vecMdc[jdc,1] <- vi ; vecMdc[jdc,2] <- vi2 ; vecMdc[jdc,3] <- vi3 ; vecMdc[jdc,4] <- vi4 } # fin de si } # fin pour jdc } # fin pour idc } # fin de si # affichage nomC <- rownames(corMat) ; for (i in 1:tvec) { f1 <- formatC(vecMdc[i,1],"f",width=9,dig=3) ; f2 <- formatC(vecMdc[i,2],"f",width=9,dig=3) ; f3 <- formatC(vecMdc[i,3],"d",width=5) ; n2 <- nomC[ as.integer(vecMdc[i,2]) ] ; n3 <- nomC[ as.integer(vecMdc[i,3]) ] ; fn2 <- sprintf("%-9s",n2) fn3 <- sprintf("%-9s",n3) #f1bis <- sprintf("%6.4f",vecMdc[i,4]) f1bis <- sprintf("%6.4f",as.numeric(vecMdc[i,4])) if (!vecMdc[i,1]==" NaN") { cat(f1," p-value ",f1bis," pour ",fn2," et ",fn3,"\n") ; } # fin de si } # fin pour i } # fin de si } ; # fin de si meilcor=true } ; # fin de si method=pearson } ; # fin de si nc>1 } # fin de fonction mdc ################################################################ mcomed <- function(dataMatrix,nomV,meilcor=TRUE) { ################################################################ # matrice des comédianes # # nomV est le nom des variables # exemple d'utilisation : # # mcomed(lesQT,c("AGE","POIDS","TAILLE","ALCOOL")) # dm <- dim(dataMatrix) ; nl <- dm[1] nc <- dm[2] resMat <- matrix(nrow=nc,ncol=9) comMat <- matrix(nrow=nc,ncol=nc) ; rownames(comMat) <- as.vector(nomV) colnames(comMat) <- nomV meiCor <- (-3) iEtj <- "-99" for (i in 1:nc) { vi <-- dataMatrix[,i] for (j in 1:nc) { if (i chi2m cond <- cnd1 | cnd2 | cnd3 # affichage éventuel des contributions if (contr) { cat("\n Contributions au chi-deux \n\n") nc <- 7 cntr <- matrix(nrow=lthe,ncol=nc) colnames(cntr) <- c("Ind.","The","Obs","Dif","Cntr","Pct","Cumul") rownames(cntr) <- rep(" ",lthe) scnt <- 0 for (i in 1:lthe) { cntr[i,1] <- i cntr[i,2] <- vth[i] cntr[i,3] <- vobs[i] dif <- vth[i] - vobs[i] cntr[i,4] <- dif cntr[i,5] <- dif**2/vth[i] cntr[i,6] <- 100.0*(cntr[i,5]/vchi) scnt <- scnt + cntr[i,5] cntr[i,7] <- scnt } ; # fin pour i print(cntr) cat("\n") } ; # fin de si # affichage éventuel des histogrammes avec stockage en image PNG #if (!graph=="") { if (graph) { png1 <- paste(fichier,"1.png",sep="") png(png1,width=1024,height=768) scr(2) barplot(vth,col="red",main="Effectifs théoriques") barplot(vobs,col="blue",main="Effectifs observés") dev.off() png2 <- paste(fichier,"2.png",sep="") png(png2,width=1024,height=768) scr(1) barplot(rbind(vobs,vth),col=c("blue","red"),beside=TRUE, main="Effectifs observés (bleu) vs théoriques (rouge)") dev.off() cat(" Vous pouvez utiliser ",png1," et ",png2,"\n") } ; # fin de si } ; # fin de fonction chi2Adeq ################################################################ chi2Indep <- function(var1,var2,modal1,modal2) { ################################################################ # # un simple appel de la fonction suivante # chi2IndepTable( table(var1,var2),modal1,modal2 ) } ; # fin de fonction chi2Indep ################################################################ chi2IndepTable <- function(tcr,mo1="",mo2="") { ################################################################ # # calcul du chi-deux d'indépendance classique # MAIS : on détaille les contributions au chi-deux cat("\n CALCUL DU CHI-DEUX D'INDEPENDANCE\n") cat("\n =================================\n") nbl <- dim(tcr)[1] nbc <- dim(tcr)[2] mtric <- matrix(nrow=nbl+1,ncol=nbc+1) for (i in (1:nbl)) { for (j in (1:nbc)) { mtric[i,j] <- tcr[i,j] } ; # fin pour sur j } ; # fin pour sur i if (length(mo1)==1) { mo1 <- 1:nbl } if (length(mo2)==1) { mo2 <- 1:nbc } rownames(mtric) <- c(mo1,"Total") colnames(mtric) <- c(mo2,"Total") valt <- matrix(nrow=nbl+1,ncol=nbc+1) rownames(valt) <- c(mo1,"Total") colnames(valt) <- c(mo2,"Total") soml <- matrix(nrow=nbl,ncol=1) somc <- matrix(nrow=nbc,ncol=1) # calcul des marges en ligne for (i in (1:nbl)) { sl <- sum( tcr[i,] ) soml[ i ] <- sl valt[ i , nbc + 1 ] <- sl mtric[ i , nbc + 1 ] <- sl } ; # fin pour sur i # total colonne for (j in (1:nbc)) { sc <- sum( tcr[,j] ) somc[ j ] <- sc valt[ nbl + 1 , j ] <- sc mtric[ nbl + 1 , j ] <- sc } ; # fin pour sur j # total général tg <- sum(tcr) valt[nbl+1,nbc+1] <- tg mtric[nbl+1,nbc+1] <- tg cat("\n TABLEAU DES DONNEES \n\n") print(mtric,right=TRUE,print.gap=3) # remplissage des valeurs sous hypothèse d'indépendance for (i in (1:nbl)) { for (j in (1:nbc)) { valt[i,j] = soml[i]*somc[j]/tg } ; # fin pour sur j } ; # fin pour sur i cat("\n VALEURS ATTENDUES et MARGES \n\n") print(valt,digits=2,right=TRUE,print.gap=3) # calcul des contributions signées cat("\n CONTRIBUTIONS SIGNEES \n") nbt <- nbl*nbc tc <- vector(length=nbt) vl <- vector(length=nbt) vc <- vector(length=nbt) ts <- vector(length=nbt) vo <- vector(length=nbt) vt <- vector(length=nbt) ch <- 0 ic <- 0 vchi <- 0 mtcr <- matrix(nrow=nbl+1,ncol=nbc+1) rownames(mtcr) <- rep(" ",nbl+1) colnames(mtcr) <- rep(" ",nbc+1) mtcr[1,1] <- "" for (j in (1:nbc)) { mtcr[1,j+1] <- paste(" ",surncarg(8,mo2[j])) } ; # fin pour sur j for (i in (1:nbl)) { mtcr[i+1,1] <- surncarg(10,mo1[i]) for (j in (1:nbc)) { x <- (tcr[i,j]-valt[i,j]) ctr <- x*x/valt[i,j] vchi <- vchi + ctr ch <- ch + ctr if (x>0) { sgn <- "+" } else { sgn <- "-" } fvc <- formatC(ctr,format="f",width=6,dig=3) mtcr[i+1,j+1] <- paste(sgn,fvc) ic <- ic + 1 tc[ ic ] <- ctr vl[ ic ] <- i vc[ ic ] <- j ts[ ic ] <- x vo[ ic ] <- tcr[i,j] vt[ ic ] <- valt[i,j] } ; # fin pour sur j } ; # fin pour sur i print(mtcr,digits=3,right=TRUE,print.gap=3,quote=FALSE) ddl <- (nbl-1)*(nbc-1) alpha <- 5 pval <- chisq.test( tcr, correct=TRUE)$p.value chi2m <- qchisq( 1 -alpha/100, ddl ) cat("\n Valeur du chi-deux ",vchi,"\n") cat("\n Le chi-deux max (table) à 5 % est ",chi2m ) cat(" ; p-value ", pval) cat(" pour ", ddl," degrés de liberte") cat("\n") # une petite phrase pour conclure le test cat("\n décision : au seuil de ",alpha,"%") if (vchi<=chi2m) { cat(" on ne peut pas rejeter l'hypothèse") } else { cat(" on peut rejeter l'hypothèse") } ; # fin de si cat("\n qu'il y a indépendance entre ces deux variables qualitatives.\n") cat("\n") # tri des contributions signées idc <- order( tc , decreasing=TRUE) # affichage trié des contributions cat("\n PLUS FORTES CONTRIBUTIONS AVEC SIGNE DE DIFFERENCE \n\n") act <- matrix(nrow=nbt,ncol=9) rownames(act) <- rep(" ",nbt) colnames(act) <- c("Signe",surncarg(9," Valeur"),surncarg(7,"Pct"), surncarg(10,"Mligne"),surncarg(10,"Mcolonne"),"Ligne","Colonne","Obs","Th") for (i in (1:nbt)) { j <- idc[i] x <- ts[ j ] if (x>0) { act[i,1] <- " + " } else { act[i,1] <- " - " } act[i,2] <- formatC(tc[j],format="f",width=8,dig=3) ct <- tc[ j ] act[i,3] <- paste(formatC(100.0*(ct/vchi),format="f",width=5,dig=2), "%") act[i,4] <- surncarg(10,mo1[ vl[ j ] ]) act[i,5] <- surncarg(10,mo2[ vc[ j ] ]) act[i,6] <- vl[ j ] act[i,7] <- vc[ j ] act[i,8] <- formatC( vo[j] ,format="d",width=5,dig=0) act[i,9] <- formatC( vt[j] ,format="f",width=6,dig=1) } ; # fin pour sur i print(act,quote=FALSE,right=TRUE,print.gap=3) ; cat("\n") } ; # fin de fonction chi2IndepTable ################################################################ fcomp <- function() { ################################################################ cat("\n Vous disposez de deux fonctions pour comparer les moyennes : \n") cat("\n compMoyData( titre,nomA, varA, nomB,varB ) ") cat("\n si vous disposez des données") cat("\n") cat("\n compMoyNoData(titre,nA,moyA,ectA,nB,moyB,ectB) ") cat("\n si vous n'avez que les résumés statistiques") cat("\n") cat("\n Pour comparer des pourcentages :") cat("\n") cat("\n compPourc(titre,ia,na,ib,nb) \n\n") } ; # fin de fonction compMoy ################################################################ compMoyData <- function( titre,nomA, varA, nomB,varB ) { ################################################################ if (missing(titre)) { cat(" syntaxe : compMoyData( titre,nomA, varA, nomB,varB )\n") cat(" exemple : compMoyData( \"AGE Elf vs AGE Titanic\",\"AGE Elf\", age_elf, \"AGE Titanic\",age_titanic )\n") } else { # modification au 23 avril 2004 : utilisation de variance estimée # et non pas variance exacte (pour suivre sas) vab <- c(varA,varB) nab <- length(vab) mab <- sum(vab)/nab vab <- (sum(vab*vab)/nab-mab*mab)*nab/(nab-1) sab <- sqrt(vab) cab <- round(100*sab/mab) wab <- vab/nab na <- length(varA) ma <- sum(varA)/na va <- (sum(varA*varA)/na-ma*ma)*na/(na-1) sa <- sqrt(va) ca <- round(100*sa/ma) wa <- va/na nb <- length(varB) mb <- sum(varB)/nb vb <- (sum(varB*varB)/nb-mb*mb)*nb/(nb-1) sb <- sqrt(vb) cb <- round(100*sb/mb) wb <- vb/nb nu <- abs(ma-mb) de <- wa+wb r <- sqrt(de) delta <- nu/r cat("\n COMPARAISON DE MOYENNES (valeurs fournies) : ",titre,"\n\n") ; cat("Variable NbVal Moyenne Variance Ecart-type Cdv\n") ; cat(nomA," ",formatC(na,format="d",width=9), formatC(ma,format="f",width=12,digits=3), formatC(va,format="f",width=12,digits=3), formatC(sa,format="f",width=12,digits=3), formatC(ca,format="d",width=9)," %\n") ; cat(nomB," ",formatC(nb,format="d",width=9), formatC(mb,format="f",width=12,digits=3), formatC(vb,format="f",width=12,digits=3), formatC(sb,format="f",width=12,digits=3), formatC(cb,format="d",width=9)," %\n") ; cat("Global",formatC(nab,format="d",width=9), formatC(mab,format="f",width=12,digits=3), formatC(vab,format="f",width=12,digits=3), formatC(sab,format="f",width=12,digits=3), formatC(cab,format="d",width=9)," %\n") ; cat("\n") tt <- t.test(varA,varB) ; cat(" différence réduite : ",formatC(delta,format="f",digits=4)," ; p-value ",tt$p.value,"\n\n") ; cat(" au seuil de 5 % soit 1.96, on ") if (delta<1.96) { cat("ne peut pas rejeter") } else { cat("peut rejeter") } ; # fin de si cat(" l'hypothèse d'égalité des moyennes.\n\n") ; cat(" En d'autres termes, ") if (delta<1.96) { cat("il n'y a pas de") } else { cat("il y a une") } ; # fin de si cat(" différence significative entre les moyennes au seuil de 5 %.\n") cat("\n Voici ce qu'affiche la fonction t.test de R :\n" ) ; print(tt) ; cat("\n") ; cat("A l'aide d'anova, on obtient \n") valeurs <- c(varA,varB) groupes <- rep( c(1,2) , c(length(varA),length(varB)) ) print (anova(lm( valeurs ~ groupes ))) } # fin de si cat("\n") ; } ; # fin de function compMoyData ################################################################ compMoyNoData <- function(titre,na,ma,sa,nb,mb,sb) { ################################################################ if (missing(titre)) { cat(" syntaxe : compMoyNoData( titre,na,ma,sa,nb,mb,sb )\n") } else { ca <- round(100*sa/ma) va <- sa*sa wa <- va/na cb <- round(100*sb/mb) vb <- sb*sb wb <- vb/nb nu <- abs(ma-mb) de <- wa+wb r <- sqrt(de) delta <- nu/r cat("\n COMPARAISON DE MOYENNES (valeurs non fournies) ",titre,"\n\n") ; cat(" Variable nbVal Moyenne Variance Ecart-type Cdv\n") ; cat(" A ",formatC(na,format="d",width=9), formatC(ma,format="f",width=12,digits=3), formatC(va,format="f",width=12,digits=3), formatC(sa,format="f",width=12,digits=3), formatC(ca,format="d",width=9)," %\n") ; cat(" B ",formatC(nb,format="d",width=9), formatC(mb,format="f",width=12,digits=3), formatC(vb,format="f",width=12,digits=3), formatC(sb,format="f",width=12,digits=3), formatC(cb,format="d",width=9)," %\n") ; cat("\n") cat(" différence réduite : ",formatC(delta,format="f",digits=4),"\n\n") ; cat(" au seuil de 5 % soit 1.96, on peut ") if (delta<1.96) { cat("accepter") } else { cat("refuser") } ; # fin de si cat(" l'hypothèse d'égalité des moyennes.\n\n") ; cat(" En d'autres termes, ") if (delta<1.96) { cat("il n'y a pas de") } else { cat("il y a une") } ; # fin de si cat(" différence significative entre les moyennes au seuil de 5 %.\n") cat("\n") cat(" Vous pourriez utiliser à titre informatif\n\n") cat(" compMoyData(\"",titre,"\", simule(",ma,",",sa,",",na,"), simule(",mb,",",sb,",",nb,") ) \n") } # fin de si cat("\n") } ; # fin de function compMoyNoData ################################################################ compPourc <- function(titre,ia,na,ib,nb) { ################################################################ if (missing(titre)) { cat(" syntaxe : compPourc( titre,ia,na,ib,nb ) \n") } else { pa <- ia / na ; pb <- ib / nb ; p <- (ia+ib)/(na+nb) dp <- pa - pb df <- abs(dp) q <- p*(1-p)*(1/na+1/nb) r <- 1 eps <- 0 r <- sqrt(q) eps <- df/r rtest <- prop.test(c(ia,ib),c(na,nb)) pval <- rtest$p.value cat("\n COMPARAISON DE POURCENTAGES ",titre,"\n\n") ; cat(" population A, ",formatC(ia,width=5)," individus marqués sur ",formatC(na,format="d",width=7), " soit une proportion de ",formatC(100.0*pa,format="f",digits=5)," % \n") cat(" population B, ",formatC(ib,width=5)," individus marqués sur ",formatC(nb,format="d",width=7), " soit une proportion de ",formatC(100.0*pb,format="f",digits=5)," % \n") cat(" globalisation, ",formatC(ia+ib,width=5)," individus marqués sur ",formatC(na+nb,format="d",width=7), " soit une proportion de ",formatC(100.0*p,format="f",digits=5)," % \n") cat("\n") ; cat(" écart-réduit : ",formatC(eps,format="f",digits=4)," ; ") ; cat("\"p-value\" associée : ",pval,"\n\n") ; cat(" au seuil de 5 % soit 1.96, on peut ") if (eps<1.96) { cat("accepter") } else { cat("refuser") } ; # fin de si cat(" l'hypothèse d'égalité des pourcentages.\n\n") ; cat(" En d'autres termes, ") if (eps<1.96) { cat("il n'y a pas de") } else { cat("il y a une") } ; # fin de si cat(" différence significative entre les pourcentages au seuil de 5 %.\n\n") cat(" L'instruction R associée est \n\n prop.test( c(",ia,",",ib,") , c(",na,",",nb,") )\n\n") } # fin de si cat("\n") } ; # fin de fonction compPourc ################################################################ ################################################################ ## ## ## ghf FONCTIONS DE LECTURE DE FICHIERS ET DE VARIABLES ## ## ## ################################################################ ################################################################ ################################################################ lit <- function() { ################################################################ # # fonction générale de présentation des fonctions de lecture et # et d'écriture de fichiers de données # cat("\n") cat("Pour la lecture des fichiers, vous pouvez utiliser :\n\n") cat(" lit.dbf( fichier .DBF ) pour Dbase\n") cat(" lit.dar( fichier .DAR ) avec nom de lignes et de colonnes \n") cat(" lit.dat( fichier .DAT ) avec nom de lignes seulement \n") cat(" lit.xls( fichier .XLS ) pour Excel [via le package gdata] \n") cat(" ecrit.dar( dataframe ) avec nom de lignes et de colonnes \n") cat(" xls2dar( fichier .XLS) conversion .XLS ==> .DAR \n") cat("\n") } ; # fin de fonction lit ################################################################ lit.dbf <- function(dbf.name) { ################################################################ # # fonction copiée par (gH) avec l'autorisation de son auteur # pour lire un fichier Dbase # #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Shapefile Format - Read/Write shapefile format within R #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Ben Stabler benjamin.stabler@odot.state.or.us # Copyright (C) 2003 Oregon Department of Transportation #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Read DBF format #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ################################################################ infile<-file(dbf.name,"rb") # Header file.version <- readBin(infile,integer(), 1, size=1, endian="little") file.year <- readBin(infile,integer(), 1, size=1, endian="little") file.month <- readBin(infile,integer(), 1, size=1, endian="little") file.day <- readBin(infile,integer(), 1, size=1, endian="little") num.records <- readBin(infile,integer(), 1, size=4, endian="little") header.length <- readBin(infile,integer(), 1, size=2, endian="little") record.length <- readBin(infile,integer(), 1, size=2, endian="little") file.temp <- readBin(infile,integer(), 20, size=1, endian="little") header <- list(file.version,file.year, file.month, file.day, num.records, header.length, record.length) names(header) <- c("file.version","file.year","file.month","file.day","num.records","header.length","record.length") rm(file.version,file.year, file.month, file.day, num.records, header.length, record.length) # Calculate the number of fields num.fields <- (header$header.length-32-1)/32 field.name <- NULL field.type <- NULL field.length <- NULL field.decimal <- NULL # Field Descriptions (32 bytes each) for (i in 1:num.fields) { field.name.test <- readBin(infile,character(), 1, size=10, endian="little") field.name <- c(field.name,field.name.test) if (nchar(field.name.test)!=10) { file.temp <- readBin(infile,integer(), 10-(nchar(field.name.test)), 1, endian="little") } # fin de si field.type <- c(field.type,readChar(infile, 1)) file.temp <- readBin(infile,integer(), 4, 1, endian="little") field.length <- c(field.length,readBin(infile,integer(), 1, 1, endian="little")) field.decimal <- c(field.decimal, readBin(infile,integer(), 1, 1, endian="little")) file.temp <- readBin(infile,integer(), 14, 1, endian="little") } # fin de pour i # Create a table of the field info fields <- data.frame(NAME=field.name,TYPE=field.type,LENGTH=field.length,DECIMAL=field.decimal) # Set all fields with length<0 equal to correct number of characters fields$LENGTH[fields$LENGTH<0]<-(256+fields$LENGTH[fields$LENGTH<0]) # Read in end of attribute descriptions terminator - should be integer value 13 file.temp <- readBin(infile,integer(), 1, 1, endian="little") # Increase the length of field 1 by one to account for the space at the beginning of each record fields$LENGTH[1]<-fields$LENGTH[1]+1 # Add fields to the header list header <- c(header,fields=NULL) header$fields <- fields # Read in each record to a list element all.records <- list() for (i in 1:header$num.records) { all.records <- c(all.records, list(readChar(infile, header$record.length))) } # fin de pour i # Close the dbf file connection close(infile) # Function to split the strings and replace all " " with "" at the end of string format.record <- function(record) { record <- substring(record, c(1,cumsum(fields$LENGTH)[1:length(cumsum(fields$LENGTH))-1]+1),cumsum(fields$LENGTH)) record <- gsub(" +$","", record) record } # fin de fonction format.record # Split each record into columns and save as data.frame dbf <- data.frame(t(data.frame(lapply(all.records, format.record)))) rm(all.records) dimnames(dbf) <- list(1:header$num.records, header$fields$NAME) # Set the numeric fields to numeric for (i in 1:ncol(dbf)) { if(fields$TYPE[i]=="C") { dbf[[i]] <- as.character(dbf[[i]]) } if(fields$TYPE[i]=="N") { dbf[[i]] <- as.numeric(as.character(dbf[[i]])) } if(fields$TYPE[i]=="F") { d bf[[i]] <- as.numeric(as.character(dbf[[i]])) warning("Possible trouble converting numeric field in the DBF\n") } # fin de si type F } # fin de pour i # If the first field is of type character then remove the first # character of each record since the DBF stores a space for a # valid record and an * for a deleted record. # If the field is numeric then R removes the white space if(fields[1,2]=="C") { dbf[[1]] <- gsub("^[ *]", "", as.character(dbf[[1]])) } colnames(dbf) <- as.character(fields$NAME) colnames(dbf) <- gsub("_",".",colnames(dbf)) # Return the dbf as a list with a data.frame and a header list list(dbf=dbf, header=header) } # fin de fonction lit.dbf ################################################################ lit.dar <- function( nomfic ) { ################################################################ # # cette fonction traite la première ligne du fichier comme le nom # des colonnes et ensuite la première colonne de chaque ligne # comme le nom des lignes # read.table(nomfic,head=TRUE,row.names=1) } ; # fin de fonction lit.dar ################################################################ ecrit.dar <- function( df , nomfic ) { ################################################################ # # cette fonction écrit un "dataframe" avec comme première ligne du fichier comme le nom # des colonnes et ensuite la première colonne de chaque ligne # est le nom des lignes # #write.table(df,nomfic,quote=FALSE) options(width=380) ; sink(nomfic) print(df) sink() cat("en principe, vous pouvez utiliser le fichier ",nomfic,"\n") } ; # fin de fonction ecrit.dar ################################################################ lit.dat <- function( nomfic ) { ################################################################ # # cette fonction lit un fichier de données # la première valeur de chaque ligne devient le nom de la ligne # la première ligne contient des données, ce n'est pas le nom des colonnes read.table(nomfic,head=FALSE,row.names=1) } ; # fin de fonction lit.dat ################################################################ lit.xls <- function( nomfic ) { ################################################################ # # lecture d'un fichier Excel library(gdata) dfxls <- read.xls( nomfic ) return( dfxls ) } ; # fin de fonction lit.xls ################################################################ nombase <- function( nomfic ) { ################################################################ # # on récupère le nom du fichier (partie avant le point) # tabstr <- strsplit(nomfic,"\\.") return( unlist(tabstr)[1]) } ; # fin de fonction nombase ################################################################ xls2dar <- function( nomficxls ) { ################################################################ # # convertit un fichier Excel en un fichier .DAR if (missing(nomficxls)) { cat("xls2dar : convertit un fichier .xls en un fichier texte .dar avec noms des colonnes \n") ; cat("syntaxe : xls2dar(fichierExcel) \n") cat("exemple : xls2dar(\"elf.xls\") \n\n") return() ; } ; # fin de si library(gdata) df <- read.xls( nomficxls ) row.names(df) <- df[,1] df <- df[,(-1)] nomdar <- paste(nombase(nomficxls),".dar",sep="") ecrit.dar(df,nomdar) } ; # fin de fonction xls2dar ################################################################ ################################################################ ################################################################ ## ## ## ghf AUTRES FONCTIONS ## ## ## ################################################################ ################################################################ ################################################################ print10 <- function( matd ) { ################################################################ # # affiche les 10 premières de la matrice des données # nbl <- dim(matd)[1] dix <- 10 if (nbl>dix) { lemsg <- paste("Voici les 10 premières lignes de données (il y en a ",nbl," en tout)",sep="") } else { lemsg <- paste("Voici l'ensemble des ",nbl," lignes de données",sep="") dix <- nbl } ; # fin de si catln(lemsg) print(matd[(1:dix),]) } ; # fin de fonction read.dar ################################################################ litcol <- function(fichier,numCol,entete=TRUE) { ################################################################ mydata <- read.table(fichier,head=entete) vdata <- mydata[,numCol] cat(" vous disposez de ",length(vdata)," valeurs \n") return( vdata ) } ; # fin de fichier litcol ################################################################ ## FONCTIONS SUR CHAINES DE CARACTERES ################################################################ catln <- function( x) { ################################################################ # affiche x suivi d'un saut de ligne cat(x,"\n") } ; # fin de fonction catln ################################################################ surncarg <- function( long,chen ) { ################################################################ # # renvoie la chaine "chen" sur "long" caractères # sans tronquer ni déborder # newchen <- chen while (nchar(newchen)prec) && (netap<20)) { netap <- netap + 1 xdata <- xdata - lmdif lmdif <- mean(xdata) - moy_s } ; # fin de tant que ledif <- ect(xdata) - ect_s # correction de l'écart-type vpc <- 1 ntour <- 1 ledif <- ect(xdata) - ect_s xdata <- sort(xdata) while ((abs(ledif)>prec) && (ntour<2000)) { netap <- 1 idvam <- 1 xdata <- sort(xdata) ledif <- ect(xdata) - ect_s while ((abs(ledif)>prec) && (netap<2000)) { netap <- netap + 1 if ((xdata[idvam]>moy_s) && (ledif>0)) { xdata[idvam] <- xdata[idvam] - ledif/vpc sens <- 1 } else { xdata[idvam] <- xdata[idvam] + ledif/vpc sens <- 0 } # fin de si jdvam <- nb_val + 1 - idvam if (sens==0) { xdata[jdvam] <- xdata[jdvam] - ledif/vpc } else { xdata[jdvam] <- xdata[jdvam] + ledif/vpc } # fin de si idvam <- idvam + 1 if (idvam>nb_val) { idvam <- 1 } ledif <- ect(xdata) - ect_s } ; # fin de tant que ntour <- ntour + 1 } ; # fin de tant que return(xdata) } # fin si cat("\n") } # fin de fonction simule ################################################################ allqtdbf <- function(nomdbf,unites="???") { # exemples : # allqtdbf("vins",unites) # allqtdbf("vins","hl") # allqtdbf("vins") # allqtdbf("chu825qt") # allqtdbf("chu825qt",chaineEnListe("G/l UI/l mmo/l mug/l % mg/dl an")) 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 de fonction allqtdbf ################################################################ col2fact <- function(tdd) { ################################################################ # transforme des colonnes en facteur # par exemple 10 15 30 # 20 21 19 # devient 10 1 # 20 1 # 15 2 # 21 2 # 30 3 # 19 3 dd <- dim(tdd) nblig <- dd[1] nbcol <- dd[2] nblc <- nblig*nbcol mr <- matrix(nrow=nblc,ncol=2) idl <- 0 for (idc in 1:nbcol) { deb <- (idc-1)*nblig + 1 fin <- deb + nblig - 1 plg <- deb:fin mr[plg,1] <- tdd[,idc] mr[plg,2] <- idc } ; # fin pour idc return( mr) } # fin de fonction col2fact ########################################### bestValAssoc <- function(x,y) { if (y>0.5) { return(x) } else { return(1-x) } } # fin de fonction bestValAssoc ########################################### ########################################### identif <- function(nomv,vec,mat,affiche=1) { ########################################### if (affiche==1) { cat("\nIDENTIFICATION PROBABILISTE\n\n") } nbl <- dim(mat)[1] nbc <- dim(mat)[2] if (affiche==1) { cat("Vecteur\n") print( vec ) cat("Matrice\n") print( mat ) } # fin si vpa <- vector(length=nbl) fpa <- vector(length=nbl) rpa <- vector(length=nbl) cvm <- vector(length=nbl) pvm <- vector(length=nbl) for (idl in (1:nbl)) { p1 <- 1 p2 <- 1 for (jdc in (1:nbc)) { x <- mat[idl,jdc] y <- vec[jdc] p1 <- p1*bestValAssoc(x,y) p2 <- p2*bestValAssoc(x,x) } ; # fin pour jdc vpa[idl] <- p1 cvm[idl] <- p2 } ; # fin pour idl vpa <- round( vpa , 4 ) fpa <- round( 100*vpa/sum(vpa) , 2 ) rpa <- round( 100*vpa/max(vpa) , 2 ) pvm <- round( 100*vpa/cvm , 2 ) matres <- cbind(vpa,fpa,rpa,cvm,pvm) row.names(matres) <- row.names(mat) if (affiche==1) { cat("\nRésultats par groupe\n") print(matres) cat("\n") } # fin si # gestion des résultats pour tests sur groupes connus tauxMs <- 80 tauxVm <- 80 de <- 0 ms <- 0 is <- 0 vm <- 0 for (idl in (1:nbl)) { if (fpa[idl]>ms) { is <- idl ms <- fpa[idl] vm <- pvm[idl] } # fin si if (fpa[idl]>tauxMs) { de <- 1 if (affiche==1) { cat("Le vecteur semble appartenir au groupe ") cat(idl) cat(" soit ") cat(row.names(mat)[idl]) cat("\n") } # fin si if (pvm[idl]>tauxVm) { if (affiche==1) { cat("En fait, le vecteur appartient vraiment au groupe ") cat(idl) cat(" soit ") cat(row.names(mat)[idl]) cat(".\n") } # fin si } else { if (affiche==1) { cat("En fait, le vecteur n'appartient pas vraiment au groupe ") cat(idl) cat(" soit ") cat(row.names(mat)[idl]) cat(".\n") } # fin si } # fin si } # fin si } ; # fin pour idl ldr <- c(nomv,de,ms,vm,is) if (affiche==1) { cat("\nLigne résumée d'affectation\n") cat(ldr) cat("\n") } # fin si return(ldr) } # fin de fonction identif ########################################### ########################################### identifgc <- function(fdac,fmdp,fngr) { ########################################### cat("\nVERIFICATION D'IDENTIFICATION PROBABILISTE SUR GROUPES CONNUS\n") dac <- read.table(fdac,row.names=1) mdp <- read.table(fmdp,row.names=1) ngr <- read.table(fngr) ngr <- ngr[,-1] grp <- dac[,1] mdac <- dac[,-1] nbl <- dim(mdac)[1] nbc <- dim(mdac)[2] colnames(mdac) <- 1:nbc colnames(mdp) <- 1:nbc cat("Colonnes\n") print(ngr) cat("Données binaires avec indication de groupe\n") print(dac) cat("Matrice des fréquences de positivité\n") print(mdp) matres <- matrix(nrow=nbl,ncol=9) row.names(matres) <- row.names(mdac) nb_de <- 0 nb_bc <- 0 nb_mc <- 0 for (idl in (1:nbl)) { nmvec <- row.names(mdac)[idl] levec <- mdac[idl,] lr <- identif(nmvec,levec,mdp,0) print(c(idl,lr)) matres[idl,1] <- lr[2] matres[idl,2] <- grp[idl] matres[idl,3] <- lr[5] matres[idl,4] <- lr[3] matres[idl,5] <- lr[4] nb_de <- nb_de + as.numeric(lr[2]) matres[idl,6] <- nb_de if (matres[idl,3]==grp[idl]) { vbc <- 1 nb_bc <- nb_bc + 1 } else { vbc <- 0 nb_mc <- nb_mc + 1 } # fin de si matres[idl,7] <- vbc matres[idl,8] <- nb_bc matres[idl,9] <- nb_mc } ; # fin pour idl print(matres) } # fin de fonction identifgc ########################################### cdr <- function(vecteurQT) { return( cdrn(vecteurQT,min(vecteurQT),max(vecteurQT)) ) } ; # fin de fonction cdr cdrn <- function(vecteurQT,a,b) { nbv <- length(vecteurQT) rng <- b-a ect <- sd(vecteurQT) smax <- rng/2 rbiais <- nbv/(nbv-1) fcdrn <- ect/(sqrt(rbiais)*smax) ; return( fcdrn ) } ; # fin de fonction cdrn ########################################### ########################################### ic <- function() { ########################################### # # rappel des diverses fonctions pour les intervalles de confiance # cat("\n Vous disposez de cinq fonctions pour les intervalles de confiance ") cat("\n et de deux fonctions pour les calculs de taille d'échantillons : \n") cat("\n icp \n") cat(" pour l'estimation par intervalle d'une proportion\n") cat(" syntaxe : icp(nbVal,pChap,nivConf=0.05,echo=FALSE) \n") cat("\n icm \n") cat(" pour l'estimation par intervalle d'une moyenne sachant m et s\n") cat(" syntaxe : icm(nbVal,moyQt,ectQt,nivConf=0.05,affichage=FALSE) \n") cat("\n icmQT \n") cat(" pour l'estimation par intervalle de la moyenne d'une QT\n") cat(" syntaxe : icmQT(varQT,nivConf=0.05,affichage=FALSE) \n") cat(" (on peut aussi utiliser t.test(varQT) \n" ) ; cat("\n ice \n") cat(" pour l'estimation par intervalle de l'écart-type \n") cat(" syntaxe : ics(nbVal,ectQt,nivConf=0.05,echo=FALSE) \n") cat("\n iceQT \n") cat(" pour l'estimation par intervalle d'une écart-type d'une QT\n") cat(" syntaxe : icsQT(varQT,nivConf=0.05,echo=FALSE) \n") cat("\n tailleEchProp \n") cat(" pour la taille d'un échantillon en vue d'estimer une proportion\n") cat(" syntaxe : tailleEchProp(pChapeau,margeErreur,nivConf=0.05) \n") cat("\n") cat("\n tailleEchMoy \n") cat(" pour la taille d'un échantillon en vue d'estimer une moyenne\n") cat(" syntaxe : tailleEchMoy(ecartType,margeErreur,nivConf=0.05) \n") cat("\n") } ; # fin de fonction ic ########################################### icp <- function(nbVal,pChap,nivConf=0.05,affichage=FALSE) { ########################################### if (missing(nbVal)) { cat("icp : intervalle de confiance d'une proportion \n") cat("syntaxe : icp(nbVal,pChap,nivConf=0.05,echo=FALSE) \n") cat("exemples : icp(580,0.262) \n") cat(" icp(580,0.262,0.05) \n") cat(" icp(580,0.262,0.05,TRUE) \n") return() } ; # fin de si alp <- nivConf/2; z <- qnorm(1-alp); f_z <- sprintf("%6.3f",z) ; margeE <- z*sqrt(pChap*(1-pChap)/nbVal) liminf <- pChap - margeE limsup <- pChap + margeE if (affichage) { cat(" Pour ",nbVal," et une proportion estimée p-chapeau ",pChap," au niveau ",nivConf,"\n") cat(" la valeur critique issue de la loi normale est ",z," soit ",f_z,"\n") cat(" l'intervalle de confiance est donc sans doute [",liminf," ; ",limsup," ]\n ") cat(" soit, en arrondi, [ ",sprintf("%0.3f",liminf)," ; ",sprintf("%0.3f",limsup)," ]\n ") } # fin de si return(c(liminf,limsup)) } ; # fin de fonction icm ########################################### icm <- function(nbVal,moyQt,ectQt,nivConf=0.05,affichage=FALSE) { ########################################### if (missing(nbVal)) { cat("icm : intervalle de confiance d'une moyenne \n") ; cat("syntaxe : icm(nbVal,moyQt,ectQt,nivConf=0.05,echo=FALSE) \n") cat("exemples : icm(50,25.03,1.341641) \n") cat(" icm(50,25.03,1.341641,0.01) \n") cat(" icm(50,25.03,1.341641,0.01,TRUE) \n") return() ; } ; # fin de si ddl <- nbVal-1; alp <- nivConf/2; valt <- qt(alp,ddl,lower.tail=FALSE) f_t <- sprintf("%6.3f",valt) ; margeE <- valt*ectQt/sqrt(nbVal) liminf <- moyQt - margeE limsup <- moyQt + margeE if (affichage) { cat(" Pour ",nbVal," valeurs de moyenne ",moyQt," et d'écart-type ",ectQt," au niveau ",nivConf,"\n") cat(" la valeur critique issue de la loi du t de Student pour ",ddl," ddl est ",valt," soit ",f_t,"\n") cat(" l'intervalle de confiance est donc sans doute [",liminf," ; ",limsup," ]\n ") cat(" soit, en arrondi, [ ",sprintf("%0.3f",liminf)," ; ",sprintf("%0.3f",limsup)," ]\n ") } # fin de si return(c(liminf,limsup)) } ; # fin de fonction icm ########################################### icmQT <- function(varQT,nivConf=0.05,affichage=FALSE) { ########################################### if (missing(varQT)) { cat("icmQT : intervalle de confiance d'une moyenne sachant ses valeurs\n") ; cat("syntaxe : icmQT(varQT,nivConf=0.05,echo=FALSE) \n") cat("exemples : icm(lng) \n") cat(" icm(AGE_ELF,0.01) \n") return() ; } ; # fin de si n <- length(varQT) m <- mean(varQT) s <- sd(varQT) return( icm(n,m,s,nivConf,affichage) ) } ; # fin de fonction icmQT ########################################### ice <- function(nbVal,ectQt,nivConf=0.05,affichage=FALSE) { ########################################### if (missing(nbVal)) { cat("ice : intervalle de confiance d'un écart-type \n") ; cat("syntaxe : ice(nbVal,ectQt,nivConf=0.05,echo=FALSE) \n") ; cat("exemples : ice(5,3.676955) \n") ; cat(" ice(5,3.676955,0.05) \n") ; cat(" ice(5,3.676955,0.01,TRUE) \n") ; return() ; } ; # fin de si s <- ectQt ddl <- nbVal-1 alpr <- nivConf/2 chir <- qchisq(alpr,ddl,lower.tail=FALSE) f_chir <- sprintf("%6.3f",chir) alpl <- 1-nivConf/2 chil <- qchisq(alpl,ddl,lower.tail=FALSE) f_chil <- sprintf("%6.3f",chil) liminf <- sqrt(ddl*s*s/chir) limsup <- sqrt(ddl*s*s/chil) if (affichage) { cat(" Pour ",nbVal," valeurs d'écart-type ",ectQt," au niveau ",nivConf,"\n") cat(" la valeur critique gauche issue de la loi du chi2 est ",chir," soit ",f_chir,"\n") cat(" la valeur critique droite issue de la loi du chi2 est ",chil," soit ",f_chil,"\n") cat(" l'intervalle de confiance est donc sans doute [",liminf," ; ",limsup," ]\n ") cat(" soit, en arrondi, [ ",sprintf("%0.3f",liminf)," ; ",sprintf("%0.3f",limsup)," ]\n ") } # fin de si return(c(liminf,limsup)) } ; # fin de fonction ice ########################################### iceQT <- function(varQT,nivConf=0.05,affichage=FALSE) { ########################################### n <- length(varQT) s <- sd(varQT) return( ics(n,s,nivConf,affichage) ) } ; # fin de fonction iceQT ########################################### tailleEchProp <- function(pChapeau,margeErreur,nivConf=0.05) { ########################################### if (missing(pChapeau)) { cat("tailleEchProp : taille d'Echantillon pour estimer une proportion\n") ; cat("syntaxe : tailleEchProp(pChapeau,margeErreur,nivConf=0.05) \n") cat("exemples : tailleEchProp(0.169,4) \n") cat(" tailleEchProp(0.169,4,0.10) \n") return() ; } ; # fin de si z <- qnorm(1-nivConf/2) f_z <- sprintf("%6.3f",z) cat(" la valeur critique issue de la loi normale est ",f_z,"\n") m <- margeErreur/100 n1 <- z*z*pChapeau*(1-pChapeau)/(m*m) n2 <- floor(n1) + 1 ; cat("La taille minimale requise est sans doute ",n1," qu'on arrondira certainement en ",n2,"\n\n") return(n2) } ; # fin de fonction tailleEchProp ########################################### tailleEchMoy <- function(ecartType,margeErreur,nivConf=0.05) { ########################################### if (missing(ecartType)) { cat("tailleEchProp : taille d'Echantillon pour estimer une moyenne\n") ; cat("syntaxe : tailleEchMoy(ecartType,margeErreur,nivConf=0.05) \n") cat("exemples : tailleEchMoy(0.64,0.25) \n") cat(" tailleEchMoy(0.64,0.25,0.10) \n") return() ; } ; # fin de si z <- qnorm(1-nivConf/2) f_z <- sprintf("%6.3f",z) cat(" la valeur critique issue de la loi normale est ",f_z,"\n") m <- margeErreur/100 rn <- z*ecartType/margeErreur n1 <- rn*rn n2 <- floor(n1) + 1 ; cat("La taille minimale requise est sans doute ",n1," qu'on arrondira certainement en ",n2,"\n\n") return(n2) } ; # fin de fonction tailleEchMoy ########################################### acpFacteur <- function(titre,dataMatrix,facteur,grFile="") { ########################################### if (missing(titre)) { cat("acpFacteur : une analyse en composante principale avec modalités coloriées\n") ; cat("syntaxe : acpFacteur(titre,dataMatrix,facteur,grFile) \n") return() ; } ; # fin de si # chargement de la librairie library(ade4) # préparation de la liste des couleurs listeCouleurs <- c("blue","red","green","black","orange") couleurs <- rep("blue",length(facteur)) moda <- sort(unique(facteur)) nbMod <- length(moda) for (idm in (1:nbMod)) { couleurs[facteur==moda[idm]] <- listeCouleurs[idm] } # fin pour idm # diagrammes de dispersion colorisés #pairs(dataMatrix,col=couleurs,main=titre) # réalisation de l'acp acp <- dudi.pca(dataMatrix,scannf = FALSE, nf = 5) # on met en forme les valeurs propres # et on trace son histogramme vp <- acp$eig vp_i <- 1:length(vp) vp_r <- round( 100*vp/sum(vp),0 ) vp_rc <- cumsum(vp_r) dsc_vp <- cbind(vp_i,vp,vp_r,vp_rc) colnames(dsc_vp) <- c("Axe","Valeur propre","Pourcent","Cumul") print(dsc_vp) ## barplot(vp_r,col=heat.colors(length(vp_r)),main=titre) # quelques tracés pour mieux comprendre # comment se situent les variables # et les groupes if (!grFile=="") { png(grFile,width=1600,height=1200) } par(mfrow=c(2,2)) barplot(vp_r,col=heat.colors(length(vp_r)),main=titre) scatter(acp,clab.row=0,posieig="none",main=paste(titre,"plan 1-2")) s.class(dfxy=acp$li,fac=facteur,col=listeCouleurs,xax=1,yax=2,add.plot=TRUE) scatter(acp,clab.row=0,posieig="none",main=paste(titre,"plan 1-3")) s.class(dfxy=acp$li,fac=facteur,col=listeCouleurs,xax=1,yax=3,add.plot=TRUE) scatter(acp,clab.row=0,posieig="none",main=paste(titre,"plan 2-3")) s.class(dfxy=acp$li,fac=facteur,col=listeCouleurs,xax=2,yax=3,add.plot=TRUE) if (!grFile=="") { dev.off() cat ("vous pouvez utiliser ",grFile,"\n") } # fin si } ; # fin de fonction acpFacteur #################################################################################### numeroteId <- function(df) { #################################################################################### nbl <- dim(df)[1] ids <- matrix(nrow=nbl,ncol=1) for (idl in (1:nbl)) { ids[idl,1] <- paste("I",sprintf("%04d",idl),sep="") } # fin pour idl ndf <- as.data.frame(cbind(ids,df)) colnames(ndf) <- c("Iden",colnames(df)) return(ndf) } ; # fin de fonction numeroteId #################################################################################### datagh <- function(dossier) { #################################################################################### if (missing(dossier)) { # on prépare la liste des dossiers idossier <- 0 nbdossier <- 8 matdossier <- matrix(nrow=nbdossier,ncol=5) colnames(matdossier) <- c(surncarg(8,"Dossier"),"individus","variables","QT","QL") rownames(matdossier) <- 1:nbdossier idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"ELF") matdossier[idossier,(2:5)] <- c(99,7,1,1) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"PBIO") matdossier[idossier,(2:5)] <- c(419,12,0,12) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"CETOP") matdossier[idossier,(2:5)] <- c(246,122,0,0) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"TITANIC") matdossier[idossier,(2:5)] <- c(2201,5,0,5) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"VINS") matdossier[idossier,(2:5)] <- c(18,8,8,0) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"HER") matdossier[idossier,(2:5)] <- c(80,14,13,1) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"LEA778") matdossier[idossier,(2:5)] <- c(778,29,28,1) idossier <- idossier + 1 matdossier[idossier,1] <- surncarg(8,"IRIS") matdossier[idossier,(2:5)] <- c(150,5,4,1) # et on affiche cat("\n LISTE DES ",nbdossier,"DOSSIERS DISPONIBLES VIA datagh() \n\n") print(matdossier,right=TRUE,print.gap=3,quote=FALSE) cat("\n Pour utiliser un dossier, taper datagh(\"DOSSIER\") ; ") ; cat("\n où DOSSIER est le nom (en minuscules) du dossier ") ; cat("\n par exemple pour elf, taper datagh(\"elf\") ; \n\n") ; return(-1) } else { dossier <- tolower(dossier) lieu1 <- "http://forge.info.univ-angers.fr/~gh/Datasets" lieu2 <- "/home/info/gh/public_html/Datasets" lieu3 <- "/home/gh/public_html/Datasets" lieux <- c(lieu1,lieu2,lieu3,".","K:/Stat_ad","~/Crs/Stat/Data","I:/Crs/Stat/Data") # on cherche les données nomfic <- paste(dossier,".dar",sep="") idl <- 1 vu <- 0 nbl <- length(lieux) while (idl<=nbl) { nf <- paste(lieux[idl],nomfic,sep="/") fc <- file.access(nf,mode=0) if (fc==0) { # le fichier est vu cat("\n DOSSIER ",dossier) cat(" (après lecture de ",nf,") :\n") idl <- nbl + 1 vu <- 1 # récupération de la matrice des données matdata <- lit.dar(nf) } ; # fin de si fc ==0 idl <- idl + 1 } ; # fin de tant que sur idl if (vu==0) { cat("\n dossier ",sQuote(dossier)," non vu...\n") } else { # on cherche le descriptif des données src1 <- "/home/info/gh/Bin/" src2 <- "/home/gh/Bin" src3 <- "http://forge.info.univ-angers.fr/~gh/" src4 <- "http://forge.info.univ-angers.fr/~gh/Datasets" srcs <- c(src1,src2,src3,src4,".","X:/","K:/Stat_ad","~/Crs/Stat/Data","I:/Crs/Stat/Data") nds <- "datagh.r" jdl <- 1 wu <- 0 nbs <- length(srcs) while (jdl<=nbs) { ndf <- paste(srcs[jdl],nds,sep="/") gc <- file.access(ndf,mode=0) if (gc==0) { # le fichier est vu jdl <- nbs + 1 wu <- 1 source(ndf,encoding="latin1") # et on exécute la fonction liée au dossier fdossier <- paste("datagh_",dossier,sep="") ghcmd <- call(fdossier,matdata) eval(ghcmd) # pour elf cela serait datagh_elf(matdata) } ; # fin de si gc ==0 jdl <- jdl + 1 } ; # fin de tant que sur idl if (wu==0) { cat("\n dossier ",sQuote(dossier)," aucun descriptif vu...\n") } ; # fin de si pas vu } ; # fin de si pas vu } ; # finsi sur paramètre manquant } ; # fin de fonction datagh ########################################### # # fonctions de "fainéant de la frappe au clavier" # ########################################### dql <- function(vql,a=FALSE,b="") { ########################################### # # un decritQL automatique tdql <- paste("tql <- t_",vql,sep="") mdql <- paste("m_",vql,sep="") decritQL(eval(tdql),eval(vql),eval(mdql),a,b) } # fin de fonction dql ########################################### dqt <- function(vqt,a=FALSE,b="") { ########################################### # # un decritQT automatique tdqt <- paste("tqt <- t_",vqt,sep="") udqt <- paste("u_",vqt,sep="") decritQT(eval(tdqt),eval(vqt),eval(udqt),a,b) } # fin de fonction dql ########################################### scr <- function(n) { ########################################### # # découpage de l'écran graphique if (n==1) { par (mfrow=c(1,1)) } if (n==2) { par (mfrow=c(1,2)) } if (n==4) { par (mfrow=c(2,2)) } } # fin de fonction scr ########################################### sinksource <- function(fName,encoding="latin1",split=TRUE) { ########################################### # # fait un sink puis un source et ferme le sink # exemple : sinksource("test") ==> exécute test.r et produit test.sor baseName <- fName lng <- nchar(fName) if (substr(fName,lng-1,lng)==".r") { baseName <- substr(fName,1,lng-2) } # fin si fs <- paste(baseName,".sor",sep="") fp <- paste(baseName,".r",sep="") unlink(fs) sink(fs,split=split) source(fp,encoding=encoding) sink() cat(" vous pouvez utiliser ",fs,"\n") } # fin de fonction sinksource ########################################### sigcodes <- "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 'NS' 1" ########################################### as.sigcode <- function(pvalue) { ########################################### # as.sigcode(p) pour mettre # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 'NS' 1 retc <- " NS" if (pvalue<=0.100) { retc <- " ." } if (pvalue<=0.050) { retc <- " *" } if (pvalue<=0.010) { retc <- " **" } if (pvalue<=0.001) { retc <- "***" } return(retc) } # fin de fonction as.sigcode ########################################### sdmax <- function(n,a,b) { ########################################### # renvoie l'écart-type maximal de n points sur [a,b] lesval <- c(rep(a,n/2),rep(b,n/2)) if ((n%%2)==1) { lesval <- c(a,lesval) ; } ; # n impair return(sd(lesval)) } # fin de fonction sdmax ########################################### sdr <- function(x,a,b) { ########################################### # calcule le rapport écart-type sur écart-type maximal # en % return(100*sd(x)/sdmax(length(x),a,b)) } # fin de fonction sdmax ########################################### ifiab <- function(x,a,b) { ########################################### # renvoie l'indice de fiabilité # en % return(100-sdr(x,a,b)) } # fin de fonction ifiab ########################################### bendist <- function(data,distcol=TRUE) { ########################################### # calcule la distance entre les profils de la matrice des données # au sens de Benzécri (distance du chi2 centré) # on utilise les profils-colonnes si distcol=TRUE, # les profils-lignes sinon. # total général du tableau totgen <- sum(data) # fréquences conditionnelles nblig <- dim(data)[1] nbcol <- dim(data)[2] frq <- matrix(nrow=nblig,ncol=nbcol) fql <- vector(length=nblig,mode="numeric") fqc <- vector(length=nbcol,mode="numeric") somLig <- apply(data,1,sum) somCol <- apply(data,2,sum) pdlig <- 1:nblig pdcol <- 1:nbcol for (idl in pdlig) { for (jdc in pdcol) { frq[idl,jdc] <- data[idl,jdc] / totgen } # fin pour jdc } # fin pour idl for (idl in pdlig) { fql[idl] <- somLig[idl] / totgen } # fin pour idl for (jdc in pdcol) { fqc[jdc] <- somCol[jdc] / totgen } # fin pour jdc if (distcol) { # Distances colonnes distB <- matrix(nrow=nbcol,ncol=nbcol) row.names(distB) <- colnames(data) colnames(distB) <- colnames(data) for (ja in pdcol) { for (jb in pdcol) { distB[ja,jb] <- 0 } # fin pour ja } # fin pour ja for (ja in pdcol) { for (jb in pdcol) { ladist <- 0 for (i in pdlig) { denx <- fqc[ja]*fql[i] if (denx==0) { x <- 0 } else { x <- frq[i,ja] / denx } # fin si deny <- fqc[jb]*fql[i] if (deny==0) { y <- 0 } else { y <- frq[i,jb] / deny } # fin si ladist <- ladist + fql[i] * (x-y) * (x-y) #cat(ja,jb,i,fqc[ja],frq[i,ja],fqc[jb],fql[i],frq[i,jb],ladist,"\n") distB[ja,jb] <- ladist } # fin pour i } # fin pour ja } # fin pour ja } else { # Distances lignes distB <- matrix(nrow=nblig,ncol=nblig) row.names(distB) <- row.names(data) colnames(distB) <- row.names(data) for (ia in pdlig) { for (ib in pdlig) { distB[ia,ib] <- 0 } # fin pour ib } # fin pour ia for (ia in pdlig) { for (ib in pdlig) { ladist <- 0 for (j in pdcol) { denx <- fql[ia]*fqc[j] if (denx==0) { x <- 0 } else { x <- frq[ia,j] / denx } # fin si deny <- fql[ib]*fqc[j] if (deny==0) { y <- 0 } else { y <- frq[ib,j] / deny } # fin si ladist <- ladist + fqc[j] * (x-y) * (x-y) } # finpour j distB[ia,ib] <- ladist } # fin pour ib } # fin pour ia } # fin si return(as.dist(distB,diag=TRUE)) } # fin de fonction bendist ######################################################################## decritClass <- function(cla) { ######################################################################## noms <- cla$labels fusi <- cla$merge haut <- cla$height nbe <- dim(fusi)[1] nomc <- rep("",nbe) nbc <- 0 cat("Classe Hauteur Contenu \n") for (ide in (1:length(noms))) { nbc <- nbc + 1 cat(sprintf("%5d",nbc)) cat(" ",noms[ide]) cat("\n") } # fin pour ide for (ide in (1:nbe)) { nbc <- nbc + 1 cat(sprintf("%5d",nbc)) cat(sprintf(" %6.0f ",haut[ide])) #cat(" ",fusi[ide,]) cdlc <- c("") for (edc in fusi[ide,]) { if (edc<0) { cdlc <- paste(cdlc,sprintf("%-8s",noms[-edc]),sep="") } else { cdlc <- paste(cdlc,sprintf("%-8s",nomc[edc]),sep="") } # fin si } # fin pour cat(cdlc) nomc[ide] <- cdlc #cat(" donc en ",ide," : ", nomc[ide]) cat("\n") } # fin pour ide } # fin de fonction decritClass ######################################################################## lit.dac <- function(fn) { ######################################################################## if (missing(fn)) { cat(" syntaxe : lit.dac( nom_fichier ) \n") cat(" cette fonction renvoie une liste dont les composantes $nomgrp, $grps et $dac \n") cat(" donnent respectivement les noms des groupes, les numéros de groupes et les données ;\n") cat(" la composante $basename contient la base du nom des fichiers à utiliser.\n ") return() } # fin si # préparation des noms de fichier cat(" Base des noms de fichier : ",fn,"\n") fngr <- paste(fn,".ngr",sep="") fgal <- paste(fn,".gal",sep="") fdac <- paste(fn,".dac",sep="") ngrp <- read.table(fngr,as.is=TRUE) nomgrp <- ngrp[,2] cat(" ",length(nomgrp)," noms de groupes lus dans ",fngr,"\n") cols <- read.table(fgal,as.is=TRUE) nomcol <- cols[,2] cat(" ",length(nomcol)," noms de colonnes lus dans ",fgal,"\n") dac <- read.table(fdac,row.names=1) grps <- dac[,1] dac <- dac[,-1] colnames(dac) <- nomcol cat(" ",dim(dac)[1]," lignes et ",dim(dac)[2],"colonnes dans ",fdac,"\n") return( list(nomgrp=nomgrp,grps=grps,dac=dac,basename=fn) ) } # fin de fonction lit.dac ######################################################################## verifSol <- function(data,numg,formule) { ######################################################################## if (missing(data)) { cat(" syntaxe : verifSol(data,numg,formule) ) \n") cat(" cette fonction vérifie si la formule et le groupe correspondent.\n") cat(" exemple : verifSol(raphv100,1,c(-30,63)) \n") return() } # fin si cat("Vérification de la formule ",formule," pour le groupe ",numg," dans les données ",data$basename,"\n") eff <- table(data$grps)[numg] cat(" nombre de lignes dans le groupe ",numg," : ",eff,"\n") pdi <- 1:length(formule) for (ind in pdi) { edf <- formule[ind] if (ind==1) { if (edf>0) { cnd <- data$dac[,edf]==1 } else { cnd <- data$dac[,-edf]==0 } # fin si } else { if (edf>0) { cnd <- cnd & data$dac[,edf]==1 } else { cnd <- cnd & data$dac[,-edf]==0 } # fin si } # fin si } # fin pour nbv <- sum(cnd) cat(" nombre de lignes vérifiant la formule : ",nbv,"\n") w1 <- which(data$grps==numg) w2 <- which(cnd) cat(" Groupe : ",w1,"\n") cat(" Formule : ",w2,"\n") oui <- 0 if ((length(w1)==length(w2))) { if (sum(abs(w1-w2))==0) { oui <- 1 } # fin si } # fin si if (oui==1) { cat(" oui, numéro de groupe et formule correspondent.\n") } else { cat(" non, numéro de groupe et formule ne correspondent pas.\n") } # fin si } # fin de fonction verifSol ############################################################################## ############################################################################## ## ## fonctions pour la statistique lexicale ## ############################################################################## ############################################################################## ############################################################################## ucfirst <- function(chaine) { # renvoie la chaine avec l'initiale en majuscule ############################################################################## return(paste(toupper(substr(chaine,1,1)),tolower(substr(chaine,2,nchar(chaine))),sep="")) } # fin de fonction ucfirst ############################################################################## analexies<- function(texte,minu=TRUE,dbg=FALSE,show=TRUE) { ############################################################################## # analyse lexicale d'un texte : on détermine les mots du texte # et leurs occurrences library(hash) if (missing(texte)) { cat("Syntaxe : res <- analexies(texte) \n") cat(" res$nmots contient le nombre de mots, res$nmotsd le nombre de mots différents, \n") cat(" res$tmots est le tableau des mots et res$toccs le tableau des mots trié par fréquence ;\n") cat(" de plus res$hmots est un tableau associatif (\"hash\") des mots. \n") return() } # fin si # nettoyage du texte tac <- texte if (minu) { tac <- tolower(tac) } tac <- gsub("[';.?:!(),]"," ",tac) tac <- gsub("--","",tac) tac <- gsub("\""," ",tac) tac <- gsub("[«»]","",tac) if (dbg) { cat("Texte avant : \n") cat(texte) cat("\n") cat("Texte après : \n") cat(tac) cat("\n") } # fin de si # découpage en mots tbdm <- strsplit(tac," ",perl=TRUE) nbm <- 0 nbmd <- 0 hdm <- hash() for (mot in tbdm[[1]]) { if (mot!="") { nbm <- nbm + 1 if (has.key(mot,hdm)) { hdm[[mot]] <- hdm[[mot]] + 1 } else { hdm[[mot]] <- 1 nbmd <- nbmd + 1 } # fin si } # fin si } # fin pour # structuration en tableau tdm <- matrix(nrow=nbmd,ncol=1) ldm <- matrix(nrow=nbmd,ncol=1) idm <- 0 for (mot in keys(hdm)) { idm <- idm + 1 ldm[idm,1] <- mot tdm[idm,1] <- hdm[[mot]] } # fin pour row.names(tdm) <- ldm colnames(tdm) <- "occ" # tableau des occurences décroissantes idx <- rev(order(tdm[,1])) tdo <- matrix(nrow=nbmd,ncol=1) tdo[,1] <- tdm[idx] row.names(tdo) <- ldm[idx] colnames(tdo) <- "occ" if (show) { cat(" ",nbm,"mots en tout dont ",nbmd," mots distincts.\n") cat(paste(" chaque mot est donc répété ",sprintf("%9.2f",nbm/nbmd)," fois en moyenne.",sep=""),"\n") cat("\n") } # fin si if (dbg) { for (mot in keys(hdm)) { cat(sprintf(" %-30s",mot),sprintf(" %5d",hdm[[mot]]),"\n") } # fin pour } # fin de si return(list(nmots=nbm,hmots=hdm,nmotsd=nbmd,tmots=tdm,toccs=tdo)) } # fin de fonction analexies ############################################################################## lit.texte <- function(fn) { ############################################################################## if (missing(fn)) { cat(" syntaxe : txt <- lit.texte(nom_de_fichier) \n") cat(" renvoie une variable chaine de caractères qui contient tout le texte ;\n") cat(" on peut ensuite l'analyser avec analexies(txt).\n") return() } # fin si if (!file.exists(fn)) { cat("désolé, le fichier ",fn," n'existe pas.\n") return() } # fin si ftxt <- readLines(fn,n=-1) leTxt <- "" for (lig in ftxt) { leTxt <- paste(leTxt,lig) } # fin pour return(leTxt) } # fin de fonction lit.texte ############################################################################## cpt <- function(ldf,ldm,fs="") { ############################################################################## library(hash) if (missing(ldm)) { cat(" syntaxe : tdo <- cpt(fichiers,mots[,fic_sor]) \n") cat(" renvoie un tableau qui dénombre les mots indiqués dans les fichiers.\n") cat(" le troisième paramètre stocke les résultats dans un fichier texte.\n") cat(" exemples : tdo <- cpt(\"bete.txt candide.txt\",\"maison temps train\") \n") cat(" tdf <- cpt(\"bete.txt candide.txt\",\"maison temps train\",\"resultat.txt\")\n") return() } # fin si # préparation des noms de fichier fics <- strsplit(ldf," ",perl=TRUE) nbf <- 0 nldf <- c() for (nfic in fics[[1]]) { if (nfic!="") { nbf <- nbf + 1 #cat(" fichier ",nbf," : ",nfic,"\n") nldf <- c(nldf,nfic) } # fin si } # fin pour # préparation des mots mots <- strsplit(ldm," ",perl=TRUE) nbm <- 0 nldm <- c() for (mot in mots[[1]]) { if (mot!="") { nbm <- nbm + 1 #cat(" mot ",nbm," : ",mot,"\n") nldm <- c(nldm,mot) } # fin si } # fin pour # initialisation du tableau résultat tres <- matrix(nrow=nbf,ncol=nbm,rep(0,nbm*nbf)) row.names(tres) <- nldf colnames(tres) <- nldm # passage en revue des fichiers via analexies idf <- 0 for (fic in nldf) { idf <- idf + 1 cat(" fichier ",idf," / ",nbf," : ",fic,"\n") adf <- analexies(lit.texte(fic)) idm <- 0 for (mot in nldm) { idm <- idm + 1 if (has.key(mot,adf$hmots)) { tres[idf,idm] <- adf$hmots[[mot]] } # fin si } # fin pour mot } # fin pour fic if (fs!="") { sink(fs) print(tres,quote=FALSE) sink() cat(" vous pouvez utiliser ",fs,"\n") } # fin si return(tres) } # fin de fonction cpt ############################################################################## pchShow <- function(extras = c("*",".", "o","O","0","+","-","|","%","#"), cex = 3, ## good for both .Device=="postscript" and "x11" col = "red3", bg = "gold", coltext = "brown", cextext = 1.2, main = paste("plot symbols : points (... pch = *, cex =", cex,")")) { ############################################################################## nex <- length(extras) np <- 26 + nex ipch <- 0:(np-1) k <- floor(sqrt(np)) dd <- c(-1,1)/2 rx <- dd + range(ix <- ipch %/% k) ry <- dd + range(iy <- 3 + (k-1)- ipch %% k) pch <- as.list(ipch) # list with integers & strings if(nex > 0) pch[26+ 1:nex] <- as.list(extras) plot(rx, ry, type="n", axes = FALSE, xlab = "", ylab = "", main = main) abline(v = ix, h = iy, col = "lightgray", lty = "dotted") for(i in 1:np) { pc <- pch[[i]] ## 'col' symbols with a 'bg'-colored interior (where available) : points(ix[i], iy[i], pch = pc, col = col, bg = bg, cex = cex) if(cextext > 0) text(ix[i] - 0.3, iy[i], pc, col = coltext, cex = cextext) } } # fin de fonction pchShow ############################################################################## ############################################################################## ## ## fonctions pour l'étude de la discordance ## ############################################################################## ############################################################################## ############################################################################## rchModeleLogistique <- function(df) { ############################################################################## # la variable à prédire est forcément en colonne 1 cible <- df[,1] library(ROCR) mdc(df,meilcor=FALSE) # étude du modèle complet, # des modèles à 1 variable # des modèles à n-1 variables nbvar <- ncol(df) - 1 # car on retire la cible fnbv <- sprintf("%02dV",nbvar) nbcdr <- 3 + nbvar # car on rajoute AIC en colonne 1, il y a la constante et AUROC en fin nbmod <- 1 + 2*nbvar + 3 # les 3 derniers modèles sont issus de step matres <- matrix(nrow=nbmod,ncol=nbcdr) colnames(matres) <- c("AIC","(Intercept)",colnames(df)[-1],"AUROC") rownames(matres) <- c(paste("M",fnbv,sprintf("%02d",0:nbvar),sep=""), paste("M",fnbv,sprintf("%03d",-(1:nbvar)),sep=""), paste("M",fnbv,c("BS","FS","MS"),sep="")) pdv <- 1:nbmod # boucle sur les modèles models <- glm( cible ~ . ,data=df[,-1],family="binomial") # modèle saturé utilisé par step model0 <- glm( cible ~ 1 ,data=df[,-1],family="binomial") # modèle minimal utilisé par step for (idm in pdv) { tmpData <- df[,-1] if (idm==1) { # modèle complet modele <- glm( cible ~ .,data=tmpData,family="binomial") } # fin si if ((idm>1) & (idm<=nbvar+1)) { # modèles à 1 variable tmpData <- df[,idm,drop=FALSE] modele <- glm( cible ~ .,data=tmpData,family="binomial") } # fin si if ((idm>nbvar+1) & (idmym) { cm <- seuil ; ym <- you ; sm <- sensi ; sp <- speci ; vppm <- vpp ; vpnm <- vpn ; vsiy <- seuil ; } ; # fin de si if ((vpp>vcip) & (na+nb>0) ) { vcip <- vpp ; vsip <- seuil ; } ; if ((vpn>vcin) & (nc+nd>0) ) { vcin <- vpn ; vsin <- seuil ; } ; if (voa>vcoa) { vcoa <- voa ; vsoa <- seuil ; } ; # stockage pour les autres boucles ts_na[iseuil] <- na ; ts_nb[iseuil] <- nb ; ts_nc[iseuil] <- nc ; ts_nd[iseuil] = nd ; ts_spec[iseuil] = speci ; ts_sens[iseuil] = sensi ; ts_vpp[iseuil] = vpp ; ts_vpn[iseuil] = vpn ; ts_you[iseuil] = you ; ts_ova[iseuil] = voa ; iseuil <- iseuil + 1 } ; # fin de tant que ### ### # deuxième boucle : on affiche ### cols <- chaineEnListe("CUTOFF CRIT VP FP FN VN TOTAL SPEC 1-SPEC SENS VPP VPN YOUDEN OA ") ; nbc <- length(cols) tabRes <- matrix(nrow=nbseuils,ncol=nbc) colnames(tabRes) <- cols #cat(" CUTOFF CRIT VP FP FN VN TOTAL SPEC 1-SPEC SENS VPP VPN YOUDEN OA\n") ; iseuil <- 1 ; pseuil <- -1 ; is_iyou <- 0 ; is_iova <- 0 ; is_ivpp <- 0 ; is_ivpn <- 0 ; while (iseuil <= nbseuils) { seuil <- vseuils[iseuil] ; # comptage des biens et mal classés # nbl est le nombre de données # en idl, sco[idl] est la valeur de la QT et fcs[idl] la valeur de la QL # pour savoir si on vaut 0 ou 1 on utilise la variable seuil na <- ts_na[iseuil] ; nb <- ts_nb[iseuil] ; nc <- ts_nc[iseuil] ; nd <- ts_nd[iseuil] ; tabRes[iseuil,"CUTOFF"] <- sprintf("%9.4f",seuil) tabRes[iseuil,"VP"] <- na tabRes[iseuil,"FP"] <- nb tabRes[iseuil,"FN"] <- nc tabRes[iseuil,"VN"] <- nd tabRes[iseuil,"TOTAL"] <- na+nb+nc+nd sensi <- ts_sens[iseuil] ; speci <- ts_spec[iseuil] ; vpp <- ts_vpp[iseuil] ; vpn <- ts_vpn[iseuil] ; you <- ts_you[iseuil] ; ova <- ts_ova[iseuil] ; umspi <- 1 - speci ; #ls <- paste(ls,sprintf("%9.3f",seuil)) cri <- "" ; if ((idsvpp0p9==0) & (vpp>=0.9)) { cri <- "VPP0.9" idsvpp0p9 <- iseuil } # fin si if ((idsvpn0p9==0) & (vpn<=0.9)) { cri <- "VPN0.9" idsvpn0p9 <- iseuil } # fin si if (seuil==vsio) { cri <- "IO " ; } ; if (seuil==vsoa) { cri <- "OA " ; is_ova <- iseuil ; } ; if (seuil==vsip) { cri <- "IP " ; is_vpp <- iseuil ; } ; if (seuil==vsin) { cri <- "IN " ; is_vpn <- iseuil ; } ; if (seuil==vsiy) { cri <- "IY " ; is_you <- iseuil ; } ; #ls <- paste(ls," ",cri) ; tabRes[iseuil,"CRIT"] <- cri #ls <- paste(ls,sprintf("%6d",na),sprintf("%6d",nb),sprintf("%6d",nc),sprintf("%6d",nd)) ; #ls <- paste(ls,sprintf(" %6d ",(na+nb+nc+nd))) ; #ls <- paste(ls,sprintf("%12.7f",speci)) ; #ls <- paste(ls,sprintf("%12.7f",umspi)) ; #ls <- paste(ls,sprintf("%12.7f",sensi)) ; #ls <- paste(ls,sprintf("%12.7f",vpp)) ; #ls <- paste(ls,sprintf("%12.7f",vpn)) ; #ls <- paste(ls,sprintf("%12.7f",you)) ; #ls <- paste(ls,sprintf("%12.7f",ova)) ; #ls <- paste(ls,"\n") ; tabRes[iseuil,"SPEC"] <- sprintf("%11.6f",speci); tabRes[iseuil,"1-SPEC"] <- sprintf("%11.6f",umspi) ; tabRes[iseuil,"SENS"] <- sprintf("%11.6f",sensi) ; tabRes[iseuil,"VPP"] <- sprintf("%11.6f",vpp) ; tabRes[iseuil,"VPN"] <- sprintf("%11.6f",vpn) ; tabRes[iseuil,"YOUDEN"] <- sprintf("%11.6f",you) ; tabRes[iseuil,"OA"] <- sprintf("%11.6f",ova) ; n_seuil <- seuil + 0.0 ; n_pseuil <- pseuil + 0.0 ; #if (n_seuil!=n_pseuil) { cat(ls) ; } ; seuil <- seuil + pas ; pseuil <- seuil ; iseuil <- iseuil + 1 } ; # fin de tant que tabRes <- as.data.frame(tabRes) colnames(tabRes) <- cols cat("\nListe des seuils et de leurs caractéristiques \n\n") print(tabRes,row.names=FALSE) cat("\nSeuils remarquables parmi les ",nbseuils," seuils vus\n\n") ; tabRes2 <- tabRes[ tabRes[,"CRIT"]!="", ] print(tabRes2,row.names=FALSE) } # fin de fonction sensSpec ############################################################################## ############################################################################## ############################################################################## ## ## autres fonctions ## ############################################################################## ############################################################################## ################################################################################################ lesColonnes <- function(df) { ################################################################################################ nbc <- ncol(df) nbl <- nrow(df) cat("Voici les ",nbc,"colonnes (stat. sur ",nbl," lignes en tout)\n") noms <- names(df) nbvs <- rep(NA,nbc) mins <- rep(NA,nbc) maxs <- rep(NA,nbc) ddc <- cbind((1:nbc),nbvs,mins,maxs) row.names(ddc) <- noms colnames(ddc) <- c("Num","NbVal","Min","Max") pdc <- 1:nbc for (idc in pdc) { laCol <- as.numeric(df[,idc]) laCol <- laCol[!is.na(laCol)] ddc[idc,2] <- length(laCol) ddc[idc,3] <- min(laCol) ddc[idc,4] <- max(laCol) } # fin pour idx <- order(noms) ddc <- ddc[idx,] print.data.frame(as.data.frame(ddc),quote=FALSE,row.names=TRUE) cat("\n") } # fin de fonction lesColonnes ############################################################################## rownames <- function(nom) { return( row.names(nom) ) } ############################################################################## gr <- function(fichier,largeur=1600,hauteur=1200) { ############################################################################## png(fichier,width=largeur,height=hauteur) } # fin de fonction gr ############################################################################## z <- function() { quit("no") } ############################################################################## versiongh <- function() { ############################################################################## cat(" statgh.r : ",version_gH," pour ",R.Version()$version.string,"\n") } # fin de fonction versiongh ############################################################################## aide <- function() { ############################################################################## cat("\n\n (gH) version ",version_gH,"\n\n") cat("fonctions d'aides : lit() fqt() fql() ic() fapprox() chi2() fcomp() datagh() \n") ; cat("taper aide() pour revoir cette liste\n\n") } ; # fin de fonction aide ############################################################################## ############################################################################## aide() ## fin de statgh.r