Cours complet d’initiation au logiciel R pas à pas

Cours complet d’initiation au logiciel R pas à pas

LOGICIEL R 

VERSION 2.6.0.

  PRISE EN MAIN

STATISTIQUES DE BASE ET

DEMARCHE EXPERIMENTALE

ANALYSE DE DONNEES

Non-supervisées

Supervisées

2008-2009

Yves Le Roux 

Version 2 (01/02/2008)

(avec utilisation de multiples supports disponibles sur Internet)

LOGICIEL R

R est un logiciel de Statistique créé par Ross Ihaka & Robert Gentleman (1996, J. Comput. Graph. Stat.,5 :299-314). R est à la fois un langage et un environnement : les commandes sont exécutées grâce à des instructions codées dans un langage relativement simple, les résultats sont affichées sous forme de texte et les graphiques sont visualisées directement dans une fenêtre qui leur est propre. C’est un clone de S-plus basé sur le langage de programmation orienté objet S développé par les AT&T Bell Laboratories en 1988.

Ce logiciel sert à manipuler des données, à tracer des graphiques et à faire des analyses statistiques sur ces données.

Pourquoi utiliser R? Tout d’abord R est gratuit !

Mais c’est aussi un outil très puissant et très complet, particulièrement bien adapté pour la mise en oeuvre informatique de méthodes statistiques. Il est plus difficile d’accès que certains logiciels car il n’est pas conçu pour être utilisé à l’aide de “clics” dans des menus. Il faut au contraire passer du temps à apprendre la syntaxe et les commandes du logiciel.

L’avantage est double : l’approche est pédagogique puisqu’il faut maîtriser la méthode statistique pour la mettre en oeuvre ; l’outil est très efficace lorsque l’on domine le langage S (on peut alors faire des manipulations très sophistiquées sur les données).

R est une suite intégrée de composants logiciels pour la manipulation de données, le calcul et l’affichage de graphiques. Parmi d’autres choses, il a :

-  des procédures efficaces de traitement des données et des capacités de stockage de ces données

-  une suite d’opérateurs pour des calculs sur des tableaux et en particulier sur des matrices

-  une vaste et cohérente collection d’outils intermédiaires pour l’analyse de données

-  des capacités graphiques pour l’analyse de données et des possibilités d’affichage à l’écran ou sur papier

-  un langage de programmation simple et efficace intégrant les conditions, les boucles, la récursivité, et des possibilités d’entrée-sortie. (En fait, la plupart des fonctions du système sont elles-mêmes écrites dans le langage S).

R est plutôt un véhicule pour le développement interactif de nouvelles méthodes d’analyse de données. La compatibilité avec les versions antérieures n’est pas toujours respectée. Par conséquent, il vaut mieux considérer les programmes R comme éphémères.

R est un logiciel de calcul statistique distribuée selon la licence GNU General Public License. La version officielle actuelle est la version R 2.6.0. Elle est disponible dans les archives du réseau CRAN (Comprehensive R Archive Network) dont les sites de base sont :

On peut utiliser le miroir de l’université 

Je remercie l’ensemble des personnes, groupes, communautés qui permettent, en autorisant l’utilisation de leurs cours-exemples-tutoriels, une approche opérationnelle de R sans pour autant en oublier les fondamentaux de la théorie statistiques.

En particulier, ce manuel s’inspire (parfois beaucoup, parfois moins) de :

Toute la suite des documents de formation de Mr. Bouchier, formateur INRA

Toute la suite des documents de formation (cours ou TD) proposés par une excellente équipe pédagogique de l’Univeristé de Lyon

D’autres contributions 

.

Forums

Un forum français d’aide

Un forum anglais d’aide

LOGICIEL R 

VERSION 2.6.0.

PRISE EN MAIN

Environnement de travail

Dans cible (clique droit) de l’icône R sur le bureau qui lance le programme : 

"C:\Program Files\R\R-2.5.1\bin\" R_PROFILE="c:\logiciel R\"

Dans le répertoire de travail (ici c:\logiciel R\), créer un fichier texte

Mettre dans ce fichier :

print("Cher utilisateur de R : BONJOUR !") # Affichage à l’ouverture de R setwd("C:/Logiciel R") # definition du répertoire de travail

cat(paste("Working directory =",getwd(),"\n")) # affichage du repertoire de travail cat("Rprofile file =",Sys.getenv("R_PROFILE"),"\n")

dirperm <- dirname(Sys.getenv("R_PROFILE"))

Des fonctions très importantes

getwd() : identification du répertoire de travail

library() : permet de lister les différentes librairies disponibles

ls() : La commande fait la liste des objets de l’espace de travail.

() : caractérisation du type de variable en mémoire

setwd() : changement du répertoire de travail en indiquant l’adressage

setwd(dirname(file.choose())) : choix du ficher pour changer le répertoire de travail

setwd("~") : répertoire de travail de l’utilisateur

R manipule des objets

Les objets ont tous deux attributs intrinsèques : le mode et la longueur.

Le mode

il en existe quatre principaux : numérique, caractère, complexe, et logique, les objets peuvent être aussi des fonction ou expression

Un vecteur est une variable.Un facteur est une variable catégorielle. Un arrayest un tableau à k dimensions. Une matrice est un cas particulier d’array avec k=2. Un data.frameest un tableau de données composé de un ou plusieurs vecteurs et/ou facteurs ayant tous la même longueur mais pouvant être de modes différents.Un ts est un jeu de données de type séries temporelles (time series) et comporte donc des attributs supplémentaires comme la fréquence et les dates. Une listepeut contenir n’importe quel type d’objet.

Manipulation de données

Quelques exemples d’utilisation de ces fonctions

2+2

[1] 4

x <- 2 + 2 10 * x

[1] 40

seq(2,12,by=2) [1] 2 4 6 8 10 12

seq(4,5,length=5)

[1] 4.00 4.25 4.50 4.75 5.00

rep(1, 30)

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

seq(length=9, from=1, to=5)

[1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0

gl(3, 5)

[1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

Lister et effacer les objets en mémoire

Lister

name <- "Laure" ; n1 <- 10 ; n2 <- 100 ; m <- 0.5 ls()

[1] "m" "n1" "n2" "name"

() : caractérisation du type de variable

ls(pat = "m") : ne liste que les variables contenantla lettre m ls(pat = "^m") : ne liste que les variables commençant par la lettre m

Effacer

Effacement sélectif : rm(m)

Effacement total : rm(list=ls())

Aide en ligne

help.start()

apropos("anova")

Expressions mathématiques de base

Gestion des indéterminations (évite ainsi les plantages…) de la division par 0 et des données manquantes

Inf + 1 [1] Inf 1/Inf [1] 0 Inf/0 [1] Inf Inf - Inf [1] NaN

Inf/Inf [1] NaN 1/0 [1] Inf -1/0 [1] -Inf 0/0 [1] NaN NaN pour Not a Number.

log(-1) [1] NaN NA + 3 [1] NA 4 * NA [1] NA NA pour Not Available.

Notions de Fonctions

De nombreuses fonctions mathématiques sont disponibles

Nous avons vu des exemples d’objets contenant des valeurs numériques. Ces objets sont des instances de classe « numeric » :

x <- 2 + 2 class(x)

[1] "numeric"

Notion de vecteurs

c(2,10,5,37,4,8) [1] 2 10 5 37 4 8 Ages<-c(10,18,25,60) ages[c(3,1)] [1] 25 10 ages[-3] [1] 10 18 60 v1 <- 1:4 v2 <- c(5,6,1,2) v1+v2 [1] 6 8 4 6 x=c(1,4,9)

y=c(x,2,3) y [1] 1 4 9 2 3 Un vecteur peut-être de type caractère letters[c(3, 5, 7)] [1] "c" "e" "g" ind <- c(3, 5, 7) letters[ind] [1] "c" "e" "g" letters[8:13] [1] "h" "i" "j" "k" "l" "m" letters[c(1, 2, 1, 2)] [1] "a" "b" "a" "b"

L’utilisation d’entiers négatifs permet d’exclure les éléments correspondants.

letters[-5]

[1] "a" "b" "c" "d" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q"

[17] "r" "s" "t" "u" "v" "w" "x" "y" "z"

ind <- c(3, 5, 7) letters[-ind]

[1] "a" "b" "d" "f" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"

[17] "t" "u" "v" "w" "x" "y" "z"

x <- 1:10 names(x) <- letters[1:10]

x

a b c d e f g h i j

1 2 3 4 5 6 7 8 9 10

x[c("a", "b", "c", "f")]

a b c f 1 2 3 6

Logiciel R version 2.4.1 (2006-12-18) –le 2007- Création de Tableaux et matrices

Les matrices (et plus généralement les tableaux de dimensions quelconques) sont stockées comme des vecteurs ayant des dimensions :

Deux solutions pour le même résultats

x <- 1:12 dim(x) <- c(3, 4) x

 ou

x = matrix(c(1 :12),ncol=4)

x 

Les matrices sont créées avec la fonction matrix() à partir d'un vecteur. On doit fixer le nombre de colonnes ncol et/ou le nombre de lignes nrow.

x = matrix(c(2,3,5,7,11,13),ncol=2) x

            [,1]      [,2]

[1,]     2           7

[2,]     3           11

[3,]     5           13

nrow(x)

[1] 3

ncol(x)

[1] 2

Par défaut la matrice est remplie colonne par colonne. Pour remplir ligne par ligne, on ajoute l'argument byrow=T

y = matrix(c(2,3,5,7,11,13),ncol=2, byrow=T) y

                [,1]      [,2]

[1,]     2           3

[2,]     5           7

[3,]    11          13

On peut donner des noms aux lignes et aux colonnes avec les fonctions rownames() et colnames()

colnames(x) <- paste("colonne", 1:ncol(x)) x

     colonne1 colonne2 colonne3 colonne4

[1,]           1     4     7    10

rownames(x) <- paste("ligne", 1:nrow(x)) x

       colonne1 colonne2 colonne3 colonne4

ligne1        1     4        7    10 ligne2        2            5             8            11 ligne3 3            6             9            12

y <- matrix(1:12, nrow = 2, ncol = 6) class(y) [1] "matrix" y

                [,1]      [,2]      [,3]      [,4]     [,5]       [,6]

[1,]      1          3          5          7         9          11

[2,]      2          4          6          8         10        12

Attention : si la dimension du vecteur n'est pas égale au produit (ncol * nrow) alors l'opération effectuée est la suivante :

matrix(c(1:3), ncol=2)

                [,1]      [,2]

[1,]      1          3

[2,]      2          1

            Extraire les composantes >8 ou <2 Extraire des éléments d’une matrice ou vect[(vect>8) | (vect<2)] d’un vecteur          [1] 1 9

vect=c(1:9)      Extraire les composantes >8 et <10 vect     vect[(vect>8) & (vect<10)]

[1] 1 2 3 4 5 6 7 8 9                                                  [1] 9

             mat=matrix(vect,ncol=3,nrow=3)       

mat                                                                            mat[2,1]

                [,1]      [,2]      [,3]                                         [1] 2

[1,]     1           4          7                                            mat[,1]

[2,]     2           5          8                                            [1] 1 2 3

[3,]     3           6          9                                            mat[, 3]

                                                                                  [1] 3 6 9

                                                                                  mat[2:3,1:2] #sous matrice

Extraire les composantes >8                        [,1]     [,2] vect[vect>8]         [1,] 2           5 [1] 9 [2,]     3     6 Manipulation des données et Calcul arithmétique et fonctions simples

Accéder à une valeur particulière d’un objet : 

Pour accéder à, par exemple, la 3ème valeur d’un vecteur x, on tape x[3].

Si x est une matrice ou un data.frame, on accédera à la valeur de la ième ligne et jème colonne par x[i,j]. Pour changer toutes les valeurs de la 3ème colonne, on peut taper : > x[,3] <- 10.2 Ce système d’indexation se généralise facilement pour les array, on aura alors autant d’indices que l’array a de dimensions (par exemple pour une array à trois dimensions : x[i,j,k], x[,,3], ...). Il faut retenir que l’indexation se fait à l’aide de crochets, les parenthèses étant réservées pour les arguments d’une fonction

L’indexation peut aussi être utilisée pour supprimer une (des) ligne(s) ou colonne(s). Par exemples, x[-1,] supprimera la 1ère ligne, ou x[-c(1,15),] fera de même avec les 1ère et 15ème lignes.

Pour les vecteurs, matrices et array il est possible d’accéder aux valeurs de ces éléments à l’aide d’une expression de comparaison en guise d’indice :

x <- 1:10 x[x >= 5] <- 20

x

[1] 1 2 3 4 20 20 20 20 20 20 x[x == 1] <- 25

x

[1] 25 2 3 4 20 20 20 20 20 20

Les six opérateurs de comparaison utilisés par R sont : < (inférieur à), > (supérieur à), <= (inférieur ou égal à), >= (supérieur ou égal à), == (égal), et != (différent de). On notera que ces opérateurs ont pour résultat de retourner une variable de mode logical (TRUE ou FALSE).

Calcul arithmétique et fonctions simples

Il existe de nombreuses fonctions dans R pour manipuler les données. La plus simple, on l’a vue plus haut, est c() qui concatène les objets énumérés entre parenthèses. Par exemple :

c(1:5, seq(10, 11, 0.2))

[1] 1.0 2.0 3.0 4.0 5.0 10.0 10.2 10.4 10.6 10.8 11.0

Les vecteurs peuvent être manipulés selon des expressions arithmétiques classiques :

x <- c(1,2,3,4) y <- c(1,1,1,1) z <- x + y z

[1] 2.0 3.0 4.0 5.0

Des vecteurs de longueurs différentes peuvent être additionnés, dans ce cas le vecteur le plus court est “recyclé”. Exemples :

x <- c(1,2,3,4) y <- c(1,2) z <- x+y

z

[1] 2 4 4 6

x <- c(1,2,3) y <- c(1,2) z <- x+y Warning message: longer object length

is not a multiple of shorter object length in: x + y z

[1] 2 4 4

On notera que R a retourné un message d’avertissement et non pas un message d’erreur, l’opération a donc été effectuée. Si l’on veut ajouter (ou multiplier) la même valeur à tous les

éléments d’un vecteur :

x <- c(1,2,3,4) a <- 10 z <- a*x

z

[1] 10 20 30 40

Opérations sur les matrices

x = matrix(c(1 :12),ncol=4)

x 

                [,1]      [,2]      [,3]       [,4]

[1,]      1          4          7             10

[2,]      2          5          8             11

[3,]      3          6          9             12

t(x)

                [,1]      [,2]      [,3]

[1,]     1           2          3

[2,]     4           5          6

[3,]     7           8          9

[4,]    10          11        12

x + x

	[,1]	[,2]	[,3]	[,4]
[1,]	2	4	6	8
[2,]	10	12	14	16
[3,]  x - x	18	20	22	24
	[,1]	[,2]	[,3]	[,4]
[1,]	0	0	0	0
[2,]	0	0	0	0
[3,]	0	0	0	0

L’opérateur de multiplication matricielle est %*%

x %*% t(x)

	[,1]	[,2]	[,3]
[1,]	30	70	110
[2,]	70	174	278
[3,]	110	278	446

Pour inverser une matrice on utilise la fonction solve() :

x <- matrix((1:9)^2, nrow = 3) x

[,1]      [,2]      [,3] [1,]      1          16       49 [2,]      4          25        64 [3,]      9          36        81 solve(x)
[,1]                  [,2]	[,3]
[1,]      1.291667         -2.166667	0.9305556
[2,]      -1.166667       1.666667	-0.6111111
[3,]      0.375000         -0.500000  zapsmall(solve(x) %*% x) [,1]  [,2]      [,3] [1,]      1          0         0 [2,]      0          1         0 [3,]      0          0          1	0.1805556

Valeurs propres et vecteurs propres eigen(x)

$values

[1] 112.9839325 -6.2879696 0.3040371

$vectors

	[,1]	[,2]	[,3]
[1,]	-0.3993327	-0.8494260	0.7612507
[2,]	-0.5511074	-0.4511993	-0.6195403
[3,]	-0.7326760	0.2736690	0.1914866

zapsmall(solve(eigen(x)$vectors) %*% x %*% eigen(x)$vectors)

	[,1]	[,2]	[,3]
[1,]	112.9839	0.00000	0.00000
[2,]	0.0000	-6.28797	0.00000
[3,]	0.0000	0.00000	0.30404

Création de Facteurs (pour variables qualitatives et quantitatives

douleur <- c("rien", "fort", "moyen", "moyen", "leger") douleur

[1] "rien"  "fort"  "moyen" "moyen" "leger"

La fonction levels sert à extraire les niveaux possibles d’un facteur : ff <- factor(c(2, 4), levels=2:5) ff

[1] 2 4 Levels: 2 3 4 5 levels(ff)

[1] "2" "3" "4" "5"

Variables qualitatives non ordonnées

fdouleur <- factor(douleur) fdouleur

[1] rien  fort  moyen moyen leger

Levels: fort leger moyen rien 

Variables qualitatives ordonnées (classement alphabéiques)

fdouleur <- factor(douleur, ordered = TRUE) fdouleur

[1] rien  fort  moyen moyen leger

Levels: fort < leger < moyen < rien

table(fdouleur)

fdouleur
fort	leger	moyen	rien
1	1	2	1

Des fonctions plus complexes

GESTION ET IMPORTATION DES DONNEES

read.table("", header=TRUE)[ , ] # header permet d’affecter la première ligne du fichier comme titre des variables

# 3 formes équivalentes d’affectation du tableau dans un objet t3var <- read.table("", h=T)  t3var = read.table("", h=T) 

read.table("", h=T) -> t3var

La fonction read.table possède de multiples options

Lecture dans un fichier

rm(list=ls())

t3var <- read.table("", h=T) t3var

Créer une nouvelle variable dans un tableau de données

IMC<-t3var$poi*10000/t3var$tai^2

IMC

IMC<-round(IMC, digits = 2) summary(IMC)

L’ajouter dans le tableau de données

t3var<-data.frame(t3var,IMC)

Renommer une variable names(t3var)[4]<-"BMI" t3var

Utilisation de attach(base) detach(base)

Comparer les trois programmes suivants

1-      rm(list=ls()) t3var <- read.table("", h=T) attach(t3var)

IMC<-poi*10000/tai^2

IMC

IMC<-round(IMC, digits = 2) summary(IMC) t3var<-data.frame(t3var,IMC) detach(t3var)

t3var

2-      rm(list=ls())

t3var <- read.table("", h=T)

IMC<-poi*10000/tai^2

IMC

IMC<-round(IMC, digits = 2) summary(IMC) t3var<-data.frame(t3var,IMC)

t3var

3-      rm(list=ls()) t3var <- read.table("", h=T) t3var

IMC<-t3var$poi*10000/t3var$tai^2

IMC

IMC<-round(IMC, digits = 2) summary(IMC) t3var<-data.frame(t3var,IMC)

t3var

Modifier une variable t3var$tai <- t3var$tai/100 t3var

Supprimer les données manquantes de tout le tableau (t3var) t3<(t3var)

Calculer un paramètre statistique pour chaque niveau d’une variable qualitative

Charger les variables quantitatives

taille<-c(2,4)

summary(t3var[,taille])

Choisir la variable qualitative de tri

aggregate(t3var[,taille], list(sexe=t3var$sexe), mean)

Découper la fenêtre graphique (pour avoir plusieurs graphiques visibles) par(mfcol=c(2,3))

Découper une variable quantitative en classes en choisissant les bornes (pour Khi-deux par exemple) histc1<-c(0,55,75,200)

x<-cut(t3var$poi, breaks=histc1, include.lowest = TRUE) tb<- table(x, t3var$sexe, dnn=c("poids","sexe")) tb (tb)

summary(tb)

Sauvegarde de fichier de données obtenus lors d’un script : exemple "t3var" write.table(t3var,file="t3var3, sep="\t",dec=",",row.names=F, col.names=T)

Copier-coller des données vers Excel ou Open-Office

Ouvrir Excel et copier le presse papier dans une feuille de calcul

Lire un fichier MS-Excel library(xlsReadWrite)

(file.choose(),colNames= TRUE, sheet= 1, from= 1, rowNames=T)

Rediriger les sorties textes vers un fichier dans un répertoire précis (plus sur la console)

setwd("c:/temp")

sink("", append=TRUE)

 à la fin du programme (tout ce suit entre les deux sink sera redirigé) sink()

LOGICIEL R 

VERSION 2.6.0.

STATISTIQUES DE BASE

ET DEMARCHE

EXPERIMENTALE

Les distributions disponibles

Opérations sur la loi Normale qnorm(0.5) [1] 0 qnorm(0.1) [1] -1.281552 qnorm(0.2) [1] -0.8416212 pnorm(8,mean=10,sd=2) [1] 0.1586553 pnorm(c(2,3),mean=2) [1] 0.5000000 0.8413447 Loi de Fisher

qf(p=0.95, df1=1, df2=2) [1] 18.51282 pf(18.51, df1=1, df2=2) [1] 0.949993 Loi de Student qt(0.05, df= 330) [1] -1.649484 qt(0.05, df= 3) [1] -2.353363 pt(1, df=330) [1] 0.8409784 pt(0, df=330) [1] 0.5

[1] 0.0004792445

Loi du Khi-2 qchisq(0.95, df=1) [1] 3.841459 pchisq(3.84, df=1) [1] 0.9499565 Loi Binomiale pbinom(1, size=1000, prob=0.01)

pbinom(5, size=1000, prob=0.01) [1] 0.06613951 pbinom(10, size=1000, prob=0.01) [1] 0.5830408 pbinom(11, size=1000, prob=0.01) [1] 0.6973501 pbinom(15, size=1000, prob=0.01) [1] 0.9521294 qbinom(0.95, size=1000, prob=0.01) [1] 15

Graphiques sous R

(Adapté de J. R. Lobry, )

On distingue : 

Les fonctions graphiques de bas niveau

Ces fonctions permettent de retoucher un graphique déjà existant (e.g. points(), abline(), arrows(), lines(), seqments(), polygon(), rect(), box(), axis(), title(), rug(), grid(), legend(), text(), mtext()).

Les fonctions graphiques de bas niveau

Ce sont celles que l’on utilise le plus souvent parce qu’elles donnent un graphique complet. Elles sont très nombreuses (e.g. plot(), hist(), dotchart(), stripchart(), pie(), barplot(), boxplot(), curve(), sunflowerplot(), symbols(), pairs(), stars(), assocplot(), mosaicplot(), coplot(), contour(), image(), persp())

Histogrammes

Simple

t3var <- read.table("", h = T)

hist(t3var$tai ,col = grey(0.9), border = grey(0.2), main = paste("Taille de", nrow(t3var),

"étudiants"), xlab = "Taille [cm]", ylab = "Effectifs", labels = TRUE, las = 1, ylim = c(0, 15))

Avec ajustement à la loi Normale t3var <- read.table("", h = T)

hist(t3var$tai , col = grey(0.9), border = grey(0.2), main = paste("Taille de", nrow(t3var),

"étudiants"), xlab = "Taille [cm]", proba = TRUE)

x <- seq(from = min(t3var$tai, = T), to = max(t3var$tai, = T), length = 100) lines(x, dnorm(x, mean(t3var$tai, = TRUE), sd(t3var$tai,  = TRUE))) mtext("Ajustement (moyen) a une loi normale")

Avec ajustement à la densté de probabilité réelle 

t3var <- read.table("", h = T) dst <- density(t3var$tai, = TRUE) hist(t3var$tai, col = grey(0.9), border = grey(0.8), main = paste("Taille de", nrow(t3var), "étudiants"), xlab = "Taille [cm]", proba = TRUE, ylim = c(0, max(dst$y))) lines(dst$x, dst$y)

t3var <- read.table("", h = T) ng <- sum(t3var$sex== "h", = TRUE) nf <- sum(t3var$sex== "f", = TRUE) n <- ng + nf

dst <- density(t3var$tai, = TRUE)

dstg <- density(t3var$tai[t3var$sex == "h"], = TRUE) dstf <- density(t3var$tai[t3var$sex == "f"], = TRUE)

hist(t3var$tai, col = grey(0.9), border = grey(0.8),main = paste("Taille de", nrow(t3var),

"étudiants"),xlab = "Taille [cm]", proba = TRUE, ylim = c(0, max(dst$y))) lines(dstg$x, ng/n * dstg$y, lwd = 3, col = "darkblue") lines(dstf$x, nf/n * dstf$y, lwd = 3, lty = 2, col = "darkred")

lines(dst$x, dst$y)

legend("topright", inset = 0.01, legend = c("Filles","Garcons", "Total"), col = c("darkred", "darkblue","black"), lty = c(2, 1, 1), lwd = 2, = 2)

Relation de type X Y

t3var <- read.table("", h = T)

plot(t3var$tai, t3var$poi, pch = ifelse(t3var$sex =="h", 1, 19)) legend("topleft", inset = 0.01, c("Homme", "Femme"), pch = c(1, 19))

Boîtes à Moustaches

t3var <- read.table("", h = T)

boxplot(t3var$tai ~ t3var$sex, col = c("lightpink", "lightblue"), main = paste("Comparaison taille Homme-Femme"), ylab = "Taille", las = 1)

Diagramme en secteur ou Camembert

t3var <- read.table("", h = T)

pie(table(t3var$sex), col = c("orangered", "yellow"), main = "Répartition Homme-Femme")

Les diagrammes en secteur sont la plupart du temps à éviter. On peut s’en convaincre facilement à l’aide de l’exemple suivant :

(1071966)

data <- rep(10, 10) + rep(2 * runif(5), rep(2, 5)) + rep(c(-2, 2), 5) data <- 100 * data/sum(data) names(data) <- letters[1:10] par(mfrow = c(1, 2))

pie(data, main = "Diagramme en secteur")

dotchart(data, xlim = c(0, 14), pch = 19, main = "Graphe de Cleveland")

par(mfrow = c(1, 1))

Graphique de haut niveau particulier

r=rnorm(100,1) z=hist(r,plot=F)

plot(z)

w=density(r)

plot(w)

plot(sin,xlim=c(0,15))

Les fonctions de format de sortie des graphiques

Il existe de nombreuses fonctions pour ouvrir un nouveau dispositif graphique (e.g. pdf(), jpeg(), postscipt(), x11(), png(), gnome(), quartz(), xfig(), bitmap(), pictex()).

Elles ne sont pas toutes disponibles pour tous les systèmes d’exploitation.

getOption("device").

pdf("")

plot(0)

()

Complément pour les Histogrammes

Les données quantitatives

hist(t3var$poi)

Choix du nombre de classes

hist(t3var$poi,nclass=5)

Bornes et les effectifs des classes

hist(t3var$poi, nclass=5, plot=F)

Bornes choisies  histc<-c(35,55,75,95)

hist(t3var$poi, breaks=histc)

Effectif dans les classes choisies histc<-c(35,55,75,95) hist(t3var$poi, breaks=histc, plot=F)

#Attention hist(t3var$sexe) barplot(summary(t3var$sexe))

barplot(summary(as.factor(t3var$sexe)))

Boite simple

boxplot(t3var$poi)

Une boite par niveau de facteur

boxplot(split(t3var$poi,t3var$sexe))

Pour rendre le graphique plus lisible, ajouter un titre title("poids", sub="Données YLR", ylab= "Effectifs", xlab= "sexe")

Représenter les individus en dehors de la boite boxplot(split(t3var$poi,t3var$sexe)) rug(t3var$poi, side=1)

Fonctions graphiques disponibles

mal <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112) fem <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111) mal fem

summary(mal, fem)# synthèse des données

Min. 1st Qu.  Median    Mean 3rd Qu.      Max. 

107.0 111.2   113.5         113.4     115.5   120.0

tot <- c(mal, fem)

tot

par(mfrow = c(1, 2))

hist(tot, nclass = 8) qqnorm(tot)

qqline(tot)

(tot)

Shapiro-Wilk normality test

data:  tot 

W = 0.9593, p-value = 0.5298

(fem)

Shapiro-Wilk normality test

data:  fem 

W = 0.8809, p-value = 0.1335

(mal)

     Shapiro-Wilk normality test

data:  mal 

W = 0.9937, p-value = 0.9995

plan <- as.factor(rep(c("mal", "fem"), c(10, 10)))

(tot, plan)

(mal, fem)

  F test to compare two variances

data:  mal and fem 

F = 2.681, num df = 9, denom df = 9, p-value = 0.1579 alternative hypothesis: true ratio of variances is not equal to 1  95 percent confidence interval:   0.665931 10.793829  sample estimates: ratio of variances 

          2.681034

y1 <- c(124, 42, 25, 45, 412, 51, 1112, 46, 103, 876, 146, 340,

396)

y2 <- c(1235, 24, 1581, 1166, 40, 727, 3808, 791, 1804, 3460, 719) par(mfrow = c(1, 2)) qqnorm(y1) qqline(y1) qqnorm(y2)

qqline(y2)

(y1) (y2)

Shapiro-Wilk normality test

data:  y1 

W = 0.7547, p-value = 0.002075

Shapiro-Wilk normality test

data:  y2 

W = 0.8686, p-value = 0.07431

Le TEST T

TMF <- t.test(mal, fem, var.equal = TRUE)

TMF

  Two Sample t-test

data:  mal and fem  t = 3.4843, df = 18, p-value = 0.002647

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:  1.905773 7.694227  sample estimates: mean of x mean of y 

    113.4     108.6

ANOVAMF <- anova(lm(tot ~ plan))

ANOVAMF

Analysis of Variance Table

Response: tot

          Df  Sum Sq Mean Sq F value   Pr(>F)    plan       1 115.200 115.200  12.140 0.002647 **

Residuals 18 170.800   9.489                    

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test de student avec test préalable de l’égalité des variances mal <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112) fem <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111) mal fem

