#Exercice 1
#a)
taille = 50
Mangue=rnorm(150/3,50,0.2)
Orange=rnorm(150/3,50,0.2)
Eau=rnorm(150/3,50,0.2)
col_mangue=rgb(red=244,green=187,blue=68,maxColorValue=255)
boxplot(Mangue,ylim=c(49,51),col=col_mangue)
plot(Mangue,pch=19,col=col_mangue,ylim=c(49,51),
cex=1,xaxt="n",xlab="")
boxplot(Orange,ylim=c(49,51),col='orange')
plot(Orange,pch=19,col='orange',ylim=c(49,51),
cex=1,xaxt="n",xlab="")
boxplot(Eau,ylim=c(49,51),col='cadetblue')
plot(Eau,pch=19,col='cadetblue',ylim=c(49,51),
cex=1,xaxt="n",xlab="")
boxplot(Mangue,Orange,Eau,
col=c(col_mangue,'orange','cadetblue'),
names=c("Mangue","Orange","Eau"),
ylim=c(49,51),
main="Distributions de Mangue, Orange, Eau")
#b)
Bouteille=Mangue+Orange+Eau
par(mfrow=c(1,3))
boxplot(Bouteille,ylim=c(148.5,151.5))
plot(Bouteille,ylim=c(148.5,151.5))
hist(Bouteille) #la distribution est gaussienne
sd(Bouteille) #estimation de l'écart-type
## [1] 0.3415017
#95% des valeurs sont entre 149,5 et 150,5 (écart de 1) donc 3*sigma donc un écart-type de 1/3
#c)
dev.off()
## null device
## 1
qqnorm(Bouteille) #nuage longiligne -> on valide l'hypothèse de normalité
#shapiro.test(Bouteille) # p-value > 5% (pas d'intérêt ici car qqnorm est bon)
#d)
plot(Bouteille,ylim=c(149,151),
col=ifelse(Bouteille<150.5&Bouteille>149.5,"darkgreen","red"),
xlab="",xaxt="n",pch=19,ylab="", #pch = type de points
main="Quantités de jus de fruits")
abline(h=mean(Bouteille),col="darkgreen",lty=3,lwd=2)
abline(h=149.5,col="blue",lty=3,lwd=2)
abline(h=150.5,col="blue",lty=3,lwd=2)
abline(h=150,col="blue",lty=3,lwd=2)
nb_def=sum(Bouteille>150.5)+sum(Bouteille<149.5) #on cherche cb de bouteille non conforme
propdef=nb_def/taille
propdef
## [1] 0.16
round(prop.test(x=nb_def,n=taille,)$conf.int,2)
## [1] 0.08 0.30
## attr(,"conf.level")
## [1] 0.95
barplot(propdef*100,ylim=c(0,100),col='black',main="proportion de bouteilles \n présentant un défaut de volume")
abline(h=2,col='red')
#var(bouteille)=var(mangue)+var(orange)+var(eau) = 3*0.2² -> sigma(bouteille)=0,2v3=0,34
#Exercice 2
#1)
#-> soit calcul de proba pour confirmer 95% ou calcul des intervalles mu-2sig et mu+2sigma
mu=75
sigma=0.3
mu-2*sigma
## [1] 74.4
mu+2*sigma
## [1] 75.6
pnorm(75.5,mean=mu,sd=sigma)-pnorm(74.5,mean=mu,sd=sigma)
## [1] 0.9044193
#-> 95% de bouteilles contiendront entre 74.4cl et 75.6cl
#2) taille de l'effet = dérive / écart-type. si taille proche de 0, écart petit par rapport à la variabilité, effet non significatif, inutile de faire un test.
#si ratio proche de 0.5, on peut faire un test. ratio >0.5 , effet fort. test sera puissant on peut minimiser la taille de l'éch
ecart=0.2
taille=ecart/sigma
taille
## [1] 0.6666667
#taille 0.66 effet assez fort -> ok pour test
#3) H0: moy=75 -> test de student
power.t.test(n=NULL, delta=ecart,sd=sigma, power=0.8,type="one.sample",alternative="two.sided")
##
## One-sample t test power calculation
##
## n = 19.66697
## delta = 0.2
## sd = 0.3
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
# 1-beta = proba de rejeter H0 quand il faut rejeter H0
#Exercice 3
tableau = matrix(c(39,52,77,61,48,23), nrow=2, byrow=TRUE)
colnames(tableau)=c("A","B","C")
rownames(tableau)=c("Bon","Très Bon")
tableau
## A B C
## Bon 39 52 77
## Très Bon 61 48 23
couleur = colorRampPalette(colors=c("pink","lightblue"))(3)
barplot(c(61,48,23),ylim=c(0,100), col=couleur,names.org=c("A","B","C"),main="Pourcentage de personnes jugeant \n les produits TB")
## Warning in plot.window(xlim, ylim, log = log, ...): "names.org" n'est pas un
## paramètre graphique
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "names.org" n'est pas un paramètre graphique
## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "names.org"
## n'est pas un paramètre graphique
barplot(tableau, col=c("pink","cadetblue"))
chisq.test(x=tableau)
##
## Pearson's Chi-squared test
##
## data: tableau
## X-squared = 30.276, df = 2, p-value = 2.665e-07
# pvalue<1/1000, H0 rejetée de manière hautement significative. ces produits n'ont pas été jugés de la même manière. l'appréciation dépend des produit
#comparaison deux à deux
tableau_AB=tableau[ ,-3]
tableau_AB
## A B
## Bon 39 52
## Très Bon 61 48
chisq.test(tableau_AB)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tableau_AB
## X-squared = 2.9035, df = 1, p-value = 0.08839
# on ne peut pas rejeter H0. on ne peut pas dire que les produits sont significativement différents
# on a pas réussi à mettre en évidence que A est meilleur que B mais pvalue proche de 0.05
tableau_BC=tableau[ ,-1]
tableau_BC
## B C
## Bon 52 77
## Très Bon 48 23
chisq.test(tableau_BC)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tableau_BC
## X-squared = 12.578, df = 1, p-value = 0.0003904
#on rejet H0. diff très hautement significative
#pas besoin de comparer AC car AB non significativement diff et BC significativement diff
tableau_AC=tableau[ ,-2]
tableau_AC
## A C
## Bon 39 77
## Très Bon 61 23
chisq.test(tableau_AC)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tableau_AC
## X-squared = 28.099, df = 1, p-value = 1.152e-07
couleur = colorRampPalette(colors=c("lightblue","cadetblue","orange"))(3)
barplot(c(61,48,23),ylim=c(0,100), col=couleur,names.org=c("A","B","C"),main="Pourcentage de personnes jugeant \n les produits TB")
## Warning in plot.window(xlim, ylim, log = log, ...): "names.org" n'est pas un
## paramètre graphique
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "names.org" n'est pas un paramètre graphique
## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "names.org"
## n'est pas un paramètre graphique
#Exercice 4
quantité = c(33.03,32.87,33.03,33.03,32.94,32.75,32.81,32.87,32.90,33.14,32.98,32.82,32.79,32.87,32.87,32.86,32.93,32.81,32.94,32.78)
qqnorm(quantité)
t.test(quantité,mu=33,alternative="less")
##
## One Sample t-test
##
## data: quantité
## t = -4.3638, df = 19, p-value = 0.0001671
## alternative hypothesis: true mean is less than 33
## 95 percent confidence interval:
## -Inf 32.94023
## sample estimates:
## mean of x
## 32.901
#Exercice 6
directe =c(27.85,28.35,28.45,27.9,27.45,27.4,28.15,27.65,28.35,28.15,27.95)
IR=c(27.9,28.4,28.4,27.95,27.4,27.5,28.25,27.5,28.5,28.35,27.7)
#représentation avec nuages de points car échantillons appariés
dif=directe-IR
plot(dif)
abline(h=c(0,mean(dif)))
boxplot(dif)
qqnorm(dif)
#shapiro.test(dif)
t.test(x=directe,y=IR,paired=TRUE,alternative="two.sided")
##
## Paired t-test
##
## data: directe and IR
## t = -0.45408, df = 10, p-value = 0.6595
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.10739940 0.07103576
## sample estimates:
## mean difference
## -0.01818182
t.test(dif,mu=0,alternative="two.sided")
##
## One Sample t-test
##
## data: dif
## t = -0.45408, df = 10, p-value = 0.6595
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.10739940 0.07103576
## sample estimates:
## mean of x
## -0.01818182
#p-value>0.5, on ne rejette pas H0.les résultats ne sont pas significativement dif, on n'a pas mettre en évidence une différences entre les deux.
#Exercice 7
mu=500
#valeurs comprises entre 495 et 505 g -> mu-3*sigma et mu+3*sigma (99,7% des échantillons compris dans cette plage)
#-> 3*sigma = 5 -> sigma = 5/3
sigma=5/3
ecart=2
#pour une doseuse
power.t.test(n=NULL, delta=ecart,sd=sigma, power=0.8,type="one.sample",alternative="two.sided")
##
## One-sample t test power calculation
##
## n = 7.587719
## delta = 2
## sd = 1.666667
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
#pour deux doseuses
power.t.test(n=NULL, delta=ecart,sd=sigma, power=0.8,type="two.sample",alternative="two.sided")
##
## Two-sample t test power calculation
##
## n = 11.94228
## delta = 2
## sd = 1.666667
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
dos1=c(499.8,500.4,500.1,500.8,499.9,500.5,501.1,499.3,498.7,500,499.8,499.5)
dos2=c(498.6,498.4,499.7,500.2,500,498.6,500.2,498.7,500.8,499.6,498.5,498.2)
boxplot(dos1,dos2,
col=c('lightgreen','darkgreen'),
names=c("Doseuse 1","Doseuse 2"),
ylim=c(497,502),
main="Distributions des quantités")
stripchart(list(dos1,dos2),
method="jitter",
vertical=TRUE,
add=TRUE,
pch=16)
#la doseuse 1 semble fonctionner conformement au cahier de charge
#la doseuse 2 semble sousdoser avec une grande variabilité
par(mfrow=c(1,2))
qqnorm(dos1)
qqnorm(dos2)
var.test(x=dos1,y=dos2,paired=FALSE,alternative="two.sided")
##
## F test to compare two variances
##
## data: dos1 and dos2
## F = 0.55629, num df = 11, denom df = 11, p-value = 0.3451
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1601425 1.9323733
## sample estimates:
## ratio of variances
## 0.5562869
var(dos1) #facultatif
## [1] 0.4353788
var(dos2)
## [1] 0.7826515
var(dos1)/var(dos2)
## [1] 0.5562869
#pvalue > 0.05 -> H0 est validée, les variances ne semblent pas significativement différentes, même variabilité
#homoscédaticité
t.test(x=dos1,y=dos2,var.equal=TRUE,alternative="two.sided")
##
## Two Sample t-test
##
## data: dos1 and dos2
## t = 2.1971, df = 22, p-value = 0.03884
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.03927522 1.36072478
## sample estimates:
## mean of x mean of y
## 499.9917 499.2917
#les moy sont significativement différentes
t.test(x=dos1,mu=500,alternative="two.sided")
##
## One Sample t-test
##
## data: dos1
## t = -0.04375, df = 11, p-value = 0.9659
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 499.5724 500.4109
## sample estimates:
## mean of x
## 499.9917
t.test(x=dos2,mu=500,alternative="two.sided")
##
## One Sample t-test
##
## data: dos2
## t = -2.7736, df = 11, p-value = 0.01811
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 498.7296 499.8538
## sample estimates:
## mean of x
## 499.2917
#mean(dos1)=499,99. non dif de 500. la doseuse 1 est conforme
#mean(dos2)=499,3. sousdosage -> faire estimation de la moy par IC
#Estimation de la variance
SCE_1=sum((dos1-mean(dos1))^2)
SCE_2=sum((dos2-mean(dos2))^2)
n1=length(dos1)
n2=length(dos2)
SCE=SCE_2+SCE_1
ddl=n1+n2-2
var_estimee=SCE/ddl
sigma_estimee=sqrt(var_estimee)
sigma_estimee
## [1] 0.7803942
((n1-1)*var(dos1)+(n2-1)*var(dos2))/(n1+n2-2)
## [1] 0.6090152
#Exercice 8
valeurs=c(25.2,24.3,26.8,25.9,22.1,23.8,21.9,22.6,18.4,19.5,18.9,19.9)
Dose=as.factor(c(rep("5ppm",4),rep("10ppm",4),rep("20ppm",4)))
data=data.frame(Dose,valeurs)
data
## Dose valeurs
## 1 5ppm 25.2
## 2 5ppm 24.3
## 3 5ppm 26.8
## 4 5ppm 25.9
## 5 10ppm 22.1
## 6 10ppm 23.8
## 7 10ppm 21.9
## 8 10ppm 22.6
## 9 20ppm 18.4
## 10 20ppm 19.5
## 11 20ppm 18.9
## 12 20ppm 19.9
str(data)
## 'data.frame': 12 obs. of 2 variables:
## $ Dose : Factor w/ 3 levels "10ppm","20ppm",..: 3 3 3 3 1 1 1 1 2 2 ...
## $ valeurs: num 25.2 24.3 26.8 25.9 22.1 23.8 21.9 22.6 18.4 19.5 ...
dev.off()
## null device
## 1
plot(valeurs~Dose,col="cadetblue")
# il semble qu'il ait un effet du facteur dose sur le brunissement. on peut conjecturer l'existence de trois groupes suivant le niveau de dose
# H0: mu1=mu2=mu3 ; Il n'y a pas d'effet du facteur dose
#Extraction des données pour les modalités 5ppm
D5=subset(valeurs,Dose=="5ppm")
moy_D5=mean(D5)
residus_D5=D5-moy_D5
residus_D5
## [1] -0.35 -1.25 1.25 0.35
#avec la fonction aov
anova=aov(valeurs~Dose)
residus=resid(anova)
residus
## 1 2 3 4 5
## -3.500000e-01 -1.250000e+00 1.250000e+00 3.500000e-01 -5.000000e-01
## 6 7 8 9 10
## 1.200000e+00 -7.000000e-01 8.604228e-16 -7.750000e-01 3.250000e-01
## 11 12
## -2.750000e-01 7.250000e-01
#anova= test unilatéral à droite. pas de parametre "alternative"
par(mfrow=c(2,2))
plot(anova)
#qq residuals : nuage longiligne, distribution des residus normale ou gaussienne, constant leverage: nuage à allure similaire, hypothese d'homoscedastidicité est validée
round(sum(residus),3)
## [1] 0
SCEres=sum(residus^2)
SCEres
## [1] 6.8575
SCEres/9 #estimation de la variance de l'essai
## [1] 0.7619444
SCEtot=sum((valeurs-mean(valeurs))^2)
tapply(valeurs,Dose,summary)
## $`10ppm`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 21.90 22.05 22.35 22.60 22.90 23.80
##
## $`20ppm`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.40 18.77 19.20 19.18 19.60 19.90
##
## $`5ppm`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 24.30 24.98 25.55 25.55 26.12 26.80
tapply(valeurs,Dose,mean)
## 10ppm 20ppm 5ppm
## 22.600 19.175 25.550
summary(anova)
## Df Sum Sq Mean Sq F value Pr(>F)
## Dose 2 81.43 40.72 53.44 1.01e-05 ***
## Residuals 9 6.86 0.76
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#rejet du H0 ; effet factoriel très hautement significatif
pairwise.t.test(valeurs,Dose,p.adjust.method="holm")
##
## Pairwise comparisons using t tests with pooled SD
##
## data: valeurs and Dose
##
## 10ppm 20ppm
## 20ppm 0.00071 -
## 5ppm 0.00100 8.2e-06
##
## P value adjustment method: holm
# donne des p-value. toutes les p-value sont < 5%. toutes les moyennes sont deux à deux différents
# si je veux augmenter le niveau de brun, je privilégie 5pmm
#exercice 9
# donnees = read.csv2("data_exo9.csv", sep=";",dec=",", header=TRUE, stringsAsFactors = TRUE)
#
# Doseuse=donnees$Doseuse
# Masse=donnees$Masse
#
# str(donnees)
#
# # 1) facteur étudié : doseuse , variable : masse
#
# plot(Masse~Doseuse,col=c("lightpink","lightblue","lightyellow"))
# points(Masse~Doseuse, col=c("blue"),pch=19)
#
# # boites de tailles différentes, on peut conjecturer que la doseuse 3 ne se comportent pas comme les autres et à une tendance à surdoser, les doseuses 1 et 2 semblent se comportent de manière similaire et permettent de répondre au cahier de charge (proche de 800 g)
# # ne pas oublier titre
#
# #ANOVA
# anova=aov(Masse~Doseuse)
# par(mfrow=c(2,2))
# plot(anova)
#
# # qqnorm: nuage longiligne -> distribution des résidus est gaussienne, residuals vs factor levels : 3 nuages à allures similaires -> on vérifie l'hypothèse d'homoscédasticité
#
# #normalité des résidus (facultatif)
# residus=resid(anova)
# shapiro.test(residus)
# residus
#
# #tableau anova
# summary(anova)
# summary(anova)[[1]][2,3] #carré moyen résiduel = 22.64 (estimation de la variance)
# ET=sqrt(22.64)
# ET
# # l'écart-type est trop grand car 95% des éch pèsent entre mu-3*sigma et mu+3*sigma
# #il y a un effet significatif du facteur doseuse, on peux considérer qu'au moins une des doseuses à un fonctionnement différent des autres
#
# #comparaison multiples
# sort(tapply(Masse,Doseuse,mean))
# pairwise.t.test(Masse,Doseuse,p.adjust.method="holm")
# # les doseuses 1 et 2 ont un fonctionnement similaire, elles forment un groupe. la doseuse 3 est différente des autres, elle met environ en moy 5 ou 6 g de trop
#
# dev.off()
# QN=800
# Comp_cible=tapply(Masse,Doseuse,mean)-QN
# Comp_cible
# barplot(Comp_cible,col=c("lightpink","lightblue","lightyellow"),main="Comparaison à QN",ylim=c(-6,6))
#
# # est ce que la doseuse 2 est loin de 800 -> test de student
# Masses_2=subset(donnees,Doseuse=="D2")$Masse
# t.test(Masses_2,mu=800,alternative="less")
# # la doseuse 2 ne sousdose pas