library(ggplot2)
library(lattice)
library(gridExtra)

Campionamento stratificato

Ipotesi a monte: esiste una correlazione tra la variabile oggetto di studio (pil pro capite) e la variabile ausiliaria (continente)

#caricare dataset dati (pilprocapite mondo)
N=142 
y=dati
mu=with(y,mean(gdpPercap))
mu=round(mu,digits=2)
mu
## [1] 11680.07
sigma=with(y,sd(gdpPercap))
sigma=round(sigma,digits=2)
sigma
## [1] 12859.94
round(with(y, tapply(gdpPercap,continent,mean)),digits=0)
##   Africa Americas     Asia   Europe  Oceania 
##     3089    11003    12473    25054    29810
bwplot(~y$gdpPercap|y$continent,xlab="Pil pro capite",main="Boxplot")

Creazione di 3 strati

I=which(y$continent=="Asia" & y$gdpPercap>30000)
y$country[I]
## [1] Hong Kong, China Japan            Kuwait           Singapore       
## 142 Levels: Afghanistan Albania Algeria Angola Argentina ... Zimbabwe
y$gdpPercap[I]
## [1] 39724.98 31656.07 47306.99 47143.18
#attenzione! il primo paese corrisponde Honk Kong in Cina, poi ci sono il Giappone, Kuwait e Singapore
# Creiamo la variabile che identifica i 3 strati
y$strato <- NA
# strato=1 se il paese si trova in Europa o Oceania
y$strato[y$continent=="Europe"] <- 1
y$strato[y$continent=="Oceania"] <- 1
# strato=2 se il paese si trova in Asia o nelle Americhe
y$strato[y$continent=="Asia"] <- 2
y$strato[y$continent=="Americas"] <- 2
# strato=33 se il paese si trova in Africa
y$strato[y$continent=="Africa"] <- 3
# collochiamo USA e Canada nello strato 1 
y$strato[y$country=="United States"] <- 1
y$strato[y$country=="Canada"] <- 1
# collochiamo nello strato 1 i 4 paesi dell'Asia precedentemente identificati
y$strato[I] <- 1

Pilprocapite negli strati

# media negli strati
round(with(y, tapply(gdpPercap,strato,mean)),digits=0)
##     1     2     3 
## 27799  8492  3089
#mediana negli strati
round(with(y, tapply(gdpPercap,strato,median)),digits=0)
##     1     2     3 
## 31063  5877  1452
with(y, tapply(gdpPercap,strato,median))
##         1         2         3 
## 31063.042  5876.864  1452.267
bwplot(~y$gdpPercap|as.factor(y$strato),layout=c(1,3),xlab="Pil pro capite",ylab="strati",main="Boxplot")

Dimensione e pesi degli strati

table(y$strato)
## 
##  1  2  3 
## 38 52 52
N1=38;N2=52;N3=52;
N=N1+N2+N3
N
## [1] 142
#Pesi strati
W1=round(N1/N,digits=3);
W1
## [1] 0.268
W2=round(N2/N,digits=3);
W2
## [1] 0.366
W3=round(N3/N,digits=3);
W3
## [1] 0.366

Immaginiamo di dover estrarre un campione stratificato di dimensione n=20.

Calcolo dimensione campionaria e frazione di campionamento negli strati

n=20
n1<-round(W1*n)
n1
## [1] 5
n2<-round(W2*n)
n2
## [1] 7
n3<-round(W3*n)
n3
## [1] 7
#la dimensione campionaria totale deve essere pari a n=20
n1+n2+n3
## [1] 19
# siccome il totale fa 19, portiamo n1 a 6 in modo che il totale sia pari a n=20
n1=6
n1+n2+n3
## [1] 20
f1=round(n1/N1,digits=2);f1
## [1] 0.16
f2=round(n2/N2,digits=2);f2
## [1] 0.13
f3=round(n3/N3,digits=2);f3
## [1] 0.13

