Statistiques en sciences humaines et sociales avec R, Jean-Herman Guay

Chapitre 8 Bâtir des fonctions en explorant l'inférence

Ce chapitre explore les principes de l'inférence pratiquement: 1) par le calcul de la marge d'erreur (CR47); 2) par le calcul de la taille de l'échantillon pour une marge d'erreur (CR48); 3) par la fabrication de tirages aléatoires et la construction d'un histogramme (CR52 et CR54) «reproduisant» la loi des grands nombres. On y trouve aussi des codes permettant l'utilisation des tests de signification (CR57 à CR59) puis l'utilisation des tests de puissance avec le package pwr, principe développé par Cohen (1988) et travaillé par Stéphane Champely.

#CR48


par(mfrow=c(2,2), oma=c(2,2,2,2),mar=c(1.9,1.5,3,1.5))

z=1
x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l",lwd=1, main="Courbe normale
 Z=1", xlab="", ylab="")
x=seq(-z,z,length=200)
y=dnorm(x,mean=0,sd=1)
polygon(c(-z,x,z),c(0,y,0),col="gray")
abline(v=-z, lty="dotted")
abline(v=z, lty="dotted")
abline(h=0)
text(0,0.2,"68 %", cex=1.6)
text(-2.5,0.2,"16% %")
text(2.5,0.2,"16 %")
arrows(-2.5, 0.18, -2.5, 0.015, lty=1)
arrows(2.5, 0.18, 2.5, 0.015, lty=1)
text(1.2,0.3, srt=90, adj = 0, labels = z, cex=.9)
text(-1.2,0.3, srt=90, adj = 0, labels = -z, cex=.9)


z=1.64
x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l",lwd=1, main="Courbe normale
 Z=1.64", xlab="", ylab="")
x=seq(-z,z,length=200)
y=dnorm(x,mean=0,sd=1)
polygon(c(-z,x,z),c(0,y,0),col="gray")
abline(v=-z, lty="dotted")
abline(v=z, lty="dotted")
abline(h=0)
text(0,0.2,"90 %", cex=1.6)
text(-2.5,0.2,"5 % %")
text(2.5,0.2,"5 %")
arrows(-2.5, 0.18, -2.5, 0.015, lty=1)
arrows(2.5, 0.18, 2.5, 0.015, lty=1)
text(1.9,0.3, srt=90, adj = 0, labels = z, cex=.9)
text(-1.9,0.3, srt=90, adj = 0, labels = -z, cex=.9)



z=1.96
x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l",lwd=1, main="Courbe normale
 Z=1.96",xlab="", ylab="")
x=seq(-z,z,length=200)
y=dnorm(x,mean=0,sd=1)
polygon(c(-z,x,z),c(0,y,0),col="gray")
abline(v=-z, lty="dotted")
abline(v=z, lty="dotted")
abline(h=0)

text(0,0.2,"95 %", cex=1.6)
text(-2.5,0.2,"2.5 %")
text(2.5,0.2,"2.5 %")
arrows(-2.5, 0.18, -2.5, 0.015, lty=1)
arrows(2.5, 0.18, 2.5, 0.015, lty=1)
text(2.2,0.3, srt=90, adj = 0, labels = z, cex=.9)
text(-2.2,0.3, srt=90, adj = 0, labels = -z, cex=.9)


z=2.58
x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l",lwd=1, main="Courbe normale
 Z=2.58", xlab="", ylab="")
x=seq(-z,z,length=200)
y=dnorm(x,mean=0,sd=1)
polygon(c(-z,x,z),c(0,y,0),col="gray")
abline(v=-z, lty="dotted")
abline(v=z, lty="dotted")
abline(h=0)
text(0,0.2,"99 %", cex=1.6)
text(-2.8,0.2,".5 % ")
text(2.8,0.2,".5 %")
arrows(-2.8, 0.18, -2.8, 0.015, lty=1)
arrows(2.8, 0.18, 2.8, 0.015, lty=1)
text(2.8,0.3, srt=90, adj = 0, labels = z, cex=.9)
text(-2.8,0.3, srt=90, adj = 0, labels = -z, cex=.9)

#CR49


pr=.5
n=1000
Z=1.96
Z*(sqrt(pr*(1-pr)/n))

#CR50


Z=1.96
pr=.5
margeErreur=.05
Z^2*pr*(1-pr)/margeErreur^2

Complément: marge d'erreur d'une moyenne avec l'écart-type


