POL6022 - Méthodes quantitatives avancées

3.1 L’idée

Bootstrap = incertitude de l’estimation

• On observe un échantillon de taille \(n\) : \((x_1, \ldots, x_n)\)
• On veut une approximation de la distribution d’une statistique \(\widehat{\theta}\) (ex. \(\bar{X}\), médiane, coefficient de régression)

Procédure (bootstrap non paramétrique) :

1. Tirer un échantillon bootstrap de taille \(n\) avec remise à partir des données observées
2. Calculer \(\widehat{\theta}^*\) sur cet échantillon
3. Répéter \(B\) fois pour obtenir \(\widehat{\theta}^*_1, \ldots, \widehat{\theta}^*_B\)

À quoi ça sert ?

• Estimer l’erreur-type d’un estimateur : \(\widehat{SE}_{boot}(\widehat{\theta}) = sd(\widehat{\theta}^*_1, \ldots, \widehat{\theta}^*_B)\)
• Construire des intervalles de confiance
• Visualiser la variabilité de \(\widehat{\theta}\)

3.2 Application : Incarcération et participation électorale

Question de recherche : Les ex-détenu·e·s votent-ils moins que les autres ?

Nous utilisons notre jeu de données real_data pour estimer cet effet, puis le bootstrap pour quantifier l’incertitude.

# Étape 1 : Régression sur les données observées
modele_obs <- lm(voted ~ incarcerated, data = real_data)
summary(modele_obs)

## 
## Call:
## lm(formula = voted ~ incarcerated, data = real_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7388 -0.5041  0.2612  0.2612  0.4959 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.738811   0.002177   339.4   <2e-16 ***
## incarcerated -0.234733   0.005629   -41.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4489 on 49998 degrees of freedom
## Multiple R-squared:  0.03362,    Adjusted R-squared:  0.0336 
## F-statistic:  1739 on 1 and 49998 DF,  p-value: < 2.2e-16

L’incarcération est associée à une baisse de la probabilité de voter de 23.5 en moyenne.

Maintenant, utilisons le bootstrap pour estimer l’incertitude de ce coefficient.

set.seed(123)

# Coefficient observé sur les données complètes
modele_obs <- lm(voted ~ incarcerated, data = real_data)
coef_obs <- coef(modele_obs)[2]

# Paramètres bootstrap
B <- 1000  
bootstrap_coefs <- numeric(B)

# Boucle bootstrap
for (i in 1:B) {
  # Ré-échantillonner les observations avec remplacement
  indices <- sample(1:nrow(real_data), size = nrow(real_data), replace = TRUE)
  data_boot <- real_data[indices, ]
  
  # Estimer le modèle sur l'échantillon bootstrap
  modele_boot <- lm(voted ~ incarcerated, data = data_boot)
  bootstrap_coefs[i] <- coef(modele_boot)[2]
}

# Calcul de l'IC à 95%
ic_95 <- quantile(bootstrap_coefs, probs = c(0.025, 0.975))
se_boot <- sd(bootstrap_coefs)

# Erreur-type bootstrap
round(se_boot, 5)

## [1] 0.00615

# IC 95% 
round(ic_95[1], 4)

##    2.5% 
## -0.2474

round(ic_95[2], 4)

##   97.5% 
## -0.2232

hist(bootstrap_coefs,
     breaks = 40,
     col = "#ec8cac",
     xlab = "Coefficient estimé (effet de l'incarcération sur le vote)",
     main = "Distribution bootstrap du coefficient",
     xlim = c(-0.25, -0.15))
abline(v = coef_obs, col = "red", lwd = 3, lty = 2)
abline(v = ic_95, col = "blue", lwd = 2, lty = 3)
legend("topright", 
       legend = c("Coefficient observé", "IC 95%"),
       col = c("red", "blue"), 
       lwd = c(3, 2), 
       lty = c(2, 3))

On bootstrap le coefficient de régression en ré-échantillonnant les observations pour approximer sa distribution d’échantillonnage, c’est-à-dire voir à quel point le coefficient de régression varierait si on refaisait l’étude, et donc quantifier son incertitude.

# Étape 1 : Calculer les écarts-types observés
sd_incarc <- sd(real_data$voted[real_data$incarcerated == 1])
sd_non_incarc <- sd(real_data$voted[real_data$incarcerated == 0])
ratio_obs <- sd_incarc / sd_non_incarc

# Étape 2 : Bootstrap pour l'IC
# Fixez le seed à 123
set.seed(123)