(fem, mal)# test de la variance

F test to compare two variances

data:  fem and mal 

F = 0.373, num df = 9, denom df = 9, p-value = 0.1579 alternative hypothesis: true ratio of variances is not equal to 1  95 percent confidence interval:  0.09264553 1.50165699  sample estimates: ratio of variances 

         0.3729904

TMF <- t.test(mal, fem, paired=FALSE, var.equal = TRUE) # test bilateral, variance égale TMF 

Two Sample t-test

data:  mal and fem 

t = 3.4843, df = 18, p-value = 0.002647

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:  1.905773 7.694227  sample estimates: mean of x mean of y 

    113.4     108.6

TMF <- t.test(mal, fem, alternative = “less”, paired=FALSE, var.equal = TRUE) # test unilateral gauche

TMF 

Two Sample t-test

data:  mal and fem  t = 3.4843, df = 18, p-value = 0.9987

alternative hypothesis: true difference in means is less than 0  95 percent confidence interval:

     -Inf 7.188844  sample estimates: mean of x mean of y 

    113.4     108.6

TMF <- t.test(mal, fem, “greater”, paired=FALSE, var.equal = TRUE) # test unilateral droit TMF 

Two Sample t-test

data:  mal and fem 

t = 3.4843, df = 18, p-value = 0.001324

alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval:

 2.411156      Inf  sample estimates: mean of x mean of y 

    113.4     108.6

Khi2

Khi2 d’ajustement à une loi de probabilité théorique

Test direct

(c(1790, 547, 548, 213), p = c(9/16, 3/16, 3/16, 1/16))

Chi-squared test for given probabilities

data:  c(1790, 547, 548, 213) 

X-squared = 7.0628, df = 3, p-value = 0.06992

Autre méthode de saisie

effobs <- c(315, 101, 108, 32) probtheo <- c(9/16, 3/16, 3/16, 1/16)

(effobs, p= probtheo)

    Chi-squared test for given probabilities

data:  effobs 

X-squared = 0.47, df = 3, p-value = 0.9254

Normalité d’une distribution

# chargement des deux variables eff : effecti et val : class eff <- c(1, 1, 4, 2, 5, 2, 7, 6, 11, 17, 13, 25, 38, 72, 102, 140,

107, 77, 52, 34, 21, 4, 6, 4, 2, 0, 0, 1, 0, 1, 0, 0, 1)

val <- seq(275, 595, by = 10)

x <- rep(val, eff)

# paramètres de la série m0 <- mean(x)

sd0 <- sd(x)

# séparation de l’écran graphique en deux parties par(mfrow = c(1, 2))

plot(val, eff/sum(eff)/10, type = "s")

lines(x0 <- seq(275, 595, le = 50), dnorm(x0, m0, sd0), lwd = 2) qqnorm(x)

qqline(x)

(x)

Shapiro-Wilk normality test data: x

W = 0.9344, p-value < 2.2e-16

Le test du Khi2 d’indépendance ou de Contingence

Khi2 : Tableau de contingence 2*2

(matrix(c(49, 7, 11, 18), nc = 2))

Pearson's Chi-squared test with Yates' continuity correction

data:  matrix(c(49, 7, 11, 18), nc = 2) 

X-squared = 20.2872, df = 1, p-value = 6.665e-06

noncampeur <- c(4, 21, 37, 78, 62, 25, 15, 15, 24) campeur <- c(7, 43, 37, 99, 79, 19, 15, 9, 2) client <- matrix(c(campeur, noncampeur), nrow = 9)

nomligne = c("16-19", "20-24", "25-29", "30-39", "40-49", "50-54",

+ "55-59", "60-65", "+65") nomcol = c("camping", "non_camping") rownames(client) <- nomligne colnames(client) <- nomcol

client

camping	non_camping
16-19	7	4
20-24	43	21
25-29	37	37
30-39	99	78
40-49	79	62
50-54	19	25
55-59	15	15
60-65	9	15
+65	2	24

(client)

Pearson's Chi-squared test

data:  client 

X-squared = 32.5107, df = 8, p-value = 7.543e-05

(client[-9, ]) # sans la dernière ligne

        Pearson's Chi-squared test

data:  client[-9, ]

X-squared = 10.7239, df = 7, p-value = 0.1511

(client[c(-8,-9), ])

data:  client[c(-8, -9), ] 

X-squared = 7.8213, df = 6, p-value = 0.2515

client[c(-1,-3,-4,-5,-6,-7,-9), ]

                        camping         non_camping

20-24               43          21

60-65               9            15

(client[c(-1,-3,-4,-5,-6,-7,-9), ])

        Pearson's Chi-squared test with Yates' continuity correction

data:  client[c(-1, -3, -4, -5, -6, -7, -9), ] 

X-squared = 5.1949, df = 1, p-value = 0.02265

Autres exemples

#saisie des données

ki2 =matrix(c(10,20,30,20,30,40,40,50,60),byrow=T,nrow=3) colnames(ki2)<-c("enfant","ado","adulte" ) rownames(ki2)<-c("oui","peut-être", "non") ki2

enfant  ado       adulte
oui	10         20      30
peut-être	20         30      40
non	40         50      60

(ki2)

        Pearson's Chi-squared test

data:  ki2 

X-squared = 2.9304, df = 4, p-value = 0.5695

#saisie des données, à comparer avec ou sans byrow ki2 =matrix(c(10,20,30,20,30,40,40,50,60),nrow=3) colnames(ki2)<-c("enfant","ado","adulte" ) rownames(ki2)<-c("oui","peut-être", "non") ki2

enfant  ado       adulte
oui	10         20      40
peut-être	20         30      50
non	30         40      60

 (ki2)

        Pearson's Chi-squared test

data:  ki2 

X-squared = 2.9304, df = 4, p-value = 0.5695

 #version rapide

 (matrix(c(10,20,30,20,30,40,40,50,60), nc = 3))

  Pearson's Chi-squared test

data:  matrix(c(10, 20, 30, 20, 30, 40, 40, 50, 60), nc = 3) X-squared = 2.9304, df = 4, p-value = 0.5695

Exemple avec découpage en classe

#Découper une variable quantitative en classes en choisissant les bornes (pour Khi-deux par exemple)

unicoolait <- read.table("", header=TRUE) () histc1<-c(0,1,2,3) histc2<-c(0,50,100,150,200,250,305,500,900) max(unicoolait$Jours)

instabilité<-cut(unicoolait$Inst4, breaks=histc1, include.lowest = TRUE) stade<-cut(unicoolait$Jours, breaks=histc2, include.lowest = TRUE)

#regarder le format de x1 et x2 ()

x<-table(stade,instabilité) x

          instabilité stade           [0,1]  (1,2]  (2,3] [0,50]              52     11     10

  (50,100]      43     10      3

  (100,150]     53     22      3

  (150,200]    100     35     18

  (200,250]    107     34     46

  (250,305]    124     48     87

  (305,500]     16     12     32

  (500,900]      9          5      11

(x)

        Pearson's Chi-squared test

data:  x 

X-squared = 101.3307, df = 14, p-value = 2.635e-15

Warning message:

Chi-squared approximation may be incorrect in: (x) 

