TP 2 : Estimation (suite)
Exercice 1 : Sondage sur les communes
Un statisticien a réalisé un sondage auprès d’un échantillon de \(n\) communes : le tirage a été réalisé à l’aide d’un tirage poissonien. Les probabilités d’inclusion d’ordre 1 ont été définies à l’aide de la taille des ménages : 0.9 si la commune dispose de plus de 100 000 ménages et 0.1 sinon. Les résultats sont décrits dans la table donnees_ech.csv.
Cette table est constituée des variables suivantes :
| Nom de la variable | Description de la variable |
|---|---|
code |
Code géographique de la commune |
lib |
Nom de la commune |
nb_ul |
Nombre d’unités légales |
nb_appart |
Nombre d’appartements dans la commune |
nb_log |
Nombre de logements dans la commune |
nb_men |
Nombre de ménages dans la commune |
nb_fam_mono |
Nombre de familles monoparentales dans la commune |
nb_fam |
Nombre de familles |
prob |
Probabilité d’inclusion d’ordre 1 |
- Donnez la population \(\mathcal{U}\). > Ensemble des communes de France (hexagonale + DROMs)
- Donnez la taille de l’échantillon.
library("sample")
library("dplyr")
donnees_ech <- read.csv("donnees_ech.csv")
cat("La taille de l'échantillon est ", nrow(donnees_ech))
- Proposez un estimateur sans biais du total de n’importe quelle variable d’intérêt. Justifiez.
Les probabilités d’inclusions d’ordre 1 valent 0.1 ou 0.9 donc sont non nulles. L’estimateur du total de Narain-Horvitz-Thompson est donc sans biais.
- Proposez une estimation du :
- nombre de logements en France,
- part d’appartement par département,
- part de familles monoparentales,
- nombre de communes avec au moins 100 unités légales,
- nombre de communes en France,
- nombre de communes dans le Nord.
#Estimation du nombre de logements en France
donnees_ech |>
summarise(sum(nb_log/prob))
#Estimation de la part d'appartements par département
donnees_ech |>
mutate("dep" = if_else(stringr::str_sub(code,1,2) == "97", stringr::str_sub(code,1,2), stringr::str_sub(code,1,3))) |>
group_by(dep) |>
summarise(sum(nb_log/prob))
#Estimation de la part de familles monoparentales
donnees_ech |>
summarise("nb_mono" = sum(nb_fam_mono/prob), "nb_famille" = sum(nb_fam/prob)) |>
summarise(nb_mono/nb_famille)
#Estimation du nombre de communes avec au moins 100 unités légales
donnees_ech |>
mutate(bcp_ul = as.integer(nb_ul > 100)) |>
summarise(sum(bcp_ul/prob))
#Estimation du nombre de communes en France
sum(1/donnees_ech$prob)
#Estimation du nombre de communes dans Nord
donnees_ech |>
filter(stringr::str_starts(code, "59")) |>
summarise(sum(1/prob))
Exercice 2 : Visualisation du biais de l’estimateur d’Horvitz-Thompson
Dans cet exercice, on considère deux populations de \(N = 10 000\) individus. Dans chaque population, une probabilité d’inclusion d’ordre 1 est associée à chaque individu : des échantillons vont être tirés à l’aide d’un plan poissonien. Le but est d’estimer le total d’une variable d’intérêt \(y\). Les bases de sondage sont décrites dans les tables pop1.csv et pop2.csv. (Remarque : le but de cet exercice est d’illustrer l’existence d’un biais pour l’estimateur d’Horvitz-Thompson - on suppose ici disposer de la variable d’intérêt sur toute la population, ce qui n’est pas le cas en pratique).
- Est-ce possible de proposer, dans les deux cas, un estimateur sans biais du total de la variable \(y\) ?
- On souhaite estimer le biais (éventuel) de l’estimateur proposé à la question précédente par simulation : l’idée étant de faire comme si on pouvait tirer un très grand nombre d’échantillons. Proposez une fonction
simulationqui prend en argument une data.framepopulation(composée de deux variablesyetpi) et un nombre de simulationsN_simet qui renvoie un vecteur de tailleN_simqui contient l’estimation du total obtenu en tirant un échantillon selon le plan poissonien associé aux probabilités d’inclusionpi. La fonction suivante permet de tirer selon un plan poissonien :
tirage_poisson <- function(pi){
if(any(pi > 1 | pi < 0)){
stop("pi doit appartenir à [0,1]")
}
sapply(X = pi, FUN = function(p){rbinom(1,1,p)})
}
simulation <- function(population, N_sim){
res <- c()
for(i in 1:N_sim){
ech <- tirage_poisson(population$pi)
res[i] <- with(population[ech == 1, ], sum(y/pi))
}
return(res)
}
- Réalisez \(1000\) tirages à l’aide de la fonction
simulationpour les deux populations (et plans de sondage). Comparez aux totaux.
pop1 <- read.csv("pop1.csv")
pop2 <- read.csv("pop2.csv")
res1 <- simulation(pop1, 1000)
res2 <- simulation(pop2, 1000)
hist(res1, breaks = 30)
sum(pop1$y)
mean(res1)
hist(res2, breaks = 30)
sum(pop2$y)
mean(res2)