inst/doc/models.R

In SamplingStrata: Optimal Stratification of Sampling Frames for Multipurpose Sampling Surveys

## ----setup, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 999)
knitr::opts_chunk$set(echo = TRUE, fig.width=6, fig.height=4)


## ---- eval = T, message=F------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(SamplingStrata)
data("swissmunicipalities")
swissmunicipalities$id <- c(1:nrow(swissmunicipalities))
swissmunicipalities$dom <- 1
head(swissmunicipalities[,c(3,4,5,6,9,22)])

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
set.seed(1234)
swiss_sample <- swissmunicipalities[sample(c(1:nrow(swissmunicipalities)),500),]

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
mod_Airbat_POPTOT <- lm(swiss_sample$Airbat ~ swiss_sample$POPTOT)
summary(mod_Airbat_POPTOT)

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
mod_Surfacesbois_HApoly <- lm(swiss_sample$Surfacesbois ~ swiss_sample$HApoly)
summary(mod_Surfacesbois_HApoly)

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Airbat <- computeGamma(mod_Airbat_POPTOT$residuals,
             swiss_sample$POPTOT,
             nbins = 10)
Airbat

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Surfacesbois <- computeGamma(mod_Surfacesbois_HApoly$residuals,
             swiss_sample$HApoly,
             nbins = 10)
Surfacesbois

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
model <- NULL
model$beta[1] <- mod_Airbat_POPTOT$coefficients[2]
model$sig2[1] <- Airbat[2]^2
model$type[1] <- "linear"
model$gamma[1] <- Airbat[1]
model$beta[2] <- mod_Surfacesbois_HApoly$coefficients[2]
model$sig2[2] <- Surfacesbois[2]^2
model$type[2] <- "linear"
model$gamma[2] <- Surfacesbois[1] 
model <- as.data.frame(model)
model

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
frame <- buildFrameDF(swissmunicipalities,
                      id = "COM",
                      domainvalue = "dom",
                      X = c("POPTOT", "HApoly"),
                      Y = c("POPTOT", "HApoly"))
frame$Airbat <- swissmunicipalities$Airbat
frame$Surfacesbois <- swissmunicipalities$Surfacesbois

## ---- eval = T-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
cv <- as.data.frame(list(DOM=rep("DOM1",1),
                         CV1=rep(0.05,1),
                         CV2=rep(0.05,1),
                         domainvalue=c(1:1)
                    ))
cv

## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1234)
#  solution <- optimStrata(
#      method = "continuous",
#      errors = cv ,
#      framesamp = frame,
#      model = model,
#      iter = 25,
#      pops = 20,
#      parallel = FALSE,
#      nStrata = 5)
#  
#  # *** Domain :  1   1
#  #  Number of strata :  2896
#  # GA Settings
#  #   Population size       = 20
#  #   Number of Generations = 25
#  #   Elitism               = 4
#  #   Mutation Chance       = 0.111111111111111
#  #
#  #
#  #
#  #  *** Sample cost:  295.1979
#  #  *** Number of strata:  5
#  #  *** Sample size :  296
#  #  *** Number of strata :  5
#  # ---------------------------
#  

## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  outstrata <- solution$aggr_strata
#  framenew <- solution$framenew
#  framenew$Y3 <- framenew$AIRBAT
#  framenew$Y4 <- framenew$SURFACESBOIS
#  results <- evalSolution(framenew, outstrata, 500, progress = FALSE)
#  results$coeff_var
#  #      CV1    CV2    CV3    CV4  dom
#  # 1 0.0369 0.0241 0.0352 0.0366 DOM1

## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  set.seed(1234)
#  solution <- optimStrata(
#      method = "continuous",
#      errors = cv ,
#      framesamp = frame,
#      model = NULL,
#      iter = 25,
#      pops = 20,
#      parallel = FALSE,
#      nStrata = 5)
#  
#  #  *** Domain :  1   1
#  #  Number of strata :  2896
#  # GA Settings
#  #   Population size       = 20
#  #   Number of Generations = 25
#  #   Elitism               = 4
#  #   Mutation Chance       = 0.111111111111111
#  #
#  #
#  #
#  #  *** Sample cost:  151.9848
#  #  *** Number of strata:  5
#  #  *** Sample size :  152
#  #  *** Number of strata :  5

## ---- eval = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  outstrata <- solution$aggr_strata
#  framenew <- solution$framenew
#  framenew$Y3 <- framenew$AIRBAT
#  framenew$Y4 <- framenew$SURFACESBOIS
#  results <- evalSolution(framenew, outstrata, 500, progress = FALSE)
#  results$coeff_var
#  #      CV1    CV2    CV3    CV4  dom
#  # 1 0.0501 0.0491 0.0502 0.0686 DOM1