#calcul du khi-2 théorique ddl<-(nrow(x)-1)*(ncol(x)-1)

ddl

[1] 14

qchisq(0.95, df=ddl)

[1] 23.68479

# calcul des distances de Khi-2 par cases library(pgirmess)

round(valchisq(x),digit=2)

instabilité
stade	[0,1]  (1,2] (2,3]
[0,50]	2.78   0.85 3.02
(50,100]	4.05   0.11 7.88

  (100,150]   1.79   2.73  12.87

  (150,200]   2.09   0.70   9.05

  (200,250]   0.01   0.27   0.08

  (250,305]   3.46   0.23  11.04

  (305,500]   9.48   0.00  22.55   (500,900]   1.87   0.00   4.43 Warning message:

In valat(matr) : 4.2 percent of cells <5

#on enlève des colonnes ou des lignes pour connaître l'effet de ces lignes ou colonnes x[c(-5,-6,-7,-8),]

(x[c(-5,-6,-7,-8),])

          instabilité

stade                [0,1]  (1,2]  (2,3]   [0,50]        52 11     10

  (50,100]      43     10      3

  (100,150]     53     22      3

  (150,200]    100     35     18

        Pearson's Chi-squared test

data:  x[c(-5, -6, -7, -8), ] 

X-squared = 10.186, df = 6, p-value = 0.1170

x[c(-6,-7,-8),]

(x[c(-6,-7,-8),])

          instabilité

stade                [0,1]  (1,2]  (2,3]   [0,50]        52 11     10

  (50,100]      43     10      3

  (100,150]     53     22      3

  (150,200]    100     35     18

  (200,250]    107     34     46

  Pearson's Chi-squared test

data:  x[c(-6, -7, -8), ] 

X-squared = 31.007, df = 8, p-value = 0.0001401

x[,-3]

(x[,-3])

          instabilité

stade                 [0,1] (1,2]

  [0,50]              52 11

  (50,100]      43     10

  (100,150]     53     22

  (150,200]    100     35

  (200,250]    107     34

  (250,305]    124     48

  (305,500]     16     12

  (500,900]      9          5

    Pearson's Chi-squared test

data:  x[, -3] 

X-squared = 9.6327, df = 7, p-value = 0.2104

Warning message:

Chi-squared approximation may be incorrect in: (x[, -3])

Comparaison entre khi-deux contingence à 2*2 et test de comparaison de fréquence

test=matrix(c(10, 20, 10, 40), ncol=2) test

(test)

  Pearson's Chi-squared test with Yates' continuity correction

data:  test 

X-squared = 1.1378, df = 1, p-value = 0.2861

(test)

2-sample test for equality of proportions with continuity correction

data:  test 

X-squared = 1.1378, df = 1, p-value = 0.2861 alternative hypothesis: two.sided  95 percent confidence interval:  -0.1161579  0.4494912  sample estimates:    prop 1    prop 2 

0.5000000 0.3333333

TEST de Fisher pour les faibles effectifs

#test exact de fisher, pour les effectifs un  trop faibles

ki2 =matrix(c(2,25,3,4), nrow=2) colnames(ki2)<-c("enfant","adulte" ) rownames(ki2)<-c("oui", "non")

ki2

(ki2)

        Fisher's Exact Test for Count Data

data:  ki2  p-value = 0.04762

alternative hypothesis: true odds ratio is not equal to 1  95 percent confidence interval:  0.007512789 1.351647922  sample estimates:

odds ratio 

 0.1178116

(ki2)

   Pearson's Chi-squared test with Yates' continuity correction

data:  ki2 

X-squared = 3.1016, df = 1, p-value = 0.07821

Warning message:

Chi-squared approximation may be incorrect in: (ki2) 

Test de Mac Nemar et test de Cochran 

()

Ce test permet de comparer une proportion ou une distribution sur des échantillons appariés.

Mac Nemar : deux échantillons appariés, une variable nominale dichotomique ou un échantillon et deux variables nominales dichotomiques non indépendantes

Cochran : plusieurs échantillons appariés, une variable nominale dichotomique ou un échantillon et plusieurs variables nominales dichotomiques non indépendantes

Exemple : Un directeur d’école souhaite sensibiliser les élèves de CM2 à l’équilibre de l’alimentation. Il réalise une enquête dont il ressort qu’à la cantine, 18 élèves sur 52 (35% des élèves) choisissent leur menu en fonction de l’équilibre alimentaire et pas du seul critère de goût. Pendant deux semaines, une information est proposée dans les classes et pendant la demi-pension. Une seconde enquête est réalisée qui montre que seulement 30 enfants (58% des élèves) choisissent maintenant leur menu en tenant compte de l’équilibre alimentaire. Que penser de l’efficacité de la campagne d’information menée auprès des élèves

L’information donnée ci-dessus n’est pas suffisante  et ne permet que remplir les marges, il nous faut connaître les enfants qui ont changé d’avis (dans les deux sens)

On sait par ailleurs que 19 élèves qui ne prenaient pas en compte l’équilibre alimentaire pour composer leur repas ont changé et 7 ont inversé leur comportement, on complète le tableau dont le nombre de degré de liberté n’est que de 1

Sur les 26 enfants ayant changé, on en a 19 de goût vers équilibre et 7 de équilibre vers goût

# test de McNemar (non indépendance des individus)

nemar =matrix(c(11,7,19,15), nrow=2) colnames(nemar)<-c("avant-équilibre","avant-goût" ) rownames(nemar)<-c("après-équilibre", "après-goût") nemar

                                    avant-équilibre            avant-goût après-équilibre               11                                19 après-goût                          7                                  15

#option False car il n’y a pas d’effectif théorique <5

(nemar, correct=FALSE)

McNemar's Chi-squared test

data:  nemar 

McNemar's chi-squared = 5.5385, df = 1, p-value = 0.01860

Test de Cochran 

C’est le même principe que le test de MacNemar, on travaille dans ce cas sur l’étude d’un pourcentage sur plusieurs échantillons appariés (la programmation R reste la même)

rm(list=ls())

nemar1 =matrix(c(22,3,7,19,3,12,0,1,11), nrow=3)

colnames(nemar1)<-c("avant-Sarko","avant-Royal", "avant-Besancenot") rownames(nemar1)<-c("après-Sarko", "après-Royal","après-Besancenot") nemar1

                                      avant-Sarko avant-Royal  avant-Besancenot

après-Sarko	22	19                 0
après-Royal	3	3                         1
après-Besancenot	7	12                     11

(nemar1)

   McNemar's Chi-squared test

data:  nemar1 

McNemar's chi-squared = 27.9441, df = 3, p-value = 3.732e-06

CORRELATION et REGRESSION

taimdc <- read.table("",header = TRUE)

summary(taimdc)

taimdc <- 2.54 * taimdc #passage en cm

summary(taimdc)

Représentation graphique appelée”sunflowerplot”. Nous avons d’une part les points correspondant au croisement des tailles ”parents, enfants”et les branches décomptent le nombre de répétitions.

sunflowerplot(taimdc, xlab = "Taille du mid-parent [cm]", ylab = "Taille de l'enfant [cm]",las = 1, main = "Les donnees de F. Galton (1886)")

On peut utiliser également une représentation avec une carte de densité :

library(MASS)

densite <- kde2d(taimdc$x, taimdc$y, n = 50)

filled.contour(densite, color = topo.colors, xlab = "Taille du mid-parent [cm]", ylab = "Taille de l'enfant [cm]", las = 1, main = "Les donnees de F. Galton (1886)")

Exemple à peu près complet…

olympic <- read.table("")

names(olympic) <- c("d100", "long", "poid", "haut", "d400", "d110",

"disq", "perc", "jave", "d1500", "score")

names(olympic)

cor(olympic$d100, olympic$d400)

[1] 0.6977326

#test de la corrélation

(olympic$d100, olympic$d400)

Pearson's product-moment correlation

data:  olympic$d100 and olympic$d400  t = 5.5098, df = 32, p-value = 4.507e-06 alternative hypothesis: true correlation is not equal to 0  95 percent confidence interval:  0.4706060 0.8381398  sample estimates:

      cor 

0.6977326

#test de la régression des coefficients a et b de y=aX+b lm.olympic<-lm(olympic$d100 ~ olympic$d400)  summary(lm.olympic)

Call:

lm(formula = olympic$d100 ~ olympic$d400)

Residuals:

     Min       1Q   Median       3Q      Max 

-0.55133 -0.08515  0.01149  0.12298  0.39289 

Coefficients:
Estimate          Std. Error        t value	Pr(>|t|)
(Intercept)    2.80741             1.52791           1.837	0.0754 .
olympic$d400  0.17048          0.03094          5.510   --- Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1	4.51e-06 ***

Residual standard error: 0.209 on 32 degrees of freedom

Multiple R-Squared: 0.4868,     Adjusted R-squared: 0.4708 F-statistic: 30.36 on 1 and 32 DF,  p-value: 4.507e-06

# matrice de corrélation 

cor(olympic)

d100       long       poid       haut       d400       d110 d100   1.0000000 -0.6905079 -0.4201671 -0.3637426  0.6977326  0.7512681 long  -0.6905079 1.0000000  0.3905859  0.4713932 -0.6355823 -0.6544919 poid  -0.4201671 0.3905859  1.0000000  0.3205137 -0.1422790 -0.4888436 haut  -0.3637426 0.4713932  0.3205137  1.0000000 -0.2754966 -0.4868472 d400   0.6977326 -0.6355823 -0.1422790 -0.2754966  1.0000000  0.6547162 d110   0.7512681 -0.6544919 -0.4888436 -0.4868472  0.6547162  1.0000000 disq  -0.3526016 0.3748494  0.8561846  0.3760049 -0.1543935 -0.4034591 perc  -0.6272388 0.6322881  0.6425315  0.4715259 -0.5211362 -0.7087536 jave  -0.3441302 0.4462843  0.7026397  0.3377261 -0.1498775 -0.3504906 d1500  0.2537741 -0.3559770  0.2020151 -0.1315514  0.5544555  0.1553157 score -0.7627374 0.8019923  0.7163215  0.6340494 -0.6527022 -0.8049267             disq perc        jave       d1500      score d100  -0.3526016 -0.6272388 -0.34413018  0.25377410 -0.7627374 long   0.3748494  0.6322881  0.44628427 -0.35597703  0.8019923 poid   0.8561846  0.6425315  0.70263971  0.20201511 0.7163215 haut   0.3760049  0.4715259  0.33772613 -0.13155143  0.6340494 d400 -0.1543935 -0.5211362 -0.14987754  0.55445554 -0.6527022 d110  -0.4034591 -0.7087536 -0.35049057  0.15531567 -0.8049267 disq   1.0000000  0.6195217 0.61782022  0.28819060  0.6746586 perc   0.6195217  1.0000000  0.55713064 -0.07024619  0.8702426 jave   0.6178202  0.5571306  1.00000000  0.04497636 0.6746978 d1500  0.2881906 -0.0702462  0.04497636  1.00000000 -0.2553643 score 0.6746586  0.8702426  0.67469777 -0.25536426  1.0000000

plot(olympic)

#encore mieux vous avez tout library(pgirmess)

cormat(olympic, method = "pearson", sep = TRUE)

$method

[1] "Pearson's product-moment correlation"

$coef.estimates

d100  long  poid  haut  d400  d110  disq  perc  jave d1500 score d100   1.00 -0.69 -0.42 -0.36  0.70  0.75 -0.35 -0.63 -0.34  0.25 -0.76 long  -0.69  1.00 0.39  0.47 -0.64 -0.65  0.37  0.63  0.45 -0.36  0.80 poid  -0.42  0.39  1.00 0.32 -0.14 -0.49  0.86  0.64  0.70  0.20  0.72 haut  -0.36  0.47  0.32  1.00 -0.28 -0.49  0.38  0.47  0.34 -0.13  0.63 d400   0.70 -0.64 -0.14 -0.28  1.00 0.65 -0.15 -0.52 -0.15  0.55 -0.65

d110 0.75 -0.65 -0.49 -0.49  0.65  1.00 -0.40 -0.71 -0.35  0.16 -0.80 disq  -0.35 0.37  0.86  0.38 -0.15 -0.40  1.00  0.62  0.62  0.29  0.67 perc  -0.63  0.63 0.64  0.47 -0.52 -0.71  0.62  1.00  0.56 -0.07  0.87 jave  -0.34  0.45  0.70 0.34 -0.15 -0.35  0.62  0.56  1.00  0.04  0.67 d1500  0.25 -0.36  0.20 -0.13 0.55  0.16  0.29 -0.07  0.04  1.00 -0.26 score -0.76  0.80  0.72  0.63 -0.65 -0.80  0.67  0.87  0.67 -0.26  1.00

$p.values

d100  long  poid  haut  d400  d110  disq  perc  jave d1500 score d100   0.000 0.000 0.013 0.034 0.000 0.000 0.041 0.000 0.046 0.148 0.000 long  0.000 0.000 0.022 0.005 0.000 0.000 0.029 0.000 0.008 0.039 0.000 poid  0.013 0.022 0.000 0.065 0.422 0.003 0.000 0.000 0.000 0.252 0.000 haut  0.034 0.005 0.065 0.000 0.115 0.004 0.028 0.005 0.051 0.458 0.000 d400  0.000 0.000 0.422 0.115 0.000 0.000 0.383 0.002 0.398 0.001 0.000 d110  0.000 0.000 0.003 0.004 0.000 0.000 0.018 0.000 0.042 0.380 0.000 disq  0.041 0.029 0.000 0.028 0.383 0.018 0.000 0.000 0.000 0.098 0.000 perc  0.000 0.000 0.000 0.005 0.002 0.000 0.000 0.000 0.001 0.693 0.000 jave  0.046 0.008 0.000 0.051 0.398 0.042 0.000 0.001 0.000 0.801 0.000 d1500 0.148 0.039 0.252 0.458 0.001 0.380 0.098 0.693 0.801 0.000 0.145 score 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.145 0.000

# représentation graphique plot(olympic$d100, olympic$d400) abline(lsfit(olympic$d100, olympic$d400))

hist(residuals(lm.olympic))

Exemple direct

annee <- 75:87

inclinaison <- c(642,644,656,667,673,688, 696,698,713,717,725,742,757) plot(annee,inclinaison) <- lm(inclinaison ~ annee)

cor(annee, inclinaison)

 [1] 0.9939717

Call: lm(formula = inclinaison ~ annee)

Coefficients:

(Intercept)        annee  

    -61.121        9.319  

abline(-61.121,9.319)

summary()

Call: lm(formula = inclinaison ~ annee)

Residuals:

    Min      1Q  Median      3Q     Max 

-5.9670 -3.0989  0.6703  2.3077  7.3956 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)     (Intercept) -61.1209    25.1298  -2.432   0.0333 *  

annee         9.3187     0.3099  30.069  6.5e-12 ***

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 4.181 on 11 degrees of freedom

Multiple R-Squared: 0.988,      Adjusted R-squared: 0.9869 

F-statistic: 904.1 on 1 and 11 DF,  p-value: 6.503e-12

Exemple de Régression multiple

Fichier Poussin poussin<-read.table(file="",header=TRUE)

lm.poussin<-lm(poussin$GMQT6 ~ poussin$GMQT1+ poussin$GMQT2+ poussin$GMQT3+ poussin$GMQT4+ poussin$GMQT5) summary(lm.poussin)

Call:

lm(formula = poussin$GMQT6 ~ poussin$GMQT1 + poussin$GMQT2 +      poussin$GMQT3 + poussin$GMQT4 + poussin$GMQT5)

Residuals:

    Min      1Q  Median      3Q     Max 

-30.704 -10.638  -2.988   7.145 110.110 

Coefficients:
Estimate          Std. Error                  t value	Pr(>|t|)
(Intercept)                   -0.6602            19.6103                       -0.034	0.9733
poussin$GMQT1        -0.1895           0.4179                         -0.453	0.6519
poussin$GMQT2        0.3359             0.4528                         0.742	0.4611
poussin$GMQT3        1.6304             0.3211                         5.078	4.1e-06 ***
poussin$GMQT4        0.3694             0.1668                         2.215	0.0307 *
poussin$GMQT5        -0.3472           0.1258                         -2.760 --- Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  Residual standard error: 21.04 on 59 degrees of freedom Multiple R-Squared: 0.4952,     Adjusted R-squared: 0.4524  F-statistic: 11.57 on 5 and 59 DF,  p-value: 8.209e-08 confint(lm.poussin) 2.5 %                                       97.5 % (Intercept)       -39.90042060                         38.57995129 poussin$GMQT1          -1.02570873                0.64668520 poussin$GMQT2            -0.57011175                1.24192390 poussin$GMQT3                0.98794133                2.27280955 poussin$GMQT4                0.03564542                0.70308689 poussin$GMQT5            -0.59892166                 -0.09544432	0.0077 **

residuals(lm.poussin)

1           2           3           4           5           6 7           8           9          10          11          12          13 

 -3.9569792 1.6390383  -8.4238402   1.7599483   7.7799530  -1.2166056  16.7519952 -18.0174593  -4.2804064   7.0074868 -11.2608538  -0.2252806  21.3250295 

14          15          16          17          18          19 20          21          22          23          24          25          26 

 -5.9104490 15.2306556  -1.4938839   3.5859821   6.4163923   7.2050599  -4.0747102 -19.8629673 -17.5120991 -27.4056233  -9.9629939  14.9806517  -6.7666494 

27          28          29          30          31          32          33 34          35          36          37          38          39 

-11.2315913 -5.6894417 -18.9421643  -4.8384983  13.8816697 110.1095436 -10.6380158 4.1714380  34.7076507 -30.7042500 14.7716424 -20.7655660  -5.4084583 

40          41          42          43          44          45 46          47          48          49          50          51          52 

 31.6333394 -20.3909248  -0.2584892  18.6972813  20.7514943 -10.0703904   7.1446354 -5.6262470  -8.7955474 -27.0668720 10.0197122   8.9769332  -2.9877063 

53          54          55          56          57          58 59          60          61          62          63          64          65 

5.6123232 -24.7084838 -11.1802845  35.9747981  -0.6225396 -25.0702911 -11.4451849  16.2734827  12.2155891   0.5373977  

5.0711217 -4.9764381  -2.8613503

require(car)

vif(lm.poussin)

poussin$GMQT1 poussin$GMQT2 poussin$GMQT3 poussin$GMQT4 poussin$GMQT5 

     2.398600      2.127667      1.468087      1.361920 1.116827

hist(residuals(lm.poussin))