Z=1.96
et=3.5   #ecart-type
n=1000
Z*(et/sqrt(n))

#CR51


calculMargeErreur=function(Z,pr,n,N)
{
    et=(sqrt(pr*(1-pr)/n))*sqrt((N-n)/(N-1))
    e=Z*et
    e=round(e,2)
    print(c("Niveau de confiance:",Z),quote=F)
    print(c("Proportion:",pr),quote=F)
    print(c("Taille de la population:",round(N)),quote=F)
    print(c("Taille de l'échantillon:",n),quote=F)
    print(c("Erreur type:",et),quote=F)
    print(c("Marge d'erreur:",e),quote=F)
    print(c("Intervalle de confiance (borne inférieure):",pr-e),quote=F)
    print(c("Intervalle de confiance (borne supérieure):",pr+e),quote=F)
}

#CR52 après exécution du CR51


calculMargeErreur(1.96,.5,1000,1000000)
calculMargeErreur(1.96,.5,60,1000000)

#CR53


x=c(30,60,100,200,500,1000)
#Création d'une matrice qui contient le même nombre de rangées que x et deux colonnes
resultat=matrix(nrow=length(x), ncol=3)
#Produit une boucle qui remplit la matrice: en première colonne, le x; en deuxième colonne, le calcul de la marge d'erreur
for(i in 1:length(x)){
resultat[i,1]=x[i]
resultat[i,2]=(1.96*sqrt(.25/(x[i]))) *100
resultat[i,3]=(2.56*sqrt(.25/(x[i]))) *100
}
#affiche le résultat
resultat
#produit un graphique
par(mfcol=c(1,1))
plot(resultat[,1], resultat[,2], yaxt="n",xaxt="n", type='l', ylab="Marge d'erreur (en point de %)", xlab="Taille des échantillons", main="La marge d'erreur  (19 fois sur 20)") #ajoute le quadrillé vertical et horizontal au graphique: facultatif


axis(1,labels=x, at=x, las=2)
for (j in 1:length(x)) {
abline(v=x[j], lty=3)
}

for (j in 1:6) {
abline(h=resultat[j,2], lty=3)
text(2.7,resultat[j,2]+.2, round(resultat[j,2]))
}

#CR54


faireEchantillons=function(k)
{
 j=100 #Nombre fixe d'échantillons tirés

 dd=seq(4.5,100,6.1)
 piece=c(0,1)
 tableau=array(0, dim=c(j,1))
 for(j in 1:j){
   tirage=sample(piece,k,replace=T)
   prop=sum(tirage)/k
   tableau[j]=round(prop*100)
 }

hist(tableau, xaxt="n", breaks=dd, xlim=c(0,100), ylim=c(0,j), xlab="Pourcentages de piles par échantillon", main=paste("100 échantillons de", k ,"pièces"), ylab="Nombre d'échantillons", col="gray")
axis(1,labels=round(dd), at=dd, las=2)

}

faireEchantillons(60)

#CR55 fonction nécessaire: faireEchantillon() du CR52


par(mfcol=c(3,2))
faireEchantillons(30)
faireEchantillons(100)
faireEchantillons(500)
faireEchantillons(60)
faireEchantillons(200)
faireEchantillons(1000)
par(mfcol=c(1,1))

Distribution binomiale avec R

En utilisant la loi binomiale (Jacques Bernoulli 1654-1705), et la fonction correspondante dans R -dbinom()-, on peut obtenir un résultat analogue. On suppose deux issues ou résultats ou valeurs: vrai et faux, succès et échec, homme et femme, pile et face, etc. On veut déterminer la proportion d'une distribution qui aura une valeur donnée ou un intervalle de valeurs données. Pour mieux comprendre, examiner le code suivant: en utilisant la fonction on obtient la proportion de la distribution dont le nombre (de piles par exemple) se situe dans l'intervalle défini dans le premier argument. Dans le premier exemple, on obtient la proportion des 100 (deuxième argument) observations ou pièces de monnaie qui auront entre 47 et 53 piles (premier argument) selon l'hypothèse d'une distribution de 0,5 (troisième argument).


#Voici quatre exemples:

sum(dbinom(47:53, 100, 0.5))
sum(dbinom(94:106, 200, 0.5))
sum(dbinom(235:265, 500, 0.5))
sum(dbinom(470:530, 1000, 0.5))


#Pour mieux comprendre le premier exemple

