1 Comparaison de deux moyennes

Test-t de Student
L’exercice vise à déterminer l’effet d’un traitement (nature de la vitamine C) sur la longueur cellulaire.

Pour cela l’objectif est de comparer les deux moyennes de longueurs cellulaires soumises à un traitement vitamine C naturelle et un traitement vitamine C synthétique.

30 minutes

1.0 Charger les données orange.csv

orange <- read.table("orange.csv", sep=";", dec=",", header=T)
edit(orange) # Penser à fermer la fenetre ensuite, sinon impossible d'exécuter des commandes dans la console

##    longueur traitement
## 1       8.2     Orange
## 2       9.4     Orange
## 3       9.6     Orange
## 4       9.7     Orange
## 5      10.0     Orange
## 6      14.5     Orange
## 7      15.2     Orange
## 8      16.1     Orange
## 9      17.6     Orange
## 10     21.5     Orange
## 11      4.2 Ascorbique
## 12      5.2 Ascorbique
## 13      5.8 Ascorbique
## 14      6.4 Ascorbique
## 15      7.0 Ascorbique
## 16      7.3 Ascorbique
## 17     10.1 Ascorbique
## 18     11.2 Ascorbique
## 19     11.3 Ascorbique
## 20     11.5 Ascorbique

head(orange)

##   longueur traitement
## 1      8.2     Orange
## 2      9.4     Orange
## 3      9.6     Orange
## 4      9.7     Orange
## 5     10.0     Orange
## 6     14.5     Orange

dim(orange)

## [1] 20  2

1.1 Décrire les données étudiées

type de variables, effectifs, représentation graphiques …

# lister les variables
names(orange)

## [1] "longueur"   "traitement"

# lister le contenu d'une variable
#longueur
# reconnaissance de la variable
attach(orange)
longueur

##  [1]  8.2  9.4  9.6  9.7 10.0 14.5 15.2 16.1 17.6 21.5  4.2  5.2  5.8  6.4  7.0
## [16]  7.3 10.1 11.2 11.3 11.5

# Nombre d'observations pour chaque modalité d'une variable qualitative
table(traitement)

## traitement
## Ascorbique     Orange 
##         10         10

# Statistiques calculées sur une variable quantitative
mean(longueur)

## [1] 10.59

var(longueur)

## [1] 20.01884

sd(longueur)

## [1] 4.474242

summary(longueur)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.200   7.225   9.850  10.590  12.250  21.500

# Statistiques par modalité de la variable qualitative
tapply(longueur, traitement, mean)

## Ascorbique     Orange 
##       8.00      13.18

tapply(longueur, traitement, sd)

## Ascorbique     Orange 
##   2.768072   4.437667

tapply(longueur, traitement, summary)

## $Ascorbique
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.20    5.95    7.15    8.00   10.93   11.50 
## 
## $Orange
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.200   9.625  12.250  13.180  15.875  21.500

Que déduisez vous des comparaisons des valeurs de moyennes, écart types et des autres statistiques ?
Pour cet échantillon, la longueur cellulaire moyenne est plus grande pour celles ayant reçu le traitement jus d’orange que celle ayant reçu le traitement acide ascorbique.

Pour cet échantillon, la longueur cellulaire de celles ayant reçu le traitement jus d’orange est plus variable que celle ayant reçu le traitement acide ascorbique.

On peut aussi visualiser ces données grâce à un box plot

data<-split(orange, orange$traitement)
data

## $Ascorbique
##    longueur traitement
## 11      4.2 Ascorbique
## 12      5.2 Ascorbique
## 13      5.8 Ascorbique
## 14      6.4 Ascorbique
## 15      7.0 Ascorbique
## 16      7.3 Ascorbique
## 17     10.1 Ascorbique
## 18     11.2 Ascorbique
## 19     11.3 Ascorbique
## 20     11.5 Ascorbique
## 
## $Orange
##    longueur traitement
## 1       8.2     Orange
## 2       9.4     Orange
## 3       9.6     Orange
## 4       9.7     Orange
## 5      10.0     Orange
## 6      14.5     Orange
## 7      15.2     Orange
## 8      16.1     Orange
## 9      17.6     Orange
## 10     21.5     Orange