Histogram of residuals(lm.poussin)

La sélection des variables (méthode du meilleur sous-ensemble) par le critère du BIC

par(mfrow=c(1,1)) x<-poussin[,5:9] y<-poussin[,4] require(leaps) leaps(x,y) res<-regsubsets(x,y)

plot(res)

ANALYSE DE VARIANCE

ali <- rep(c("t1", "t2", "t3"), c(5, 5, 5)) ali <- as.factor(ali)

gain <- c(37.7, 44.6, 42.1, 45.1, 43.2, 45.2, 54.2, 38.1, 48.3,

55.1, 71.3, 63, 67, 65, 69)

anova(lm(gain ~ ali))

Analysis of Variance Table

Response: gain

          Df Sum Sq Mean Sq F value Pr(>F) ali        2 126.15 63.07  1.9529 0.1844

Residuals 12 387.58   32.30       

TukeyHSD(aov(gain~ ali))

$ali diff       lwr           upr     p adj

t2-t1  5.64 -3.949229 15.22923 0.2956934 t3-t1  6.56 -3.029229 16.14923 0.2031780 t3-t2  0.92 -8.669229 10.50923 0.9646347

(ali, gain)

  Kruskal-Wallis rank sum test

data:  ali and gain 

Kruskal-Wallis chi-squared = 13.3, df = 13, p-value = 0.4249

library(pgirmess)

kruskalmc(gain, ali, probs=0.05)

Multiple comparison test after Kruskal-Wallis  p.value: 0.05 

Comparisons

      difference t1-t2     4.9 6.771197      FALSE t1-t3     4.7     6.771197      FALSE t2-t3     0.2 6.771197      FALSE

Tableau graphique des moyennes

library(agricolae)

#le test LSD

<- aov(gain ~ ali) df<-df.residual() MSerror<-deviance()/df

compar<- (gain, ali, df, MSerror, alpha=0.01,

group=TRUE)

Study:

LSD t Test for gain 

                            ......

Alpha                     0.010000

Error Degrees of Freedom 12.000000

Error Mean Square        22.662667

Critical Value of t       3.054540

Treatment Means

  ali  gain  replication 1  t1 42.54 1.320076           5

2    t2 48.18 3.119519           5

3    t3 67.06 1.457258           5

Least Significant Difference 9.196684

Means with the same letter are not significantly different.

Groups, Treatments and means

a      t3      67.06  b       t2      48.18 

b     t1      42.54

#le graphique

bar.group(compar, horiz=TRUE,density=8,col="blue", border="red", xlim=c(0,100)) title(=0.8,main="Comparaison entre traitements",xlab="mesure", ylab="traitements")

CALCUL DE LA PUISSANCE

Test T

Calcul de la puissance : évaluation de n

(power = .80, delta = 3, sig.level = 0.05, sd = 3, type=”two.sample“,alternative = "two.sided") 

# option pour type : two.sample, paired, option pour alternative : “one sided ou two.sided

Two-sample t test power calculation 

              n = 16.71477           delta = 3              sd = 3       sig.level = 0.05           power = 0.8

    alternative = two.sided

 NOTE: n is number in *each* group

Calcul de la puissance : évaluation de la puissance (power=NULL,n=15, delta = 3, sig.level = 0.05, sd = 3, type=”two.sample“, alternative = "two.sided")

Two-sample t test power calculation 

              n = 15           delta = 3              sd = 3 sig.level = 0.05           power = 0.752921

    alternative = two.sided

 NOTE: n is number in *each* group

Calcul de la puissance : évaluation de l’écart-type résiduel

(power=0.8, n=15, delta = 3, sd=NULL, sig.level = 0.05, type=”two.sample“, alternative = "two.sided") 

Two-sample t test power calculation 

              n = 15           delta = 3              sd = 2.830725       sig.level = 0.05           power = 0.8

    alternative = two.sided

 NOTE: n is number in *each* group

Calcul de la puissance : évaluation de l’écart attendu avant traitement (power=0.8, n=15, delta=NULL, sd = 3, sig.level = 0.05, type=”two.sample“, alternative = "two.sided") 

Two-sample t test power calculation 

              n = 15           delta = 3.179403              sd = 3       sig.level = 0.05           power = 0.8

    alternative = two.sided

 NOTE: n is number in *each* group

Calcul de Puissance pour ANOVA

(groups=4, n=5, =3, =3)

Balanced one-way analysis of variance power calculation 

         groups = 4               n = 5

    = 3      = 3       sig.level = 0.05

          power = 0.8303491

 NOTE: n is number in each group

(groups=4, =1, =3, power=.80)

Balanced one-way analysis of variance power calculation 

         groups = 4               n = 11.92613     = 1      = 3       sig.level = 0.05

          power = 0.8

 NOTE: n is number in each group

groupmeans <- c(120, 130, 140, 150)      (groups = length(groupmeans), =var(groupmeans),                       =500, power=.90) Balanced one-way analysis of variance power calculation 

         groups = 4

              n = 15.18834     = 166.6667 = 500       sig.level = 0.05

          power = 0.9

 NOTE: n is number in each group

EXEMPLECOMPLET (A PEU PRES ….)

Cas de l’influence d’un facteur à 4 traitements sur le poids de Cerf

# Importation du fichier

cerf<-read.table(“”,header=TRUE) [ , ] cerf     poi    reg

1   60.5  rA

2   62.1  rA

3   57.3  rA

4   55.0  rA

5   64.2  rA

6   61.1  rA

7   60.0  rA

8   59.7  rA

9   60.2  rA

10  59.9  rA

11  72.1  rB

12  70.7  rB

13  72.5  rB

14  68.0  rB

15  67.4  rB

16  72.6  rB

17  67.2  rB

18  68.9  rB

19  74.2  rB

20  71.4  rB 21 62.0  rC

22  60.3  rC

23  57.5  rC

24  61.8  rC

25  62.5  rC

26  61.2  rC

27  64.5  rC

28  56.3  rC

29  63.1  rC

30  60.8  rC

31  40.1  rD

32  36.5  rD

33  39.7  rD

34  42.3  rD

35  45.7  rD

36  41.4  rD

37  40.6  rD

38  39.8  rD

39  42.0  rD

40  42.9  rD

# Test de l’homogénéité des variances

(cerf$poi ~ cerf$reg)

Bartlett test of homogeneity of variances

data:  cerf$poi by cerf$reg 

Bartlett's K-squared = 0.0083, df = 3, p-value = 0.9998

# pour comparer deux variances (test de Fisher, donc ne marche pas ici)  (cerf$poi ~ cerf$reg)

Erreur dans .formula(cerf$poi ~ cerf$reg) : grouping factor must have exactly 2 levels

# Test de la normalité des données

(cerf$poi) # A ne pas faire car test la normalité sur l’ensemble des données et pas par traitement

(cerf$poi[cerf$reg=="rA"])

Shapiro-Wilk normality test

data:  cerf$poi[cerf$reg == "rA"] 

W = 0.9368, p-value = 0.5177

(cerf$poi[cerf$reg=="rB"]) (cerf$poi[cerf$reg=="rC"])

(cerf$poi[cerf$reg=="rD"])

# Analyse de variance classique avec Test de la normalité des résidus # 2 méthodes # Méthode 1

ana=anova(lm(cerf$poi ~ cerf$reg)) ana

sum=summary(ana)

Analysis of Variance Table

Response: cerf$poi

          Df Sum Sq Mean Sq F value    Pr(>F)    

cerf$reg   3 4547.7  1515.9  249.17 < 2.2e-16 ***

Residuals 36  219.0     6.1                      

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Méthode 2

ana=aov(poi ~ reg, data = cerf)

sum=summary(ana) sum

Df         Sum Sq          Mean Sq         F value		Pr(>F)
reg         3          4547.7             1515.9             249.17  Residuals        36        219.0               6.1 --- Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1		< 2.2e-16 ***
resid=residuals(ana) res=c((resid)) res $statistic         W  0.9701124  $p.value [1] 0.3628720 $method [1] "Shapiro-Wilk normality test" $ [1] "resid" # test de comparaison multiples Tukey TukeyHSD(aov(cerf$poi ~ cerf$reg)) Tukey multiple comparisons of means     95% family-wise confidence level Fit: aov(formula = cerf$poi ~ cerf$reg) $`cerf$reg`
diff        lwr        upr	p adj
rB-rA  10.5   7.529164  13.470836	0.0000000
rC-rA   1.0  -1.970836   3.970836	0.8014418
rD-rA -18.9 -21.870836 -15.929164	0.0000000
rC-rB  -9.5 -12.470836  -6.529164	0.0000000
rD-rB -29.4 -32.370836 -26.429164	0.0000000
rD-rC -19.9 -22.870836 -16.929164	0.0000000

# Analyse de variance non-paramétrique (poi ~ reg, data = cerf)

Kruskal-Wallis rank sum test

data:  poi by reg 

Kruskal-Wallis chi-squared = 33.401, df = 3, p-value = 2.651e-07

# Test de comparaison de moyennes multiples non-paramétriques

kruskalmc(cerf$poi, cerf$reg, probs=0.05, cont="one-tailed") # unilatéral

Multiple comparison test after Kruskal-Wallis, treatments vs control (one-tailed)  p.value: 0.05  Comparisons

      difference rA-rB    16.8 7.526056       TRUE rA-rC     3.6     7.526056      FALSE

rA-rD    13.2     7.526056       TRUE

kruskalmc(cerf$poi, cerf$reg, probs=0.05, cont="two-tailed") # bilatéral

Multiple comparison test after Kruskal-Wallis, treatment vs control (two-tailed)  p.value: 0.05  Comparisons

      difference rA-rB    16.8 11.12570       TRUE rA-rC     3.6     11.12570      FALSE

rA-rD    13.2     11.12570       TRUE

Analyse de Variance à deux facteurs

adv2F <- read.table("", h=T) names(adv2F) <- c("fong", "regul", "rep", "parc", "rdt")

adv2F$fong <- as.factor(adv2F$fong) adv2F$regul <- as.factor(adv2F$regul)

summary(adv2F)

lm1 <- lm(rdt ~ fong * regul , data = adv2F)

anova(lm1)

Analysis of Variance Table

Response: rdt
Df	Sum Sq	Mean Sq          F value	Pr(>F)
fong                1	8.333	8.333               1.2821	0.30073
regul               2	88.667	44.333             6.8205	0.02851 *
fong:regul       2	80.667	40.333             6.2051	0.03462 *
Residuals        6	39.000	6.500

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model.tables(aov(rdt ~ fong * regul , data = adv2F), type = "means")

Tables of means

Grand mean

70.66667 

 fong  fong

    1     2 

71.50 69.83 

 regul  regul

   1    2    3 

68.5 69.0 74.5 

 fong:regul      regul fong 1     2       3   

   1 67.5 68.0 79.0    2 69.5 70.0 70.0

resid=residuals(lm1) res=c((resid))

res

$statistic

        W 

0.8793824 

$p.value

[1] 0.08606476

$method

[1] "Shapiro-Wilk normality test"

$

[1] "resid"

hist(resid)

plot(resid)

abline(0,0)

TukeyHSD(aov(adv2F$rdt ~ adv2F$fong * adv2F$regul))

Tukey multiple comparisons of means

    95% family-wise confidence level

Fit: aov(formula = adv2F$rdt ~ adv2F$fong * adv2F$regul)

$àdv2F$fong`

          diff lwr                                upr                  p adj

2-1       -1.666667 -5.268423              1.935090         0.3007299

$àdv2F$regul`

            diff         lwr                            upr       p adj 2-1         0.5       -5.03141039    6.03141           0.9587700

3-1       6.0     �� 0.46858961               11.53141         0.0364438

3-2       5.5       -0.03141039             11.03141         0.0510850

$àdv2F$fong:adv2F$regul`

        diff         lwr       upr                   p adj 2:1-1:1 2.0  -8.1466581 12.146658       0.9605063

1:2-1:1  0.5  -9.6466581 10.646658              0.9999344

2:2-1:1  2.5  -7.6466581 12.646658              0.9088593

1:3-1:1 11.5   1.3533419 21.646658              0.0289198

2:3-1:1  2.5  -7.6466581 12.646658              0.9088593

1:2-2:1 -1.5 -11.6466581 8.646658              0.9881282

2:2-2:1  0.5  -9.6466581 10.646658              0.9999344

1:3-2:1  9.5  -0.6466581 19.646658              0.0655999

2:3-2:1  0.5  -9.6466581 10.646658              0.9999344

2:2-1:2  2.0  -8.1466581 12.146658              0.9605063

1:3-1:2 11.0   0.8533419 21.146658              0.0352816

2:3-1:2  2.0  -8.1466581 12.146658              0.9605063

1:3-2:2  9.0  -1.1466581 19.146658              0.0812384

2:3-2:2  0.0 -10.1466581 10.146658             1.0000000

2:3-1:3 -9.0 -19.1466581 1.146658              0.0812384

library(lattice)

my.m <- matrix(lm1$residuals,nrow=4)

levelplot(my.m, cuts=4)

attach(adv2F)

(fong, regul, rdt,col = 2:20, lwd=2) fong, regul, rdt,col = 2:20, lwd=2)

detach(adv2F)

(adv2F$rdt, adv2F$regul,p.adjust="bonf")

Pairwise comparisons using t tests with pooled SD 

data:  adv2F$rdt and adv2F$regul 

1     2   

2     1.00 -    3 0.15 0.21

P value adjustment method: bonferroni

Quelques compléments

#cas intéressant rm(list=ls())

adv2F <- read.table("", h=T)

names(adv2F) <- c("fong", "regul", "rep", "parc", "rdt")

adv2F$fong <- as.factor(adv2F$fong) adv2F$regul <- as.factor(adv2F$regul)

summary(adv2F)

library(agricolae)

<- aov(adv2F$rdt ~ adv2F$fong * adv2F$regul)

DFerror<-df.residual()

MSerror<-deviance()/DFerror

(adv2F$rdt, adv2F$regul,adv2F$fong,DFerror, MSerror,  alpha = 0.05, group=TRUE, main = NULL)

(adv2F$rdt, adv2F$regul, DFerror, MSerror,  alpha = 0.05, group=TRUE, main =

NULL)

Study:

HSD Test for adv2F$rdt 

                                      ......

Alpha                               0.050000

Error Degrees of Freedom            6.000000

Error Mean Square                   6.500000

Critical Value of Studentized Range 4.339195

Treatment Means

  adv2F.regul   replication

1       1      68.5 1.190238           4

2       2      69.0 1.290994           4 3           3      74.5 2.753785           4

Honestly Significant Difference 5.53141

Means with the same letter are not significantly different.

Groups, Treatments and means

a 3       74.5  ab       2       69  b       1       68.5 

trt means	M	N
1   3	74.5   a	4	2.753785
2   2	69.0  ab	4	1.290994
3   1	68.5   b	4	1.190238

compar<- (adv2F$rdt, adv2F$regul,DFerror, MSerror,  alpha = 0.05, group=TRUE, main = NULL)

#le graphique

bar.group(compar, horiz=FALSE,density=8,col="blue", border="red", ylim=c(0,90))

title(=0.8,main="Comparaison entre traitements",ylab="mesure", xlab="traitements")

Analyse de Variance en BLOC

advbloc <- read.table("", h=T) names(advbloc) <- c("var", "bloc", "rdt") advbloc$var <- as.factor(advbloc$var) advbloc$bloc <- as.factor(advbloc$bloc)

summary(advbloc)

var   bloc       rdt       

 1:4   1:4   Min.   :59.00  

 2:4   2:4   1st Qu.:63.75  

 3:4   3:4   Median :65.50  

 4:4   4:4   Mean   :66.31  

             3rd Qu.:69.25  

             Max.   :75.00

lm1 <- lm(rdt ~ var+bloc , data = advbloc)

anova(lm1)

Analysis of Variance Table

Response: rdt
Df	Sum Sq	Mean Sq          F value	Pr(>F)
var                   3	122.688	40.896             90.60	4.823e-07 ***
bloc                 3	138.687	46.229             102.42	2.829e-07 ***
Residuals        9	4.062	0.451

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(aov(advbloc$rdt ~ advbloc$var + advbloc$bloc))

TukeyHSD(aov(advbloc$rdt ~ advbloc$var + advbloc$bloc))

  Tukey multiple comparisons of means

    95% family-wise confidence level

Fit: aov(formula = advbloc$rdt ~ advbloc$var + advbloc$bloc)

$àdvbloc$var`

      diff        lwr       upr        p adj 2-1   4.00  2.5169174 5.483083    0.0000705

3-1       -3.75 -5.2330826 -2.266917                          0.0001178

4-1       1.00 -0.4830826 2.483083                            0.2224631

3-2       -7.75 -9.2330826 -6.266917                          0.0000003

4-2       -3.00 -4.4830826 -1.516917                          0.0006502

4-3       4.75  3.2669174 6.233083                            0.0000174

$àdvbloc$bloc`

     diff       lwr      upr         p adj 2-1  2.50 1.0169174 3.983083    0.0023825

3-1       4.75 3.2669174 6.233083                              0.0000174

4-1       8.00 6.5169174 9.483083                              0.0000002

3-2       2.25 0.7669174 3.733083                              0.0048056

4-2       5.50 4.0169174 6.983083                              0.0000051

4-3       3.25 1.7669174 4.733083                              0.0003571

LA MÊME, EN RANDOMISATION TOTALE

advbloc <- read.table("", h=T) names(advbloc) <- c("var", "bloc", "rdt") advbloc$var <- as.factor(advbloc$var) advbloc$bloc <- as.factor(advbloc$bloc)

summary(advbloc)

var   bloc       rdt       

 1:4   1:4   Min.   :59.00  

 2:4   2:4   1st Qu.:63.75  

 3:4   3:4   Median :65.50  

 4:4   4:4   Mean   :66.31  

             3rd Qu.:69.25  

             Max.   :75.00  

lm1 <- lm(rdt ~ var , data = advbloc)

anova(lm1)

Analysis of Variance Table

Response: rdt
Df	Sum Sq	Mean Sq          F value	Pr(>F)
var                   3	122.688	40.896             3.4378	0.052 .
Residuals        12	142.750	11.896

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

TukeyHSD(aov(advbloc$rdt ~ advbloc$var)) # A comparer avec le modèle en bloc

Tukey multiple comparisons of means

    95% family-wise confidence level

Fit: aov(formula = advbloc$rdt ~ advbloc$var)

$àdvbloc$var`

