Session 2 : Statistiques pratiques sur données continues

Pour les utilisateurs de R, on commence par charger un certain nombre de modules intéressants.

# le package moderne standard pour manipuler ses données
library(tidyverse)

# quelques couleurs manuelles
bleufonce <- "#3d5468"
bleuclair <- "#5b7c98"
rose <- "#ff5555"

Préparation des données

Simulation

Pour se focaliser sur les stats et uniquement les stats, on travaille sur un jeu de données idéalisé. Et puisqu’on bosse dans les essais cliniques, on va simuler le RCT (Randomized Controlled Trial) ayant les caractéristiques suivantes :

• l’étude possède deux bras: control et treatment.
• on enregistre des variables démographiques à l’inclusion, type sex, age, age_group, height, weight_t0.
• une variable d’outcome liés à la prise de poids : diff_weight.
• une variable d’outcome à trois modalités : feels_sleepy.
• une variable de temps de suivi du traitement : t_compliance.
• une variable d’outcome de nombre de chutes pendant le suivi : n_falls.

Pour ceux que ça intéresse, voici la loi dans laquelle les données sont simulées :

# data that do not depend on the group
d <- data.frame(id = seq(1,500),
                group = as.factor(c(rep("control",250), rep("treatment",250))),
                sex = ifelse(rbinom(n=500, size=1, prob=0.4) == 0, "F", "M"),
                age = round(runif(n=500, min=19, max=77)),
                height = round(rnorm(n=500, mean=1.7, sd=0.1), digits=2),
                weight_t0 = round(rnorm(n=500, mean=70, sd=10)),
                t_compliance = round(rexp(n=500, rate=0.05))  )

# possibilities for the "feels_sleepy" variable
modalites <- c("morning", "afternoon", "always")

# data that do depend on the group
d <- d %>% mutate(diff_weight = c(rnorm(n=250, mean=0, sd=2.5), rnorm(n=250, mean=0.025*d$t_compliance[251:500], sd=2.5)),
                  feels_sleepy = c( sample(modalites, 250, prob = c(0.3, 0.4, 0.3), replace = TRUE),
                                    sample(modalites, 250, prob = c(0.2, 0.3, 0.5), replace = TRUE) ),
                  n_falls = c(rpois(n=250, lambda=6), rpois(n=250, lambda=5.5)) )

d <- d %>% mutate(weight_tf = weight_t0 + diff_weight,
                  age_group = as.factor(ifelse(age < 40, "moins de 40", ifelse(age < 60, "entre 40 et 60", "plus de 60"))) ) 

Ecriture/Lecture du dataset

Pour être sûr de travailler sur le fichier tel qu’il sera lu par tout le monde, on enregistre le dataset simulé puis on écrase le dataset avec l’import.

write.csv(d, "etude_fictive.csv")
d <- read.csv("etude_fictive.csv")
head(d)

##   X id   group sex age height weight_t0 t_compliance diff_weight feels_sleepy
## 1 1  1 control   F  77   1.59        72           16   8.1030007      morning
## 2 2  2 control   M  62   1.72        72           14   1.8256008      morning
## 3 3  3 control   F  74   1.83        79            4   2.0608949       always
## 4 4  4 control   F  42   1.56        65           51   2.4725727    afternoon
## 5 5  5 control   F  75   1.51        65           22   0.2801064       always
## 6 6  6 control   F  50   1.80        64           48   2.8282484      morning
##   n_falls weight_tf      age_group
## 1       3  80.10300     plus de 60
## 2       5  73.82560     plus de 60
## 3       5  81.06089     plus de 60
## 4       7  67.47257 entre 40 et 60
## 5       7  65.28011     plus de 60
## 6       8  66.82825 entre 40 et 60

Nous sommes prêts à travailler sur ce jeu de données, et, pour ce qui nous intéresse aujourd’hui, à analyser les données continues qu’il renferme.

T test : adéquation à une norme

On observe des différences de poids \((X_i)\) chez nos individus entre le début et la fin de l’expérience. Est-ce que la différence de poids moyenne est différente de zéro ?