boxplot(orange$longueur~as.factor(orange$traitement))

Notre question est donc : est-ce-qu’en moyenne les cellules ayant reçu de la vitamine C naturelle sont plus longue que celle ayant reçu de la vitamine C synthétique ?

1.2 Selon la question écologique quelles sont nos hypothèses ?

On pose nos hypothèses avec des mots :

Les cellules reçevant de la vitamine C synthétique sont en moyenne moins longue que les cellules reçevant de la vitamine C naturelle

On pose aussi les hypthèses du test :
H0: \(\mu_o = \mu_a\)

H1: \(\mu_o > \mu_a\) (test unilatéral)

H1: \(\mu_o != \mu_a\) (test bilatéral)

1.3 Définir le test approprié pour répondre à la question

On effectue une comparaison de moyennes, c’est pourquoi on fait un test de Student

1.4 Vérification des conditions de validité

A. Evaluer la normalité des données

Hypothèse nulle: les données sont normalement distribué

H0 : Wobs = Wth(alpha-1,n)

On effectue pour tester la normalité des données un test de Shapiro

shapiro.test(data$Ascorbique$longueur)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$Ascorbique$longueur
## W = 0.88648, p-value = 0.1547

Que déduisez vous concernant la normalité des données ?
On ne rejette pas H0 au risque alpha = 5%, car la pvalue est supérieur à 0.05. On considère donc que les données sont normalement distribuées

On peut aussi visualiser grâce à la fonction qqnrom() . Dans ce cas les données doivent être alignée sur la première bisectrice (ligne rouge)

qqnorm(data$Ascorbique$longueur)
qqline(data$Ascorbique$longueur, col="red")

Quelle est donc l’étape suivante ? Vérifions l’homoscédasticité des variables (homogénéité des variances)

B. Evaluer l’homogénéité des variances (homoscédasticité)

On effecture un test F qui compare 2 variances

Hypothèse nulle : les variances sont égales

H0: F = 1 où F est le rapport des variances

H1 : F != 1

var.test(data$Ascorbique$longueur, data$Orange$longueur)

## 
##  F test to compare two variances
## 
## data:  data$Ascorbique$longueur and data$Orange$longueur
## F = 0.38909, num df = 9, denom df = 9, p-value = 0.1759
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.09664339 1.56645691
## sample estimates:
## ratio of variances 
##          0.3890857

Que déduisez vous concernant l’homogénéité des variances ?
On ne rejette pas H0, car pvalue > 0.05. Les variances sont donc considérées comme égales

1.5 Choissir le risque de première espèce

On prend alpha = 5%

1.6. Réaliser le test en calculant la statistique observée et en la comparant à la statistique attendue sous l’hypothèse H0

On réalise un test de student sur échantillons non appariés avec hypothèse d’homoscédasticité et de normailité vérifiées

t.test(data$Orange$longueur, data$Ascorbique$longueur,  alternative = c("two.sided"), var.equal = TRUE, paired=F) # Test bilatéral

## 
##  Two Sample t-test
## 
## data:  data$Orange$longueur and data$Ascorbique$longueur
## t = 3.1319, df = 18, p-value = 0.005762
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.705205 8.654795
## sample estimates:
## mean of x mean of y 
##     13.18      8.00

t.test(data$Orange$longueur, data$Ascorbique$longueur,  alternative = c("greater"), var.equal = TRUE, paired=F) # test unilatéral àdroite avec H1: mu(Orange)>mu(ascorbique)

## 
##  Two Sample t-test
## 
## data:  data$Orange$longueur and data$Ascorbique$longueur
## t = 3.1319, df = 18, p-value = 0.002881
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  2.311966      Inf
## sample estimates:
## mean of x mean of y 
##     13.18      8.00