La teoria ci dice a quanto ammonta l’errore standard della media campionaria nel campionamento stratificato

#ci calcoliamo le varianze dei singoli strati
sig2st1=var(y$gdpPercap[y$strato==1])
sig2st2=var(y$gdpPercap[y$strato==2])
sig2st3=var(y$gdpPercap[y$strato==3])
#errore standard media campionaria nel campionamento stratificato
ES_Tst<-(W1^2*(1-f1)*sig2st1/n1+W2^2*(1-f2)*sig2st2/n2+W3^2*(1-f3)*sig2st3/n3)^0.5
ES_Tst=round(ES_Tst,digits=2)
ES_Tst
## [1] 1631.39

Vi ricordate quanto era l’errore standard nel campionamento casuale semplice per n=20?

Con il campionamento stratificato le nostre stime sono tali da garantire una maggiore precisione?

Estrazione di un campione casuale stratificato di n=20 elementi

Costruzione delle 3 urne - 3 strati - dove andremo ad effettuare il CCS

strato1<-y$gdpPercap[y$strato==1]
strato2<-y$gdpPercap[y$strato==2]
strato3<-y$gdpPercap[y$strato==3]

Estrazione di un campione casuale in ogni strato

set.seed(1)
ccs_s1<-sample(strato1,n1)
ccs_s2<-sample(strato2,n2)
ccs_s3<-sample(strato3,n3)

Calcolo media in ogni strato

xbars1<-mean(ccs_s1)
xbars2<-mean(ccs_s2)
xbars3<-mean(ccs_s3)

Calcolo varianza in ogni strato

s1<-var(ccs_s1)
s2<-var(ccs_s2)
s3<-var(ccs_s3)

Calcolo media nel campionamento stratificato

media_st<-(W1*xbars1)+(W2*xbars2)+(W3*xbars3)
media_st
## [1] 11645.34
# indichiamo l'errore di stima con ec 
ec=media_st-mu
round(ec,digits=2)
## [1] -34.73

Calcolo errore standard della stima ottenuta

es_st<-(W1^2*(1-f1)*s1/n1+W2^2*(1-f2)*s2/n2+W3^2*(1-f2)*s3/n3)^0.5
es_st
## [1] 1560.266

Simuliamo la distribuzione campionaria ripetendo l’estrazione del campione stratificato (n=20) per 500 volte

M=500
set.seed(2)
s1=replicate(M,sample(strato1,n1))
s2=replicate(M,sample(strato2,n2))
s3=replicate(M,sample(strato3,n3))
xbarstrato1<-apply(s1,2,FUN=mean)
xbarstrato2<-apply(s2,2,FUN=mean)
xbarstrato3<-apply(s3,2,FUN=mean)
xbarst<-(W1*xbarstrato1)+(W2*xbarstrato2)+(W3*xbarstrato3)

Valore atteso dello stimatore

round(mean(xbarst),digits=2)
## [1] 11666.4

parametro

mu
## [1] 11680.07

Si nota come le stime in media sono pari al parametro incognito mu

Quanto sono precise le nostre stime?

Mediamente, l’errore di stima, corrisponde a quanto postulato dalla teoria ovvero a ES_Tst?

round(sd(xbarst),digits=2)
## [1] 1717.71

Ovviamente si!

Rappresentiamo le prime nostre 50 stime e vediamo come “ballano” intorno al valore del parametro mu

Alcune stime sono vicine al parametro mu altre sono lontane. Si nota come le stime oscillano casualmente attorno alla media della popolazione, ovvero l’errore di stima in media vale zero: stimatore corretto

Distribuzione campionaria dello stimatore media nel campionamento stratificato

Confronto fra le distribuzioni campionarie per n=20 del campionamento casuale semplice e del campionamento stratificato

Si nota chiaramente che, attraverso la stratificazione, le stime sono meno variabili e quindi maggiormente precise