Mise en place

## Chargement des données et des packages ----

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.0.4     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(labelled)
library(questionr)
library(gtsummary)
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(effects)
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(ggeffects)

theme_gtsummary_language("fr", decimal.mark = ",", big.mark = " ")
## Setting `language: fr` theme
theme_gtsummary_mean_sd()

data(hdv2003)

## Étiquettes de variables ------

hdv2003 <- hdv2003 %>%
  set_variable_labels(
    sport = "Pratique du sport ?",
    age = "Âge",
    sexe = "Sexe",
    nivetud = "Niveau d'études",
    relig = "Rapport à la religion",
    heures.tv = "Heures quotidiennes de TV"
  )

## Recodages -------

hdv2003$groupe_age <- hdv2003$age %>% 
  cut(
    include.lowest = TRUE,
    right = FALSE,
    dig.lab = 4,
    breaks = c(18, 25, 45, 60, 97)
  ) %>% 
  fct_recode(
    "18-24" = "[18,25)",
    "25-44" = "[25,45)",
    "45-59" = "[45,60)",
    "60 et +" = "[60,97]"
  )
var_label(hdv2003$groupe_age) <- "Groupe d'âges"

hdv2003$sexe <- hdv2003$sexe %>%
  fct_relevel("Femme")

hdv2003$etudes <- hdv2003$nivetud %>% 
  fct_recode(
    "Primaire" = "N'a jamais fait d'etudes",
    "Primaire" = "A arrete ses etudes, avant la derniere annee d'etudes primaires",
    "Primaire" = "Derniere annee d'etudes primaires",
    "Secondaire" = "1er cycle",
    "Secondaire" = "2eme cycle",
    "Technique/Professionnel" = "Enseignement technique ou professionnel court",
    "Technique/Professionnel" = "Enseignement technique ou professionnel long",
    "Supérieur" = "Enseignement superieur y compris technique superieur"
  ) %>%
  fct_explicit_na("Manquant")
var_label(hdv2003$etudes) <- "Niveau d'études"

Analyse uni et bivariée

## Analyse univariée --------

hdv2003 %>%
  select(sport, groupe_age, sexe, etudes, relig, heures.tv) %>%
  tbl_summary(
    statistic = all_categorical() ~ "{p}% [{n}]",
    digits = all_categorical() ~ c(1, 0)
  )
Caractéristique N = 2 0001
Pratique du sport ?
Non 63,8% [1 277]
Oui 36,1% [723]
Groupe d'âges
18-24 8,5% [169]
25-44 35,3% [706]
45-59 30,4% [609]
60 et + 25,8% [516]
Sexe
Femme 55,0% [1 101]
Homme 45,0% [899]
Niveau d'études
Primaire 23,3% [466]
Secondaire 19,4% [387]
Technique/Professionnel 29,7% [594]
Supérieur 22,0% [441]
Manquant 5,6% [112]
Rapport à la religion
Pratiquant regulier 13,3% [266]
Pratiquant occasionnel 22,1% [442]
Appartenance sans pratique 38,0% [760]
Ni croyance ni appartenance 20,0% [399]
Rejet 4,7% [93]
NSP ou NVPR 2,0% [40]
Heures quotidiennes de TV 2,25 (1,78)
Manquant 5

1 % [n]; Moyenne (ET)

## Analyse bivariée --------

hdv2003 %>%
  select(sport, groupe_age, sexe, etudes, relig, heures.tv) %>%
  tbl_summary(
    by = "sport",
    percent = "row",
    statistic = all_categorical() ~ "{p}% [{n}]",
    digits = all_categorical() ~ c(1, 0)
  ) %>%
  add_overall(last = TRUE) %>%
  add_p()