t.test(data$Orange$longueur, data$Ascorbique$longueur,  alternative = c("less"), var.equal = TRUE,  paired=F) # test unilatéral à gauche avec H1: mu(Orange)<mu(ascorbique)

## 
##  Two Sample t-test
## 
## data:  data$Orange$longueur and data$Ascorbique$longueur
## t = 3.1319, df = 18, p-value = 0.9971
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##      -Inf 8.048034
## sample estimates:
## mean of x mean of y 
##     13.18      8.00

1.7 Rejeter ou non l’hypothèse nulle en accord avec le risque alpha

On rejette H0 au risque alpha = 5%.

1.8 Interpréter les résultats au niveau biologique

Les moyennes des longueurs cellulaires sont différentes en fonction du traitement reçu

Dans le cas où la normalité n’aurait pas été respectée, on aurrait appliqué le test non paramétrique de : Wilcoxon si échantillons sont appariés :

#wilcox.test(data$Orange$longueur, data$Ascorbique$longueur, paired=T)

Wilcoxon Mann-Whitney si ils n’étaient pas non appariés

#wilcox.test(data$Orange$longueur, data$Ascorbique$longueur, paired=F)

2 Etude de la relation linéaire entre deux variables quantitatives

Variance, covariance, correlation, draftsman plot
Regression lineaire simple
L’objectif de l’excercice est d’étudier la relation linéaire entre la taille de becs d’oiseaux femelles et males
1 heure

2.0 Charger les données size.txt

#do<-read.table(file.choose(), header=T, sep="\t", dec=",")
do<-read.table("Size.txt", header=T, sep="\t", dec=",")

2.1 Décrire les données

names(do)

## [1] "id"          "size"        "survival"    "mother_id"   "father_id"  
## [6] "mother_size" "father_size"

do$size

##  [1] 10.81  8.74 10.04  9.75  9.76  9.55  8.70  9.62  9.22  8.79  9.61  9.02
## [13]  7.85  9.01  8.26  8.25  9.25  9.45  8.52  8.61  8.12  7.81

head(do)

##    id  size survival mother_id father_id mother_size father_size
## 1 669 10.81        1       737       476       11.71        9.90
## 2 128  8.74        0       683       688       10.95       10.86
## 3 481 10.04        0       668       410       10.77        9.81
## 4 364  9.75        0       667       353       10.76        9.64
## 5 365  9.76        0       617       529       10.53       10.27
## 6 335  9.55        0       524       269       10.22        9.38

Moyenne, variance et ecart type de chaque variable

mean(do$father_size)

## [1] 9.301364

mean(do$mother_size)

## [1] 9.675

var(do$father_size)

## [1] 0.3866314

var(do$mother_size)

## [1] 0.7566452

sd(do$father_size)

## [1] 0.6217969

sd(do$mother_size)

## [1] 0.8698536

Que déduisez vous de la comparaison de ces statistiques entre femelles et males ?
La taille des becs est similaire entre mâles et femelles en moyenne

Covariance

cov(do$father_size, do$mother_size)

## [1] 0.4507976

Relation entre deux variables

plot(father_size ~ mother_size, data=do, col="blue")

Que déduisez vous de l’étude graphique de cette relation ?

Ces 2 variables semblent corrélées

2.2 Mesure du coefficient de correlation de Pearson r et Hypothèses

Le coefficient de Pearson est utilisé, car la relation est linéaire monotone sans points extremes. Sinon on ajoute l’argument method=“spearman” pour coefficient non parametrique. S’il y avait des données manquantes, on ajouterais use=“complete” pour le calcul uniquement sur les individus disponibles pour les deux variables.

cor(do$father_size, do$mother_size)

## [1] 0.8334641

2.3 Test de corrélation entre deux échantillons appariés

Etape 2 : Poser les hypothèses
H0: r = 0 (les 2 variables ne sont pas corrélées)
H1: r != 0

