3 Inference statistique sur une et deux variables

Cette section présente différentes fonctions R utiles pour réaliser de l’inférence statistique sur une ou deux variables.

3.1 Données Chocolat

Les fonctions de cette section sont illustrées sur base du jeu de données : chocolat.xlsx (disponible sur moodle). Ce data.frame de taille 20x2 reprend des mesures de la teneur en cacao de 20 plaquettes de chocolat de marque “Cabeau” et “Gallet”.

Voici le code pour lire et préparer ces données

choco <- readxl::read_xlsx("data/chocolat.xlsx", col_names = TRUE)
choco$Marque <- as.factor(choco$Marque)
cabeau <- subset(choco, Marque =="Cabeau") 
gallet <- subset(choco, Marque =="Gallet")

Voici quelques résumés descriptifs de ces données

Visualiser via un box-plot

boxplot(Cacao ~ Marque, data = choco, main = "Box plot des données Chocolat", ylab = "Cacao")

Tableau résumé par marque

N <- table(choco$Marque)
Moyenne<-tapply(choco$Cacao, choco$Marque, mean)
EcartType<-tapply(choco$Cacao, choco$Marque, sd)
cbind(N,Moyenne,EcartType)

	N	Moyenne	EcartType
Cabeau	10	6.063	0.2568
Gallet	10	6.258	0.108

3.2 Arguments utiles dans les fonctions d’inférence

Argument alternative = dans une fonction de test

$H_0: \theta = \theta_0 \textrm{ contre } H_1: \theta \neq \theta_0$ ~> alternative = "two.sided"
$H_0: \theta \ge \theta_0 \textrm{ contre } H_1: \theta < \theta_0$ ~> alternative = "less"
$H_0: \theta \le \theta_0 \textrm{ contre } H_1: \theta > \theta_0$ ~> alternative = "greater"

Argument conf.level = dans une fonction de test

Par défaut, le niveau de confiance des intervalles de confiance est 95%. On peut le modifier via l’argument conf.level = 0.9 pour, par exemple, un niveau de confiance de 90%.

3.3 Tests sur une ou deux moyennes

3.3.1 Test t sur une moyenne

$H_0: \mu = 6 \textrm{ contre } H_1:\mu \neq 6$

t.test(choco$Cacao, mu = 6, alternative = "two.sided")


	One Sample t-test

data:  choco$Cacao
t = 3.3192, df = 19, p-value = 0.003606
alternative hypothesis: true mean is not equal to 6
95 percent confidence interval:
 6.059293 6.261707
sample estimates:
mean of x 
   6.1605

la fonction pander dans RMarkdown permet de mieux présenter les résultats mais ne présente pas tout. Ceci est aussi vrai pour les autres fonctions de test ou de modélisation.

res= t.test(choco$Cacao, mu = 6, alternative = "two.sided")
pander::pander(res)

One Sample t-test: `choco$Cacao`
Test statistic	df	P value	Alternative hypothesis	mean of x
3.319	19	0.003606 * *	two.sided	6.16

3.3.2 Test t de comparaison de deux moyennes

$H_0: \mu_1 = \mu_2 \textrm{ contre } H_1:\mu_1 \neq \mu_2$

Echantillons indépendants et variances égales

t.test(choco$Cacao~choco$Marque, alternative = "two.sided", var.equal = TRUE)

Echantillons indépendants et variances différentes

t.test(choco$Cacao~choco$Marque, alternative = "two.sided")


	Welch Two Sample t-test

data:  choco$Cacao by choco$Marque
t = -2.2137, df = 12.087, p-value = 0.04682
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.386775949 -0.003224051
sample estimates:
mean in group Cabeau mean in group Gallet 
               6.063                6.258

Données pairées (non applicable aux données chocolat)

# Observations pairées
t.test(y~x, alternative = "two.sided", paired = TRUE)

3.4 Tests sur une ou plusieurs variances

3.4.1 Test $\chi^2$ sur une variance

$H_0: \sigma^2 = 0.04 \textrm{ contre } H_1: \sigma^2 \neq 0.04$

EnvStats::varTest(x = cabeau$Cacao, sigma.squared = 0.04, alternative = "two.sided")


	Chi-Squared Test on Variance

data:  cabeau$Cacao
Chi-Squared = 14.835, df = 9, p-value = 0.1911
alternative hypothesis: true variance is not equal to 0.04
95 percent confidence interval:
 0.03119472 0.21974978
sample estimates:
  variance 
0.06593444

3.4.2 Test $F$ de comparaison de 2 variances

$H_0: \sigma_1^2 = \sigma_2^2 \textrm{ contre } H_1:\sigma_1^2 \neq \sigma_2^2$

var.test(choco$Cacao~choco$Marque, ratio = 1, alternative = "two.sided")


	F test to compare two variances

data:  choco$Cacao by choco$Marque
F = 5.6537, num df = 9, denom df = 9, p-value = 0.0166
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
  1.404294 22.761673