Caractéristique Non, N = 1 2771 Oui, N = 7231 Total, N = 2 0001 p-value2
Groupe d'âges <0,001
18-24 34,3% [58] 65,7% [111] 100,0% [169]
25-44 50,8% [359] 49,2% [347] 100,0% [706]
45-59 72,2% [440] 27,8% [169] 100,0% [609]
60 et + 81,4% [420] 18,6% [96] 100,0% [516]
Sexe <0,001
Femme 67,8% [747] 32,2% [354] 100,0% [1 101]
Homme 59,0% [530] 41,0% [369] 100,0% [899]
Niveau d'études <0,001
Primaire 89,3% [416] 10,7% [50] 100,0% [466]
Secondaire 69,8% [270] 30,2% [117] 100,0% [387]
Technique/Professionnel 63,6% [378] 36,4% [216] 100,0% [594]
Supérieur 42,2% [186] 57,8% [255] 100,0% [441]
Manquant 24,1% [27] 75,9% [85] 100,0% [112]
Rapport à la religion 0,14
Pratiquant regulier 68,4% [182] 31,6% [84] 100,0% [266]
Pratiquant occasionnel 66,7% [295] 33,3% [147] 100,0% [442]
Appartenance sans pratique 62,2% [473] 37,8% [287] 100,0% [760]
Ni croyance ni appartenance 59,9% [239] 40,1% [160] 100,0% [399]
Rejet 64,5% [60] 35,5% [33] 100,0% [93]
NSP ou NVPR 70,0% [28] 30,0% [12] 100,0% [40]
Heures quotidiennes de TV 2,48 (1,90) 1,83 (1,43) 2,25 (1,78) <0,001
Manquant 2 3 5

1 % [n]; Moyenne (ET)

2 Pearson's Chi-squared test; Welch Two Sample t-test

hdv2003 %>%
  ggbivariate(
    outcome = "sport",
    explanatory = c("groupe_age", "sexe", "etudes", "relig", "heures.tv")
  )
## Warning: Removed 5 rows containing non-finite values (stat_boxplot).

Régression logistique

## Calcul de la régression logistique -----

mod <- glm(
  sport ~ groupe_age + sexe + etudes + relig + heures.tv,
  family = binomial(),
  data = hdv2003
)

mod %>% 
  tbl_regression(
    intercept = TRUE,
    exponentiate = TRUE,
    add_estimate_to_reference_row = TRUE
  )
Caractéristique OR1 95% CI1 p-value
(Intercept) 0,44 0,23 – 0,82 0,011
Groupe d'âges
18-24 1,00
25-44 0,66 0,42 – 1,03 0,069
45-59 0,35 0,21 – 0,55 <0,001
60 et + 0,27 0,16 – 0,45 <0,001
Sexe
Femme 1,00
Homme 1,56 1,27 – 1,92 <0,001
Niveau d'études
Primaire 1,00
Secondaire 2,62 1,79 – 3,87 <0,001
Technique/Professionnel 2,90 2,01 – 4,23 <0,001
Supérieur 6,76 4,65 – 9,96 <0,001
Manquant 8,84 4,67 – 17,0 <0,001
Rapport à la religion
Pratiquant regulier 1,00
Pratiquant occasionnel 0,98 0,68 – 1,42 >0,9
Appartenance sans pratique 0,99 0,71 – 1,40 >0,9
Ni croyance ni appartenance 0,80 0,55 – 1,18 0,3
Rejet 0,68 0,39 – 1,19 0,2
NSP ou NVPR 0,92 0,40 – 2,03 0,8
Heures quotidiennes de TV 0,89 0,83 – 0,95 <0,001

1 OR = rapport de cotes, CI = intervalle de confiance

ggcoef_model(mod, exponentiate = TRUE)

## Effets marginaux ------

mod %>% allEffects() %>% plot()

# plot(allEffects(mod))

ggeffect(mod) %>% 
  plot() %>%
  cowplot::plot_grid(plotlist = .)

## Effet global de chaque variable ------