Etape 3: Définir le test approprié
Test de corrélation entre deux échantillons appariés

Etape 4 : Déterminer les conditions de validité
La relation est linéaire monotone sans points extremes

Etape 5 : Choisir le risque de première espèce
On prend un risque alpha = 5%

Etape 6 : Réaliser le test

cor.test(do$father_size, do$mother_size)

## 
##  Pearson's product-moment correlation
## 
## data:  do$father_size and do$mother_size
## t = 6.7455, df = 20, p-value = 1.461e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6349874 0.9287233
## sample estimates:
##       cor 
## 0.8334641

Etape 7 : Rejetter ou non l’hypothèse nulle
p-value = 1.461e-06 < alpha. On rejette H0, les 2 variables sont linéairement corrélées

Etape 8 : Interpréter les résultats
Le bec des mâles est linéairement corrélé au bec des femelles

2.4 Modèle linéaire : regression lineaire simple

On souhaite quantifier la force de la relation. On effectue pour cela un modèle linéaire.

model <- lm(father_size ~ mother_size, data=do)
?lm
names(model)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

2.4.1 Verification des conditions de validite du modele sur les résidus:

A. Independance des residus : Test de Durbin Watson

model.res <- residuals(model)

Hypothèse nulle : les résidus sont indépendants (i.e., absence d’autocorrelation)

Hypothèse alternative : Autocorrelation : les résidus successifs se ressemblent

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(model)

## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.6265, p-value = 0.903
## alternative hypothesis: true autocorrelation is greater than 0

Qu’en déduisez vous ?

p>0.05, on ne rejette pas H0, les résidus sont indépendants

par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(model)

A quoi correspond chacun de ces graphiques ?
residuals vs fitted : Homegeneite des residus (doivent être autour de 0 sans structure particulière)
Normal QQ : Normalité des résidus (doivent suivre la droite y=x)
Scale location : normalité des résidus standardisé (pas de structure)
Residual vs leverage : influence d’une observation yi sur sa valeur estimée ŷi (pas de structure) : identifier des outliers

Qu’en déduisez vous ? Tout semble ok, mais on réalise des tests pour en être sur

B. Normalité des résidus

On effectue un Test de Shapiro H0 : Normalité des données

# Droite de Henry (graph ci-dessus) ainsi que test de Shapiro :
qqnorm(model.res)
qqline(model.res)

shapiro.test(model.res)

## 
##  Shapiro-Wilk normality test
## 
## data:  model.res
## W = 0.97199, p-value = 0.7566

Qu’en déduisez vous ? On ne rejette pas H0 -> normalité des données

C. Homogeneité des variances des résidus : Homoscédasiticité

Deux graphiques ci-dessus des résidus en fonction des valeurs prédites et test defpaganch-Pagan H0 : Homogeneité des variances des résidus

bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 4.5853, df = 1, p-value = 0.03225

Qu’en déduisez vous ? On rejette H0 : pas d’homoscédasticité

D. Outliers

Détection d’individus influençants les paramètres du modèle : graphique distance de Cook et effet levier.

L’individu n°1 (ie. premiere ligne du tableau) semble influant. On le retire et refaisons le modele. Siil correspond à une réalité biologique, écologique (ie. pas d’erreur de mesure, de saisie de données etc.), on l’indique dans les résultats.

model1 <- lm(father_size[-1] ~ mother_size[-1], data=do[-1,])

On verifie les conditions sur le nouveau modèle :

model.res1 <- residuals(model1)

A. Independance des residus : Test de Durbin Watson

dwtest(model1)

## 
##  Durbin-Watson test
## 
## data:  model1
## DW = 2.4835, p-value = 0.8073
## alternative hypothesis: true autocorrelation is greater than 0

Qu’en déduisez vous ?
p>0.05, on ne rejette pas H0, les résidus sont indépendants

par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(model1)

B. Normalité des résidus