sample estimates:
ratio of variances 
          5.653678

3.4.3 Test de Levene d’homogénéité de variances

Ce test permet de tester si 2 ou plus de 2 variances sont égales.

$H_0: \sigma_1^2 = \sigma_2^2=...=\sigma_k^2 \textrm{ contre } H_1:\textrm{ Au moins deux variances sont différentes }$

car::leveneTest(choco$Cacao, group = choco$Marque)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)  
group  1  6.1977 0.0228 *
      18                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.5 Tests sur une ou plusieurs proportions

Pour les tests ci-dessous, utilisons le simple tableau de contingence pour 2 variables qualitatives X et Y de niveaux respectifs A1,A2,A3 et B1,B2 :

data=matrix(c(10,8,17,5,5,10),nrow=3,ncol=2)
dimnames(data)=list(c("A1","A2","A3"),c("B1","B2"))
pander(data)

	B1	B2
A1	10	5
A2	8	5
A3	17	10

3.5.1 Test $\chi^2$ sur une proportion

Teste si la proportion lié au niveau B1 de la variable Y est plus grande que 0.5. Ce test permet des hypothèses alternatives unilatérales >ou <. Il effectue aussi une correction de continuité.

$H_0: \pi_{B1} \le 0.5 \textrm{ contre } H_1:\pi_{B1} > 0.5$

nB1=sum(data[,1]) # Nb de données telles que Y=B1
n=sum(data) # Nb total de données
prop.test(nB1,n, p=0.5, alternative = "greater")


	1-sample proportions test with continuity correction

data:  nB1 out of n, null probability 0.5
X-squared = 3.5636, df = 1, p-value = 0.02953
alternative hypothesis: true p is greater than 0.5
95 percent confidence interval:
 0.5164367 1.0000000
sample estimates:
        p 
0.6363636

3.5.2 Test d’ajustement $\chi^2$

On teste pour la variable X si la proportion d’individus dans chaque niveau est identique $H_0: \pi_{A1} = \pi_{A2}= \pi_{A3} \textrm{ contre } H_1: \textrm{ non } H_0$

dataY=apply(data,1,sum)
chisq.test(dataY,p=c(1/3,1/3,1/3))


	Chi-squared test for given probabilities

data:  dataY
X-squared = 6.2545, df = 2, p-value = 0.04384

3.6 Tests de corrélation / indépendance

Test sur un coefficient de corrélation de Pearson

cor.test(x = .., y = .., alternative = "two.sided") #x,y quantitatives

Test $\chi^2$ d’indépendance

chisq.test(data)

Tests de corrélation 2 à 2

res <- Hmisc::rcorr(df_quanti)
                # df_quanti = data frame composé uniquement de variables quantitatives

~> res$r = matrice de corrélation
~> res$p = matrice des p-valeurs associées

3.7 Intervalles de confiance

Il y a 2 possibilités pour obtenir un intervalle de confiance dans R:

1) Calcul à la main
Si on dispose de formules, nous pouvons calculer à la main les bornes inférieures et supérieures d’IC. Pour cela on a besoin d’estimations des paramètres et de quantiles de lois de probabilité bien connues.

Les quantiles sont obtenus via les fonctions suivantes:

Normale ~> qnorm(p = .., mean = .., sd = ..)
Student ~> qt(p = .., df = ..)
Chi_Carré ~> qchisq(p = .., df = ..)

Un intervalle de confiance à 95% sur une moyenne des données chocolat de la société cabeau sera par exemple obtenu comme suit:

n=length(cabeau$Cacao)
mean(cabeau$Cacao)+c(-1,1)*qt(0.975,n-1)*sd(cabeau$Cacao)/sqrt(n)

[1] 5.879313 6.246687

2) Utiliser les fonctions de test
Comme nous pouvons voir plus haut, certains tests calculent automatiquement les intervalles de confiance 95%. Il suffit donc de voir quel test donne l’IC recherché et d’utiliser des hypothèses nulles quelconques:

Moyenne ~> test t sur 1 échantillon t.test()
Différence de moyennes ~> test t sur 2 échantillons t.test()
Variance ~> test Chi_Carré sur la variance EnvStats::varTest()
Ratio de variances ~> test F de comparaison de variances var.test()
Proportion ~> test de proportion à 1 échantillon prop.test()
Odds ratio ~> test exact de Fisher fisher.test()
Corrélation ~> test sur le coefficient de corrélation de Pearson cor.test()

Pour ces tests on peut ajouter l’argument conf.level = .. dans la fonction pour choisir soi-même le niveau de confiance de l’intervalle - par défaut 95%.

Extraire l’IC

res <- t.test(cabeau$Cacao, mu = 6, alternative = "two.sided")
ic <- data.frame(Estimate = mean(cabeau$Cacao),
                 Lwr = res$conf.int[1],
                 Upr = res$conf.int[2])
pander::pander(ic)

Estimate	Lwr	Upr
6.063	5.879	6.247