On souhaite réaliser le test d’hypothèse suivant :

• H0 : la moyenne des \((X_i)\) vaut zero.
• H1 : la moyenne des \((X_i)\) est différente de zéro.

Conditions d’application :

• les variables aléatoires sont indépendantes et de même loi.
• soit les variables aléatoires sont distribuées normalement, soit l’échantillon est “grand”.

On commence toujours par représenter nos données.

plot_diff_weight <- d %>% ggplot(aes(x=diff_weight)) +
    geom_histogram(alpha=0.7, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois")
plot_diff_weight

Ici les conditions d’application semblent respectées.

Un résultat théorique certifie qu’une statistique de test calculée à partir de l’effectif, de la moyenne empirique et de la variance empirique, est distribuée suivant une loi de Student avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

t.test(x=d$diff_weight)

## 
##  One Sample t-test
## 
## data:  d$diff_weight
## t = 3.7315, df = 499, p-value = 0.0002122
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.1980453 0.6385225
## sample estimates:
## mean of x 
## 0.4182839

L’output de la fonction nous rappelle :

1. en première ligne l’échantillon sur lequel on a appliqué le test,
2. en seconde ligne on récupère la valeur observée de la statistique de test, ainsi que le nombre de “degrés de liberté” de la loi théorique.
3. en fin de seconde ligne, super importante : la p-value !!! Probabilité que la valeur de la statistique soit au moins aussi extrême qu’observée.
4. troisième ligne : un rappel de l’hypothèse alternative du test (modifiable avec l’argument alternative).
5. un intervalle de confiance à 95% de la moyenne (bornes modifiables via l’argument conf.level).
6. tout à la fin, l’estimation ponctuelle de la moyenne.

Que concluons nous ???

T test : comparaison de 2 moyennes

On observe des valeurs \((X_{i,j})\) de nos individus \(i\) dans deux groupes \(j=1\) et \(j=2\). Ici les deux groupes correspondent au control vs. treatment, et les observations correspondent à une différence de poids entre le début et la fin de l’expérience.

Est-ce que la différence de poids moyenne est différente dans les deux groupes ?

On souhaite réaliser le test d’hypothèse suivant :

• H0 : la moyenne des \((X_{i,1})\) est identique à celle des \((X_{i,2})\).
• H1 : les deux moyennes sont différentes.

Conditions d’application :

• les variables aléatoires sont indépendantes et de même loi.
• soit les variables aléatoires sont distribuées normalement, soit l’échantillon est “grand”.

On commence toujours par représenter nos données et quantités d’intérêt.

m <- d %>% group_by(group) %>% summarise(mean=mean(diff_weight)) %>% select(mean)

plot_diff_weight_by_group <- d %>% ggplot(aes(x=diff_weight, fill=group)) +
    geom_histogram(alpha=0.5, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois") +
    scale_fill_manual(values=c(bleuclair, rose)) +
    geom_vline(xintercept = c(m$mean[1], m$mean[2]), color=c(bleuclair, rose))
plot_diff_weight_by_group

Ici les conditions d’application semblent respectées.

Un résultat théorique certifie qu’une statistique de test calculée à partir de l’effectif, des moyennes empiriques des deux groupes et de la variance empirique des deux groupes, est distribuée suivant une loi de Student avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

t.test(x=d$diff_weight[d$group == "control"], y=d$diff_weight[d$group == "treatment"])

## 
##  Welch Two Sample t-test
## 
## data:  d$diff_weight[d$group == "control"] and d$diff_weight[d$group == "treatment"]
## t = -1.7256, df = 495.02, p-value = 0.08505
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.82571065  0.05352049
## sample estimates:
## mean of x mean of y 
## 0.2252363 0.6113314

Une syntaxe alternative pratique existe également :

t.test(diff_weight ~ group, data=d)

## 
##  Welch Two Sample t-test
## 
## data:  diff_weight by group
## t = -1.7256, df = 495.02, p-value = 0.08505
## alternative hypothesis: true difference in means between group control and group treatment is not equal to 0
## 95 percent confidence interval:
##  -0.82571065  0.05352049
## sample estimates:
##   mean in group control mean in group treatment 
##               0.2252363               0.6113314

Ca se lit de façon très similaire à ce qu’on avait précédemment, les différences étant :

1. L’IC à 95% correspond à celui de la différence des deux moyennes.
2. Les estimations ponctuelles correspondent aux deux moyennes.

Que concluons-nous ???

T test : option “apparié”

Je rajoute une section à propos de cette option d’appariement, qui bien souvent embrouille plus qu’elle n’aide.

On possède des observations appariées, typiquement la même quantité mesurée “avant” \((X_i)\) et “après” \((Y_i)\) intervention sur chaque individu. Par exemple, le poids du patient avant intervention et après intervention. Y a t’il un changement de poids du patient avant et après intervention ?

Cela revient à tester : * H0 : la moyenne des \((Y_i - X_i)\) vaut zero. * H1 : la moyenne des \((Y_i - X_i)\) est différente de zéro.

C’est donc totalement similaire au premier test qu’on avait appliqué directement sur la différence de poids.

La façon de faire ça en R est la suivante (et on vérifie que ça donne bien le même résultat) :

t.test(x=d$weight_tf, y=d$weight_t0, paired=TRUE)

## 
##  Paired t-test
## 
## data:  d$weight_tf and d$weight_t0
## t = 3.7315, df = 499, p-value = 0.0002122
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.1980453 0.6385225
## sample estimates:
## mean difference 
##       0.4182839

t.test(d$diff_weight)

## 
##  One Sample t-test
## 
## data:  d$diff_weight
## t = 3.7315, df = 499, p-value = 0.0002122
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.1980453 0.6385225
## sample estimates:
## mean of x 
## 0.4182839

ANOVA : comparaison de moyennes

On souhaite plus précisément tester ici la différence d’au moins une moyenne parmi plusieurs groupes.

On dispose cette fois d’observations de la même quantité $(X_{i,j}) chez l’individu \(i\) dans le groupe \(j\), dans un nombre de groupes potentiellement supérieur à 2. Est-ce que la moyenne d’au moins un groupe diffère des autres groupes ?

Par exemple, ici, la différence de poids moyenne change-t-elle dans les trois classes d’âge ? (< 40, 40-60, > 60)

On teste formellement :

• H0 : les moyennes des \((X_{i,1})\), des \((X_{i,2})\) et des \((X_{i,3})\) sont identiques.
• H1 : au moins une des moyennes diffère.

Conditions d’application :

• les variables aléatoires sont indépendantes suivent une loi normale.
• il y a égalité de variance entre les différents groupes.

Une fois n’est pas coutûme, commençons par représenter les données :

m <- d %>% group_by(age_group) %>% summarise(mean=mean(diff_weight)) %>% select(mean)

plot_diff_weight_by_age_group <- d %>% ggplot(aes(x=diff_weight, fill=as.factor(age_group))) +
    geom_histogram(alpha=0.5, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois") +
    scale_fill_manual(values=c(bleufonce, bleuclair, rose)) +
    geom_vline(xintercept = c(m$mean[1], m$mean[2], m$mean[3]), color=c(bleufonce, bleuclair, rose))
plot_diff_weight_by_age_group

Les conditions d’application semblent respectées (au moins graphiquement).

Un résultat théorique certifie cette fois qu’une certaine statistique de test, calculée à partir des effectifs des différents groupes et des variances intra-groupes / inter-groupes, est distribuée suivant une loi de Fisher avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

res <- aov(diff_weight ~ age_group, data=d)
summary(res)

##              Df Sum Sq Mean Sq F value Pr(>F)
## age_group     2    0.4   0.177   0.028  0.972
## Residuals   497 3134.7   6.307

Les deux points les plus intéressants de l’output sont :

1. la “F value” qui est la valeur de la statistique de test,
2. et la p-value : la probabilité sous H0 d’obtenir une statistique de test au moins aussi extrême qu’observée.

Que concluons-nous ???

Finissons par un petit test, permettant de se convaincre qu’une ANOVA avec deux groupes correspond à un T test.

res <- aov(diff_weight ~ group, data=d)
summary(res)

##              Df Sum Sq Mean Sq F value Pr(>F)  
## group         1   18.6  18.634   2.978  0.085 .
## Residuals   498 3116.5   6.258                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

t.test(diff_weight ~ group, data=d)

## 
##  Welch Two Sample t-test
## 
## data:  diff_weight by group
## t = -1.7256, df = 495.02, p-value = 0.08505
## alternative hypothesis: true difference in means between group control and group treatment is not equal to 0
## 95 percent confidence interval:
##  -0.82571065  0.05352049
## sample estimates:
##   mean in group control mean in group treatment 
##               0.2252363               0.6113314

Test de corrélation de Pearson

On dispose cette fois d’observations d’une variables \((Y_i)\) d’intérêt, par exemple encore une fois la différence de poids. Et d’une autre variable quantitative \((X_i)\) correspondant à quelque chose qui pourrait être corrélé à nos \((Y_i)\), par exemple ici le temps pendant lequel le patient a strictement observé son traitement (compliance).

Est-ce que le gain de poids dépend de la compliance, dans le groupe traité ?

Le coefficient de corrélation \(\rho\) correspond à la covariance des deux quantités normalisée par les écarts-types des deux quantités. Il est compris entre -1 et 1 et renseigne sur la dispersion des données en 2D. En cas de corrélation linéaire parfaite, il vaut 1 ou -1.

Formellement on peut tester, avec le test de Bravais-Pearson :

• H0 : \(\rho = 0\),
• H1 : \(\rho \neq 0\).

Conditions d’application : les \((Y_i, X_i)\) doivent être appariées et distribuées normalement.

Commençons par représenter les données :

plot_diff_weight_vs_compliance <- d %>% ggplot(aes(y=diff_weight, x=t_compliance, color = group)) +
    geom_point() +
    geom_smooth(method="lm", se=F, formula="y ~ x") +
    theme_bw() +
    xlab("t_compliance") + 
    ylab("différence de poids à 6 mois") +
    scale_color_manual(values=c(bleuclair, rose))
plot_diff_weight_vs_compliance

On peut à présent effectuer un test de corrélation en prenant tous les points ensemble :

cor.test(d$diff_weight, d$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d$diff_weight and d$t_compliance
## t = 1.5369, df = 498, p-value = 0.125
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01909837  0.15546122
## sample estimates:
##        cor 
## 0.06870729

Puis un test de corrélation en prenant uniquement le groupe control :

d_control <- d %>% filter(group == "control")
cor.test(d_control$diff_weight, d_control$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d_control$diff_weight and d_control$t_compliance
## t = -1.2745, df = 248, p-value = 0.2037
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.20270677  0.04383722
## sample estimates:
##         cor 
## -0.08066852

Et enfin, le même test de corrélation en prenant uniquement le groupe treatment :

d_treatment <- d %>% filter(group == "treatment")
cor.test(d_treatment$diff_weight, d_treatment$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d_treatment$diff_weight and d_treatment$t_compliance
## t = 2.923, df = 248, p-value = 0.003787
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.05977989 0.29977312
## sample estimates:
##       cor 
## 0.1824934

Que concluons-nous ???

Le modèle linéaire

Le modèle linéaire a le bon goût d’unifier tous ces tests dans un unique framework.

Il consiste à modéliser une variable de réponse \(Y_i\) chez l’individu \(i\) comme une fonction linéaire de nos variables explicatrices \(X_{i,j}\), à laquelle on rajoute un bruit Gaussien centré. Une fonction linéaire n’est rien de plus qu’une somme pondérée de nos variables \(X_{i,j}\), soit :

\[ Y_i = a_1 X_{i,1} + a_2 X_{i,2} + ... + a_J X_{i,J} + \epsilon_i, ~~~\text{ avec } \epsilon_i \sim \mathcal{N}(0, \sigma^2) \]

On suppose que tous les \((Y_i)\) sont indépendantes et suivent la même loi. Le but du jeu, pour le statisticien, consiste ensuite à déterminer les valeurs des coefficients \(a_j\) grâce aux observations.

Bonne nouvelle : la syntaxe de R pour fitter un modèle linéaire est assez rapide et intuitive.

Un simple intercept

Premier cas simple : on aimerait modéliser la différence de poids comme étant uniquement une moyenne + un petit bruit Gaussien, et calculer cette moyenne comme on l’avait fait lors du premier test de Student.

res1 <- lm(formula="diff_weight ~ 1", data=d)
summary(res1)

## 
## Call:
## lm(formula = "diff_weight ~ 1", data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8317 -1.7631 -0.0587  1.6520  7.6847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.4183     0.1121   3.731 0.000212 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.507 on 499 degrees of freedom

On vérifie qu’on retombe bien sur la même chose que précédemment !

Un effet de variable discrète

Deuxième cas simple : on aimerait modéliser la différence de poids comme étant la somme d’un coefficient particulier chez les individus du groupe “control” et un autre coefficient chez les individus du groupe “treatment”. Soit exactement ce qu’on avait fait avec un test de Student comparant les deux moyennes.

Il se trouve que R, lorsqu’on lui donne à manger une variable à plusieurs modalités, va automatiquement la “scinder” en sous-variables indicatrices de chaque modalité. La syntaxe reste donc très simple :

# première version avec moyenne commune aux deux groupes
# et un petit ajout pour le groupe traitement
res2 <- lm(formula="diff_weight ~ group", data=d)
summary(res2)

## 
## Call:
## lm(formula = "diff_weight ~ group", data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.639 -1.705 -0.110  1.602  7.878 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)      0.2252     0.1582   1.424    0.155  
## grouptreatment   0.3861     0.2237   1.726    0.085 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.502 on 498 degrees of freedom
## Multiple R-squared:  0.005944,   Adjusted R-squared:  0.003947 
## F-statistic: 2.978 on 1 and 498 DF,  p-value: 0.08505

# deuxième version avec moyenne dans chaque groupe
res3 <- lm(formula="diff_weight ~ 0 + group", data=d)
summary(res3)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group", data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.639 -1.705 -0.110  1.602  7.878 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## groupcontrol     0.2252     0.1582   1.424 0.155185    
## grouptreatment   0.6113     0.1582   3.864 0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.502 on 498 degrees of freedom
## Multiple R-squared:  0.03293,    Adjusted R-squared:  0.02904 
## F-statistic: 8.478 on 2 and 498 DF,  p-value: 0.0002394

On a donc exactement le même modèle sous-jacent, avec les même estimations de paramètres. Cependant, ici, les p-value indiquées correspondent à un test d’adéquation du paramètre à 0. Si on souhaite comparer les deux moyennes et retomber sur le même résultat qu’au début du document, il faut donc passer par la première version et regarder le test d’adéquation du paramètre à 0.

Troisième possibilité, notre différence de poids dépend d’une variable prenant 3 modalités ou plus (comme on avait proposé lors de la réalisation du test ANOVA ci-dessus). La même syntaxe fonctionne toujours aussi bien :

res4 <- lm(formula="diff_weight ~ 0 + age_group", data=d)
summary(res4)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + age_group", data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8008 -1.7469 -0.0475  1.6451  7.6861 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## age_groupentre 40 et 60   0.4499     0.1862   2.417   0.0160 *
## age_groupmoins de 40      0.3873     0.1872   2.069   0.0390 *
## age_groupplus de 60       0.4169     0.2138   1.950   0.0517 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.511 on 497 degrees of freedom
## Multiple R-squared:  0.02726,    Adjusted R-squared:  0.02138 
## F-statistic: 4.642 on 3 and 497 DF,  p-value: 0.003269

L’output standard du modèle linéaire ne nous calcule des p-values que pour les tests d’adéquation de chaque paramètre à 0. Bien que l’ANOVA soit effectivement un test réalisé au sein du même framework de modèle linéaire, la p-value associée n’est récupérable, en pratique avec R, que via la fonction aov vue précédemment.

Un effet de variable continue

Dernier exemple similaire à ce qu’on avait fait précédemment : expliquer notre variable de réponse “différence de poids” par le “temps de compliance”.

res5 <- lm(formula="diff_weight ~ t_compliance", data=d)
summary(res5)

## 
## Call:
## lm(formula = "diff_weight ~ t_compliance", data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.048 -1.772 -0.077  1.651  7.715 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.251948   0.155707   1.618    0.106
## t_compliance 0.008513   0.005539   1.537    0.125
## 
## Residual standard error: 2.503 on 498 degrees of freedom
## Multiple R-squared:  0.004721,   Adjusted R-squared:  0.002722 
## F-statistic: 2.362 on 1 and 498 DF,  p-value: 0.125

Vers plus de complexification

Et si on souhaite un effet additif du groupe et du temps de compliance :

res6 <- lm(formula="diff_weight ~ 0 + group + t_compliance", data=d)
summary(res6)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group + t_compliance", data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8461 -1.6987 -0.1009  1.6003  7.8840 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## groupcontrol   0.101159   0.183925   0.550   0.5826  
## grouptreatment 0.447466   0.201003   2.226   0.0265 *
## t_compliance   0.007368   0.005581   1.320   0.1874  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.5 on 497 degrees of freedom
## Multiple R-squared:  0.03631,    Adjusted R-squared:  0.03049 
## F-statistic: 6.242 on 3 and 497 DF,  p-value: 0.0003639

Ce n’est pas exactement ce qu’on avait fait tout à l’heure sur le graphe : on avait alors deux régressions linéaires, une dans chaque groupe. R possède encore une fois une syntaxe pratique pour ça :

res7 <- lm(formula="diff_weight ~ 0 + group * t_compliance", data=d)
summary(res7)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group * t_compliance", data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3082 -1.6814 -0.1062  1.5706  7.8679 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)   
## groupcontrol                 0.422840   0.216077   1.957  0.05092 . 
## grouptreatment               0.169405   0.223178   0.759  0.44818   
## t_compliance                -0.011734   0.008814  -1.331  0.18370   
## grouptreatment:t_compliance  0.031605   0.011337   2.788  0.00551 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.483 on 496 degrees of freedom
## Multiple R-squared:  0.05117,    Adjusted R-squared:  0.04352 
## F-statistic: 6.688 on 4 and 496 DF,  p-value: 3.011e-05

On a introduit super facilement ce qu’on appelle des interactions entre nos variables. Le modèle linéaire, et la syntaxe de R, sont donc des outils de modélisation extrêmement puissants/versatiles. Ils demandent à bien réfléchir en amont à ce qu’on souhaite étudier pour éviter une trop grande multiplication de tests statistiques.

A vous de proposer : quel modèle linéaire vous semble le plus intéressant pour “diff_weight” ?

Conclusion

Nous avons vu aujourd’hui un certain nombre de méthodes permettant d’analyser des variables continues. Les points principaux à retenir sont les suivants :

1. Pour comparer deux moyennes, ou une moyenne à une norme, le T test est la bonne option, à condition que les données soient en grand nombre, ou qu’elles soient distribuées normalement.
2. Pour comparer plus que deux moyennes, une ANOVA fait l’affaire, sous des conditions identiques.
3. Pour quantifier la corrélation entre deux variables continues distribuées normalement, un test de corrélation de Pearson peut être utilisé.
4. Le modèle linéaire unifie tous les tests précédents dans un unique framework, où une variable de réponse continue \(Y\) est modélisée comme une combinaison linéaire de variables explicatives \(X\) (discrètes ou continues).

Notez pour finir qu’il existe d’autres types de régressions bien adaptées à des types de données différentes. Parmi les plus célèbres, vous rencontrerez très certainement la régression logistique, pour des données binaires. Vous rencontrerez éventuellement également la régression softmax, pour des données catégorielles, ou la régression de Poisson, pour des données de comptage d’événements.