diff        lwr        upr       p adj 2-1   4.00  -3.240660 11.2406603 0.3943038 3-1  -3.75 -10.990660  3.4906603  0.4469881

4-1       1.00  -6.240660  8.2406603   0.9756916 3-2             -7.75 -14.990660 -0.5093397            0.0347946

4-2       -3.00 -10.240660 4.2406603            0.6207763

4-3       4.75  -2.490660 11.9906603              0.2601686

Un peu plus difficile…

Création de fonctions : il est possible de créer une fonction synthétique qui permet de regrouper plusieurs commandes à la suite .

La structure générale d’une fonction est la suivante :

noms <- function( arguments ) { corps de la fonction }

Un exemple: créer une fonction qui calcule le coefficient de variation (écartype/moyenne) sur la série 1 :100

x<- c(1 :100) cv <- function(x)

{ <- sqrt(var(x))

100*(x)

} cv(x)

Génération d’histogramme 

par(mfrow=c(2,2)) for(k in 1:4)

{

hist(rnorm(k*50),main="",xlab="")

 title(paste("n =",k*50))

 }

Encore plus fort….

Programme qui réalise une régression (Fastoche)

# lecture des données :sélectionner le fichier unicoolait unico <- read.table(file.choose(), header=T)

# régression linéaire z<-aov(unico$TB ~ unico$TP)

# graphe XY

plot(unico$TP , unico$TB , ylab="TB",xlab="TP")

# ajout de la droite de régression abline(z)

Créer une fonction qui permettra de répéter les régressions

rm(list=ls())

unico <- read.table(file.choose(), header=T)

grapheXY<-function(X,Y,TAB)

{ attach(TAB) z<-aov(Y ~ X)

titre<-paste("Relation entre ",substitute(X), " et ",substitute(Y)) plot(X, Y, ylab= substitute (Y),xlab= substitute (X),main=titre) abline(z) detach(TAB)

}

Sauvegarder la fonction

sauvegarder la fonction dans un fichier : grapheXY.R (attention à l'extension)

Charger la fonction

pour l’utiliser, lancer le fichier avec le menu « Fichier/ouvrir script »

Utiliser la fonction grapheXY(TB,TP,unico)

grapheXY(logcel,TP,unico)

Un exemple qui utilise presque tout….A vous de jouer et de comprendre

#Analyse de variance

#effacement de toutes les variables en mémoire rm(list=ls())

#importation fichier txt (tabulation) avec la première ligne correspondant aux noms de variable

unicoolait <- read.table("", header=TRUE,check.names=TRUE)

#indique que ces deux variables sont à considérer comme des variables qualitative

unicoolait$STA2006 <- as.factor(unicoolait$STA2006) unicoolait$STAB2007 <- as.factor(unicoolait$STAB2007)

#donne le type des variables class(unicoolait$STA2006)

class(unicoolait$TP)

# description des modalités de la variable instabilité

levels(unicoolait$STA2006)

# est-ce pertinent ? levels(unicoolait$TP)

#synthèse summary(unicoolait)

#analyse de variance paramétrique

lm1 <- lm(TP ~ STA2006 , data = unicoolait) anova(lm1)

#affichage des moyennes par traitement

model.tables(aov(TP ~ STA2006 , data = unicoolait), type = "means")

#analyse des résidus resid=residuals(lm1) res=c((resid))

res plot(resid)

hist(resid)

#visualisation des résidus (est-ce pertinent ?) library(lattice)

my.m <- matrix(lm1$residuals,nrow=50)

levelplot(my.m, cuts=4)

# test de comparaison multiples paramétrique

TukeyHSD(aov(unicoolait$TP ~ unicoolait$STA2006))

#analyse de variance non paramétrique

(TP ~ STA2006, data = unicoolait)

#test de comparaison non paramétrique library(pgirmess)

kruskalmc(unicoolait$TP, unicoolait$STA2006, probs=0.05, cont="one-tailed") # unilatéral kruskalmc(unicoolait$TP, unicoolait$STA2006, probs=0.05, cont="two-tailed") # bilatéral

#regression simple

rm(list=ls())

unicoolait <- read.table("", header=TRUE,check.names=TRUE) unicoolait$STA2006 <- as.factor(unicoolait$STA2006) unicoolait$STAB2007 <- as.factor(unicoolait$STAB2007)

#affichage graphique

plot(unicoolait$TP,unicoolait$TB)

<- lm(unicoolait$TP ~ unicoolait$TB)

#synthèse summary()

#intervalle de confiance

confint()

#graphique

abline(lsfit(unicoolait$TP,unicoolait$TB))

resid1=residuals() res1=c((resid1)) res1 plot(resid1)

hist(resid1)

#regression multiple

rm(list=ls())

unicoolait <- read.table("", header=TRUE,check.names=TRUE) unicoolait$STA2006 <- as.factor(unicoolait$STA2006) unicoolait$STAB2007 <- as.factor(unicoolait$STAB2007)

lm.unicoolait<-lm(unicoolait$Inst4 ~ unicoolait$Lait+ unicoolait$TB + unicoolait$TP + unicoolait$logcel + unicoolait$Nlac + unicoolait$Jours) summary(lm.unicoolait)

confint(lm.unicoolait)

hist(residuals(lm.unicoolait)) resid3=residuals(lm.unicoolait) res1=c((resid3)) res1

# ACP

rm(list=ls())

unicoolait <- read.table("", header=TRUE,check.names=TRUE) unicoolait$STA2006 <- as.factor(unicoolait$STA2006)

unicoolait$STAB2007 <- as.factor(unicoolait$STAB2007)

#sélection des variables dans l'ACP uni<-unicoolait[,c(3,6,8,16,14,13)] uni

#matrice de corrélation

round(cor(uni), digits = 2)

#calcul des composantes principales, valeur propre, vecteur propre prcomp1 <- prcomp(uni, scale = T)

summary(prcomp1)

print(round(var(prcomp1$x), digits = 3))

print(round(t(prcomp1$rotation) %*% prcomp1$rotation, digits = 3))

prcomp1$sdev^2

#valeurs propres

eigen(cor(uni))$values

plot(prcomp1)

#graphique corrélation et individus prin1<-prcomp(uni, scale = T)

biplot(prin1)

#cercle des corrélations (axe 1 et 2) library(ade4) par(mfrow = c(1, 3))

s.corcircle(.frame(prcomp1$rotation[, 1] * prcomp1$sdev[1], prcomp1$rotation[, 2] * prcomp1$sdev[2]))

# Classification Ascendante Hiérarchique en 15 groupes

don<-scale(uni, center = TRUE, scale = TRUE)

dc<-dist(don, method ="euclidean", diag=FALSE, upper=FALSE) hier<-hclust(dc,"ward") plot(hier,hang=-1)

cl15<-cutree(hier,15)

# moyenne par groupe pour matrice uni doncl<-merge(uni, cl15,by= "row.names") moyennes<-aggregate(doncl[,-1], list(doncl$y), mean) moyennes[,-1]<-round(moyennes[,-1],2)

moyennes

#ecartype par groupe

doncl<-merge(uni, cl15,by= "row.names") ecartype<-aggregate(doncl[,-1], list(doncl$y), sd) ecartype[,-1]<-round(ecartype[,-1],2) ecartype

#fusion des tableaux doncluni<-doncl[,-1]

unitotal<-merge(doncluni, unicoolait)

#moyenne par groupe

moyennesbis<-aggregate(unitotal[,c(-18,-19)],list(unitotal$y), mean) round(moyennesbis[,-1],2)

#corrélation entre stade de lactation et instabilité

cor(moyennesbis$Jour,moyennesbis$Inst4)

#matrice de corrélation générale

round(cor(moyennesbis[,-1]),2)

plot(moyennesbis$Jour,moyennesbis$Inst4)

abline(lsfit(moyennesbis$Jour,moyennesbis$Inst4))

#comptage du nombre d'animaux par classe doncl$y<-as.factor(doncl$y)

summary(doncl$y)

# très joli résumé (mais nécessite 3 minutes)

library(cluster)

clusplot(dc,cl15,diss=T,shade=T,color=T,labels=5,main="") abline(v=0,h=0)

#extraction des données de stade de lactation supérieur à 250 jours et instabilité >1 unitotal

soustotal<-(unitotal[unitotal[4]>250 & unitotal[unitotal[8]>1, c(1,2,3,4,5,6,7,8)]) soustotal nrow(soustotal)

ncol(soustotal)

lm.soustotal<-lm(soustotal$Inst4 ~ soustotal$Lait+ soustotal$TB + soustotal$TP + soustotal$logcel + soustotal$Nlac + soustotal$Jours)

summary(lm.soustotal)

hist(residuals(lm.soustotal)) resid3=residuals(lm.soustotal) res1=c((resid3)) res1

#AFD unicoolait sur instabilité 2 classes

() rm(list=ls())

unicoolait <- read.table("", header=TRUE,check.names=TRUE)

uniafd<-unicoolait[,c(3,6,8,16,14,13,18)] uniafd

uniafd$STAB2007<-as.factor(uniafd$STAB2007)

#caractéristiques de toutes les variables en mémoire ()

#découpe l'écran en 4 (lignes)* 3( colonnes) pour les graphiques par(mfcol = c(4, 3))

# création d'une fonction boucle pour création de graphiques  for (k in 1:6) { j0 <- names(uniafd)[k]

br0 <- seq(min(uniafd[, k]), max(uniafd[, k]), le = 12) x0 <- seq(min(uniafd[, k]), max(uniafd[, k]), le = 891) for (i in 1:2) {

i0 <- levels(uniafd$STAB2007)[i] x <- uniafd[uniafd$STAB2007 == i0, j0] hist(x, br = br0, proba = T, col = grey(0.8), main = i0, xlab = j0)

lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)

}

}

#représentation graphique des deux groupes de l'AFD par rapport aux variables étudiées library(ade4)

pan1 <- function(x, y, ...) { xy <- .frame(x, y)

s.class(xy, uniafd$STAB2007, = F, add.p = T, clab = 1.5,

col = c("blue", "black"), cpoi = 2, csta = 0.5)

}

pairs(uniafd[, 1:6], panel = pan1)

summary(xy)

# représentation 3D library(scatterplot3d) par(mfrow = c(1, 1)) mar0 = c(2, 3, 2, 3)

scatterplot3d(uniafd[, 2], uniafd[, 3], uniafd[, 4], mar = mar0, color = c("black",

"red")[uniafd$STAB2007], pch = 19)

Modèle en Mesures Répétées

()

Petit rappel théorique et contextuel

On parle de mesures répétées :

*        que ce soit des mesures contemporaines, récupérées par exemple à un instant donné, en différents emplacements de la peau ou de la carcasse dudit animal,

*        ou qu’elles soient espacées plus ou moins régulièrement, au cours du temps, mais alors dans des conditions en tout point comparables.

En conséquence, les diverses mesures de chaque unité ne sont pas indépendantes. Elles peuvent même être très fortement corrélées. La conséquence principale de la non indépendance est que le passage d’une mesure par animal à plusieurs mesures ne va, parfois, renforcer que faiblement les conclusions d’un essai. Si les résultats de ce dernier sont peu significatifs quand on n’utilise que la première mesure de chaque animal, il est fréquent qu’ils ne le soient pas davantage quand on les utilise toutes.

En fait l’apport informatif d’une nouvelle mesure, sur chaque animal d’une étude, va dépendre du degré de corrélation des mesures entre elles. Si celle-ci est égale à l’unité, la nouvelle mesure sera inutile. Si au contraire elle est proche de zéro, on pourra presque admettre qu’il s’agit d’un nouvel animal.

Dans ce cadre tout se passe comme si avec chaque mesure nouvelle, on gagnait soit un degré de liberté entier, soit une part de degré de liberté seulement. C’est cette information qu’il convient  

Lorsque la répétition est temporelle, et que les écarts entre mesures successives sont réguliers, les corrélations entre n’importe quel couple de mesures recueillies à des temps différents, ne sont pas obligatoirement du même ordre. Elles peuvent être du même ordre entre des couples de mesures successives, mais être plus faibles pour des mesures plus éloignées. C’est ce qu’on remarque en particulier dans le cas de pesées à intervalle régulier  

Exemple de programme : ne fonctionne que si il n’existe pas de données manquantes

La gestion des données manquantes est centrale en analyse de données longitudinales, il n’y a pas de règles pour résoudre ce point, cependant il est fondamental de savoir si les données manquantes surviennent de manière complètement aléatoire (par rapport au phénomène étudié).  Pour le nombre de ddl, le calcul peut être très complexe voire sans solution quand on travaille avec des structures aléatoires complexes et des plans déséquilibrés. Une manière de se tirer du guêpier est d'estimer les modèles avec une méthode MCMC (Markov Chain Monte Carlo) et d'utiliser les chaînes de valeurs pour calculer les statistiques d'intérêt et construire les intervalles de confiance

Première approche : avec AOV (sans données manquantes)

demo1<-read.table("", header=T, sep=",") demo1 attach(demo1) par(cex=.6)

(time, factor(group), pulse, ylim=c(5, 20), lty=c(1, 12), ylab="mean of pulse", xlab="time", lwd=3, trace.label="group")

<- aov(pulse ~ factor(group)*factor(time) + Error(factor(id))) 

summary() detach(demo1)

id group pulse time

1        1     1    10    1

2        1     1    10    2

3        1     1    10    3

4        2     1    10    1

5        2     1    10    2

6        2     1    10    3

7        3     1    10    1

8        3     1    10    2

9        3     1    10    3 10        4     1    10    1

11       4     1    10    2

12       4     1    10    3

13       5     2    15    1

14       5     2    15    2

15       5     2    15    3

16       6     2    15    1

17       6     2    15    2

18       6     2    15    3

19       7     2    16    1

20       7     2    15    2

21       7     2    15    3

22       8     2    15    1

23       8     2    15    2

24       8     2    15    3

Error: factor(id)
Df	Sum Sq	Mean Sq         F value    Pr(>F)
factor(group)	1	155.042	155.042              3721  1.305e-09 ***
Residuals	6	0.250	0.042

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Within
Df	Sum Sq	Mean Sq          F value Pr(>F)
factor(time)                            2	0.08333	0.04167        1              0.3966
factor(group):factor(time)   2	0.08333	0.04167        1              0.3966
Residuals                                12	0.50000	0.04167

Deuxième approche : utilisation de la library nlme

library(nlme) attach(demo1)

fm <- lme(pulse ~ factor(group)*factor(time), random = ~ 1 | id)  anova(fm)

detach(demo1)

Error: factor(id)
Df	Sum Sq	Mean Sq         F value    Pr(>F)
factor(group)	1	155.042	155.042              3721  1.305e-09 ***
Residuals	6	0.250	0.042

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Within
Df	Sum Sq	Mean Sq          F value Pr(>F)
factor(time)                            2	0.08333	0.04167        1              0.3966
factor(group):factor(time)   2	0.08333	0.04167        1              0.3966

les résultats sont similaires (encore heureux…)

Test avec une donnée manquante : les résultats ne sont plus similaires

demo1bis<-read.table("", header=T) demo1bis attach(demo1bis)

<- aov(pulse ~ factor(group)*factor(time) + Error(factor(id))) summary() 

detach(demo1bis)

Error: factor(id)                                     Df     Sum Sq              Mean Sq          F value            Pr(>F)     factor(group)               1      148.30             148.30             2966   3.947e-08 *** factor(time)                 1          3.298e-29        3.298e-29        6.595e-28 1
Residuals                     5       0.25                  0.05 --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  Error: Within
Df        Sum Sq           Mean Sq         F value	Pr(>F)
factor(time)                 2         0.09524           0.04762         1.0476	0.3833
factor(group):factor(time)  2  0.07143            0.03571           0.7857 Residuals                  11           0.50000           0.04545 library(nlme) attach(demo1bis) fm <- lme(pulse ~ factor(group)*factor(time), random = ~ 1 | id)  anova(fm) detach(demo1bis)	0.4798
numDF          denDF           F-value	p-value
(Intercept)                              1                      11                   83453.02	<.0001
factor(group)                          1                      6                      3361.45	<.0001
factor(time)                            2                      11              1.04	0.3861
factor(group):factor(time)      2                      11                    0.85	0.4537

Autres exemples

demo3<-read.table("", header=T, sep=",")

demo3

id group pulse time 1     1     1    35    1 2     1     1    25    2 3     1     1    16    3 4     2     1    32    1 5     2     1    23    2 6     2     1    12    3 7     3     1    36    1 8     3     1    22    2 9     3     1    14    3 10    4     1    34    1 11    4     1    21    2 12    4     1    13    3

13    5     2    57    1 14    5     2    43    2 15    5     2    22    3 16    6     2    54    1 17    6     2    46    2 18    6     2    26    3 19    7     2    55    1 20    7     2    46    2 21    7     2    23    3 22    8     2    60    1 23    8     2    47    2 24    8     2    25    3

par(cex=.6) attach(demo3)

(time, factor(group), pulse, ylim=c(10, 60), lty=c(1, 12),                  lwd=3, ylab="mean of pulse", xlab="time", trace.label="group")

<- aov(pulse ~ factor(group)*factor(time) + Error(factor(id))) summary()

detach(demo3)

Error: factor(id)

              Df  Sum Sq Mean Sq F value    Pr(>F) factor(group)  1 2035.04 2035.04  343.15 1.596e-06 *** Residuals      6 35.58    5.93                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Within

                           Df  Sum Sq Mean Sq F value Pr(>F)     factor(time)                2 2830.33 1415.17 553.761 1.517e-12 ***

factor(group):factor(time)  2  200.33  100.17  39.196 5.474e-06 ***

Residuals                  12   30.67    2.56                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

demo4<-read.table("", header=T, sep=",")

demo4

par(cex=.6) attach(demo4)

(time, factor(group), pulse, ylim=c(10, 60), lty=c(1, 12), lwd=3,ylab="mean of pulse", xlab="time", trace.label="group") <- aov(pulse ~ factor(group)*factor(time) + Error(factor(id)))

summary() detach(demo4)

      id group pulse time

1     1     1    35    1

2     1     1    25    2

3     1     1    12    3

4     2     1    34    1

5     2     1    22    2

6     2     1    13    3

7     3     1    36    1

8     3     1    21    2

9     3     1    18    3

10    4     1    35    1

11    4     1    23    2

12    4     1    15    3

13    5     2    31    1

14    5     2    43    2

15    5     2    57    3

16    6     2    35    1

17    6     2    46    2

18    6     2    58    3

19    7     2    37    1

20    7     2    48    2

