Elément analyse
- Tools
- importation de données
- Fonctions basiques
- Calculs et fonctions
- ACP
- AFC
- AFM (A venir )
- Aide
library(ggplot2)
library(magrittr)
library(rmarkdown)
library(class)
library(readr)
library(pwr)
library(readxl)
library(survival)
library(ISwR)
library(magrittr)
library(tidyr)
library(dplyr)
library(survminer)
library(survsim)
library(ade4)
library(rgl)
library(gplots)
library(adegraphics)
library(factoextra)Tools
importation de données
read_csv est plus recent que read.table, plus rapide et détecte les éléments du jeu de données.
# Pour tout ce qui est manipulation de repertoire (working dir)
# setwd()
# getwd()
# dir()
# Pour les fichier de type de csv
f="http://pbil.univ-lyon1.fr/R/donnees/CrimeStateDate.csv"
CSD=read_csv(f ) # remarque elle détecte les types
csd1=read.csv(f) # souvent meilleur, mais prend du temps sur les grosses base
head(CSD,n=3 # par defaut 6
)head(csd1,n=3) # Pour les fichier de types txt
f="http://pbil.univ-lyon1.fr/R/donnees/essence.txt"
data=read.table(f,header = T,dec = ',' # pour separer les colonnes
) # On peut aussi utiliser read_tsv# Importation Via un package
library(ade4)
data(banque)
head(banque)Fonctions basiques
# Information sur fonction
help("ggplot")
# Ou
?gglot
# trouver le minimum, le maximum, la moyenne
min() max() mean()
# Pour connnaitre la classe ou les types
class() #different de
typeof()
# Pour connaitre les structures de functions, de données
str( ) , glimpse( ) # package dplyr
# nombre de modalités / facteur
levels(csd1$Etat)
# nom des colonnes
names(CSD)
colnames(CSD)
# dim de jeu de donnée
dim(banque)
# Pour recupere la taille / nbre d'observations
length(levels(csd1$Etat))
length(banque$csp)La fonction quantile permet de repartition des individus compte tenu d’une distribution
quantile(CSD$Population)## 0% 25% 50% 75% 100%
## 226167 1111000 3115000 5439692 36132147# center, reduire
CSDN=scale(CSD$Population ) # ou scalewt()
summary(CSDN)## V1
## Min. :-0.8595
## 1st Qu.:-0.6840
## Median :-0.2865
## Mean : 0.0000
## 3rd Qu.: 0.1747
## Max. : 6.2632Calculs et fonctions
Fréquences relatives
\[f_k=\frac{n_k}{n}\]
summary(banque$csp)/length(banque$csp)## agric artis cadsu inter emplo ouvri
## 0.03580247 0.05925926 0.12716049 0.12592593 0.18641975 0.22592593
## retra inact etudi
## 0.06419753 0.10493827 0.07037037# Version arrondie
round(summary(banque$csp)/length(banque$csp),4)## agric artis cadsu inter emplo ouvri retra inact etudi
## 0.0358 0.0593 0.1272 0.1259 0.1864 0.2259 0.0642 0.1049 0.0704#ou
table(banque$csp)/length(banque$csp)##
## agric artis cadsu inter emplo ouvri
## 0.03580247 0.05925926 0.12716049 0.12592593 0.18641975 0.22592593
## retra inact etudi
## 0.06419753 0.10493827 0.07037037Variance écart-types
moyenne=function(x){
sum(x)/length(x)
}
variance=function(x){
v=0
for(i in 1:length(x)){
v=v+(x[i]-moyenne(x))^2
}
v/length(x)
}
ecart=function(x){
sd=variance(x)
return(sqrt(sd))
}
moyenne(data$Gazole)## [1] 1.157mean(data$Gazole)## [1] 1.157variance(data$Gazole) # Variance avec 1/n## [1] 0.007341var(data$Gazole) # Variance avec 1/(n-1) (version estimé)## [1] 0.008156667sd(data$Gazole) ## [1] 0.09031427ecart(data$Gazole)## [1] 0.08567964 # moyenne et variance pondérée
ip=c(1,3,7) # les indices à pondérer
ph=rep(1,10) # le vecteur des poids
ph[ip]=2 # pondération des indices
# moyenne pondérés
moyenneN=mean(data$Gazole*ph)
# Variance pondérée
var=variance(data$Gazole*ph)Calcul sur les Vecteurs
• produit scalaire
sum(data$Gazole*data$SP95)## [1] 15.4087• Norme
#norme
norme=function(x,poid=1){
sqrt(sum(poid*x*x))
}
norme(data$Gazole,ph)## [1] 4.180957norme(data$SP95,ph)## [1] 4.769319norme(data$SP98,ph)## [1] 5.035742# distance entre deux vecteur (ecart)
distance=function(x,y,omega) norme(x-y,omega)
distance(data$Gazole,data$SP95,ph)## [1] 0.6136774• Coefficient de projections
# coefficient de projection
prog=function(omega,x,y){
sum(omega*x*y)/sum(omega*x*x)
}
z=prog(rep(1/10,10),y=data$Gazole,x=rep(1,length(data$Gazole)))
z## [1] 1.157#norme de la difference
norme(data$Gazole-z,rep(1/10,length(data$Gazole)))## [1] 0.08567964• Cosinus d’angle
#
cosangle=function(poid,x,y){
sum(poid*x*y)/(norme(x,poid)*norme(y,poid))
}
# proche de 1
cos=cosangle(ph,data$Gazole,data$SP95)
cos## [1] 0.999237Relation entre variables
Quantitative
\[cov(X,Y)=\frac{1}{n}\sum^n_1(x_i-\bar{x})(y_i-\bar{y})\]
cov=function(X,Y){
cov=0
for(i in 1:length(X)){
cov=cov+ (X[i]-mean(X))*(Y[i]-mean(Y))
}
cov=cov/length(X)
cov
}
cov=cov(CSD$Population,CSD$Crime_Violent)
cov## [1] 175896310094• Coefficient de Corrélation (r) et détermination (r^2)
\[s_X , s_Y \rightarrow \] écart type resp. de X et Y
\[r=\frac{cov(X,Y)}{s_Xs_Y}\]
cov/(sd(CSD$Population)*sd(CSD$Crime_Violent) )## [1] 0.9012103# equivaut à
cor(CSD$Population,CSD$Crime_Violent) ## [1] 0.9015955# Le difference vient du fait que les écarts types du premier sont des écarts-type estimés Qualitative
• Calcul du khi Deux
Soit \(EO\rightarrow=n_{ij}\) effectif observer sur une table de contingence pour les modalités i et j resp des variables X,Y.
Soit \(ET \rightarrow \frac{n_i*n_j}{n}\)
la formules du khi-deux est la suivantes:
\[\chi^2=\sum\frac{(EO-ET)^2}{ET}\] Plus le khi-deux est grand plus les variables sont liés.
# Calcul du khi et test
chisq.test(banque$age,banque$interdit)##
## Pearson's Chi-squared test
##
## data: banque$age and banque$interdit
## X-squared = 11.423, df = 4, p-value = 0.0222Rejet de \(H_0\) dependance des variable age et interdiction.
Relation Croisée
On definit la somme des carrés Totaux comme étant: \[SCE_T=\sum^n_{i=1}(x_i-\bar{x})^2=n*\sigma^2\]
SCT=function(x){
length(x)*var(x)
}• Carré des écarts intergroupe. \[SCIG=\sum^p_{i=1}n_i(x_i-\bar{x})^2\]
$n_k $ effectif sous groupe \(k\), et \(p\) groupe.
• Rapport de corrélation
\[\eta^2=\frac{SCIG}{SCE_T}\] C’est l’équivalent du \(R^2\) ## visualisation
Histogramme
crime=CSD$Crime_Violent
# histogram (frequence)
hist(crime,col = "lightblue")
# format Densité
hist(crime,col = "lightblue",freq = F) 
# Avec un log
hist(log(crime),col = "lightblue",freq = F) 
hist(log(crime),col = "lightblue") 
Dotchart
par(mfrow=c(1,2)) # aligner les deux plots la logique est matricielle
dotchart(crime,main = "Titre", pch = 20 # type de point
)
dotchart(sort(crime), # la fonction sort range dans l'ordre croissant
main = "Titre", pch = 5 # type de point
)
Boxplot
# boxplot avec une sous partie
boxplot( sample(crime,300), col = "lightblue")
Pour trouver la seuil \(val_max=Q_3+1.5x(Q_3-Q_1)\) , intéressant en présence de valeur abérante.
quantile(crime)[[4]]+1.5*(quantile(crime)[[4]]-quantile(crime)[[2]])## [1] 61055\(Val_min=Q_1-1.5*(Q_3-Q1)\)
quantile(crime)[[2]]+1.5*(quantile(crime)[[4]]-quantile(crime)[[2]])## [1] 37747Diagramme camembert (pie chart)
pie(summary(banque$csp), main="titre")
Diagramme en Batons
barplot(table(banque$csp),main="Title", col="lightblue",las=2
# las pour la lecture des paramétres
)
Plot (x,y)
plot(x=CSD$Population,y=CSD$Crime_Violent,
main="Titre", ylab="Nom", xlab="Nom")
Balloplot (Croisement variable qualitative)
balloonplot(table(banque$age,banque$interdit), main="Titre")
mosiac plots (Croisement variable qualitative)
mosaicplot(table(banque$age,banque$interdit), main="Titre",shade = TRUE)
ACP
Notion clés Et élément
notion
- Les données doivent être quantitatives.
- Faire un peu de data visualization au préalable.
- S’il y a une trop grande différence entre les ordres de grandeurs (par visualisation boxplot), il faudra centrer et réduire mes données .
- Traitement les problémes de NA (missing value).
- L’objectif étant de sonder les strutures de la population par l’intermédiaire des variables.
- le caractére de l’acp est purement descriptif.
$$ Elle ne sert à valider quelconque hypothéses.
\(\rightaarow\) Le résultat est purement descriptif.
* Elle sert à la reduction de dimension quand on observe une certaine redondance au niveau des variables.
* L’analyse permet de ce fait de comprendre la relation entre les variables.
$$ Intuition, intensité,relation entre les variables.
Les individus sont associés à un poid (par défaut 1/N).
Il sagit aussi de rechercher la petite sous espace qui restituer l’image de N la meilleur façon (cf inertie, choix des axes)
Inertie
\[Inertie=\sum v(X^k)\] \(k\) étant le nombre d’axes.
- Le choix de l’axes se fait par le vecteur propre qui a la plus grande valeur.
Elément
Choix axes
Visualisation et analyse
# package ad4
acp <- dudi.pca(mesures, center=TRUE, scale=TRUE, scannf = F, nf = 3)
# affiche la fonction de l'on à fait
acp$call## dudi.pca(df = mesures, center = TRUE, scale = TRUE, scannf = F,
## nf = 3)# permet d'évaluer la fonction sortie
eval(acp$call)## Duality diagramm
## class: pca dudi
## $call: dudi.pca(df = mesures, center = TRUE, scale = TRUE, scannf = F,
## nf = 3)
##
## $nf: 3 axis-components saved
## $rank: 3
## eigen values: 1.737 0.9016 0.3613
## vector length mode content
## 1 $cw 3 numeric column weights
## 2 $lw 20 numeric row weights
## 3 $eig 3 numeric eigen values
##
## data.frame nrow ncol content
## 1 $tab 20 3 modified array
## 2 $li 20 3 row coordinates
## 3 $l1 20 3 row normed scores
## 4 $co 3 3 column coordinates
## 5 $c1 3 3 column normed scores
## other elements: cent norm# comparé au premier
identical(eval(acp$call), acp)## [1] TRUE# Avec summary() nous avons toutes ces informations
summary(acp)## Class: pca dudi
## Call: dudi.pca(df = mesures, center = TRUE, scale = TRUE, scannf = F,
## nf = 3)
##
## Total inertia: 3
##
## Eigenvalues:
## Ax1 Ax2 Ax3
## 1.7372 0.9016 0.3613
##
## Projected inertia (%):
## Ax1 Ax2 Ax3
## 57.91 30.05 12.04
##
## Cumulative projected inertia (%):
## Ax1 Ax1:2 Ax1:3
## 57.91 87.96 100.00# les éléménts de l'acp
acp ## Duality diagramm
## class: pca dudi
## $call: dudi.pca(df = mesures, center = TRUE, scale = TRUE, scannf = F,
## nf = 3)
##
## $nf: 3 axis-components saved
## $rank: 3
## eigen values: 1.737 0.9016 0.3613
## vector length mode content
## 1 $cw 3 numeric column weights
## 2 $lw 20 numeric row weights
## 3 $eig 3 numeric eigen values
##
## data.frame nrow ncol content
## 1 $tab 20 3 modified array
## 2 $li 20 3 row coordinates
## 3 $l1 20 3 row normed scores
## 4 $co 3 3 column coordinates
## 5 $c1 3 3 column normed scores
## other elements: cent norm# attach permet d'y acceder directement par leurs noms sans utilisr acp$
attach(acp)
eig ## [1] 1.7371612 0.9015541 0.3612847# poid des colonnes (par defaut 1)
acp$cw## [1] 1 1 1#poid des lignes (par defaut 1/n)
acp$lw## [1] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
## [15] 0.05 0.05 0.05 0.05 0.05 0.05# base *n
head(acp$lw)*nrow(mesures)## [1] 1 1 1 1 1 1# valeurs propres
acp$eig## [1] 1.7371612 0.9015541 0.3612847# somme des valeurs
sum(acp$eig)## [1] 3pve <- 100*acp$eig/sum(acp$eig)
pve## [1] 57.90537 30.05180 12.04282 #somme cumulative des pourcentages
cumsum(pve)## [1] 57.90537 87.95718 100.00000# donnne le rank de la matrice
acp$rank## [1] 3 # donne les coordonnées des variables (cosinus) version normé
acp$c1# calcule des coordonnée des variables sur le premier
acp$c1$CS1 * sqrt(acp$eig[1])## [1] 0.8676353 -0.4555425 -0.8813916acp$co$Comp1 ## [1] 0.8676353 -0.4555425 -0.8813916# coordonnée pour toute les autres variables
t(t(acp$c1) * sqrt(acp$eig))## CS1 CS2 CS3
## maladie 0.8676353 -0.2662627 -0.41989654
## dentaire -0.4555425 -0.8896976 0.03031894
## optique -0.8813916 0.1977284 -0.42901321 # donne les coordonnées des individus
acp$l1 # donne les coordonnées des individus( version normé à la racine carré de la valeur propre )
head(acp$li)acp$eig## [1] 1.7371612 0.9015541 0.3612847# lien entre L1 et li pour l'axes 1
head(acp$l1$RS1 * sqrt(acp$eig[1]))## [1] 0.1095310 -1.0926721 -2.3449799 -0.7392648 -1.2495023 -0.7687505# lien entre L1 et Li pour tous les axes
head(t(t(acp$l1) * sqrt(acp$eig)))## RS1 RS2 RS3
## 1 0.1095310 0.97223653 1.1502870
## 2 -1.0926721 -0.09929729 0.2840670
## 3 -2.3449799 -3.31729801 0.3766787
## 4 -0.7392648 1.07477518 0.4809507
## 5 -1.2495023 1.03732155 -0.2082590
## 6 -0.7687505 0.85057328 -0.4949950#Ce vecteur donne les moyennes (cent pour centrage)
acp$cent## maladie dentaire optique
## 116.950 49.395 29.885colMeans(mesures)## maladie dentaire optique
## 116.950 49.395 29.885# Pour calculer variance pondéré.
var.n <- function(x) sum((x-mean(x))^2)/length(x)
sd.n <- function(x) sqrt(var.n(x))
# Ce vecteur donne les ecarts-types par rapport à squrt(n)
acp$norm ## maladie dentaire optique
## 52.53436 53.20291 12.84330apply(mesures, 2, sd.n) ## maladie dentaire optique
## 52.53436 53.20291 12.84330# décentrage et déreduction (on retrouve les données initiales)
(mesures.cr*acp$norm )+acp$cent## maladie dentaire optique
## [1,] 64.200000 18.44374 12.58920
## [2,] 2.050059 142.43073 36.40000
## [3,] 18.626961 258.90000 148.91656
## [4,] 57.900000 20.54394 69.34097
## [5,] -2.102118 88.02312 45.30000
## [6,] 24.763265 19.80000 171.82286
## [7,] 82.200000 23.34421 35.37276
## [8,] 10.050594 119.42352 25.00000
## [9,] 25.447792 31.80000 180.82176
## [10,] 110.400000 27.90671 131.89219
## [11,] 41.850191 123.96572 40.70000
## [12,] 25.472240 47.70000 166.50532
## [13,] 110.400000 31.96227 36.61550
## [14,] 73.143425 132.85262 24.60000
## [15,] 34.762261 88.80000 111.28478
## [16,] 158.300000 27.35148 33.71577
## [17,] 88.638133 112.41274 21.10000
## [18,] 49.504058 8.90000 55.65520
## [19,] 277.100000 26.31346 -42.09134
## [20,] 92.790310 78.34627 0.00000
## attr(,"scaled:center")
## maladie dentaire optique
## 116.950 49.395 29.885
## attr(,"scaled:scale")
## maladie dentaire optique
## 52.53436 53.20291 12.84330# graph des individus
par(mfrow=c(1,3))
s.label(acp$li, xax = 1, yax = 2, clabel=1.5)## Warning: Unused parameters: clabel
s.label(acp$li, xax = 1, yax = 2,
label=as.character(depsante$age) # pour utiliser leur noms
, clabel=1.5)## Warning: Unused parameters: clabel
# s.class permets d'avoir des informations supplémentaire par rapport au variable
gcol <- c("red1", "red4","orange")
s.class(dfxy = acp$li, fac = depsante$groupe, col = gcol, xax = 1, yax = 2) 
scatter(acp,posieig="bottomright")
# pondération individus
acp$lw## [1] 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
## [15] 0.05 0.05 0.05 0.05 0.05 0.05poidsU <- acp$lw
# changement de pondération
poidsD <- depsante$effectif
poidsD <- poidsD/sum(poidsD)
round(poidsD,3)## [1] 0.059 0.058 0.067 0.083 0.050 0.054 0.036 0.032 0.040 0.063 0.096
## [12] 0.095 0.063 0.046 0.039 0.056 0.037 0.021 0.004 0.001 sum(poidsD)## [1] 1acpD <- dudi.pca(mesures, row.w=poidsD, # pour predefinir le poid
scannf=FALSE, nf=3)Version Propre
inertie
# inertie
fviz_eig(acp,addlabels = T)
Variable
# Cercle de corrélation (axes 1 et 2)
fviz_pca_var(acp,axes=1:2)
# Qualité represention variable (cos2)
d=fviz_cos2(acp,choice = "var",axes=1:2)
print(d)
d$data# table cos2
d$dataIndividus
# représentation individu
fviz_pca_ind(acp)
# Qualité representation individus
d=fviz_cos2(acp,choice="ind",axes=1:2)
print(d)
# table des somme de cos2
d$dataAFC
Notion clés Et élément
Visualisation et analyse
snee (TD5)
snee74 <- read.table("http://pbil.univ-lyon1.fr/R/donnees/snee74.txt",
header = TRUE)
names(snee74)## [1] "cheveux" "yeux" "sexe" cheveux <- snee74$cheveux
yeux <- snee74$yeux
(couleurs <- table(yeux, cheveux))## cheveux
## yeux Blond Marron Noir Roux
## Bleu 94 84 20 17
## Marron 7 119 68 26
## Noisette 10 54 15 14
## Vert 16 29 5 14(dfcouleurs <- data.frame(unclass(couleurs))) scorecheveux <- c(1, -1, -1, 1)
names(scorecheveux) <- colnames(couleurs)
scorecheveux## Blond Marron Noir Roux
## 1 -1 -1 1 #frequence
dfcouleurs <- data.frame(unclass(couleurs))
dfcouleurs["Bleu", ]/sum(dfcouleurs["Bleu", ]) # score moyen yeux
yeux.bleu <- dfcouleurs["Bleu", ]/sum(dfcouleurs["Bleu", ])
yeux.bleu * scorecheveux sum(yeux.bleu * scorecheveux)## [1] 0.03255814 freqyeux <- apply(dfcouleurs, 1, function(x) x/sum(x))
freqyeux## Bleu Marron Noisette Vert
## Blond 0.43720930 0.03181818 0.1075269 0.250000
## Marron 0.39069767 0.54090909 0.5806452 0.453125
## Noir 0.09302326 0.30909091 0.1612903 0.078125
## Roux 0.07906977 0.11818182 0.1505376 0.218750 t(freqyeux)## Blond Marron Noir Roux
## Bleu 0.43720930 0.3906977 0.09302326 0.07906977
## Marron 0.03181818 0.5409091 0.30909091 0.11818182
## Noisette 0.10752688 0.5806452 0.16129032 0.15053763
## Vert 0.25000000 0.4531250 0.07812500 0.21875000 scoreyeux <- apply(t(freqyeux), 1, function(x) sum(x * scorecheveux))
scoreyeux## Bleu Marron Noisette Vert
## 0.03255814 -0.70000000 -0.48387097 -0.06250000 ac <- dudi.coa(dfcouleurs, scannf = F, nf = 3)
rownames(ac$c1)## [1] "Blond" "Marron" "Noir" "Roux" #score optimal premier axes
ac$c1[, 1]## [1] 1.8282287 -0.3244635 -1.1042772 -0.2834725# score moyen yeux à partir de score optimaux cheveux
scoreyeux <- apply(t(freqyeux), 1, function(x) sum(x * ac$c1[,1]))
scoreyeux## Bleu Marron Noisette Vert
## 0.5474139 -0.4921577 -0.2125969 0.1617534 # meme
ac$li[,1] ## [1] 0.5474139 -0.4921577 -0.2125969 0.1617534sitpay(TD5)
sitpay <- matrix(c(209, 1483, 41, 320, 60, 34, 151, 1, 70, 10, 535,
2448, 33, 897, 135, 77, 245, 4, 139, 9), byrow = T, ncol = 5)
colnames(sitpay) <- c("celibataire", "concubin","divorcee", "marie", "veuf")
rownames(sitpay) <- c("annuel", "mensuel", "semestriel", "trimestriel")
sitpay## celibataire concubin divorcee marie veuf
## annuel 209 1483 41 320 60
## mensuel 34 151 1 70 10
## semestriel 535 2448 33 897 135
## trimestriel 77 245 4 139 9 n <- sum(sitpay)
I <- nrow(sitpay)
J <- ncol(sitpay)
freqsitpay <- sitpay/n
round(freqsitpay, digits = 4)## celibataire concubin divorcee marie veuf
## annuel 0.0303 0.2149 0.0059 0.0464 0.0087
## mensuel 0.0049 0.0219 0.0001 0.0101 0.0014
## semestriel 0.0775 0.3547 0.0048 0.1300 0.0196
## trimestriel 0.0112 0.0355 0.0006 0.0201 0.0013library(gplots)
#representation graphique de tableau de contigence
balloonplot(as.table(sitpay))
# plot de table de contigence
table.cont(sitpay, row.labels = rownames(sitpay), col.labels = colnames(sitpay),
csize = 2)
#mosaic plot
mosaicplot(sitpay, shade = TRUE)
profLignes <- prop.table(sitpay, 1 ) # definit frequence ligne
profLignes## celibataire concubin divorcee marie veuf
## annuel 0.0989115 0.7018457 0.019403691 0.1514434 0.02839565
## mensuel 0.1278195 0.5676692 0.003759398 0.2631579 0.03759398
## semestriel 0.1321640 0.6047431 0.008152174 0.2215909 0.03334980
## trimestriel 0.1624473 0.5168776 0.008438819 0.2932489 0.01898734profColonnes <- prop.table(sitpay, 2) # Definit frequence colonnes
profColonnes## celibataire concubin divorcee marie veuf
## annuel 0.24444444 0.34273168 0.51898734 0.22440393 0.28037383
## mensuel 0.03976608 0.03489716 0.01265823 0.04908836 0.04672897
## semestriel 0.62573099 0.56574994 0.41772152 0.62903226 0.63084112
## trimestriel 0.09005848 0.05662122 0.05063291 0.09747546 0.04205607# Test khideux
reschi <- chisq.test(sitpay)## Warning in chisq.test(sitpay): Chi-squared approximation may be incorrect reschi$expected ## celibataire concubin divorcee marie veuf
## annuel 261.79032 1324.8734 24.188813 436.62339 65.52413
## mensuel 32.95609 166.7848 3.045066 54.96537 8.24866
## semestriel 501.52731 2538.1388 46.339951 836.46544 125.52847
## trimestriel 58.72627 297.2030 5.426170 97.94580 14.69874 reschi##
## Pearson's Chi-squared test
##
## data: sitpay
## X-squared = 129.22, df = 12, p-value < 2.2e-16 #p-value trop petite on rejette h_0 (independance) contre hyp h_1
# on parle de liaison statistiquement significatifs
# effectifs théorique
outer(margin.table(sitpay, 1 # somme ligne
), margin.table(sitpay, 2 #somme colonnes
))/sum(sitpay) # total## celibataire concubin divorcee marie veuf
## annuel 261.79032 1324.8734 24.188813 436.62339 65.52413
## mensuel 32.95609 166.7848 3.045066 54.96537 8.24866
## semestriel 501.52731 2538.1388 46.339951 836.46544 125.52847
## trimestriel 58.72627 297.2030 5.426170 97.94580 14.69874 dfsitpay <- as.data.frame(sitpay)
afc <- dudi.coa(dfsitpay, scannf = F, nf = 3)
names(afc)## [1] "tab" "cw" "lw" "eig" "rank" "nf" "l1" "co" "li" "c1"
## [11] "call" "N" # tableau analysée
afc$tab # meme chose donc c'esdt le lien entre effectifs observée et effectifs théorique
(dfsitpay - reschi$expected)/reschi$expected# lien entre les pondérations cw de l'afc et les fréquences marginale
afc$cw ## celibataire concubin divorcee marie veuf
## 0.12389509 0.62701058 0.01144762 0.20663672 0.03101000apply(dfsitpay, 2, function(x) sum(x)/n)## celibataire concubin divorcee marie veuf
## 0.12389509 0.62701058 0.01144762 0.20663672 0.03101000afc$lw## annuel mensuel semestriel trimestriel
## 0.30618751 0.03854514 0.58658165 0.06868570apply(dfsitpay, 1, function(x) sum(x)/n)## annuel mensuel semestriel trimestriel
## 0.30618751 0.03854514 0.58658165 0.06868570afc$rank## [1] 3# InertieTOt=X^2/n
reschi$statistic/sum(sitpay)## X-squared
## 0.018725sum(afc$eig)## [1] 0.018725aides <- inertia.dudi(afc, row.inertia = TRUE, col.inertia = TRUE)
names(aides)## [1] "tot.inertia" "row.contrib" "row.abs" "row.rel" "row.cum"
## [6] "col.contrib" "col.abs" "col.rel" "col.cum" "nf"
## [11] "call" # calcule de pourcentag d'inertie
IT <- sum(afc$eig)
IT## [1] 0.018725afc$eig/IT## [1] 0.942755638 0.051058450 0.006185912aides$tot.inertiasitpay[3, ]## celibataire concubin divorcee marie veuf
## 535 2448 33 897 135# multiplié par mille pour facilité la lecture
aides$row.abs# calcule de contribution
afc$li[, 1] * afc$li[, 1] * afc$lw/afc$eig[1]## annuel mensuel semestriel trimestriel
## 0.59338040 0.05139243 0.10770248 0.24752470# contribution relatives
# multiplié par mille et le signes pour situer la modalité sur l'axe
aides$row.rel(afc$li[3, ] * afc$li[3, ])/(sum(afc$tab[3, ] * afc$tab[3, ] * afc$cw))temp <- (afc$tab[3, ] * sqrt(afc$cw)) * sqrt(afc$lw[3])
sum(temp * temp)/sum(afc$eig)## [1] 0.1112115 aides$row.cumEntfum (TD6)
entfum <- read.table("http://pbil.univ-lyon1.fr/R/donnees/entfum.txt",h=T,
row.names=1)
#entfum<-read_tsv("http://pbil.univ-lyon1.fr/R/donnees/entfum.txt")
entfumafc=dudi.coa(entfum,scannf = F)
barplot(afc$eig/sum(afc$eig))
moynat <- c(42,29,20,9)
moynat <- data.frame(t(moynat))
colnames(moynat) <- colnames(entfum)
sup=suprow(afc,moynat)
row.names(sup$lisup)='moyenne'
afc$li=rbind(afc$li,sup$lisup)
scatter(afc,posieig = "none")
Version Propre
casA <- read.table("http://pbil.univ-lyon1.fr/R/donnees/qualprodA.txt", h=TRUE, row.names=1)
apply(casA,1,sum)## nettoyage resistance rincage odeur ongles coins
## 4452 3985 3886 4236 3454 4149
## tenu rayure brillance recurage
## 4593 3568 2871 4607 apply(casA,2,sum)## P1 P2 P3 P4 P5 P8 P9 P10
## 3973 5537 4609 5702 3817 4123 5856 6184 afc1=dudi.coa(casA,scannf = F,nf=6)
# Test khi deux
chisq.test(casA)##
## Pearson's Chi-squared test
##
## data: casA
## X-squared = 1770, df = 63, p-value < 2.2e-16summary(afc1)## Class: coa dudi
## Call: dudi.coa(df = casA, scannf = F, nf = 6)
##
## Total inertia: 0.04447
##
## Eigenvalues:
## Ax1 Ax2 Ax3 Ax4 Ax5
## 0.0247031 0.0156372 0.0023617 0.0012057 0.0003968
##
## Projected inertia (%):
## Ax1 Ax2 Ax3 Ax4 Ax5
## 55.5469 35.1615 5.3104 2.7112 0.8923
##
## Cumulative projected inertia (%):
## Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
## 55.55 90.71 96.02 98.73 99.62
##
## (Only 5 dimensions (out of 7) are shown)# Inertie
fviz_screeplot(afc1, addlabels = TRUE, ylim = c(0, 70))
# plot variable modalité
fviz_ca(afc1,
col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
) ## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
fviz_ca_biplot(afc1,
map ="rowprincipal", arrow = c(TRUE, TRUE),
repel = TRUE)
fviz_ca_biplot(afc1,
map ="rowprincipal", arrow = c(TRUE, TRUE),
repel = TRUE,axes = c(2,3) )
data=matrix(c("B","M","V","M","F","M"),ncol = 2)
data=tbl_df(data)
colnames(data)=c("yeux","sexe")
data$yeux=as.factor(data$yeux)
data$sexe=as.factor(data$sexe)