# Votre code bootstrap ici :
# 1. Créez 2000 échantillons bootstrap
# 2. Pour chaque échantillon, calculez le ratio des écarts-types
# 3. Calculez l'IC à 95%

# Réponse

set.seed(789)

B <- 2000
bootstrap_ratios <- numeric(B)

for (i in 1:B) {
  # Ré-échantillonner
  indices <- sample(1:nrow(real_data), size = nrow(real_data), replace = TRUE)
  data_boot <- real_data[indices, ]
  
  # Calculer les écarts-types
  sd_incarc_boot <- sd(data_boot$voted[data_boot$incarcerated == 1])
  sd_non_incarc_boot <- sd(data_boot$voted[data_boot$incarcerated == 0])
  bootstrap_ratios[i] <- sd_incarc_boot / sd_non_incarc_boot
}

# IC à 95%
ic_95 <- quantile(bootstrap_ratios, probs = c(0.025, 0.975))
cat("IC 95% pour le ratio des écarts-types :", round(ic_95[1], 4), "-", round(ic_95[2], 4), "\n")

## IC 95% pour le ratio des écarts-types : 1.1322 - 1.1443

# Visualisation
hist(bootstrap_ratios,
     breaks = 30,
     col = "#d5e8f0",
     xlab = "Ratio des écarts-types (incarcérés / non-incarcérés)",
     main = "Distribution bootstrap du ratio")
abline(v = ratio_obs, col = "red", lwd = 3, lty = 2)
abline(v = ic_95, col = "blue", lwd = 2, lty = 3)
legend("topright", 
       legend = c("Ratio observé", "IC 95%"),
       col = c("red", "blue"), 
       lwd = c(3, 2), 
       lty = c(2, 3),
       cex = 0.8)

4.1 L’idée

Monte Carlo = performance de la méthode dans un monde simulé

Question : « Si le vrai monde ressemble à XYZ, est-ce que ma procédure d’inférence est fiable ? »

• Point de départ : un DGP supposé (vraie valeur de \(\beta\), taux de participation, structure par province, etc.)
• Procédure : simuler beaucoup de jeux de données, ré-estimer \(\widehat{\beta}\), vérifier :
  • Biais : \(E[\widehat{\beta}] - \beta\)
  • Couverture : les IC à 95% contiennent-ils \(\beta\) ≈ 95% du temps ?

4.2 Application : Combien d’observations faut-il ?

Vous planifiez une étude sur l’effet de l’incarcération sur le vote. Problème : chaque observation coûte cher (accès aux données, entrevues, etc.).

Vous devez décider : Collecter 500, 1000 ou 2000 observations ?

• Trop peu → Vous ne détecterez peut-être pas l’effet même s’il existe (c’est aussi ce qu’on appelle une erreur de type II)
• Trop → Vous gaspillerez de l’argent pour rien

Monte Carlo permet de répondre avant de collecter les données.

Basé sur des études antérieures, vous estimez que l’effet de l’incarcération est environ -0.20 (20 points de % de baisse de participation dans l’aggrégé pour les personnes qui ont fait de la prison). Vous simulez des études de différentes tailles pour voir laquelle détecte cet effet de manière fiable.

set.seed(2026)

# Paramètres basés sur la littérature
effet_attendu <- -0.20
taux_incarc <- 0.15

# Tailles d'échantillon à tester
tailles <- c(500, 1000, 2000)
n_simulations <- 500

# Stockage des résultats
resultats_puissance <- data.frame(
  n = numeric(),
  puissance = numeric(),
  largeur_ic_moyenne = numeric()
)

for (n in tailles) {
  
  significatif <- numeric(n_simulations)
  largeurs_ic <- numeric(n_simulations)
  
  for (i in 1:n_simulations) {
    
    # Simuler une étude de taille n
    incarcerated <- rbinom(n, 1, taux_incarc)
    prob_vote <- ifelse(incarcerated == 1, 0.50, 0.70)
    voted <- rbinom(n, 1, prob_vote)
    
    # Analyser
    modele <- lm(voted ~ incarcerated)
    coef_est <- coef(modele)[2]
    se_est <- summary(modele)$coefficients[2, 2]
    p_value <- summary(modele)$coefficients[2, 4]
    
    # L'effet est-il significatif à 5% ?
    significatif[i] <- (p_value < 0.05)
    
    # Quelle est la largeur de l'IC ?
    ic <- coef_est + c(-1.96, 1.96) * se_est
    largeurs_ic[i] <- ic[2] - ic[1]
  }
  
  # Puissance = % d'études qui détectent l'effet
  puissance <- mean(significatif)
  largeur_moyenne <- mean(largeurs_ic)
  
  resultats_puissance <- rbind(resultats_puissance, 
                               data.frame(n = n, 
                                         puissance = puissance,
                                         largeur_ic_moyenne = largeur_moyenne))
}