21    7     2    51    3

22    8     2    32    1

23    8     2    45    2

24    8     2    53    3

Error: factor(id)

              Df  Sum Sq Mean Sq F value    Pr(>F) factor(group)  1 2542.04 2542.04  628.96 2.646e-07 ***

Residuals      6   24.25    4.04                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Error: Within

                           Df  Sum Sq Mean Sq  F value Pr(>F)    

factor(time)                2    1.00    0.50   0.0789    0.9246

factor(group):factor(time)  2 1736.33  868.17 137.0789 5.438e-09 ***

Residuals                  12   76.00    6.33

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Réflexion sur les structures de variance-covariance

Variance-Covariance Structures

Independence

As though analyzed using between subjects analysis.

s2 

0   s2

0   0   s2

Compound Symmetry

Assumes that the variance-covariance structure has a single variance (represented by s2) for all 3 of the time points and a single covariance (represented by s1) for each of the pairs of trials.  This structure is illustrated by the half matrix below.

s2 s1 s2 s1 s1 s2 

Unstructured

Assumes that each variance and covariance is unique.   Each trial has its own variance (e.g. s12 is the variance of trial 1) and each pair of trials has its own covariance (e.g. s21 is the covariance of trial 1 and trial2).  This structure is illustrated by the half matrix below.

s12 s21 s22 s31 s32 s32

Autoregressive

Another common covariance structure which is frequently observed in repeated measures data is an autoregressive structure, which recognizes that observations which are more proximate are more correlated than measures that are more distant. This structure is illustrated by the half matrix below.

s2 sr s2 sr2 srs2

Autoregressive Heterogeneous Variances

If the variances change over time, then the covariance would look like this.

s12 sr     s22 sr2   srs32

Exercice complet

exer<-read.table("", header=T) exer

> exer<-read.table("", header=T)

> exer

   id diet exertype pulse time 31 11    1        2    86    1         62 21    1        3    98    2 1   1    1        1    85    1                   32 11    1        2    86    2         63 21    1        3   110    3 2   1 1        1    85    2              33 11    1        2    84    3         64 22    1        3    98    1

3   1    1        1    88    3           34 12    1        2    93    1         65 22    1        3   104    2

4   2    1        1    90    1           35 12    1        2   103    2        66 22    1        3   112    3 5   2    1        1    92    2                  36 12    1        2   104    3        67 23    1        3    98    1

6   2    1        1    93    3           37 13    1        2    90    1         68 23    1        3   105    2

7   3    1        1    97    1           38 13    1        2    92    2         69 23    1        3    99    3

8   3    1        1    97    2           39 13    1        2    93    3         70 24    1        3    87    1

9   3    1        1    94    3           40 14    1        2    95    1         71 24    1        3   132    2

10 4    1        1    80    1           41 14    1        2    96    2         72 24    1        3   120    3 11  4    1        1    82    2                 42 14    1        2   100    3        73 25    1        3    94    1

12  4    1        1    83    3          43 15    1        2    89    1         74 25    1        3   110    2

13  5    1        1    91    1          44 15    1        2    96    2         75 25    1        3   116    3 14  5    1        1    92    2                 45 15    1        2    95    3         76 26    2        3    95    1

15  5    1        1    91    3          46 16    2        2    84    1         77 26    2        3   126    2

16  6    2        1    83    1          47 16    2        2    86    2         78 26    2        3   143    3

17  6    2        1    83    2          48 16    2        2    89    3         79 27    2        3   100    1

18  6    2        1    84    3          49 17    2        2   103    1        80 27    2        3   126    2

19  7    2        1    87    1          50 17    2        2   109    2        81 27    2        3   140    3

20  7    2        1    88    2          51 17    2        2    90    3         82 28    2        3   103    1

21  7    2        1    90    3          52 18    2        2    92    1         83 28    2        3   124    2

22  8    2        1    92    1          53 18    2        2    96    2         84 28    2        3   140    3 23  8    2        1    94    2                 54 18    2        2   101    3        85 29    2        3    94    1

24  8    2        1    95    3          55 19    2        2    97    1         86 29    2        3   135    2

25  9    2        1    97    1          56 19    2        2    98    2         87 29    2        3   130    3 26  9    2        1    99    2                 57 19    2        2   100    3        88 30    2        3    99    1

27  9    2        1    96    3          58 20    2        2   102    1        89 30    2        3   111    2

28  10    2        1   100    1       59 20    2        2   104    2        90 30    2        3   150    3

29  10    2        1    97    2        60 20    2        2   103    3

30  10    2        1   100    3       61 21    1        3    93    1

(time, factor(diet), pulse, ylim=c(90, 110), lty=c(1, 12),

+ lwd=3,ylab="mean of pulse", xlab="time", trace.label="group") > 

mat <- matrix(c(pulse[time==1], pulse[time==2], pulse[time==3]), ncol=3) var(mat) 

[,1]	[,2]	[,3]
[1,]      37.84368	48.78851	60.28506
[2,]      48.78851	212.11954	233.76092
[3,]      60.28506	233.76092	356.32299

cor(mat)

[,1]	[,2]	[,3]
[1,]      1.0000000	0.5445409	0.5191479
[2,]      0.5445409	1.0000000	0.8502755
[3,]      0.5191479	0.8502755	1.0000000

library(nlme) 

longg <- groupedData(pulse ~ exertype*time | id, data=exer)

<- gls(pulse ~ factor(exertype)*factor(time), data=longg, corr=corCompSymm(, form=

~ 1 | id) )

summary() 

Generalized least squares fit by REML

  Model: pulse ~ factor(exertype) * factor(time) 

  Data: longg 

       AIC      BIC    logLik

  612.8316 639.1706 -295.4158

Correlation Structure: Compound symmetry

 Formula: ~1 | id   Parameter estimate(s):

      Rho 

0.4558160 

anova()

Denom. DF: 81 

numDF  F-value p-value (Intercept)                       1 5802.399  <.0001 factor(exertype)                  2   27.001  <.0001 factor(time)                      2   23.543  <.0001

factor(exertype):factor(time)     4   15.512  <.0001

<- gls(pulse ~ factor(exertype)*factor(time), data=longg, corr=corSymm(form = ~ 1 | id), weights = varIdent(form = ~ 1 | time))

summary()

Generalized least squares fit by REML

  Model: pulse ~ factor(exertype) * factor(time) 

  Data: longg 

       AIC      BIC    logLik

  607.7365 643.6532 -288.8682

Correlation Structure: General  Formula: ~1 | id   Parameter estimate(s):

 Correlation: 

1    2    

2    0.434      

3    0.417 0.583 Variance function:

 Structure: Different standard deviations per stratum  Formula: ~1 | time   Parameter estimates:

       1        2        3 

1.000000 1.596720 1.877599 

anova()

Denom. DF: 81 

numDF  F-value p-value (Intercept)                       1 8184.123  <.0001 factor(exertype)                  2    6.426  0.0026 factor(time)                      2   22.324  <.0001

factor(exertype):factor(time)     4   14.387  <.0001

fit.ar1 <- gls(pulse ~ factor(exertype)*factor(time), data=longg, corr=corAR1(, form= ~ 1 | id)

)

summary(fit.ar1) 

Generalized least squares fit by REML

  Model: pulse ~ factor(exertype) * factor(time) 

  Data: longg 

       AIC      BIC    logLik

  612.1163 638.4553 -295.0582

Correlation Structure: AR(1)  Formula: ~1 | id   Parameter estimate(s):

      Phi 

0.4992423 

anova(fit.ar1)

numDF  F-value p-value (Intercept)                       1 6167.352  <.0001 factor(exertype)                  2   26.990  <.0001 factor(time)                      2   18.196  <.0001

factor(exertype):factor(time)     4   11.733  <.0001

fit.arh1 <- gls(pulse ~ factor(exertype)*factor(time), data=longg, corr=corAR1(, form= ~ 1 | id), weight=varIdent(form = ~ 1 | time)) summary(fit.arh1) 

Generalized least squares fit by REML

  Model: pulse ~ factor(exertype) * factor(time) 

  Data: longg 

       AIC      BIC    logLik

  605.7693 636.8971 -289.8846

Correlation Structure: AR(1)  Formula: ~1 | id   Parameter estimate(s):

      Phi 

0.5100781 

Variance function:

 Structure: Different standard deviations per stratum  Formula: ~1 | time   Parameter estimates:

       1        2        3 

1.000000 1.561315 1.796993 

anova(fit.arh1)

Denom. DF: 81 

numDF  F-value p-value (Intercept)                       1 8284.813  <.0001 factor(exertype)                  2    9.134   3e-04 factor(time)                      2   21.918  <.0001

factor(exertype):factor(time)     4   13.805  <.0001

We can use the anova function to compare competing models to see which model fits the data best. 

anova(, )

Model df       AIC	BIC	logLik              Test  L.Ratio  p-value
1         11        612.8316	639.1706	-295.4158
2        15        607.7365  anova(fit.ar1, )	643.6532	-288.8682          1 vs 2 13.095 12  0.0108
Model df       AIC	BIC	logLik              Test  L.Ratio  p-value
fit.ar1     1       11        612.1163	638.4553	-295.0582
2       15        607.7365   anova(fit.arh1, )	643.6532	-288.8682        1 vs 2 12.37984  0.0147
Model df       AIC	BIC	logLik              Test L.Ratio  p-value
fit.arh1    1 13             605.7693	636.8971	-289.8846
2 15             607.7365  anova(fit.ar1, )	643.6532	-288.8682        1 vs 2 2.03276  0.3619
Model df             AIC	BIC	logLik
fit.ar1    1 11               612.1163	638.4553	-295.0582
2 11               612.8316  anova(fit.ar1, fit.arh1)	639.1706	-295.4158
Model df            AIC	BIC	logLik             Test  L.Ratio p-value

fit.ar1      1 11             612.1163        638.4553         -295.0582                        

fit.arh1     2 13            605.7693        636.8971         -289.8846       1 vs 2 10.34708  0.0057

The two most promising structures are Autoregressive Heterogeneous Variances and Unstructured since these two models have the smallest AIC values and the -2 Log

Likelihood scores are significantly smaller than the -2 Log Likelihood scores of other models

LOGICIEL R 

VERSION 2.6.0.

 ANALYSE DE DONNEES

Non-supervisées

                              Supervisées

ANALYSE DE DONNEES NON SUPERVISEE

L’Analyse en Composantes Principales

Méthode 1

olympic <- read.table("")

names(olympic)

names(olympic) <- c("100", "long", "poid", "haut", "400", "110",

"disq", "perc", "jave", "1500", "score")

names(olympic)

Olympic

100     long  poid  haut  400  110    disq    perc  jave  1500 score

1   11.25 7.43 15.48 2.27 48.90 15.13 49.28  4.7 61.32 268.95  8488

2   10.87 7.45 14.97 1.97 47.71 14.46 44.36  5.1 61.76 273.02  8399

3   11.18 7.44 14.20 1.97 48.29 14.81 43.66  5.2 64.16 263.20  8328

4   10.62 7.38 15.02 2.03 49.06 14.72 44.80  4.9 64.04 285.11  8306

5   11.02 7.43 12.92 1.97 47.44 14.40 41.20  5.2 57.46 256.64  8286

6   10.83 7.72 13.58 2.12 48.34 14.18 43.06  4.9 52.18 274.07  8272

7   11.18 7.05 14.12 2.06 49.34 14.39 41.68  5.7 61.60 291.20  8216

8   11.05 6.95 15.34 2.00 48.21 14.36 41.32  4.8 63.00 265.86  8189

9   11.15 7.12 14.52 2.03 49.15 14.66 42.36  4.9 66.46 269.62  8180

10  11.23 7.28 15.25 1.97 48.60 14.76 48.02  5.2 59.48 292.24  8167

11  10.94 7.45 15.34 1.97 49.94 14.25 41.86  4.8 66.64 295.89  8143

12  11.18 7.34 14.48 1.94 49.02 15.11 42.76  4.7 65.84 256.74  8114

13  11.02 7.29 12.92 2.06 48.23 14.94 39.54  5.0 56.80 257.85  8093

14  10.99 7.37 13.61 1.97 47.83 14.70 43.88  4.3 66.54 268.97  8083

15  11.03 7.45 14.20 1.97 48.94 15.44 41.66  4.7 64.00 267.48  8036

16  11.09 7.08 14.51 2.03 49.89 14.78 43.20  4.9 57.18 268.54  8021

17  11.46 6.75 16.07 2.00 51.28 16.06 50.66  4.8 72.60 302.42  7869

18  11.57 7.00 16.60 1.94 49.84 15.00 46.66  4.9 60.20 286.04  7860

19  11.07 7.04 13.41 1.94 47.97 14.96 40.38  4.5 51.50 262.41  7859

20  10.89 7.07 15.84 1.79 49.68 15.38 45.32  4.9 60.48 277.84  7781

21  11.52 7.36 13.93 1.94 49.99 15.64 38.82  4.6 67.04 266.42  7753

22  11.49 7.02 13.80 2.03 50.60 15.22 39.08  4.7 60.92 262.93  7745

23  11.38 7.08 14.31 2.00 50.24 14.97 46.34  4.4 55.68 272.68  7743

24  11.30 6.97 13.23 2.15 49.98 15.38 38.72  4.6 54.34 277.84  7623

25  11.00 7.23 13.15 2.03 49.73 14.96 38.06  4.5 52.82 285.57  7579

26  11.33 6.83 11.63 2.06 48.37 15.39 37.52  4.6 55.42 270.07  7517

27  11.10 6.98 12.69 1.82 48.63 15.13 38.04  4.7 49.52 261.90  7505

28  11.51 7.01 14.17 1.94 51.16 15.18 45.84  4.6 56.28 303.17  7422

29  11.26 6.90 12.41 1.88 48.24 15.61 38.02  4.4 52.68 272.06  7310

30  11.50 7.09 12.94 1.82 49.27 15.56 42.32  4.5 53.50 293.85  7237 31 11.43 6.22 13.98 1.91 51.25 15.88 46.18  4.6 57.84 294.99  7231

32  11.47 6.43 12.33 1.94 50.30 15.00 38.72  4.0 57.26 293.72  7016

33  11.57 7.19 10.27 1.91 50.71 16.20 34.36  4.1 54.94 269.98  6907

34  12.12 5.83  9.71 1.70 52.32 17.05 27.10  2.6 39.10 281.24  5339

# recherche des individus aberrants ou des variables redondantes apply(olympic, 2, function(x) (1:34)[x == min(x)])

100   long   poid   haut    400    110   disq   perc   jave   1500 score 

   4     34     34 34      5                   6              34     34     34      5              34

apply(olympic, 2, function(x) (1:34)[x == max(x)])

100   long   poid   haut 400    110   disq                    perc   jave   1500  score 

   34     6           18 1               34     34     17      7                17     28      1

boxplot(.frame(scale(olympic)))

# on retire ainsi l’individu ayant le plus petit score et la variable score olympic <- olympic[-34, ] olyres <- olympic$score olympic <- olympic[, -11] names(olympic)

olympic

prcomp1 <- prcomp(olympic, scale = T)

summary(prcomp1)

Importance of components:

                                    PC1 PC2     PC3    PC4     PC5     PC6     PC7    PC8

Standard deviation     1.849 1.614 0.9712 0.9370 0.7461 0.7009 0.6562 0.5539

Proportion of Variance 0.342 0.261 0.0943 0.0878 0.0557 0.0491 0.0431 0.0307

Cumulative Proportion  0.342 0.602 0.6968 0.7846 0.8403 0.8894 0.9324 0.9631

                                     PC9 PC10

Standard deviation     0.5167 0.3191

Proportion of Variance 0.0267 0.0102

Cumulative Proportion  0.9898 1.0000

print(round(var(prcomp1$x), digits = 3))

PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10

PC1  3.418 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

PC2  0.000 2.606 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

PC3  0.000 0.000 0.943 0.000 0.000 0.000 0.000 0.000 0.000 0.000

PC4  0.000 0.000 0.000 0.878 0.000 0.000 0.000 0.000 0.000 0.000

PC5  0.000 0.000 0.000 0.000 0.557 0.000 0.000 0.000 0.000 0.000

PC6  0.000 0.000 0.000 0.000 0.000 0.491 0.000 0.000 0.000 0.000

PC7  0.000 0.000 0.000 0.000 0.000 0.000 0.431 0.000 0.000 0.000

PC8  0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.307 0.000 0.000

PC9  0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.267 0.000

PC10 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.102

print(round(t(prcomp1$rotation) %*% prcomp1$rotation, digits = 3)) prcomp1$sdev^2 eigen(cor(olympic))$values

[1] 3.4182381 2.6063931 0.9432964 0.8780212 0.5566267 0.4912275 0.4305952

 [8] 0.3067981 0.2669494 0.1018542

plot(prcomp1)

round(cor(olympic), digits = 2)

100 long  poid  haut   400   110  disq  perc  jave  1500 100   1.00 -0.54 -0.21 -0.15  0.61  0.64 -0.05 -0.39 -0.06  0.26 long -0.54  1.00  0.14  0.27 -0.52 -0.48  0.04  0.35  0.18 -0.40 poid -0.21  0.14  1.00  0.12  0.09 -0.30  0.81 0.48  0.60  0.27 haut -0.15  0.27  0.12  1.00 -0.09 -0.31  0.15  0.21  0.12 -0.11 400   0.61 -0.52  0.09 -0.09  1.00  0.55  0.14 -0.32  0.12  0.59 110 0.64 -0.48 -0.30 -0.31  0.55  1.00 -0.11 -0.52 -0.06  0.14 disq -0.05  0.04 0.81  0.15  0.14 -0.11  1.00  0.34  0.44  0.40 perc -0.39  0.35  0.48  0.21 -0.32 -0.52  0.34  1.00  0.27 -0.03 jave -0.06  0.18  0.60  0.12  0.12 -0.06 0.44  0.27  1.00  0.10

1500  0.26 -0.40  0.27 -0.11  0.59  0.14  0.40 -0.03  0.10  1.00

Méthode 2

prin1 <- princomp(olympic, cor = T) names(prin1)

[1] "sdev"     "loadings" "center" "scale"    "n.obs"    "scores" "call"   

biplot(prin1)

round(prin1$sdev^2, digit=3)

Comp.1  Comp.2  Comp.3  Comp.4  Comp.5  Comp.6  Comp.7  Comp.8  Comp.9 Comp.10 

  3.418     2.606      0.943      0.878      0.557      0.491 0.431      0.307      0.267      0.102