mod %>% 
  tbl_regression(
    intercept = TRUE,
    exponentiate = TRUE,
    add_estimate_to_reference_row = TRUE
  ) %>%
  add_global_p(keep = TRUE)
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## add_global_p: Global p-values for variable(s) `add_global_p(include =
## c("(Intercept)", "groupe_age", "sexe", "etudes", "relig", "heures.tv"))` were
## calculated with
##   `car::Anova(x$model_obj, type = "III")`
Caractéristique OR1 95% CI1 p-value
(Intercept) 0,44 0,23 – 0,82 0,011
Groupe d'âges <0,001
18-24 1,00
25-44 0,66 0,42 – 1,03 0,069
45-59 0,35 0,21 – 0,55 <0,001
60 et + 0,27 0,16 – 0,45 <0,001
Sexe <0,001
Femme 1,00
Homme 1,56 1,27 – 1,92 <0,001
Niveau d'études <0,001
Primaire 1,00
Secondaire 2,62 1,79 – 3,87 <0,001
Technique/Professionnel 2,90 2,01 – 4,23 <0,001
Supérieur 6,76 4,65 – 9,96 <0,001
Manquant 8,84 4,67 – 17,0 <0,001
Rapport à la religion 0,5
Pratiquant regulier 1,00
Pratiquant occasionnel 0,98 0,68 – 1,42 >0,9
Appartenance sans pratique 0,99 0,71 – 1,40 >0,9
Ni croyance ni appartenance 0,80 0,55 – 1,18 0,3
Rejet 0,68 0,39 – 1,19 0,2
NSP ou NVPR 0,92 0,40 – 2,03 0,8
Heures quotidiennes de TV 0,89 0,83 – 0,95 <0,001

1 OR = rapport de cotes, CI = intervalle de confiance

## Sélection pas à pas descendante ------

mod2 <- step(mod)
## Start:  AIC=2236.81
## sport ~ groupe_age + sexe + etudes + relig + heures.tv
## 
##              Df Deviance    AIC
## - relig       5   2211.1 2231.1
## <none>            2206.8 2236.8
## - heures.tv   1   2219.9 2247.9
## - sexe        1   2224.4 2252.4
## - groupe_age  3   2259.0 2283.0
## - etudes      4   2334.7 2356.7
## 
## Step:  AIC=2231.13
## sport ~ groupe_age + sexe + etudes + heures.tv
## 
##              Df Deviance    AIC
## <none>            2211.1 2231.1
## - heures.tv   1   2224.5 2242.5
## - sexe        1   2227.3 2245.3
## - groupe_age  3   2260.6 2274.6
## - etudes      4   2339.0 2351.0
mod2 %>% ggcoef_model()

ggcoef_compare(
  list("modèle complet" = mod, "modèle simplifié" = mod2),
  exponentiate = TRUE,
  type = "f"
)

t1 <- mod %>% 
  tbl_regression(exponentiate = TRUE, add_estimate_to_reference_row = TRUE) %>% 
  add_global_p()
## add_global_p: Global p-values for variable(s) `add_global_p(include =
## c("groupe_age", "sexe", "etudes", "relig", "heures.tv"))` were calculated with
##   `car::Anova(x$model_obj, type = "III")`
t2 <- mod2 %>% 
  tbl_regression(exponentiate = TRUE, add_estimate_to_reference_row = TRUE) %>% 
  add_global_p()
## add_global_p: Global p-values for variable(s) `add_global_p(include =
## c("groupe_age", "sexe", "etudes", "heures.tv"))` were calculated with
##   `car::Anova(x$model_obj, type = "III")`
tbl_merge(
  tbls = list(t1, t2),
  tab_spanner = c("**Modèle complet**", "**Modèle simplifié**")
)
Caractéristique Modèle complet Modèle simplifié
OR1 95% CI1 p-value OR1 95% CI1 p-value
Groupe d'âges <0,001 <0,001
18-24 1,00 1,00
25-44 0,66 0,42 – 1,03 0,68 0,44 – 1,06
45-59 0,35 0,21 – 0,55 0,36 0,23 – 0,58
60 et + 0,27 0,16 – 0,45 0,29 0,18 – 0,48
Sexe <0,001 <0,001
Femme 1,00 1,00
Homme 1,56 1,27 – 1,92 1,52 1,24 – 1,87
Niveau d'études <0,001 <0,001
Primaire 1,00 1,00
Secondaire 2,62 1,79 – 3,87 2,57 1,75 – 3,79
Technique/Professionnel 2,90 2,01 – 4,23 2,86 1,99 – 4,16
Supérieur 6,76 4,65 – 9,96 6,68 4,60 – 9,83
Manquant 8,84 4,67 – 17,0 8,78 4,65 – 16,9
Rapport à la religion 0,5
Pratiquant regulier 1,00
Pratiquant occasionnel 0,98 0,68 – 1,42
Appartenance sans pratique 0,99 0,71 – 1,40
Ni croyance ni appartenance 0,80 0,55 – 1,18
Rejet 0,68 0,39 – 1,19
NSP ou NVPR 0,92 0,40 – 2,03
Heures quotidiennes de TV 0,89 0,83 – 0,95 <0,001 0,89 0,83 – 0,95 <0,001