Résultats :

print(resultats_puissance)

##      n puissance largeur_ic_moyenne
## 1  500     0.882          0.2292209
## 2 1000     0.998          0.1617353
## 3 2000     1.000          0.1142288

par(mfrow = c(1, 2))

# Graphique 1 : Puissance selon la taille d'échantillon
plot(resultats_puissance$n, resultats_puissance$puissance * 100,
     type = "b", pch = 19, cex = 2, lwd = 2,
     col = "#de1d5b",
     xlab = "Taille d'échantillon",
     ylab = "Puissance (%)",
     main = "Probabilité de détecter l'effet",
     ylim = c(0, 100))
abline(h = 80, col = "blue", lty = 2, lwd = 2)
text(1500, 85, "Standard : 80%", col = "blue")
grid()

# Graphique 2 : Largeur des IC
plot(resultats_puissance$n, resultats_puissance$largeur_ic_moyenne,
     type = "b", pch = 19, cex = 2, lwd = 2,
     col = "#de1d5b",
     xlab = "Taille d'échantillon",
     ylab = "Largeur moyenne de l'IC à 95%",
     main = "Précision de l'estimation")
grid()

Décision :

• n = 500 : 0.882 de puissance
• n = 1000 : 0.998 de puissance
• n = 2000 : 1.000 de puissance

On devrait collecter au moins 1000 observations pour avoir presque 100% de chances de détecter l’effet attendu, mais inutile de débourser deux fois le montant pour aller jusqu’à 2000 observations (qui n’offre qu’un bénéfice minime).

set.seed(123)

# Paramètres du DGP
n_size <- 1000
B <- 500 # Nombre d'études

taux_incarc <- 0.15
p1 <- 0.50 # P(vote | incarceré(e)=1)
p0 <- 0.70 # P(vote | incarceré(e)=0)
beta_vrai <- p1 - p0  # -0.20

# Paramètres de non-réponse (biais de sélection)
p_rep_1 <- 0.35  # P(répond | incarceré(e)=1)
p_rep_0 <- 0.95  # P(répond | incarceré(e)=0)

# Stockage des estimations dans des vecteurs vides
beta_hat_complet <- numeric(B)
beta_hat_biaise  <- numeric(B)

for (b in 1:B) { # On utilise une loop (500 itérations)

  # 1) Le vrai monde
  
  incarcerated <- rbinom(n_size, 1, taux_incarc) # Vecteur dummy pour l'incarcération
  prob_vote <- ifelse(incarcerated == 1, p1, p0) # Vecteur avec la probabilité théorique de voter basé sur la dummy incarcération
  voted <- rbinom(n_size, 1, prob_vote) # Vecteur d'une distribution binomiale avec probabilité de voter définie pour chaque sous-échantillon (incarcérés ou pas)

  # 2) Estimation sur échantillon complet (scénario idéal)
  m_full <- lm(voted ~ incarcerated)
  beta_hat_complet[b] <- coef(m_full)[2] # On sauvegarde le coefficient

  # 3) Non-réponse dépendante de l'incarcération
  prob_rep <- ifelse(incarcerated == 1,
                     p_rep_1,
                     p_rep_0) # On crée un vecteur avec prop. de répondre
  responded <- rbinom(n_size, 1, prob_rep) # On crée un vecteur avec dummy : ont-ils répondu (1) ou pas (0)

  # 4) Estimation sur échantillon observé seulement
  m_obs <- lm(voted[responded == 1] ~ incarcerated[responded == 1]) # Modèle seulement en utilisant les gens qui ont répondu
  beta_hat_biaise[b] <- coef(m_obs)[2] # On sauvegarde les résultats encore une fois
}

# Réponse

mean(beta_hat_complet) # Coefficient dans l'échantillon complet

## [1] -0.2011338

mean(beta_hat_biaise) # Coefficient dans l'échantillon biaisé

## [1] -0.1959945

mean(beta_hat_complet) - beta_vrai # Différence entre le coefficient dans l'échantillon complet et dans la "vraie vie"

## [1] -0.001133758

mean(beta_hat_biaise)  - beta_vrai  # Différence entre le coefficient dans l'échantillon biaisé et dans la "vraie vie"

## [1] 0.004005529