# Droite de Henry (graph ci-dessus) ainsi que test de Shapiro :
qqnorm(model.res1)
qqline(model.res1)

shapiro.test(model.res1)

## 
##  Shapiro-Wilk normality test
## 
## data:  model.res1
## W = 0.95042, p-value = 0.3735

Qu’en déduisez vous ? Normalité

C. Homogeneité des variances des résidus :

Deux graphiques ci-dessus des résidus en fonction des valeurs prédites et test de Breusch-Pagan

bptest(model1)

## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 0.0032095, df = 1, p-value = 0.9548

Qu’en déduisez vous ?

On ne rejette pas H0, Homoscédasticité des résidus

2.4.2 Etude des paramètres du modèles une fois les conditions ci-dessus validées

summary(model1)

## 
## Call:
## lm(formula = father_size[-1] ~ mother_size[-1], data = do[-1, 
##     ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53784 -0.19940 -0.03606  0.17598  0.44722 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.32959    0.87699   3.797  0.00132 ** 
## mother_size[-1]  0.61664    0.09198   6.704 2.75e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2846 on 18 degrees of freedom
## Multiple R-squared:  0.714,  Adjusted R-squared:  0.6981 
## F-statistic: 44.95 on 1 and 18 DF,  p-value: 2.753e-06

plot(father_size[-1] ~ mother_size[-1], data=do[-1,], col="blue")
abline(model1, col="red")

Interprétez les différents résultats indiqués par le modèle, ceci couplé à la représentation graphique.
Call nous donne le modèle que nous avons ajusté

Residuals nous donne la distribution des résidus. Les résidus correspondent à la différence entre ce que le modèle prédit et la valeur observée des y.

Coefficients nous donne les paramètres estimés pour le modèle. L’intercept est l’ordonné à l’origine et mother_size[-1] est le coefficient directeur de la droite de régression (la pente). Pour chaque paramètre on a peut trouver la valeur estimée, l’erreur type associée (écart type divisé par la racine carré de la taille de l’échantillon). Puis on trouve la variable de decision (Estimate/Std. Error) et la pvalue associé au test. Ce test t regarde si le paramètre estimé est différent de 0.

Residual standard error:

k=length(model1$coefficients)-1 #Subtract one to ignore intercept
SSE=sum(model1$residuals^2)
n=length(model1$residuals)
sqrt(SSE/(n-(1+k)))

## [1] 0.2846133

Multiple R-squared : Une mesure de la “distance” entre le modèle des données

#Multiple R-Squared (Coefficient of Determination)
SSyy=sum((do[-1,]$father_size[-1]-mean(do[-1,]$father_size[-1]))^2)
SSE=sum(model1$residuals^2)
(SSyy-SSE)/SSyy

## [1] 0.7140368

#Alternatively
1-SSE/SSyy

## [1] 0.7140368

Adjusted R-squared : normalise Multiple R-squared en tenant compte de la taille de l’échantillon et du nombre de variables utilisé dans le modèle. Plus on ajoute de variables au modèle plus on va pouvoir expliquer de variance.

1-(SSE/SSyy)*(n-1)/(n-(k+1))

## [1] 0.6981499

F-statistic : Finalement la fonction summary() réalise un autre test. C’est un test “global” qui vérifie si au moin une variable est différente de 0.

#F-Statistic
#H0: Tous les coefficiants valent 0
#H1: Au moins un des coefficient est différent de 0
#Compare test statistic to F Distribution table
((SSyy-SSE)/k) / (SSE/(n-(k+1)))

## [1] 44.94516

Que concluez vous concernant la relation entre la taille des femelles et des males ?
La relation entre les deux variable est positive. Plus le bec des mâles est grand plus la taille du bec des femelles l’est aussi. D’après les résultats du modèle linéaire, on peut dire que la taille des becs des femelles est 0.6 fois plus petite que la taille des becs des mâles.