print(round(t(prin1$loadings) %*% prin1$loadings, digits = 3))

Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10

Comp.1       1      0      0      0      0      0      0      0 0       0

Comp.2       0      1      0      0      0      0      0      0 0       0

Comp.3       0      0      1      0      0      0      0      0 0       0

Comp.4       0      0      0      1      0      0      0      0 0       0

Comp.5       0      0      0      0      1      0      0      0 0       0

Comp.6       0      0      0      0      0      1      0      0 0       0

Comp.7       0      0      0      0      0      0      1      0 0       0

Comp.8       0      0      0      0      0      0      0      1 0       0

Comp.9       0      0      0      0      0      0      0      0 1       0

Comp.10      0      0      0      0      0      0      0      0 0       1

print(round(var(prin1$scores), digits = 3))

Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10

Comp.1   3.525  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000   0.000

Comp.2   0.000  2.688  0.000  0.000  0.000  0.000  0.000  0.000 0.000   0.000

Comp.3   0.000  0.000  0.973  0.000  0.000  0.000  0.000  0.000 0.000   0.000

Comp.4   0.000  0.000  0.000  0.905  0.000  0.000  0.000  0.000 0.000   0.000

Comp.5   0.000  0.000  0.000  0.000  0.574  0.000  0.000  0.000 0.000   0.000

Comp.6   0.000  0.000  0.000  0.000  0.000  0.507  0.000  0.000 0.000   0.000

Comp.7   0.000  0.000  0.000  0.000  0.000  0.000  0.444  0.000 0.000   0.000

Comp.8   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.316 0.000   0.000

Comp.9   0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.275   0.000

Comp.10  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000   0.105

round(prin1$center, digit=2)

100     long   poid     haut    400    110     disq     perc jave    1500 

 11.20   7.13  13.98   1.98  49.28  15.05  42.35   4.74  59.44 276.04

round(prin1$scale, digit=2)

100  long  poid  haut   400   110  disq  perc  jave  1500 

 0.24  0.30  1.31  0.09  1.05  0.50  3.66  0.33  5.41 13.45

round(prin1$n.obs, digit=3)

prin1$call

Analyse Factorielle des Correspondances

