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 000¹
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
¹ % [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 277¹	Oui, N = 723¹	Total, N = 2 000¹	p-value²
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
¹ % [n]; Moyenne (ET) ² 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	OR¹	95% CI¹	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
¹ 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	OR¹	95% CI¹	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
¹ 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é
Caractéristique	OR¹	95% CI¹	p-value	OR¹	95% CI¹	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
¹ 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)¹	95% CI¹	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
¹ 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

Analyse de la pratique du sport (hdv2003)

Joseph Larmarange

2021-02-05

Mise en place

Analyse uni et bivariée

Régression logistique