Si les conditions sur les résidus ci-dessus ne sont pas respectées, en particulier la normalité (mais également égalité des variances), plusieurs possibilités de transformation des données y (variable réponse) existent avant de refaire le modele :

transformation racine carrée sqrt() : permet de normaliser les données qui suivent une loi de Poisson (distribution asymétrique et variance égale à la moyenne) :

# Exemple visuel sur des données simulées : 
y <- rpois(1000, 3)
y1 <- sqrt(y) 
hist(y)

hist(y1)

transformation logarithmique log() : permet de normaliser les données qui suivent une loi lognormale ou binomiale négative. Dans le cas où y peut être égal à zéro, utiliser log(y+1). Si la transformation log ne suffit pas, essayer le log au carré.
transformation boxcox : Lorsqu’il n’y a pas de raison de choisir l’une ou l’autre de ces transformations, on peut utiliser la fonction boxcox() du package MASS ou boxCox() du package car.
transformation arc-sinus de la racine carrée asin(sqrt(y)) : dans le cas de données en proportions ou pourcentage, où la variance dépend de la moyenne et dont la distribution est souvent trop étalée. Les données doivent être entre 0 et 1

Pour des cas de relations non linéaires entre deux variables quantitatives, voir la fonction ?nls

3 Représentation de la relation entre variables deux à deux : draftsman plot

library(ade4)
library(psych)

Les données Iris concernent des variables morphométriques de plantes issues de trois espèces différentes

pairs.panels(iris, 
             method = "spearman", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  # show density plots
             ellipses = FALSE # show correlation ellipses
             )

Ci-dessus on utilise la correlation de Spearman (basee sur le rang des valeurs) car on identifie sur les nuages de points des tendances non lineaires et/ou des points avec des valeurs relativements fortes/faibles Les correlations avec la valeur “Species” ne sont pas à utiliser car c’est une variable qualitative. Ce qui est interessant ici est d’etudier les valeurs des mesures morphologiques intra-inter specifiques qui donnent deja des informations sur des variations.

Comment interpretez vous cette figure ?

On a que la longueur et la largueur des pétales sont très corrélés (r = 0.94). Ces 2 variables sont aussi corrélés avec la longueur des sépales (resp r= 0.88 et r =0.83).

Une espèce semble avoirs des valeurs différentes des 2 autres.

Objet : TP 3 - HAB711B Description et inference DESINF - UE Master 1

Enseignants TP : Camille Magneville, Fabien Leprieur, Julien Claude, Bastien Merigot

26 aout 2021

1 Comparaison de deux moyennes

1.0 Charger les données orange.csv

1.1 Décrire les données étudiées

1.2 Selon la question écologique quelles sont nos hypothèses ?

1.3 Définir le test approprié pour répondre à la question

1.4 Vérification des conditions de validité

A. Evaluer la normalité des données

B. Evaluer l’homogénéité des variances (homoscédasticité)

1.5 Choissir le risque de première espèce

1.6. Réaliser le test en calculant la statistique observée et en la comparant à la statistique attendue sous l’hypothèse H0

1.7 Rejeter ou non l’hypothèse nulle en accord avec le risque alpha

1.8 Interpréter les résultats au niveau biologique

2 Etude de la relation linéaire entre deux variables quantitatives

2.0 Charger les données size.txt

2.1 Décrire les données

2.2 Mesure du coefficient de correlation de Pearson r et Hypothèses

2.3 Test de corrélation entre deux échantillons appariés

2.4 Modèle linéaire : regression lineaire simple

2.4.1 Verification des conditions de validite du modele sur les résidus:

A. Independance des residus : Test de Durbin Watson

B. Normalité des résidus

C. Homogeneité des variances des résidus : Homoscédasiticité

D. Outliers

A. Independance des residus : Test de Durbin Watson

B. Normalité des résidus

C. Homogeneité des variances des résidus :

2.4.2 Etude des paramètres du modèles une fois les conditions ci-dessus validées

3 Représentation de la relation entre variables deux à deux : draftsman plot