a=dbinom(47, 100, 0.5)
b=dbinom(48, 100, 0.5)
c=dbinom(49, 100, 0.5)
d=dbinom(50, 100, 0.5)
e=dbinom(51, 100, 0.5)
f=dbinom(52, 100, 0.5)
g=dbinom(53, 100, 0.5)

a+b+c+d+e+f+g

Dans le quatrième exemple, le résultat permet de comprendre autrement la phrase qui dit que pour un échantillon de 1000 pièces de monnaie, on devrait obtenir 50% de piles avec une marge d'erreur de plus ou moins 3%, 19 fois sur vingt. Dit autrement, 95% des échantillons de 1000 pièces devraient avoir entre 470 et 530 piles.

#CR56 package nécessaire:vcd


library(vcd)
revenu=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0)
vote=c(1, 1, 1, 1, 1,1 ,1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1)

table(vote)
t=table(vote,revenu)
t
chisq.test(t)
assocstats(t)

#CR57 package nécessaire:vcd


revenu=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
vote=c(1, 1, 1, 1, 1,1 ,1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1)
revenu=rep(revenu,2)
vote=rep(vote,2)
table(vote)
t=table(vote,revenu)
t
chisq.test(t)
assocstats(t)

#CR58 package nécessaire: gplots


library(gplots)
deuxProportions=c(33,66)
margeErreur=18

par(mfcol=c(1,2))
plotCI(barplot(deuxProportions,col="gray",ylim=c(0,100),main="Marge d'erreur= +/- 17, n=30", sub="p=.14"),deuxProportions,margeErreur,add=TRUE, lwd=3)
abline(h=33+margeErreur, lty="dotted")
abline(h=66-margeErreur, lty="dotted")

margeErreur=13
plotCI(barplot(deuxProportions,col="gray",ylim=c(0,100),main="Marge d'erreur= +/- 12, n=60", sub="p=.02"),deuxProportions,margeErreur,add=TRUE, lwd=3)
abline(h=33+margeErreur, lty="dotted")
abline(h=66-margeErreur, lty="dotted")

#CR59


revenu=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
vote=c(1, 1, 1, 1, 1,1 ,1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1)
revenu=rep(revenu,2)
vote=rep(vote,2)
table(vote)
t=table(vote,revenu)
prop.test(t, alternative="two.sided")
prop.test(t, alternative="greater")

#CR60


prop.test(c(20,10),c(30,30), alternative="two.sided")
prop.test(c(20,10),c(30,30), alternative="greater")

#CR61


g1=c(1,2,2,3,3,3,4,4,5)
g2=c(3,4,4,5,5,5,6,6,7)
v3=c("A","B","A","B","A","B","A","B","A")
v4=c("A","A","A","A","A","B","B","B","B")
base=data.frame(g1,g2,v3,v4)

boxplot(base$g1, base$g2, main="Boxplot sur les deux groupes ",
names=c("g1","g2"))

t.test(base$g2,base$g1, alternative="greater")
t.test(base$g2~base$v3, alternative="two.sided")
t.test(base$g2~base$v4, alternative="less")


#Imaginons que g1 et g2 soient deux mesures sur les mêmes individus. On ajoutera l'argument suivant.
t.test(g1,g2, paired=TRUE, alternative="greater")

Complément: faire un t test sans les données mais uniquement avec les deux moyennes, les deux écarts-types et les deux n

Supposons qu'il s'agit du groupe1 (g1) selon la distribution de la variable v4 (cf le code précédent CR59). Le groupe des A a un n de 5, un écart-type de .83666 et une moyenne de 2,2. Ce sont les seules informations dont vous disposez pour ce groupe, sans avoir la distribution elle-même. Les valeurs du deuxième groupe sont les suivantes: le n est de 4, l'écart-type de 0.8164966 et la moyenne est de 4.


library(BSDA)
tsum.test(mean.x=2.2,s.x=.8366, n.x=5, mean.y=2.2,s.y=.8366, n.y=4, alternative="greater")

#CR62 package nécessaire: pwr


library(pwr)
tailleEffet=ES.h(.30,.10)
tailleEffet
pwr.2p.test(n=30, h=tailleEffet, sig.level=.05, alternative="greater")

tailleEffet=ES.h(.30,.10)
pwr.2p.test(n=60, h=tailleEffet, sig.level=.05, alternative="greater")

pwr.2p.test( h=.4, power=.8, sig.level=.05, alternative="greater")