Pour ces différents TP on pourra se référer aux exemples traités sur cette page
Ce premier TP est corrigé, afin de comprendre quelles questions il est important de se poser.
Le niveau annuel des précipitations dans la Beauce suit une loi Normale de moyenne \(\mu=700 \; mm\) et d’écart type \(\sigma=100 \; mm\)
Des industriels prétendent pouvoir augmenter le niveau moyen des pluies de manière significative sans changer l’écart type grâce à un procédé d’insémination de nuages à l’iodure d’argent.
Au cours d’une période de test de 9 ans, le niveau moyen des pluies pour chacune des 9 années a été le suivant (mm) :
Yobs <- c(610, 714, 880, 612, 601, 634, 703, 888, 750)Le procédé est-il efficace ?
Voici les questions qu’il faut se poser :
Que cherche t-on à démontrer ? Du point de vue de l’industriel, on veux prouver que le niveau moyen des pluies M sur la période de test est supérieur à 700 mm
Quel est le critère étudié ? Le critère étudié est donc la moyenne \(\mu\)
Quelle est votre hypothèse H0 ? H0 : le procédé est inefficace et \(\mu\) = 700
H1 est-elle unilatérale ou bilatérale ? Unilatérale H1 : le procédé est efficace et \(\mu > 700\)
Procédons au test Z (car \(\sigma\) est connu)
library(TeachingDemos) #librairie pour effectuer un test Z
z.test(Yobs,mu = 700,stdev = 100,alternative = "greater")##
## One Sample z-test
##
## data: Yobs
## z = 0.30667, n = 9.000, Std. Dev. = 100.000, Std. Dev. of the sample
## mean = 33.333, p-value = 0.3795
## alternative hypothesis: true mean is greater than 700
## 95 percent confidence interval:
## 655.3938 Inf
## sample estimates:
## mean of Yobs
## 710.2222La p-value aurait pu se trouver par :
1-pnorm(mean(Yobs), mean=700, sd=33)## [1] 0.3783701Il y a 38% de chance d’observer un tel niveau annuel des précipitations sous H0. On ne rejette pas H0 car l’écart par rapport à la valeur théorique \(\mu\) = 700mm sous H0, qui est dans le bon sens, n’est cependant pas suffisant pour conclure à l’efficacité du procédé. Autrement dit, l’écart observé ne peut être attribué qu’aux seules fluctuations aléatoires.
Dans ces circonstances, nous n’avons pas assez de preuves pour rejeter H0 et donc nous ne sommes pas soumis au risque \(\alpha\) de première espèce mais au risque \(\beta\) de deuxième espèce ou risque de ne pas croire le procédé efficace alors qu’il l’est en réalité mais que le bruit lié au hasard ne permet pas de déceler cet effet.
En fait nous aurons rejeté H0 avec le risque \(\alpha\) si
\(Pr(\overline{X} > \mu + h)= 0.05\) ce qui donne :
h<-qnorm(0.95)*100/sqrt(length(Yobs))Il aurait donc fallu trouver une moyenne au moins égale à :
700+h## [1] 754.8285La probabilité de voir un effet significatif est très peu probable. Car si nous supposons qu’effectivement le niveau des pluies a augmenté avec une moyenne de \(\mu_{1}=710\) alors la probabilité de rejeter H0 à raison est \(1-\beta\) la puissance du test :
1-pnorm(754.8285,mean = 710,sd = 100/sqrt(length(Yobs)))## [1] 0.08933603Soit moins de 1 chance sur 10…
En fait nous aurions pu très bien faire un test t (\(\sigma\) inconnu) et faire la même conclusion.
t.test(x =Yobs,alternative = "g",mu = 700 )##
## One Sample t-test
##
## data: Yobs
## t = 0.27497, df = 8, p-value = 0.3952
## alternative hypothesis: true mean is greater than 700
## 95 percent confidence interval:
## 641.0911 Inf
## sample estimates:
## mean of x
## 710.2222La statistique et la p-value (toujours comprendre d’où viennent les résultats !) sont obtenues par :
q<-(mean(Yobs)-700)*sqrt(length(Yobs))/sd(Yobs) #la statistique du t.test
p_value<-pt(q = q,lower.tail = F,df = length(Yobs)-1)
print(c(q,p_value))## [1] 0.2749661 0.3951550Si on refait la même chose en supposant que la taille de l’échantillon est de 1000 avec la même moyenne observée 710.
On aura une p-value de :
1-pnorm(710, mean=700, sd=100/sqrt(1000))## [1] 0.0007827011Nous rejetterons alors H0 et nous aurions une puissance de test de :
1-pnorm((700+qnorm(0.95)*100/sqrt(1000)), mean=710, sd=100/sqrt(1000))## [1] 0.9354202Si on s’en donne les moyens, même une petite différence de 10 mm devient hautement significative et se pose alors la question de l’effet réel pour l’agriculteur…
Une société s’approvisionne en pièces brutes, qui conformément en cahier des charges doivent avoir une masse moyenne de 782gr. Au moment où 500 pièces sont réceptionnées, on en prélève au hasard 50 dont on mesure les masses qui sont enregistrées ici :
et qu’on peut intégrer directement dans R via la commande :
donnees<-read.csv2(file ="https://sjaubert.github.io/SPCR/spc/TP2_TEST.csv",dec = ",",header = T )A partir de cet échantillon, effectuer un test d’hypothèse pour déterminer si on peut considérer que les 500 pièces ont une masse moyenne de 782gr.
Correction :
L’écart-type \(\sigma\) étant inconnu nous effectuerons un test t
Vérifions nos données :
str(donnees)## 'data.frame': 50 obs. of 2 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ mi: num 783 794 784 777 783 ...Nous avons 50 observations de deux variables.
X est l’indice de l’observation
mi est la variable étudiée, elle est bien de type numérique.
Nous posons pour hypohtèse nulle \(H_{0} : \mu=782\)
L’hypothèse alternative est alors \(H_{1} : \mu \neq 782\)
Sous \(H_{0}\), effectuons le test :
t.test(x = donnees$mi,mu = 782,conf.level = 0.95) #le niveau de confiance est par défaut 95% mais je l'ai rajouté pour des raisons pédagogiques##
## One Sample t-test
##
## data: donnees$mi
## t = -2.1655, df = 49, p-value = 0.03525
## alternative hypothesis: true mean is not equal to 782
## 95 percent confidence interval:
## 776.6066 781.7986
## sample estimates:
## mean of x
## 779.2026t=-2.1655 est la statistique du test, que nous pouvions obtenir par :
(t<-(mean(donnees$mi)-782)*sqrt(length(donnees$mi))/sd(donnees$mi))## [1] -2.165458df = 49 est le nombre de degrés de liberté (degree of freedom) ici 50 -1
La p-value = 0.03525 < 5%, probabilité d’observer une valeur aussi extrême, se détermine par :
pt(q = t,df = 49,lower.tail = T)*2## [1] 0.03524883Nous pouvons donc rejeter l’hypothèse \(H_{0}\) au risque de 5%
Soient deux ateliers A et B comprenant respectivement 50 et 40 personnes.
Nous trouverons dans ce fichier Données TP3
Deux listes de nombres na et nb
na (resp. nb) représente le nombre de jours d’arrêt de travail pour chaque personne de l’atelier A sur deux ans (resp. de l’atelier B).
En moyenne peut-on considérer qu’il y a une différence significative entre les ateliers A et B ?
Correction :
Nous allons effectuer un Test d’homogénéité, en posant comme hypothèse nulle \(H_{0}\) qu’il n’y a pas de différence entre les deux ateliers , c’est à dire quand moyenne, la moyenne des différences des arrêts de travail entre les deux ateliers est nulle.
Allons chercher les données :
donnees_TP3<-read.csv2(file ="https://sjaubert.github.io/SPCR/spc/TP3_TEST.csv",dec = ",",header = T )Examinons l’objet donnees_TP3 :
str(donnees_TP3)## 'data.frame': 50 obs. of 3 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ na: int 15 12 12 10 5 9 10 6 0 9 ...
## $ nb: int 10 10 9 2 3 0 3 10 10 0 ...Comme prévu nous obtenons les deux variables na et nb, nombre de jours d’arrêt de travail pour chaque personne de l’atelier A et B sur deux ans.
Pour simplifier renommons et examinons plus en détail :
na<-donnees_TP3$na
na## [1] 15 12 12 10 5 9 10 6 0 9 13 9 3 11 7 4 12 6 12 10 10 14 10 9 11
## [26] 14 0 0 15 9 6 4 12 1 14 12 13 6 4 15 12 2 6 4 6 14 15 12 5 13nb<-donnees_TP3$nb
nb## [1] 10 10 9 2 3 0 3 10 10 0 13 0 11 9 6 8 6 2 10 12 8 6 8 12 10
## [26] 0 10 9 2 4 1 11 3 13 5 11 5 1 12 5 NA NA NA NA NA NA NA NA NA NAPour éviter tout problème, prenons que les 40 premières valeurs (en fait cela s’avère inutile, mais a priori nous n’en savions rien…) :
nb<-nb[1:40]
nb## [1] 10 10 9 2 3 0 3 10 10 0 13 0 11 9 6 8 6 2 10 12 8 6 8 12 10
## [26] 0 10 9 2 4 1 11 3 13 5 11 5 1 12 5Effectuons le t test d’homogeneité :
t.test(na,nb,mu = 0,alternative = "two.sided")#les dernièrs item ne sont pas utiles##
## Welch Two Sample t-test
##
## data: na and nb
## t = 2.3154, df = 85.442, p-value = 0.02299
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.2982442 3.9217558
## sample estimates:
## mean of x mean of y
## 8.86 6.75La probabilité sous \(H_{0}\) de faire une telle observation est p-value = 0.02299<5%, nous rejeterons donc \(H_{0}\) au risque de 5%
Résultats que nous obtenons avec les lignes suivantes :
(t<-(mean(na)-mean(nb))/sqrt(sd(na)^2/50+sd(nb)^2/40))#la statistique du test## [1] 2.315396(1-pt(q = t,df = 50+40-2))*2#p-value légèrement différente que celle que nous obtenons avec le t.test, cela est dû à des corrections sur les estimations des écarts types population...## [1] 0.0229173Les principaux motifs d’arrêt de travail dans 3 sites différents d’une même entreprise sont :
Nous avons le tableau de contingence suivant :
## A B C D
## Site 1 54 14 6 51
## Site 2 45 14 8 31
## Site 3 33 34 12 33La distribution des arrêts de travail dépend elle du site concerné ?
Correction :
Nous allons effectuer un test du \(\chi^{2}\) d’homogénéité, le \(\chi^{2}\) d’homogénéité est un \(\chi^{2}\) d’indépendance. Il est seulement réalisé dans un but différent.
Le \(\chi^{2}\) d’homogénéité permet de vérifier que les répartitions de différents effectifs sont équivalentes.
Ce test repose ainsi sur 2 hypothèses :
\(H_{0}\) : il n’y a pas de différence significative dans la répartition des groupes étudiés
\(H_{1}\) : il y a une différence - cette hypothèse est a affiner en fonction du cas étudié
L’hypothèse \(H_{0}\) est donc : que les causes d’arrêts de travail sont indépendantes des sites (indépendance entre les lignes et les colonnes)
chisq.test(mat)##
## Pearson's Chi-squared test
##
## data: mat
## X-squared = 22.564, df = 6, p-value = 0.0009564Retrouvons et commentons ces résultats :
Calcul de la matrice théorique
V_th<-c()
for (i in 1:3){
for (j in 1:4){
V_th<-c(V_th,sum(mat[i,])*sum(mat[,j])/sum(mat))}}
(mat_th<-matrix(V_th,nrow = 3,byrow = T))## [,1] [,2] [,3] [,4]
## [1,] 49.25373 23.13433 9.701493 42.91045
## [2,] 38.61493 18.13731 7.605970 33.64179
## [3,] 44.13134 20.72836 8.692537 38.44776Calcul de la statistique Chi²
X_squared<-sum((c(54,14,6,51,45,14,8,31,33,34,12,33)-V_th)^2/V_th)
X_squared## [1] 22.56412Calcul de la p-value, probabilité d’observer la valeur du chi² sous \(H_{0}\) :
1-pchisq(q = X_squared,df = 6)## [1] 0.0009564136Il est bien évident que nous rejetons \(H_{0}\) au risque de 5%, appartenir à un site aura une influence sur la cause de l’arrêt de travail.