# Analyse de la variance en Splus. La fonction read.table crée un "data
# frame".
> kenton_read.table("exemples/CH16TA01.DAT")
> kenton
V1 V2 V3
1 11 1 1
2 17 1 2
3 16 1 3
4 14 1 4
5 15 1 5
6 12 2 1
7 10 2 2
8 15 2 3
9 19 2 4
10 11 2 5
11 23 3 1
12 20 3 2
13 18 3 3
14 17 3 4
15 27 4 1
16 33 4 2
17 22 4 3
18 26 4 4
19 28 4 5
# Comme nous n'avons pas besoin de la dernière colonne, on ne conserve
# que les deux premières et nous donnons un nom plus significatif que Vi à
# chacune des colonnes (le [[2]] plutôt que [[1]] indique colonne au lieu
# de rangées.
> kenton_kenton[,1:2]
> dimnames(kenton)[[2]]_c("Caisse","Design")
> kenton
Caisse Design
1 11 1
2 17 1
3 16 1
4 14 1
5 15 1
6 12 2
7 10 2
8 15 2
9 19 2
10 11 2
11 23 3
12 20 3
13 18 3
14 17 3
15 27 4
16 33 4
17 22 4
18 26 4
19 28 4
# La variable Design est considérée comme une variable continue et non pas
# une variable qualitative (factor). Pour bien vérifier que c'est le cas,
# on utilise la fonction is.factor. Puis on la force à devenir qualitative
# avec la fonction as.factor. Si on ne s'assure pas que la variable est
# qualitative, on obtiendra une régression de Y sur X (avec 1 degré de liberté)
# plutôt qu'un plan à classification simple pour quatre groupes (3 degrés de
# liberté).
> is.factor(kenton[,2])
[1] F
> kenton[,2]_as.factor(kenton[,2])
> kenton
Caisse Design
1 11 1
2 17 1
3 16 1
4 14 1
5 15 1
6 12 2
7 10 2
8 15 2
9 19 2
10 11 2
11 23 3
12 20 3
13 18 3
14 17 3
15 27 4
16 33 4
17 22 4
18 26 4
19 28 4
# On modélise la variable Caisse par la variable Design à l'aide de la
# fonction aov (analysis of variance). Le jeu de données est dans le data
# frame kenton. Pour voir le tableau anova, on utilise la fonction summary.
> kenton.aov_aov(Caisse ~ Design, data=kenton)
> summary(kenton.aov)
Df Sum of Sq Mean Sq F Value Pr(F)
Design 3 588.2211 196.0737 18.59106 2.584961e-05
Residuals 15 158.2000 10.5467
# Le résultat de la fonction summary est un objet à partir duquel on peut
# accéder aux différents éléments. En particulier, on peut accéder à la
# colonne "F value" avec le signe de dollar ($). Comme la seconde composante
# est manquante (NA), on obtient la valeur de la statistique F en ne prenant
# que la première composante.
> summary(kenton.aov)$F
[1] 18.59106 NA
> summary(kenton.aov)$F[1]
[1] 18.59106
####################################################################
# Simulations en Splus pour étudier la robustesse du niveau du test F
# lorsque les conditions de normalité ou d'égalité de variance ne sont
# pas satisfaites.
#####################################################################
# Cas 1000 échantillons de taille 30 d'une t à 4 degrés de liberté
# On génère une matrice de 30000 variables aléatoires dans une matrice de
# 30 lignes (et donc 1000 colonnes)
> temp_matrix(rt(30000,4),nrow=30)
# On utilise la fonction apply. Celle-ci "applique" une fonction à un
# tableau à k dimesions, dans ce cas-ci une matrice (donc k=2). Elle est
# appliquée à la deuxième dimension (donc les colonnes), c'est donc dire
# qu'on fait comme une boucle où chaque colonne devient tour à tour le jeu
# de données. La fonction qui est appliquée est définie explicitement.
# Elle prend un argument (x). On fait une analyse de variance où la variable
# indépendante est x et le facteur est le vecteur qui contient la valeur 1
# 10 fois (rep(1,10)), puis la valeur 2, 10 fois et la valeur 3, 10 fois.
# C'est donc comme si on avait trois groupes de 10 observations chacun. Il
# ne faut toutefois pas oublier de spécifier que ce vecteur doit être traité
# comme un facteur (variable qualitative) à l'aide de la fonction factor.
# À l'objet de type analyse de variance qui est créé, nous appliquons la
# fonction summary. Comme nous sommes intéressés qu'à la statistique F,
# nous allons chercher le vecteur F ($F) du sommaire. La statistique F est
# la première composante de ce vecteur ($F[1]), les autres étant des NA.
# Le résultat de la fonction apply est donc le vecteur tempo de longueur 1000
# (le nombre de colonnes de temp). Il contient la statistique F des 1000
# jeux de données.
> tempo_apply(temp,2,function(x){summary(aov(x~
+ factor(c(rep(1,10),rep(2,10),rep(3,10)))))$F[1]})
# Comme il y a trois groupes de 10 observations chacun, les degrés de liberté
# sont 2 et 27. La commande tempo>qf(.95,2,27) est un vecteur logique qui
# dit si la statistique F dépasse la borne critique pour un test de niveau
# 5%. En prenant la somme de ce vecteur, on obtient le nombre de jeux de
# données pour lesquels on rejette H0. En divisant par la longueur du vecteur
# on obtient la proportion des jeux de données pour lesquels on a rejeté H0.
# Ceci est un estimé du niveau du test et doit être comparé à 5%.
> sum(tempo>qf(.95,2,27))/length(tempo)
[1] 0.044
# On recommence avec des données exponentielle de moyenne 1.
> temp_matrix(rexp(30000),nrow=30)
> tempo_apply(temp,2,function(x){summary(aov(x~
+ factor(c(rep(1,10),rep(2,10),rep(3,10)))))$F[1]})
> sum(tempo>qf(.95,2,27))/length(tempo)
[1] 0.035
# Cette fois-ci, on considère des échantillons d'observations normales, mais
# avec des écart type de 1, 1.5 et 3, respectivement, soit des variances de 1,
# 2.25 et 9. Afin de s'assurer que les 10 premières rangées ont un écart type
# de 1, les 10 suivantes , un écart type de 1.5, etc., il faut s'assurer que
# la matrice construite à partir du vecteur de longueur 30000 soit formée
# rangée par rangée, plutôt que colonne par colonne. Ceci est assuré par
# l'option byrow=T.
> temp_matrix(c(rnorm(10000),1.5*rnorm(10000),3*rnorm(10000)),nrow=30,byrow=T)
> tempo_apply(temp,2,function(x){summary(aov(x~
+ factor(c(rep(1,10),rep(2,10),rep(3,10)))))$F[1]})
> sum(tempo>qf(.95,2,27))/length(tempo)
[1] 0.06
# Même chose, mais avec des variances de 1, 2 et 3 (plutôt que des écart type
# de 1, 2, 3).
> temp_matrix(c(rnorm(10000),sqrt(2)*rnorm(10000),sqrt(3)*rnorm(10000)),
+ nrow=30,byrow=T)
> tempo_apply(temp,2,function(x){summary(aov(x~
+ factor(c(rep(1,10),rep(2,10),rep(3,10)))))$F[1]})
> sum(tempo>qf(.95,2,27))/length(tempo)
[1] 0.052
> temp_matrix(c(rexp(10000),sqrt(2)*rexp(10000)-(sqrt(2)-1),
+ sqrt(3)*rexp(10000)-(sqrt(3)-1)),nrow=30,byrow=T)
> tempo_apply(temp,2,function(x){summary(aov(x~factor(c(rep(1,10),rep(2,10),rep(Interrupt
Interrupt>
> tempo_apply(temp,2,function(x){summary(aov(x~
+ factor(c(rep(1,10),rep(2,10),rep(3,10)))))$F[1]})
> sum(tempo>qf(.95,2,27))/length(tempo)
[1] 0.059