#### STIME DEI DISEGNI CAMPIONARI ##

load("campioni.rdata")
library(sampling)
data(swissmunicipalities)
library(survey)
n=200
N=nrow(swissmunicipalities)

### ESTIMATES OF THE SRSWOR ##

pik<-rep(n/N,n)

d_SRS<-svydesign(ids=~1,probs=~pik,data=s_SRS) #disegno con reinserimento (SRSWR)
svytotal(~H00P01+H00P02+H00P03+H00P04,d_SRS)   #tratta il campione d_SRS come se fosse con reinserimento 
                                               #ids: serve per il campionamento a + stadi (sia probs che fpc definite per ciascuno stadio)
d_SRS<-svydesign(ids=~1,data=s_SRS,fpc=~pik)   #disegno senza reinserimento (SRSWOR)
                                               #fpc=finite population correction!
totest_SRS<-svytotal(~H00P01+H00P02+H00P03+H00P04,d_SRS)
totest_SRS                                     #stime uguali al precedente ma SE più bassi (dovuto a fpc!)
        total    SE            
H00P01 725318 91997            #N.B. i valori ottenuti sono delle stime campionarie!
H00P02 734976 80081            #Stime diverse a seconda del campione estratto!
H00P03 313492 33467
H00P04 496288 44494

cv(totest_SRS)                                 #in ?surveysummary
   H00P01    H00P02    H00P03    H00P04        #coefficienti di variazione
0.1268369 0.1089577 0.1067546 0.0896546 

### ESTIMATES OF THE PPS ##

pik<-inclusionprobabilities(swissmunicipalities$Pop020,n) #pik per l'intera pop

set.seed(121212)          #per la riproducibilità del campione 0,1 PPS
PPS<-UPtille(pik)
set.seed(121212) 
PPS2<-UPtille(pik)
PPS==PPS2                 #TRUE; ho riprodotto lo stesso campione
              
s_PPS<-swissmunicipalities[PPS == 1,]
pik<-pik[PPS==1]          #pik a livello campionario

#d_PPS<-svydesign(ids=~1,data=s_PPS,fpc=~pik,pps=ppsmat(??)) #ppsmat for HT estimator
d_PPS<-svydesign(ids=~1,data=s_PPS,fpc=~pik,pps="brewer")    #Brewer's approximation per prob. 2° ordine
totest_PPS<-svytotal(~H00P01+H00P02+H00P03+H00P04,d_PPS,deff=TRUE) #inserire variabile qualitativa
totest_PPS
           total        SE   DEff   #DEff=V(PPS)/V(SRS)
H00P01 1069141.6   26203.8 0.0038   #maggiore precisione per il PPS:
H00P02  947002.4   15987.8 0.0041   #SE più bassi di quelli del SRS!
H00P03  401308.3    5685.9 0.0044   
H00P04  604850.0    3477.2 0.0012

cv(totest_PPS)
     H00P01      H00P02      H00P03      H00P04 
0.024509189 0.016882499 0.014168300 0.005748904 

### ESTIMATES OF THE STRATIFIED SAMPLING ##

d_STRAT<-svydesign(ids=~1,strata=~REG,data=s_STRAT,fpc=~Prob)  #disegno stratificato+srswor
totest_STRAT<-svytotal(~H00P01+H00P02+H00P03+H00P04+P00BMTOT+P00BWTOT,d_STRAT,deff=TRUE)
totest_STRAT
           total      SE   DEff   #DEff=V(STRAT)/V(SRS)
H00P01   1106718  278671 1.2486   #in questo caso maggiore precisione 
H00P02   1026953  221095 1.5861   #per il SRS!
H00P03    422072   75610 1.3498
H00P04    651483  100711 1.3643
P00BMTOT 3753114  700680 1.4025
P00BWTOT 3872528  751843 1.4039

cv(totest_STRAT)            
   H00P01    H00P02    H00P03    H00P04 
0.2517996 0.2152920 0.1791400 0.1545872 

#deff=TRUE design effect (confronto disegno con SRSWOR!!)
#under stratified random sampling the design effect in a subset 
#consisting of a single stratum will be 1.0.

tot<-matrix(colSums(swissmunicipalities[,c(18,19,20,21,11,12)]),6,1)  
rownames(tot)<-c("H00P01","H00P02","H00P03","H00P04","P00BMTOT","P00BWTOT")
colnames(tot)<-"totpop"                          
tot                              #totali veri da confrontare con le stime!  
          totpop
H00P01   1120878
H00P02    985971
H00P03    403032
H00P04    605518
P00BMTOT 3567567
P00BWTOT 3720443

#tabelle a doppia entrata (interazione tra 2 o più variabili qualitative)

women<-rep(0,nrow(s_SRS))
women[s_SRS$P00BWTOT>200]<-1    #women è diventata una variabile dicotomica a partire da P00BWTOT
pop020<-rep(0,nrow(s_SRS))
pop020[s_SRS$Pop020>200]<-1     #pop020 è diventata una variabile dicotomica a partire da Pop020
table(women)
table(pop020)

d_SRS2<-update(d_SRS,~women+pop020)            #aggiorna il disegno
svytotal(~interaction(women,pop020),d_SRS2)    #stima dei totali su tabella di contingenza

#stima del totale per domini di studio (es. le regioni)

svyby(~Pop020,by=~REG,d_PPS,FUN=svytotal)
svyby(~P00BMTOT+P00BWTOT,by=~REG,d_STRAT,FUN=svytotal)

### STIME BOOTSTRAP, JACKNIFE ##

#nell'inferenza da popolazioni finite il metodo del ricampionamento 
#introduce il concetto di "pesi replicati"

# nuovi disegni campionari basati sul SRS

d_SRSboot<-as.svrepdesign(d_SRS,type="bootstrap",rep=200)   #rep=n
d_SRSjkn<-as.svrepdesign(d_SRS,type="JK1")                  #n=#replicates

totest_SRSboot<-svytotal(~H00P01+H00P02+H00P03+H00P04,d_SRSboot)
totest_SRSjkn<-svytotal(~H00P01+H00P02+H00P03+H00P04,d_SRSjkn)
totest_SRSboot
totest_SRSjkn