Nous allons utiliser l’estimateur de Horvitz-Thompson du total pour calculer le total d’une variable dans une population à partir d’un échantillonnage. Nous allons pour cela commencer par creer une population de 2500 unités statistiques dans une grille (donc de 50x50 cases) et créer une variable z - qui sera la variable que nous allons chercher a estimer par échantillonnage. Comme cela est le plus répandu dans la nature, la variable z aura une dispersion suivant une loi normale N(mu,x).
#Grille de 50 x 50 cases
side=50
my.mat <- matrix(NA, nrow=side, ncol=side)
x.coord <- rep(1:side, each=side)
y.coord <- rep(1:side, times=side)
xy <- data.frame(x.coord, y.coord)
#On crée une variable réponse N(100,5), pour pas avoir de valeurs négatives
z<-rnorm(length(xy[,1]),100,5)
population<-cbind(x.coord,y.coord,z)
# plot de la population
library(raster)
my.mat[] <- z
r <- raster(my.mat)
plot(r, axes=F)
hist(r)On summary z
summary(z)
sum(z)Nous allons tester plusieurs protocoles d’échantillonnages pour le calcul de l’estimateur de Horvitz-Thompson: l’aléatoire simple (SRS - simple ransom sampling), l’échantillonnage systématique (SS) et le Generalized tesselation sampling (GRTS). Nous allons utiliser le package SDraw; qui a besoin d’un objet SpatialPointsDataFrame pour fonctionner. Pour s’affranchir de l’erreur d’échantillonnage due au caractére aléetoire de chacuns de ces protocoles, nous allons réaliser 500 simulations d’un tirage d’unités statistiques.
détermination des variables R
library(sampling)
library(SDraw)
pts=250 # nombre d'échantillons que l'on veut
nb_us<-2500
pop<-SpatialPointsDataFrame(population[,1:2], as.data.frame(z))
nsimul=500 # nombre de simulationsres_srs <-rep(list(),nsimul) # on crée une liste vide où on va mettre les résultats des 1000 simulations
for (i in 1:nsimul){
SRSpts <- SDraw::srs.point(pop, pts)
res_srs[i] <-list(cbind(SRSpts@coords,SRSpts$z))
i<-i+1
}
head(res_srs)
plot(res_srs[[2]])res_sss <-rep(list(),nsimul) # on crée une liste vide où on va mettre les résultats des 1000 simulations
for (i in 1:nsimul){
SSSpts <- SDraw::sss.point(pop, pts)
res_sss[i] <-list(cbind(SSSpts@coords,SSSpts$z))
i<-i+1
}
head(res_sss)
plot(res_sss[[2]])res_grts <-rep(list(),nsimul) # on crée une liste vide où on va mettre les résultats des 1000 simulations
for (i in 1:nsimul){
GRTSpts <- SDraw::grts.point(pop, pts)
res_grts[i] <-list(cbind(GRTSpts@coords,GRTSpts$z))
i<-i+1
}
head(res_grts)
plot(res_grts[[2]])Ici, nous allons utiliser le package sampling pour calculer l’estimateur du total de Horvitz-Thompson ainsi que sa variance. Cet estimateur a besoin de connaitre les probabilités d’inclusion de premier ordre de chaque unité d’échantillonnage selectionnées parmi les unités statistiques de la population de base. Comme nous n’utilisons que des protocoles d’échantillonnages dits “simples” et que les tirages des u.s sont avec remise; toutes les unitées statistiquse de la population ont la même probabilité d’être tirées = ont la même probabilitée d’inclusion. Dans ce cas là, le calcul de la probabilité d’inclusion d’une unisté statistiques est facile; c’est le nombre d’unités échantillonnés (n) divisé par le nombre total d’unités statistiques dans la population (N).
inclprob<- rep((pts/nb_us), pts)HT_srs<- rep(0,nsimul) # on crée un vecteur vide qui prendra les valeurs de HT pour chaque simulation
for (i in 1:nsimul){
HT_srs[i] <- HTestimator(res_srs[[i]][,3],inclprob)
i<-i+1
}
boxplot(HT_srs)HT_sss<- rep(0,nsimul) # on crée un vecteur vide qui prendra les valeurs de HT pour chaque simulation
for (i in 1:nsimul){
HT_sss[i] <- HTestimator(res_sss[[i]][,3],inclprob)
i<-i+1
}
boxplot(HT_sss)HT_grts<- rep(0,nsimul) # on crée un vecteur vide qui prendra les valeurs de HT pour chaque simulation
for (i in 1:nsimul){
HT_grts[i] <- HTestimator(res_grts[[i]][,3],inclprob)
i<-i+1
}
boxplot(HT_grts)HT_tot<-cbind(HT_srs,HT_sss,HT_grts)
boxplot(HT_tot)
abline(h = sum(z), col = "red")Refaites cet exercice avec des tailles d’échantillons différentes. Que remarquez-vous? Quel protocole d’échantillonnage est’il préférable d’utiliser?
La population précedente a une distribution assez homogène sur le site d’étude. Dans cette partie, vous allez créer des populations plus ou moins hétérogènes et mettre en évidence lequel des trois protocoles testés est le plus optimal pour ce genre de populations.
# créer z, une variable qui representera la distribution de l'espèce dans la zone d'étude (en présence absence)
a<-seq(from= 50, to=1, length.out = 2500)
pik=inclusionprobabilities(a, 2000)
z<-rep_len(UPpoisson(pik),2500)
# application d'un variable réponse de distribution normale sur z
rnorm<-rnorm(length(z),100,10)
for (i in 1:length(z)) {
if (z[i]==1) z[i]<- rnorm[i]
else i<-i+1
}
population<-cbind(x.coord,y.coord,z)autre exemple de population heterogène
a1<-seq(from= 1, to=25, length.out = 1275)
a2<-seq(from= 25, to=1, length.out = 1275)
a<-c(a1,a2,a1)