Montrer le code
df1 <- read.csv("dat.csv")
library(ggplot2)
library(DescTools)
library(MASS)
library(viridis)
library(knitr)
library(gridExtra)Projet de Statistiques
Modélisation de sinistres et de montant de police d’assurance sur données réelles
Partie I
Question 1)
Importation des données
On importe tout d’abord les données en csv dans un dataframe df1, ainsi que les packages utilisés ultérieurement.
Elimination des points aberrants
Nous effectuons tout d’abord un scatter plot de Sinistre1 (le montant des dommages en dizaines d’euros pour les sinistres de type 1) sur Police1 (cotisation au titre du complément de police de type 1).
ggplot(df1, aes(x = Sinistre1, y = Police1)) +
geom_point(color = "skyblue", size = 1)
Nous remarquons que trois points semblent aberrants. Pour confirmer cela, nous effectuons une régression linéaire de Sinistre1 sur toutes les autres variables, puis étudions les graphiques “Residuals vs Fitted”, “Q-Q Residuals”, “Scale-Location” et “Residuals vs Leverage” de cette régression.
m1 = lm(Sinistre1 ~ ., data=df1)
plot(m1)
Nous remarquons que les points 17, 75 et 159 sont bien aberrants, nous choisissons alors de les enlever, et étudierons dans la suite du projet la dataframe ainsi mise à jour.
df1 <- df1[-c(17, 75, 159), ]Nous refaisons même régression linéaire et observons le graphique “Residuals vs Fitted”.
m1 = lm(Sinistre1 ~ ., data=df1)
plot(m1, which = 1)
Nous remarquons que les valeurs semblent être regroupées en deux blocs, pouvant alors indiquer une bimodalité de Sinistre1.
Statistiques descriptives
Nous effectuons un boxplot de Sinistre1 sur enfants, qui est une variable qualitative qui indique si le ménage a au moins un enfant ou non.
ggplot(data = df1, aes(x = enfants, y = Sinistre1)) +
geom_boxplot()
Nous remarquons des valeurs de montant du sinistre 1 plus importantes chez les ménages ayant un ou plusieurs enfants, ce qui pourrait alors donner l’hypothèse qu’avoir des enfants fait accroître le montant du sinitre 1.
Nous regardons aussi les moyennes, variances et quantiles des montants des sinistres 1, 2 et 3.
Montrer le code
stats <- data.frame(
Montant = c("Sinistre 1", "Sinistre2", "Sinistre3"),
Moyenne = c(mean(df1$Sinistre1), mean(df1$Sinistre2), mean(df1$Sinistre3)),
Variance = c(var(df1$Sinistre1), var(df1$Sinistre2), var(df1$Sinistre3)),
Quantile25 = c(quantile(df1$Sinistre1, probs = 0.25), quantile(df1$Sinistre2, probs = 0.25), quantile(df1$Sinistre3, probs = 0.25)),
Quantile50 = c(quantile(df1$Sinistre1, probs = 0.5), quantile(df1$Sinistre2, probs = 0.5), quantile(df1$Sinistre3, probs = 0.5)),
Quantile75 = c(quantile(df1$Sinistre1, probs = 0.75), quantile(df1$Sinistre2, probs = 0.75), quantile(df1$Sinistre3, probs = 0.75))
)
kable(stats, format = "markdown")| Montant | Moyenne | Variance | Quantile25 | Quantile50 | Quantile75 |
|---|---|---|---|---|---|
| Sinistre 1 | 16.1708729 | 18.452722 | 13.26607 | 17.25522 | 19.2695 |
| Sinistre2 | 0.1615954 | 1.323781 | 0.00000 | 0.00000 | 0.0000 |
| Sinistre3 | 1.8378220 | 7.476037 | 0.00000 | 0.71000 | 2.5600 |
Nous remarquons qu’une grande majorité des ménages présentent un montant nul pour le sinistre 2 (et à moindre degré le 3), indiquant qu’ils n’ont pas eu ce sinistre.
Question 2)
Nous traçons un histogramme de Sinistre1 et Sinistre2, avec les proportions en ordonnées.
plot1 <- ggplot(data = df1, aes(x = Sinistre1)) +
geom_histogram(fill = "skyblue", color = "black", aes(y = after_stat(density))) +
labs(title = "Histogramme du sinistre 1", x = "Valeurs", y = "Proportions") +
geom_density(color = "red")
plot2 <- ggplot(data = df1, aes(x = Sinistre2)) +
geom_histogram(fill = "lightgreen", color = "black", aes(y = after_stat(density))) +
labs(title = "Histogramme du sinistre 2", x = "Valeurs", y = "Proportions") +
geom_density(color = "red")
grid.arrange(plot1, plot2, nrow = 1)
Nous traçons ensuite les fonctions de répartition empiriques de Sinistre1 et Sinistre2, ainsi que les lois normales qui leur correspondent (c’est à dire la loi normale qui a la même moyenne et variance que la variable étudiée).
Montrer le code
cdf_data1 <- data.frame(x = sort(df1$Sinistre1), y = ecdf(df1$Sinistre1)(sort(df1$Sinistre1)))
normal_cdf_data1 <- data.frame(x = seq(min(df1$Sinistre1), max(df1$Sinistre1), length.out = 1000),
y = pnorm(seq(min(df1$Sinistre1), max(df1$Sinistre1), length.out = 1000), mean = mean(df1$Sinistre1), sd = sd(df1$Sinistre1)))
cdf_data2 <- data.frame(x = sort(df1$Sinistre2), y = ecdf(df1$Sinistre2)(sort(df1$Sinistre2)))
normal_cdf_data2 <- data.frame(x = seq(min(df1$Sinistre2), max(df1$Sinistre2), length.out = 1000),
y = pnorm(seq(min(df1$Sinistre2), max(df1$Sinistre2), length.out = 1000), mean = mean(df1$Sinistre2), sd = sd(df1$Sinistre2)))
plot1 <- ggplot() +
geom_line(data = normal_cdf_data1, aes(x = x, y = y, color = "Normal CDF"), linewidth = 1) +
xlab("Valeurs") +
geom_line(data = cdf_data1, aes(x = x, y = y, color = "Empirical CDF"), linewidth = 1) +
ylab("Probabilité cumulée") +
scale_color_manual(values = c("blue", "red"),
labels = c("FDR Sinistre 1", "FDR Loi Normale")) +
guides(color = guide_legend(title = "Courbes"))
plot2 <- ggplot() +
geom_line(data = normal_cdf_data2, aes(x = x, y = y, color = "Normal CDF"), linewidth = 1) +
xlab("Valeurs") +
geom_line(data = cdf_data2, aes(x = x, y = y, color = "Empirical CDF"), linewidth = 1) +
ylab("Probabilité cumulée") +
scale_color_manual(values = c("blue", "red"),
labels = c("FDR Sinistre 2", "FDR Loi Normale")) +
guides(color = guide_legend(title = "Courbes"))
grid.arrange(plot1, plot2, widths = 10)
Nous effectuons ensuite un test de normalité de Kolmogorov-Smirnov sur ces deux variables.
LillieTest(df1$Sinistre1)
Lilliefors (Kolmogorov-Smirnov) normality test
data: df1$Sinistre1
D = 0.11372, p-value < 2.2e-16LillieTest(df1$Sinistre2)
Lilliefors (Kolmogorov-Smirnov) normality test
data: df1$Sinistre2
D = 0.45209, p-value < 2.2e-16Donc nous rejettons l’hypothèse nulle du test de Kolmogorov-Smirnov, les montants des sinistres 1 et 2 ne suivent pas une répartition de loi normale d’après ce test. Nous remarquons cependant que la statistique de test D est bien plus grande pour Sinistre2 que pour Sinistre1, ce qui est concordant avec les graphiques ci-dessus, le maximum de la différence entre la répartition empirique et la répartition de la loi normale est bien plus grande pour le sinistre 2.
Partie II
Question 1)
Nous effectuons une régression linéaire de Sinistre1 sur log_RUC.
df1$log_RUC <- log(df1$RUC)
reg1 <- lm(Sinistre1 ~ log_RUC, data = df1)
summary(reg1)
Call:
lm(formula = Sinistre1 ~ log_RUC, data = df1)
Residuals:
Min 1Q Median 3Q Max
-15.9670 -2.7134 0.9209 3.0180 9.7192
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.9626 0.8853 -5.606 2.18e-08 ***
log_RUC 2.4581 0.1028 23.920 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.083 on 5347 degrees of freedom
Multiple R-squared: 0.09666, Adjusted R-squared: 0.0965
F-statistic: 572.2 on 1 and 5347 DF, p-value: < 2.2e-16Nous pouvons dire grâce au test de Fisher et à l’indication de sa p-value que le modèle est globalement significatif. Les tests de Student permettent de dire que la constante et la variable log_RUC sont significatives. Le coefficient de log_RUC est positif, une augmentation de log_RUC de 1 (et donc une augmentation de RUC d’environ 1%) fera donc augmenter Sinistre1 d’environ 2,458 selon le modèle. Le \(R^2\) vaut cependant environ 10%, la variable log_RUC n’apporte donc que peu d’information sur Sinistre1, et il existe probablement d’autres variables explicatives permettant d’apporter des informations sur Sinistre1.
Question 2)
Nous traçons les données et la droite de cette régression.
ggplot(df1, aes(x = log_RUC, y = Sinistre1)) +
geom_point(alpha = 0.7, color = "skyblue") +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(x = "log_RUC", y = "Sinistre1", title = "Régression linéaire de Sinistre1 en fonction de log_RUC")
Question 3)
Nous régressons maintenant Sinistre1 sur log_RUC et pcs, qui sont les catégories socio-professionnelles.
reg2 <- lm(Sinistre1 ~ log_RUC + pcs, data = df1)
summary(reg2)
Call:
lm(formula = Sinistre1 ~ log_RUC + pcs, data = df1)
Residuals:
Min 1Q Median 3Q Max
-15.3650 -2.8176 0.5416 2.9529 10.7965
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.1036 0.9414 -1.172 0.24115
log_RUC 2.0413 0.1101 18.549 < 2e-16 ***
pcsArtisans, comm., chefs d'ent. -0.1420 0.4682 -0.303 0.76163
pcsAutres pers. sans activite prof. 1.5255 0.4647 3.283 0.00103 **
pcsCadres et prof. intellectuelles sup. -0.5629 0.4229 -1.331 0.18324
pcsEmployes -0.7492 0.3932 -1.905 0.05679 .
pcsOuvriers -1.5922 0.3801 -4.189 2.85e-05 ***
pcsProfessions intermediaires -1.0372 0.3951 -2.625 0.00868 **
pcsRetraites 1.5851 0.3885 4.080 4.57e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.896 on 5340 degrees of freedom
Multiple R-squared: 0.1789, Adjusted R-squared: 0.1776
F-statistic: 145.4 on 8 and 5340 DF, p-value: < 2.2e-16Nous remarquons que pcsétant une variable qualitative prenant plusieurs modalités, le logiciel a automatiquement séparé cette variable en plusieures variables indicatrices. La variable log_RUC reste significative, elle apporte donc toujours des informations permettant d’expliquer Sinistre1, qui ne sont pas contenues dans la variable pcs. A noter qu’il faut cependant faire attention aux problèmes de multicolinéarité imparfaite.
Question 4)
Nous effectuons la régression de Sinistre1 sur pcs.
reg3 <- lm(Sinistre1 ~ pcs, data = df1)
summary(reg3)
Call:
lm(formula = Sinistre1 ~ pcs, data = df1)
Residuals:
Min 1Q Median 3Q Max
-14.6094 -2.9407 0.5878 2.9608 10.8390
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.0185 0.3731 40.251 < 2e-16 ***
pcsArtisans, comm., chefs d'ent. 1.1061 0.4780 2.314 0.0207 *
pcsAutres pers. sans activite prof. 2.4007 0.4769 5.034 4.97e-07 ***
pcsCadres et prof. intellectuelles sup. 1.8129 0.4158 4.360 1.32e-05 ***
pcsEmployes 0.5562 0.3991 1.394 0.1635
pcsOuvriers -0.7181 0.3891 -1.845 0.0650 .
pcsProfessions intermediaires 0.6383 0.3968 1.609 0.1078
pcsRetraites 3.3499 0.3886 8.621 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.019 on 5341 degrees of freedom
Multiple R-squared: 0.126, Adjusted R-squared: 0.1248
F-statistic: 109.9 on 7 and 5341 DF, p-value: < 2.2e-16Nous remarquons alors que des trois régressions effectuées, la régression avec log_RUC et pcs a le meilleur \(R^2\) (environ 18%). Sous l’hypothèse d’absence de multicolinéarité imparfaite entre log_RUC et pcs, ce modèle semble être le meilleur. Nous effectuons aussi un stepAIC sur cette régression.
stepAIC(reg3)Start: AIC=14888.39
Sinistre1 ~ pcs
Df Sum of Sq RSS AIC
<none> 86256 14888
- pcs 7 12430 98685 15594
Call:
lm(formula = Sinistre1 ~ pcs, data = df1)
Coefficients:
(Intercept) pcsArtisans, comm., chefs d'ent. pcsAutres pers. sans activite prof.
15.0185 1.1061 2.4007
pcsCadres et prof. intellectuelles sup. pcsEmployes pcsOuvriers
1.8129 0.5562 -0.7181
pcsProfessions intermediaires pcsRetraites
0.6383 3.3499 L’algorithme de stepAIC garde bien les deux variables, ce qui indique bien que la régression avec les deux variables explicatives est ici le meilleur modèle.
Partie III
Question 1)
Nous mettons les 10 derniers individus de df1 dans df0 et on les enlève de df1.
n <- nrow(df1)
df0 <- df1[(n - 9):n, ]
df1 <- df1[1:(n - 10), ]Nous choisissons d’étudier séparément Sinistre1 et Police1, nous créont alors les dataframes dfsin et dfpol qui excluent respectivement Police1 et Sinistre1.
df_sin <- df1[, !(names(df1) == "Police1")]
df_pol <- df1[, !(names(df1) == "Sinistre1")]Nous opérons l’algorithme stepAIC sur la régression Sinistre1 sur toutes les autres variables.
Le but de l’algorithme est de simplifier le modèle, en enlevant les variables apportant le moins d’information, et donc d’éviter les problèmes d’overfitting. L’algorithme va calculer l’AIC de chaque variable explicative, et va éliminer la variable avec le plus petit AIC. Il va ensuite recalculer tous les AIC. Si la variable avec le plus petit AIC a un AIC inférieur au minimum de l’étape précédente, il l’élimine et continue l’algorithme. Si la valeur est supérieure, l’algorithme s’arrête et on choisit l’ensemble des variables de l’étape précédentes.
La formule générale de l’AIC est la suivante : \(AIC = 2k - 2 \ln (L)\), avec \(k\) le nombre de paramètres estimés qui va pénaliser l’AIC plus il est élevé (et donc pénalise la complexité), et \(L\) le maximum de vraisemblance du modèle, qui favorise l’AIC plus il est élevé (et donc favorise la qualité d’ajustement du modèle).
reg_sin1 = lm(Sinistre1 ~ ., data = df_sin)
stepAIC(reg_sin1, trace = 0)
Call:
lm(formula = Sinistre1 ~ RUC + Acompm + NSin, data = df_sin)
Coefficients:
(Intercept) RUC AcompmCouple avec enfant(s) AcompmCouple sans enfant
16.8471907 0.0002095 -8.0108592 0.0317552
AcompmPersonne seule NSin
-0.0358729 0.0164336 Nous régressons le modèle ainsi déterminé.
reg_sin1_AIC = lm(formula = Sinistre1 ~ RUC + Acompm + NSin, data = df_sin)
summary(reg_sin1_AIC)
Call:
lm(formula = Sinistre1 ~ RUC + Acompm + NSin, data = df_sin)
Residuals:
Min 1Q Median 3Q Max
-9.4840 -1.4233 -0.0391 1.4898 7.8200
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.685e+01 8.062e-02 208.982 <2e-16 ***
RUC 2.095e-04 8.730e-06 23.995 <2e-16 ***
AcompmCouple avec enfant(s) -8.011e+00 7.593e-02 -105.509 <2e-16 ***
AcompmCouple sans enfant 3.176e-02 8.357e-02 0.380 0.7040
AcompmPersonne seule -3.587e-02 1.023e-01 -0.351 0.7259
NSin 1.643e-02 8.517e-03 1.929 0.0537 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.144 on 5333 degrees of freedom
Multiple R-squared: 0.751, Adjusted R-squared: 0.7508
F-statistic: 3217 on 5 and 5333 DF, p-value: < 2.2e-16Nous notons que le modèle est d’ailleurs bien globalement significatif, grâce au test de Fisher et sa p-value; et comporte un \(R^2\) correct d’environ 75%.
Nous faisons de même pour Police1.
reg_pol1 = lm(Police1 ~ ., data = df_pol)
stepAIC(reg_pol1, trace = 0)
Call:
lm(formula = Police1 ~ RUC + region + agecat + Acompm + nbpers +
Bauto + Police2, data = df_pol)
Coefficients:
(Intercept) RUC region agecat41-50
20.758246 0.001214 0.019178 -1.006152
agecat51-60 agecat61-96 AcompmCouple avec enfant(s) AcompmCouple sans enfant
-0.560081 0.466588 -8.005988 0.082492
AcompmPersonne seule nbpers BautoPas de vehicule Police2
0.122622 3.038834 -1.386994 0.005741 reg_pol1_AIC = lm(formula = Police1 ~ RUC + region + agecat + Acompm + nbpers +
Bauto + Police2, data = df_pol)
summary(reg_pol1_AIC)
Call:
lm(formula = Police1 ~ RUC + region + agecat + Acompm + nbpers +
Bauto + Police2, data = df_pol)
Residuals:
Min 1Q Median 3Q Max
-9.4027 -1.6064 0.0019 1.6312 10.7888
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.076e+01 2.091e-01 99.266 < 2e-16 ***
RUC 1.214e-03 1.019e-05 119.077 < 2e-16 ***
region 1.918e-02 1.212e-02 1.582 0.1138
agecat41-50 -1.006e+00 9.882e-02 -10.182 < 2e-16 ***
agecat51-60 -5.601e-01 1.163e-01 -4.814 1.52e-06 ***
agecat61-96 4.666e-01 1.097e-01 4.252 2.15e-05 ***
AcompmCouple avec enfant(s) -8.006e+00 1.025e-01 -78.125 < 2e-16 ***
AcompmCouple sans enfant 8.249e-02 1.193e-01 0.692 0.4892
AcompmPersonne seule 1.226e-01 1.584e-01 0.774 0.4389
nbpers 3.039e+00 4.499e-02 67.539 < 2e-16 ***
BautoPas de vehicule -1.387e+00 1.289e-01 -10.760 < 2e-16 ***
Police2 5.741e-03 2.952e-03 1.945 0.0518 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.398 on 5327 degrees of freedom
Multiple R-squared: 0.8223, Adjusted R-squared: 0.8219
F-statistic: 2240 on 11 and 5327 DF, p-value: < 2.2e-16Le modèle est bien globalement significatif, et comporte un \(R^2\) d’environ 82%.
Par curiosité, nous refaisons la régression de Sinistre1, cette fois en gardant Police1 dans la dataframe.
stepAIC(lm(Sinistre1 ~ ., data = df1), trace = 0)
Call:
lm(formula = Sinistre1 ~ RUC + agecat + Acompm + nbpers + Bauto +
Police1, data = df1)
Coefficients:
(Intercept) RUC agecat41-50 agecat51-60
0.1806165 -0.0007522 0.7884658 0.3917138
agecat61-96 AcompmCouple avec enfant(s) AcompmCouple sans enfant AcompmPersonne seule
-0.4354412 -1.6830522 0.0208051 -0.0684577
nbpers BautoPas de vehicule Police1
-2.3894167 1.1140196 0.7955038 reg_tot = lm(formula = Sinistre1 ~ RUC + agecat + Acompm + nbpers + Bauto + Police1, data = df1)
summary(reg_tot)$adj.r.squared[1] 0.9478506Nous notons qu’inclure Police1 dans la régression de Sinistre1 améliore nettement le \(R^2\) (qui passe à environ 95%), mais introduit un problème d’endogénéité.
Question 2)
Nous avions déjà nettoyé les points aberrants lors de la Partie I Question 1), nous pouvons cependant refaire les mêmes graphiques diagnostics.
plot(reg_sin1_AIC)
Nous oservons une bimodalité sur le graphique “Residuals vs Fitted” et “Scale-location”. Nous n’observons pas de points aberrants, ceux-ci ayant été éliminés auparavant.
Le graphique “Q-Q Residuals” nous permet d’observer la normalité des résidus, qui est ici globalement respectée. La normalité des résidus est partiellement utile : d’une part, elle est nécessaire pour respecter les hypothèses de l’estimateur des MCO, ainsi qu’à la validité des tests statistiques, d’autre part, le grand nombre d’observations peut potentiellement permettre de transgresser l’hypothèse de normalité des résidus, et donc d’avoir un modèle valide asymptotiquement.
Nous pouvons aussi visualiser la normalité des résidus avec un simple histogramme de ceux-ci.
hist(resid(reg_sin1_AIC), main="Histogramme des résidus", col="skyblue", border="white")
Nous faisons le même diagnostic pour la régression de Police1 sur les autres variables.
plot(reg_pol1_AIC)
Nous observons aussi une normalité des résidus, avec pas de points aberrants.
Question 3)
Le modèle de la Question 1) a déjà été estimé sans les points aberrants.
Question 4)
Grâce aux deux modèles précédents, nous effectuons les prédictions de Sinistre1 et Police1 des 10 individus de la dataframe df0.
Montrer le code
df0_sin <- df0[, !(names(df0) == "Sinistre1")]
df0_sin <- df0_sin[, !(names(df0_sin) == "Police1")]
df0_pol <- df0[, !(names(df0) == "Police1")]
df0_pol <- df0_pol[, !(names(df0_pol) == "Sinistre1")]
stats <- data.frame(
Individu = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10),
Sinistre1_réel = df0$Sinistre1,
Sinistre1_estimé = predict(reg_sin1_AIC, df0_sin),
Police1_réel = df0$Police1,
Police1_estimé = predict(reg_pol1_AIC, df0_pol)
)
kable(stats, format = "markdown")| Individu | Sinistre1_réel | Sinistre1_estimé | Police1_réel | Police1_estimé | |
|---|---|---|---|---|---|
| 5343 | 1 | 22.76060 | 18.591777 | 36.76939 | 33.55132 |
| 5344 | 2 | 16.06724 | 17.607349 | 34.00404 | 36.91910 |
| 5345 | 3 | 21.20237 | 18.265117 | 37.97438 | 35.59925 |
| 5346 | 4 | 21.57460 | 17.689517 | 42.55233 | 37.93391 |
| 5347 | 5 | 21.96043 | 18.615112 | 44.09972 | 40.17874 |
| 5348 | 6 | 15.03810 | 18.265117 | 31.35385 | 35.57511 |
| 5349 | 7 | 18.63927 | 17.567506 | 32.34460 | 32.87165 |
| 5350 | 8 | 17.59627 | 22.083837 | 55.75561 | 59.79546 |
| 5351 | 9 | 10.30899 | 9.671337 | 30.04415 | 31.44112 |
| 5352 | 10 | 20.99097 | 19.313054 | 45.23877 | 42.94950 |
Nous notons que nous aurions pu avoir de meilleurs résultats en gardant Sinistre1 et Police1 en variable explicative du modèle opposé, au prix d’un problème d’endogénéité.
Partie IV
Question 1)
On a :
\[ \delta_i=\mathbb{1}_{\{U_i=a+bZ_i+\epsilon_i>0\}}\text{, avec } \epsilon_i \sim \mathcal{N}(0,1) \]
On remarque que \(\delta_i\) est à valeurs dans \(\{0;1\}\), avec :
\[ \begin{aligned} p_i&=\mathbb{P}(\delta_i=1) \\&=\mathbb{P}(a+bZ_i+\epsilon_i>0) \\&=\mathbb{P}(\epsilon_i>-(a+bZ_i))\text{, avec } \epsilon_i\sim\mathcal{N}(0,1) \\&=1-\Phi(-(a+bZ_i)) \\&=\Phi(a+bZ_i)\text{, par symétrie de la loi normale.} \end{aligned} \]
Ainsi \(\delta_i\sim\mathcal{B}(1,p_i)=\mathcal{B}(1,\Phi(a+bZ_i))\).
Question 2)
Nous effectuons une simulation pour n=1000 et traçons en vert la droite des données exactes, en rouge la droite de régression du modèle complet, et en bleu la droite de régression du modèle positif.
Montrer le code
#Simulation des données
set.seed(4)
n <- 1000
a <- -0.5
b <- 1
alpha <- 2*3
beta <- 2/5
Z <- rlnorm(n, meanlog = 0, sdlog = 1)
X <- rlnorm(n, meanlog = 1, sdlog = 2)
delta <- rbinom(n, size = 1, prob = pnorm(a + b * Z))
Y <- (alpha + X * beta + rnorm(n,mean = 0,sd = 1)) * delta
donnees <- data.frame(Z, X, delta, Y)
#Variables correspondant aux droites à afficher
droite <- alpha + beta * (1:n)
regcomp = lm(Y ~ X, data = donnees)
donnees_pos <- donnees[donnees$Y > 0, ]
regpos = lm(Y ~ X, data=donnees_pos)
#Affichage des données
graphique<-ggplot(donnees, aes(x = X, y = Y)) +
geom_point() +
geom_smooth(data = donnees, method = "lm", se = FALSE, color = "red")+
geom_smooth(data = donnees_pos, aes(x = X, y = Y), method = "lm", se = FALSE, color = "blue")+
geom_line(data = donnees, aes(x = 1:n, y = droite),color = "green")
print(graphique)
Nous remarquons que la droite bleue du modèle positif est une très bonne estimation de la droite verte tandis que la droite rouge du modèle complet n’est pas très bonne.
Le modèle complet ne semble pas convergent car les valeurs de \(a\) et \(b\) estimées sont très eloignées des valeurs voulues (\(a=6\) et \(b=0.4\)) même pour \(n=1000\). En revanche, la régression du modèle positif donne des estimateurs convergents. Ils sont en effet très proches des valeurs attendues. Nous instaurons néanmoins un biais lors de la suppression des valeurs nulles.
Question 3)
On estime le modèle \(\delta_i\) sur les données produites en question 2).
prob = glm(delta ~ Z, family = binomial(link = "probit"))
summary(prob)
Call:
glm(formula = delta ~ Z, family = binomial(link = "probit"))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.66630 0.08732 -7.63 2.34e-14 ***
Z 1.15138 0.09236 12.47 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1249.17 on 999 degrees of freedom
Residual deviance: 937.06 on 998 degrees of freedom
AIC: 941.06
Number of Fisher Scoring iterations: 7Nous devons trouver \(a=-0.5\) et \(b=1\). Nous remarquons que les valeurs associées aux coefficients estimés par la régression en sont très proches et elles sont toutes les deux très significatives quand \(n=1000\). On peut notamment tester avec un n plus grand, par exemple n = 100000.
n = 100000
Z <- rlnorm(n, meanlog = 0, sdlog = 1)
delta <- rbinom(n, size = 1, prob = pnorm(a + b * Z))
prob = glm(delta ~ Z, family = binomial(link = "probit"))
summary(prob)
Call:
glm(formula = delta ~ Z, family = binomial(link = "probit"))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.493433 0.008420 -58.6 <2e-16 ***
Z 0.990417 0.008492 116.6 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 121772 on 99999 degrees of freedom
Residual deviance: 94499 on 99998 degrees of freedom
AIC: 94503
Number of Fisher Scoring iterations: 7Les estimateurs sont maintenant encore plus proches des vraies valeurs \(a=-0.5\) et \(b=1\), nous concluons donc que les estimateurs sont convergents.
Question 4)
Soit \(Z\sim \mathcal{N}(0,1)\)
\[ \begin{aligned} \mathbb{E}[ Z|Z>c ] &= \mathbb{E}[ Z\,\mathbb{1}_{\{ Z>c \} } | Z>c ] \\ &= \frac{\mathbb{E}[ Z\,\mathbb{1}_{\{ Z>c \} } ] }{\mathbb{P}(Z>c)} \\ &= \frac{\int_{c}^{+\infty} z\varphi(z) \, dz }{1-\mathbb{P}(Z<c) } \\ &=\frac{\varphi(c)}{\Phi(-c)} \\ &=\frac{\varphi(-c)}{\Phi(-c)} \text{, par symétrie de la loi normale} \end{aligned} \]
Avec notamment \(\frac{\varphi(-c)}{\Phi(-c)}\) l’inverse du ratio de Mills.