sond <- read.table(« ")

sond library(MASS) sond.corresp <- corresp(sond) names(sond.corresp)

dotchart(c(sort(sond.corresp$cscore), sort(sond.corresp$rscore)))

sond <- sond[-8, ] sond.corresp <- corresp(sond)

dotchart(c(sort(sond.corresp$cscore), sort(sond.corresp$rscore)))

round(sond.corresp$cscore, digit=3)

ExtGauche                  PC            PS      Ecolo                Droite  ExtDroite     Divers 

    1.460                       0.880          1.070     0.738    -1.102    -     0.431               0.663

library(ade4)

<- (sond, scannf = F) t($c1[, 1, drop = F] ExtGauche         PC        PS      Ecolo   Droite ExtDroite     Divers

CS1 -1.460312 -0.8802587 -1.069967 -0.7382985 1.102263 0.4307402 -0.6627689

scatter()

Exemple ACP et Classification Ascendante Hiérarchique

ACP

rm(list=ls())

unicoolait <- read.table("", header=TRUE,check.names=TRUE) unicoolait$STA2006 <- as.factor(unicoolait$STA2006)

unicoolait$STAB2007 <- as.factor(unicoolait$STAB2007)

#sélection des variables dans l'ACP uni<-unicoolait[,c(3,6,8,16,14,13)] uni

#matrice de corrélation

round(cor(uni), digits = 2)

         Lait TB    TP Jours  Nlac logcel

Lait    1.00 -0.34 -0.53 -0.56  0.22  -0.18

TB     -0.34 1.00  0.53  0.24 -0.07   0.09

TP     -0.53 0.53  1.00  0.59 -0.05   0.23

Jours -0.56  0.24  0.59  1.00 -0.10   0.19 Nlac    0.22 -0.07 -0.05 -0.10  1.00 0.27

logcel -0.18 0.09  0.23  0.19  0.27   1.00

#calcul des composantes principales, valeur propre, vecteur propre prcomp1 <- prcomp(uni, scale = T)

summary(prcomp1)

Importance of components:

PC1   PC2   PC3   PC4    PC5   PC6
Standard deviation	1.586 1.131 0.911 0.793 0.6597 0.558
Proportion of Variance	0.419 0.213 0.138 0.105 0.0725 0.052
Cumulative Proportion	0.419 0.632 0.771 0.875 0.9480 1.000

print(round(var(prcomp1$x), digits = 3))

            PC1   PC2 PC3  PC4   PC5   PC6

PC1  2.515 0.000 0.000 0.00 0.000 0.000

PC2  0.000 1.279 0.000 0.00 0.000 0.000

PC3  0.000 0.000 0.829 0.00 0.000 0.000

PC4  0.000 0.000 0.000 0.63 0.000 0.000

PC5  0.000 0.000 0.000 0.00 0.435 0.000

PC6  0.000 0.000 0.000 0.00 0.000 0.312

print(round(t(prcomp1$rotation) %*% prcomp1$rotation, digits = 3))

            PC1 PC2 PC3 PC4 PC5 PC6

PC1    1   0   0   0   0   0

PC2    0   1   0   0   0   0

PC3    0   0   1   0   0   0

PC4    0   0   0   1   0   0

PC5    0   0   0   0   1   0

PC6    0   0   0   0   0   1

prcomp1$sdev^2

[1] 2.5150014 1.2792946 0.8290995 0.6295761 0.4352603 0.3117680

#valeurs propres

eigen(cor(uni))$values

[1] 2.5150014 1.2792946 0.8290995 0.6295761 0.4352603 0.3117680

plot(prcomp1)

#graphique corrélation et individus prin1<-prcomp(uni, scale = T)

biplot(prin1)

#cercle des corrélations (axe 1 et 2) library(ade4)

s.corcircle(.frame(prcomp1$rotation[, 1] * prcomp1$sdev[1], prcomp1$rotation[, 2] * prcomp1$sdev[2]))

Classification Ascendante Hiérarchique en 15 groupes

don<-scale(uni, center = TRUE, scale = TRUE) dc<-dist(don, method ="euclidean", diag=FALSE, upper=FALSE) hier<-hclust(dc,"ward") plot(hier,hang=-1)

cl15<-cutree(hier,15)

Moyenne par groupe pour matrice uni doncl<-merge(uni, cl15,by= "row.names") moyennes<-aggregate(doncl[,-1], list(doncl$y), mean) moyennes[,-1]<-round(moyennes[,-1],2)

moyennes

   Group.1  Lait    TB    TP  Jours Nlac logcel  y

1       1 22.20 36.35 30.46 203.67 3.79   5.13  1

2       2 18.41 37.36 31.21 246.67 6.40   4.47  2

3       3 11.92 50.40 37.74 308.67 2.33   4.39  3

4       4 21.33 36.18 30.33 210.27 2.11   4.46  4

5       5 18.14 42.88 33.61 272.39 3.96   4.35  5

6       6 16.35 44.37 36.83 532.83 2.03   4.30  6

7       7 25.18 43.23 31.28 118.04 3.60   4.02  7

8       8 14.55 37.83 33.94 276.40 1.98   4.50  8

9       9 22.26 35.44 28.61  86.67 1.43   3.92  9

10      10 31.16 34.43 27.49  79.23 3.05   3.83 10

11      11 17.44 35.23 30.87 269.86 1.28   4.10 11

12      12 18.19 41.28 31.33 209.97 1.59   4.03 12

13      13 20.20 36.32 30.43 235.74 1.42   3.60 13

14      14 26.83 31.68 28.57 139.65 5.35   3.85 14

15      15 17.45 44.83 34.78 264.71 1.32   3.72 15

Ecar-type par groupe doncl<-merge(uni, cl15,by= "row.names") ecartype<-aggregate(doncl[,-1], list(doncl$y), sd) ecartype[,-1]<-round(ecartype[,-1],2) ecartype

   Group.1 Lait   TB   TP Jours Nlac logcel y

1       1 4.13 5.35 2.56 72.17 1.46   0.34 0

2       2 4.62 5.74 2.46 91.29 1.44   0.36 0

3       3 3.98 5.58 3.80 72.23 1.11   0.36 0

4       4 3.22 3.05 1.82 47.96 0.96   0.23 0

5       5 3.02 3.42 1.53 36.11 0.91   0.31 0

6       6 3.51 6.92 2.98 91.84 1.03   0.29 0

7       7 4.31 6.83 2.71 76.56 1.60   0.44 0

8       8 3.15 4.19 1.89 43.99 1.00   0.25 0

9       9 3.55 3.94 2.26 53.45 0.64   0.31 0

10      10 3.73 3.82 2.12 37.73 1.08   0.45 0

11      11 3.31 3.47 1.61 45.69 0.45   0.22 0

12      12 2.67 3.60 1.65 55.04 0.69   0.17 0

13      13 2.88 3.64 1.81 54.56 0.62   0.19 0

14      14 5.42 3.71 1.32 70.45 1.14   0.21 0

15      15 4.56 2.64 1.48 64.57 0.48   0.20 0

Que donne la suite, à vous de comprendre…

#fusion des tableaux doncluni<-doncl[,-1]

unitotal<-merge(doncluni, unicoolait)

#moyenne par groupe

moyennesbis<-aggregate(unitotal[,c(-18,-19)],list(unitotal$y), mean) round(moyennesbis[,-1],2)

#corrélation entre stade de lactation et instabilité

cor(moyennesbis$Jour,moyennesbis$Inst4)

#matrice de corrélation générale

round(cor(moyennesbis[,-1]),2)

plot(moyennesbis$Jour,moyennesbis$Inst4)

abline(lsfit(moyennesbis$Jour,moyennesbis$Inst4))

#comptage du nombre d'animaux par classe doncl$y<-as.factor(doncl$y)

summary(doncl$y)

# très joli résumé library(cluster)

clusplot(dc,cl15,diss=T,shade=T,color=T,labels=5,main="") abline(v=0,h=0)

AFCM et Classification Ascendante Hiérarchique

Exemple tiré de la formation INRA

Copyright © André Bouchier.

© 2006, André Bouchier (20 Janvier 2006)

Dans un premier temps, les données quantitatives sont découpées en classes. L'objectif est d'obtenir un tableau de données qualitatives. Puis on transforme les données de façon à obtenir un tableau disjonctif. Ce tableau disjonctif peut-être considéré comme un tableau de fréquence. On peut alors le soumettre à l'analyse factorielle des correspondances. L'analyse factorielle des correspondances multiples permet d'obtenir un résumé pertinent d'un tableau de données disjonctif. Ce résumé peut-être projeté sur un plan. On visualise ainsi les relations les plus marquantes du tableau de données.

Le tableau des données contient 50 observations et 11 variables. Il contient les résultats d'un suivi agronomique sur 50 parcelles de blé dur.

RDT Rendement en grains

PLM Nb de plantes par m²

ZON Zone géographique

ARG Taux d'argile de la parcelle

LIM Taux de limon de la parcelle

SAB Taux de sable de la parcelle

VRT Variété cultivée

PGM Poids de 1000 grains

MST Matière sèche totale à la récolte

AZP Azote dans la plante à la récolte

VRTC Variété cultivée (codée en 3 classes)

Listing

rm(list=ls())

#choisir son fichier TXT (choisir le fichier ) don<-read.table(file.choose(), sep=" ", header=T, dec=",")

#découpage en 3 classes d'effectifs égaux (fonction codage)

codage<-function(nom)

{

#calcul des bornes

bornes<-quantile(nom, probs = c(0, 1/3,2/3,1), = TRUE,names = TRUE)

#description des bornes et effectifs

Amax<-aggregate(nom,list(Nom=cut(nom,bornes,include.lowest=T,label=F)),max)

Amin<-aggregate(nom,list(Nom=cut(nom,bornes,include.lowest=T,label=F)),min) Afreq<-as.matrix(summary(as.factor(cut((nom),bornes, include.lowest=T, label=F))))

limites<.frame(cbind(Amin[,1],Amin[,2],Amax[,2],Afreq)) names(limites)<-c("Classe","Mini","Maxi","Effectif")

#calcul du nombre de valeurs manquantes

manques<-length(nom)-length((nom))

#impression des bornes

cat(paste("Découpage de la variable ",deparse(substitute(nom))," - Nb de valeurs manquantes : ", manques, "\n"))

print(limites)

#découpage de la variable

varfac<-cut(nom,bornes,include.lowest=T,label=F)

#transformation en facteur as.factor(varfac)

}

RDT<-codage(don$RDT)

Découpage de la variable  don$RDT  - Nb de valeurs manquantes :  0 

  Classe  Mini   Maxi Effectif

1      1  0.12  8.540       17

2      2  8.70 12.106       16

3      3 12.23 21.670       17

PLM<-codage(don$PLM)

Découpage de la variable  don$PLM  - Nb de valeurs manquantes :  0 

  Classe Mini Maxi Effectif

1      1   46   92       17

2      2   95  130       16

3      3  132  252       17

ARG<-codage(don$ARG)

Découpage de la variable  don$ARG  - Nb de valeurs manquantes :  0    Classe Mini Maxi Effectif

1      1  9.6 19.2       17

2      2 19.4 26.7       16

3      3 27.2 48.4       17

LIM<-codage(don$LIM)

Découpage de la variable  don$LIM  - Nb de valeurs manquantes :  0    Classe Mini Maxi Effectif

1      1  4.0 29.5       17

2      2 29.7 40.9       16 3      3 42.0 68.0       17 SAB<-codage(don$SAB)

Découpage de la variable  don$SAB  - Nb de valeurs manquantes :  0 

  Classe Mini Maxi Effectif

1      1 11.7 26.2       17

2      2 28.3 49.3       16

3      3 49.4 86.4       17

PGM<-codage(don$PGM)

Découpage de la variable  don$PGM  - Nb de valeurs manquantes :  0 

  Classe  Mini  Maxi Effectif

1      1 24.00 33.91       17

2      2 34.20 41.26       16

3      3 41.53 53.90       17

MST<-codage(don$MST)

Découpage de la variable  don$MST  - Nb de valeurs manquantes :  0 

  Classe  Mini  Maxi Effectif

1      1 13.72 26.89       17

2      2 27.08 34.49       16

3      3 35.51 59.50       17

AZP<-codage(don$AZP)

Découpage de la variable  don$AZP  - Nb de valeurs manquantes :  0 

  Classe Mini Maxi Effectif

1      1 1.82 2.99       17

2      2 3.04 3.51       16

3      3 3.57 4.27       17

# vous pouvez découper les variables quantitatives en qualitatives en choisissant vous- #même les bornes de classes :  

histc2<-c(0,5,10,15,20,25) max(don$RDT)

x1<-cut(don$RDT, breaks=histc2, include.lowest = TRUE) x2<-data.frame(x1, don)

x1 Numero    RDT PLM ZON  ARG  LIM  SAB VRT   PGM   MST  AZP VRTC 1 (5,10]      1  6.490  84   1 21.5 60.6 17.9   3 43.10 34.49 3.82    2

2    (15,20] 2 15.580 112   1 21.0 58.3 20.7   3 38.30 39.18 3.78    2

3    (5,10] 3  7.290  68   1 26.2 47.6 26.2   3 45.30 26.89 2.61    2

4    [0,5] 4  1.090  88   1 29.7 54.5 15.8   3 29.09 23.09 3.78    2

5    (5,10] 5  5.100 174   1 22.8 59.0 18.2   3 42.80 18.10 3.41    2 6    [0,5]      6 2.030  63   1 19.6 68.0 12.4   3 41.26 20.43 3.04    2

7    (5,10] 7  6.330  92   1 26.7 53.7 19.6   3 38.57 20.93 2.26    2

8    (15,20] 8 17.300 117   1 34.0 44.9 21.1   1 31.80 40.70 3.85    1

9    (5,10] 9  6.970  58   1 16.7 57.6 25.7   4 42.40 29.97 2.69    3

10  (20,25] 10 20.350 116   1 37.4 49.0 13.6   6 42.90 40.75 3.62    3

11  (5,10] 11  7.888  80   1 44.2 44.1 11.7   1 39.10 22.29 3.61    1

12  (5,10]  12  9.380  46   1 27.5 44.2 28.3   1 41.06 29.78 3.32    1

13  (10,15] 13 11.000 189   1 39.4 34.5 36.1   2 41.53 29.60 3.51    2

14  (5,10] 14  9.852 112   1 48.4 29.5 22.1   1 33.56 32.35 2.78    1

15  (15,20] 15 17.330 180   1 32.4 47.5 20.1   2 40.50 38.03 4.02    2

16  (10,15] 16 12.230 234   1 37.6 40.9 21.5   2 34.70 35.63 4.27    2

17  (10,15] 17 12.570 162   1 38.5 44.8 16.7   1 44.50 53.97 3.22    1

18  (5,10] 18  6.323  72   2 15.4 27.0 57.6   1 44.20 29.32 2.99    1

19  [0,5] 19  1.416  86   2 13.0 26.7 60.3   1 33.54 18.42 3.35    1

20  [0,5] 20  0.520  67   2 14.8 23.5 61.7   1 33.91 13.72 3.06    1

21  (15,20] 21 15.980  95   2 18.2 31.6 50.2   1 50.80 34.38 4.17    1

22  (10,15] 22 12.850  88   2 19.7 32.5 47.8   1 39.20 36.65 2.04    1

23  (10,15] 23 12.040  95   2 19.1 44.3 36.6   5 42.70 34.04 3.10    3

24  (10,15] 24 14.110 140   2 22.8 34.9 42.3   5 44.90 28.31 2.72    3

25  (20,25] 25 21.670 126   2 26.6 39.8 33.6   1 47.20 57.67 4.17    1

26  [0,5] 26  3.236  96   2 18.9 40.8 40.3   1 30.22 22.64 1.82    1

27  (20,25] 27 20.210 130   2 31.9 42.0 26.1   1 40.00 45.21 3.18    1

28  [0,5] 28  0.120  90   2 20.0 36.7 43.3   1 24.00 16.52 3.36    1

29  (10,15] 29 12.106 129   2 21.0 29.7 49.3   2 34.20 35.51 3.21    2

30  (10,15] 30 11.010  87   2 22.3 28.3 49.4   1 41.60 32.01 3.26    1

31  (10,15] 31 13.420 150   2 27.5 36.7 35.8   2 29.80 40.42 4.11    2

32  (15,20] 32 15.868 174   2 32.6 33.7 33.7   1 53.90 41.37 3.57    1

33  [0,5] 33  4.800 132   3 15.6 12.1 72.3   3 36.51 16.20 2.86    2

34  (5,10] 34  9.180 185   3 12.4 10.9 76.7   1 30.40 23.18 4.19    1

35  (5,10] 35  8.700 166   3 16.3 29.9 53.8   2 37.10 20.30 3.07    2

36  (5,10] 36  7.920 102   3 16.9 27.4 55.7   2 27.10 19.32 2.67    2

37  (5,10] 37  8.540 120   3 16.1 26.4 57.5   2 25.30 30.54 2.02    2

38  (5,10] 38  9.420 119   3 19.4 40.2 40.4   2 27.50 27.82 2.89    2

39  (10,15] 39 12.460 122   3 19.8 28.3 51.9   2 33.19 25.66 3.41    2

40  [0,5] 40  4.580 194   3 27.2 40.9 31.9   3 31.80 27.38 2.71    2

41  (15,20] 41 15.100 117   3 31.1 47.7 21.2   4 30.70 41.30 3.84    3

42  (10,15] 42 12.350 108   3 19.2 30.7 50.1   1 42.50 36.15 2.26    1

43  (10,15] 43 11.560  72   3 17.6 25.5 56.9   1 47.30 33.76 2.42    1

44  (10,15] 44 10.550  72   3 13.2 14.6 72.2   1 39.80 24.95 2.97    1

45  (15,20] 45 16.460 216   3 13.8 16.9 69.3   1 49.20 44.66 3.37    1 46  (5,10]     46 9.080  67   3  9.6  4.0 86.4   1 41.20 27.08 3.07    1

47  (10,15]     47 12.072 252   3 28.3 24.2 47.5   2 31.19 57.07 3.64    2

48  (10,15] 48 10.003 198   3 27.2 25.1 47.7   1 39.80 59.50 3.64    1

49  (5,10] 49  9.015 240   3 24.4 35.5 40.1   2 27.20 31.52 2.55    2 50 (10,15]     50 11.122 180   3 23.8 20.6 55.6   1 34.44 22.62 4.23    1

#transformation des données qualitatives en facteur

ZON<-as.factor(don$ZON)

VRTC<-as.factor(don$VRTC)

#création du tableau codé

doncd<-data.frame(RDT,PLM,ARG,LIM,SAB,PGM,MST,AZP,ZON,VRTC) row.names(doncd)<-don$Numero summary(doncd)

RDT    PLM    ARG    LIM    SAB    PGM    MST    AZP    ZON VRTC  

 1:17   1:17    1:17    1:17    1:17    1:17    1:17    1:17 1:17     1:24  

 2:16   2:16    2:16    2:16    2:16    2:16    2:16    2:16 2:15     2:21  

 3:17   3:17    3:17    3:17    3:17    3:17    3:17    3:17 3:18     3: 5  

#appel de la bibliothèque de fonction ade4

#Calcul AFCM

library(ade4)

#13 axes factoriels sont retenus 

z<(df = doncd, scannf = FALSE, nf = 13)

#calcul et plot des valeurs propres

round(z$eig/sum(z$eig)*100,2)

[1] 18.45 13.86 10.18  9.26  7.96  7.29  5.40  5.05  4.57  3.65 3.38  2.39

[13]  2.09  1.64  1.48  1.06  0.85  0.65  0.63  0.19

inertie<-z$eig/sum(z$eig)*100

round(inertie,2)

barplot(inertie,ylab="% d'inertie",=round(inertie,2)) title("Evolution des valeurs propres en %")

#plot plan des modalités

plot(z$co[,1],z$co[,2],type="n",xlab="Axe 1",ylab="Axe 2", xlim=c(-1.4,1.4)) text(z$co[,1], z$co[,2], label= rownames(z$co)) title("Blé dur-Plan des modalités")

abline(h=0,v=0)

# plan des individus

plot(z$li[,1],z$li[,2],type="n",xlab="Axe 1",ylab="Axe 2", xlim=c(-1.4,1.4)) text(z$li[,1], z$li[,2]) title("Blé dur - Plan des individus") abline(h=0,v=0)

#Contributions des variables à la construction des axes (z,col.inertia = T)$

       Comp1 Comp2 Comp3 Comp4 Comp5 Comp6 Comp7 Comp8 Comp9 Comp10 Comp11 Comp12 Comp13

RDT.1     52  1178   122   185 51   625   236   509    24     56     37     11    166

RDT.2    265   128   125    99    52 2097    73   960     7      1    161    211    551

RDT.3    528   544     0   542 0   377   558    56     5     69     39    301     97

PLM.1     76   966   305    10 515    36   340    20   604      7    166     82    555

PLM.2     29   188    80   155 678   636   157  2645     4    143    247    411    251

PLM.3     12   317   683   230 7   341    39  2067   707     79      6    112     67

ARG.1    878    13    25   643 170    23    39   145    22    456    562    847      2

ARG.2     12    35   251   945 317    50  1343   322    25    547   1992    387      2

ARG.3    687    89   413    20 18   137  1748    29     0      2    384    100      7 LIM.1    915    11 425   368     8    35    47     5     1    203    363    474    648

LIM.2      1   985   214   845 60   198    45   242    70     20    353    557    420

LIM.3    869   740    41    81 21    61     0   168    52     98      1      1     31 SAB.1   1001   726 10    59    90     0     5    12    21     38     18     47     53

SAB.2      4   937   205  1147 100    20   475    12     2    331     24     50      8

SAB.3   1121     8   116   635 0    19   355    47    37    131      0    188    103 PGM.1     87    17 129   221   673   959  1714   327   284    173     52    226      3

PGM.2     25    83   558    14 961   141   535   260  1062   1632    268    179    308 PGM.3     19    21 1175   342    17   379   360  1139   218    678     75      4    354 MST.1 186   815   180    87   141   763     3    17   313    434    256     51   1759

MST.2    132    12   825    63   388 1546   156    19    82    529     91    309    588

MST.3    614   632   209   290 52   111   108     0    79      2     45    586    339

AZP.1    264    20   226   252 993    12   100   154  1855    109    312   1109    138

AZP.2     41    28    57    17 1312   114     6   576  2992     18    911    266     14

AZP.3    503    93   498   141 13    46   155   118   100     40   2203   2414     66

ZON.1    803   783     1     8 6   243    16    25    12    291    183     62    352

ZON.2     25   367  1101    27 844   695   199     1    19    100     31    121    652

ZON.3    528    94   873     4 581    79   280    35     0     55     65    311     26

VRTC.1   123    52     3   819 712    11   456     4   615    153      0    155   1231

VRTC.2    30    96   270  1486 160    15   433     5   222    228    217    329   1157

VRTC.3   169    19   882   265 1059   231    17    82   566   3374    937     98     52

#Contributions des lignes à la construction des axes (z,row.inertia = T)$

   Axis1 Axis2 Axis3 Axis4 Axis5 Axis6 Axis7 Axis8 Axis9 Axis10 Axis11 Axis12 Axis13

1     151   444   159    53    18    88    12    90    18    758    388 435    330

2     511    25    67     2     0    58   702   364   157     98     20 179    181

3     51   774   136   154    12    13   148   142    34    102     68 142    208

4     173   555    76    49    44   113   662    23    20     16     57 213      5

5     98   434     9   146    35     3   262   165  1052    158      1 45    420

6     72   900     0   164   281     5   134    85    52     24 11      4     28

7     53   891     3   269    12    26   177     0   326     83      0 61      0

8     476     2    88   272    83   152   252   190    50    301 15      6    188

9     31   413   832    90   624   285     1   161     6    596 0      8     17

10    770     3    55   553   432    16    18    32    78    153 174      1     29

11    174   546    66    30   122     3   364    25   290    155     18 235    200

12    23    32    12     7   374   955   539   353    68     31 322      0      9

13    19   100    22   375    20   915    72   108   471    204    523 187    144

14    0     4     4     8   325   134   641   459   300    575    152 156    394

15    635     5   416     9     1    26    41   109    15     11 89     35    325

16    465    42   358    16    13     0    62   304    63      3    222 221    149

17    385     0     6   259    74   106     0   193   272    208    443 943    401

18    347    43   401   220     0     0    57   454   247     89 11      1    552

19    375   131    14   151   214   410   400     4   286      7 0      8    143

20    375   131    14   151   214   410   400     4   286      7 0      8    143

21    18   246   270   374     2    88    46    61    34    310     18 2452      7

22    6   164   169    59   537   162    81     1  1337    182      1 641    130

23    0    71   994    87    68   401     6   586   863    842      5 350      0

24    12   364  1030   108    77    59   202   313    25    494    445 319     17

25    82   551   232     9   146   337    99     1   104    218    357 100     99

26    108    30   195    82    16  1017   301     1    49    273    209 113    977

27    247    19     0   346   266    78     0   526     1     44    802 201     64

28    46     3   258   301   614   531   363    13   101     14    115 130    108

29    5   236    32   355   475     9   243   968     9     45     97 137    146

30    208     0   285    79   461   310     0    67    26    807 357     39    240

31    138   578    57   144     7   252   246   262    59      2 0     35    622

32    136   690     9    19   179    31   103   817     3      2 25      6     29

33    308    98   335    19    69    22   346   382     1    497 95      0    117

34    307     9   553    93    23     0    79    31    96      1    675 187    591

35    186     3   300    81    96    74   467     3   405    466    375 688    217

36    366    72   102     2   630   500    24    21    22      1 66      1      3

37    357    14     8     1   978    75     2    44     4    142 247     25    718

38    57   162    33   914   336     9    16   639   119    163 1      0     90

39    98     0   128     2    20   354   386   603   659    435     72 312     41

40    14    78     7   863   393     2   312   303     0      5    349 121     81

41    385    51    26   313   841    57    27   433   101    526    364 197     57

42    44   230    19   308   139   135   615    58   239     44    478 19    606

43    484     2    28   198    68   509     2    42   360    151 8     24      6

44    489    67   143    98    24    55    34    13   595    384 8     21    189

45    92   108    79   620    25    19   334   326   380    104     81 665      3

46    427    13    35   187   142   577     0   214    36      8 75     26    296

47    0   242   698   112    72    81   367     0    74      6    335 81    309

48    0   248   669     3    84   317   157     1   166    195    227 68      1

49    60   176     0  1241   161   144     6     2    13     46     40 110     41

50    132     0   566     2   154    77   187     2    30     13 1556     42    328

#pour mieux interpréter l'importance des modalités modal<.frame(z$co) modal<-modal[(modal$Comp1),] dotchart(modal[,1],labels = row.names(modal),cex=0.8) title(sub="Répartition des modalités sur l'axe 1") abline(v=0)

#dendrogramme

hier<-hclust(dist(z$li), method="ward") plot(hier,hang=-1)

# découpage en 3 classes et 5 classes cl3<-cutree(hier,3) cl5<-cutree(hier,5)

#fusion de tableau de données avec les variables à 3 et 5 groupes arvalis<-data.frame(don,cl3,cl5) (pat = "^a")

arvalis : 'data.frame': 50 obs. of  14 variables:

 $ Numero: int  1 2 3 4 5 6 7 8 9 10 ...

 $ RDT   : num   6.49 15.58  7.29  1.09  5.10 ...

 $ PLM   : int  84 112 68 88 174 63 92 117 58 116 ...

 $ ZON   : int  1 1 1 1 1 1 1 1 1 1 ...

 $ ARG   : num  21.5 21 26.2 29.7 22.8 19.6 26.7 34 16.7 37.4 ...

 $ LIM   : num  60.6 58.3 47.6 54.5 59 68 53.7 44.9 57.6 49 ...

 $ SAB   : num  17.9 20.7 26.2 15.8 18.2 12.4 19.6 21.1 25.7 13.6 ...

 $ VRT   : int  3 3 3 3 3 3 3 1 4 6 ...

 $ PGM   : num  43.1 38.3 45.3 29.1 42.8 ...

 $ MST   : num  34.5 39.2 26.9 23.1 18.1 ...

 $ AZP   : num  3.82 3.78 2.61 3.78 3.41 3.04 2.26 3.85 2.69 3.62 ...

 $ VRTC  : int  2 2 2 2 2 2 2 1 3 3 ...

 $ cl3   : int  1 1 1 1 1 1 1 1 2 1 ...

 $ cl5   : int  1 2 1 1 1 1 1 2 3 2 ...

#moyenne par groupe pour les variables quantitatives

moyennes5<-aggregate(arvalis[,c(-1,-12,-13,-14)],list(arvalis$cl5), mean) round(moyennes5,2)

Group.1   RDT    PLM  ZON   ARG   LIM   SAB  VRT   PGM   MST  AZP

1     1      5.17       92.71 1.00 27.24  55.36  17.40  2.71    39.89 23.75 3.22

2     2     16.33    146.00 1.38 32.99  46.89  20.12  2.50   37.92  41.85 3.72

3     3     11.59    132.86 1.93 26.59  38.14  35.98  2.29   38.35  34.46 3.10

4     4     10.48    132.12 2.82 17.98  22.56  59.46  1.41   38.33  31.57 3.19 5       5     1.32        84.75 2.00 16.68  31.93  51.40  1.00   30.42 17.82 2.90

#moyenne par groupe pour les variables quantitatives (sans les variables qualitatives) moyennes3<-aggregate(arvalis[,c(-1,-12,-13,-14)],list(arvalis$cl3), mean) round(moyennes3,2)

Group.1   RDT    PLM  ZON   ARG   LIM   SAB  VRT   PGM   MST  AZP

1     1     11.13   121.13 1.20   30.31   50.84 18.85 2.60   38.84   33.40 3.49

2     2    11.59    132.86 1.93   26.59   38.14 35.98 2.29   38.35   34.46 3.10

3     3    8.73      123.10 2.67   17.73   24.34 57.92 1.33   36.82 28.95 3.13

#comptage des individus par classe arvalis$cl3<-as.factor(arvalis$cl3) arvalis$cl5<-as.factor(arvalis$cl5)

summary(arvalis$cl3)

1    2    3 

15 14  21

summary(arvalis$cl5)

1   2   3  4   5 

7   8 14 17  4

# très joli résumé (faire « échappe » après apparition du graphique sur le clavier sinon le système se bloque)

par(mfrow=c(1,1)) library(cluster)

clusplot(dist(z$li),cl3,diss=T,shade=T,color=T,labels=5, span =T, main="représentation graphique des clusters")

abline(v=0,h=0)

par(mfrow=c(1,1)) library(cluster)

clusplot(dist(z$li),cl5,diss=T,shade=T,color=T,labels=5, span =T, main="représentation graphique des clusters")

abline(v=0,h=0)

# exportation du fichier de données avec les partitions 3 et 5 groupes récupérables avec sur excel

write.table(arvalis, file="arvalis", sep="\t",dec=",",row.names=F, col.names=T)

ANALYSE DE DONNEES SUPERVISEE

Analyse Factorielle Discriminante Linéaire

On génère 2 populations par la loi normale (X1 et X2) chacune de ces deux populations est décrite par deux variables

# théorie de l'AFD sous forme graphique

# d'après

(24122006)

library(MASS)

s <- matrix(c(1, 0.8, 0.8, 1), 2) x1 <- mvrnorm(50, c(0.3, -0.3), s) x2 <- mvrnorm(50, c(-0.3, 0.3), s) x <- .frame(x1, x2) x <- scalewt(x, scale = F) fac <- factor(rep(1:2, rep(50, 2)))

s.class(x, fac, col = c("red", "blue")) arrows(0, 0, 2, 0, lwd = 2) text(2, 0, "A", pos = 1, cex = 2) arrows(0, 0, sqrt(2), sqrt(2), lwd = 2) text(sqrt(2), sqrt(2), "B", pos = 4, cex = 2) arrows(0, 0, 0, 2, lwd = 2) text(0, 2, "C", pos = 2, cex = 2) arrows(0, 0, -sqrt(2), sqrt(2), lwd = 2)

text(-sqrt(2), sqrt(2), "D", pos = 2, cex = 2)

On représente le plan des variables centrées. Dans ce plan, choisir une direction et projeter les points sur l’axe ainsi défini, c’est faire une combinaison linéaire des variables centrées. Suivant la direction, les deux populations apparaîtront un peu, beaucoup ou pas du tout différentes.

par(mfrow = c(2, 2)) f1 <- function(a, b, cha) { z <- a * x[, 1] + b * x[, 2] z0 <- seq(-4, 4, le = 50) z1 <- z[fac == 1] z2 <- z[fac == 2]

hist(z, proba = TRUE, col = grey(0.9), xlim = c(-4, 4), main = cha) lines(z0, dnorm(z0, mean(z), sd(z)), lwd = 2)

lines(z0, 0.5 * dnorm(z0, mean(z1), sd(z1)), col = "red", lwd = 2) lines(z0, 0.5 * dnorm(z0, mean(z2), sd(z2)), col = "blue", lwd = 2)

} f1(1, 0, "A")

f1(1/sqrt(2), 1/sqrt(2), "B") f1(0, 1, "C")

f1(-1/sqrt(2), 1/sqrt(2), "D")

Cette figure indique aussi un élément essentiel. Entre les positions la capacité discriminante varie mais varie également la variabilité. On peut faire beaucoup de variance (B), beaucoup de variance inter-classe ou beaucoup de variance interclasse en proportion (D). On peut discriminer avec peu de variabilité (qui est alors d’abord inter-classe), on peut faire beaucoup de variance (c’est le propre de l’ACP) sans rien discriminer.

La position B est presque celle qui maximise la variance (elle est en fait l’axe principal théorique) : toute la variabilité est intra-population. La position D est presque celle qui minimise la variance (elle est en fait le second axe principal théorique). On y trouve quatre fois moins de variance mais la variabilité y est pour moitié inter-groupe.

f2 <- function(a, b) { z <- a * x[, 1] + b * x[, 2] a1 <- var(z) * 99/100

a2 <- var(predict(lm(z ~ fac))) * 99/100 a3 <- a2/a1 round(c(a1, a2, a3), 3) } f2(1, 0)

f2(1/sqrt(2), 1/sqrt(2))

f2(-1/sqrt(2), 1/sqrt(2))

w <- (x, scal = F, scannf = F)

w$eig

w$c1

wbet <- between(w, fac, scannf = F) wbet$eig

wbet$c1

wdis <- discrimin(w, fac, scannf = F)

wdis$fa/sqrt(sum(wdis$fa^2))

wdis$eig

Télécharger Cours complet d’initiation au logiciel R pas à pas

Ads