1 OR = rapport de cotes, CI = intervalle de confiance

## Y a-t-il un effet d'interaction entre l'âge et le sexe ?

mod3 <- glm(
  sport ~ sexe * groupe_age + etudes + heures.tv, 
  # sport ~ sexe + groupe_age + sexe:groupe_age + etudes + heures.tv
  family = binomial, data = hdv2003
)

mod3 %>% tbl_regression() %>% add_global_p()
## add_global_p: Global p-values for variable(s) `add_global_p(include = c("sexe",
## "groupe_age", "etudes", "heures.tv", "sexe:groupe_age"))` were calculated with
##   `car::Anova(x$model_obj, type = "III")`
Caractéristique log(OR)1 95% CI1 p-value
Sexe <0,001
Femme
Homme 1,6 0,88 – 2,4
Groupe d'âges 0,002
18-24
25-44 0,06 -0,49 – 0,62
45-59 -0,43 -1,0 – 0,17
60 et + -0,63 -1,3 – 0,00
Niveau d'études <0,001
Primaire
Secondaire 0,95 0,57 – 1,3
Technique/Professionnel 1,1 0,70 – 1,4
Supérieur 1,9 1,5 – 2,3
Manquant 2,2 1,6 – 2,9
Heures quotidiennes de TV -0,12 -0,18 – -0,05 <0,001
Sexe * Groupe d'âges 0,005
Homme * 25-44 -1,2 -2,0 – -0,36
Homme * 45-59 -1,4 -2,3 – -0,60
Homme * 60 et + -1,4 -2,4 – -0,56

1 OR = rapport de cotes, CI = intervalle de confiance

ggcoef_model(mod3)

plot(allEffects(mod3))

mod3 %>%
  ggeffect(c("groupe_age", "sexe")) %>%
  plot()

mod4 <- glm(
  sport ~ sexe:groupe_age + etudes + heures.tv, 
  family = binomial, data = hdv2003
)
ggcoef_model(mod4)

mod4 %>%
  ggeffect(c("groupe_age", "sexe")) %>%
  plot()

mod5 <- glm(
  sport ~ groupe_age + sexe:groupe_age + etudes + heures.tv, 
  family = binomial, data = hdv2003
)
ggcoef_model(mod5)

mod5 %>%
  ggeffect(c("groupe_age", "sexe")) %>%
  plot()

## Interaction sexe et niveau d'étude -----


mod6 <- glm(
  sport ~ groupe_age + sexe * etudes + heures.tv, 
  family = binomial, data = hdv2003
)
mod6 %>%
  ggeffect(c("etudes", "sexe")) %>%
  plot()

car::Anova(mod6)
## Analysis of Deviance Table (Type II tests)
## 
## Response: sport
##             LR Chisq Df Pr(>Chisq)    
## groupe_age    48.744  3  1.479e-10 ***
## sexe          16.217  1  5.649e-05 ***
## etudes       127.899  4  < 2.2e-16 ***
## heures.tv     12.970  1  0.0003165 ***
## sexe:etudes   10.894  4  0.0277790 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Y a-t-il un risque de multicolinéarité ?

car::vif(mod)
##                GVIF Df GVIF^(1/(2*Df))
## groupe_age 1.844390  3        1.107411
## sexe       1.050186  1        1.024786
## etudes     1.832560  4        1.078654
## relig      1.124563  5        1.011809
## heures.tv  1.068119  1        1.033498
# add_vif() sera disponible dans lr pcohaine version de gtsummary