inst/doc/MicSim_vignette.R

In MicSim: Performing Continuous-Time Microsimulation

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# Load library
library(MicSim)

## -----------------------------------------------------------------------------
# Defining simulation horizon
startDate <- 20140101 # yyyymmdd
endDate   <- 20181231 # yyyymmdd
simHorizon <- c(startDate=startDate, endDate=endDate)

# Seed for random number generator
set.seed(12)

# Definition of maximal age (sharp, i.e. max age is 100.0)
maxAge <- 100

# Definition of state space
sex <- c("m","f")                     
country_R <- c("NL","ES","SE")
fert <- c("0","1","2","3","4+","999")           
stateSpace <- expand.grid(sex=sex,country_R=country_R,fert=fert)

# Definition of nonabsorbing and absorbing states
absStates <- c("dead","rest")   

## ----setup--------------------------------------------------------------------
# initial population included in the MicSim package
initPop <- initPopMigrExp
head(initPop)

# population of immigrants entering to virtual population over simulation time
immigrPop <- immigrPopMigrExp # included in the MicSim package
head(immigrPop)

# Definition of initial states for newborns
fixInitStates <- 2 # give indices for attribute/substate that will be taken over from the mother, here: nat (nationality)
varInitStates <- rbind(c("m","ES","0"), c("f","ES","0"), # to have possibility to define distinct sex ratios in distinct countries, 
                       c("m","NL","0"), c("f","NL","0"), # hold substate that are inherited by mother in the init state (i.e. nationality)
                       c("m","SE","0"), c("f","SE","0")) 
initStatesProb <- c(0.5151,0.4849, # probabilities for female / male newborns for each nationality separately
                    0.5124,0.4876,
                    0.5146,0.4854)

## -----------------------------------------------------------------------------
# Load rates from data included in the MicSim package (one column for one transition)
rates <- migrExpRates

# Define function to easily transform rates in function format required by MicSim 
require(glue)
for (i in 1:length(names(rates))) {
  script_var = names(rates[i])
  eval(parse(text = glue("{script_var} <- unlist(as.numeric(rates[,{i}]))")))
}

# --------------------------------------------------------------------------
# Fertility Rates
# --------------------------------------------------------------------------

fertRates_NL_0_1 <- function(age,calTime, duration){
  rate <- fert_NL_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_1_2 <- function(age,calTime, duration){
  rate <- fert_NL_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_2_3 <- function(age,calTime, duration){
  rate <- fert_NL_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_3_4 <- function(age,calTime, duration){
  rate <- fert_NL_3_4[as.integer(age)-+1]
  return(rate)  
}

fertRates_ES_0_1 <- function(age,calTime, duration){
  rate <- fert_ES_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_1_2 <- function(age,calTime, duration){
  rate <- fert_ES_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_2_3 <- function(age,calTime, duration){
  rate <- fert_ES_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_3_4 <- function(age,calTime, duration){
  rate <- fert_ES_3_4[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_0_1 <- function(age,calTime, duration){
  rate <- fert_SE_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_1_2 <- function(age,calTime, duration){
  rate <- fert_SE_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_2_3 <- function(age,calTime, duration){
  rate <- fert_SE_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_3_4 <- function(age,calTime, duration){
  rate <- fert_SE_3_4[as.integer(age)+1]
  return(rate)  
}

# --------------------------------------------------------------------------
# Internal Migration Rates
# --------------------------------------------------------------------------

`%!in%` <- Negate(`%in%`)

for(i in 1:length(country_R)) {
  other_provinces = country_R[which(country_R %!in% country_R[i])]
  for(k in 1:length(other_provinces)) {
    eq = paste(sprintf('%s_%s_rates', glue("{country_R[i]}"),  glue("{other_provinces[k]}")), 
               '<- function(age,calTime, duration)', '{',
               sprintf('rate <- rate_%s_%s[as.integer(age)+1]', glue("{country_R[i]}"), 
               glue("{other_provinces[k]}")), "\n ", 'return(rate)','}')
    eval(parse(text = eq))
  }
}

# --------------------------------------------------------------------------
# Mortality Rates
# --------------------------------------------------------------------------

# Female mortality
for(i in 1:length(country_R)) {
  eq = paste(sprintf('mortRates_f_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- mort_f_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# Male mortality
for(i in 1:length(country_R)) {
  eq = paste(sprintf('mortRates_m_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- mort_m_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# ---------------------------------------------------------------------------
# Emigration rates
# ---------------------------------------------------------------------------

# Emigration rates for females
for(i in 1:length(country_R)) {
  eq = paste(sprintf('emigrRates_f_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- emig_f_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# Emigration rates for males
for(i in 1:length(country_R)) {
  eq = paste(sprintf('emigrRates_m_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- emig_m_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}


## -----------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# Transition matrix for fertility
# ---------------------------------------------------------------------------
fertTrMatrix <- cbind(c("f/ES/0->f/ES/1","f/ES/1->f/ES/2","f/ES/2->f/ES/3","f/ES/3->f/ES/4+",
                        "f/SE/0->f/SE/1","f/SE/1->f/SE/2","f/SE/2->f/SE/3","f/SE/3->f/SE/4+",
                        "f/NL/0->f/NL/1","f/NL/1->f/NL/2","f/NL/2->f/NL/3","f/NL/3->f/NL/4+"),
                      c("fertRates_ES_0_1", "fertRates_ES_1_2", "fertRates_ES_2_3", "fertRates_ES_3_4",
                        "fertRates_SE_0_1", "fertRates_SE_1_2", "fertRates_SE_2_3", "fertRates_SE_3_4",
                        "fertRates_NL_0_1", "fertRates_NL_1_2", "fertRates_NL_2_3", "fertRates_NL_3_4"))

# ---------------------------------------------------------------------------
# Transition matrix for migration
# ---------------------------------------------------------------------------
# First: make names for transition matrix, i.e. stateOfOrigin->stateOfDestination 
testo1 <- "intmigTrMatrix <- cbind(c("
for(i in 1:length(country_R)) {
  for(m in 1:length(country_R)) {
    if(m != i) {
      eq1 = paste(sprintf("'%s->%s'", glue("{country_R[i]}"), glue("{country_R[m]}")))
      if (i == length(country_R) & m == (i-1)) {
        testo1 = paste(testo1,eq1)
      } else {
        testo1 = paste(testo1 ,eq1, ",")
      }
    }
    if (m == i & m == length(country_R)){
      testo1 = paste(testo1, "),")
      
    }
  }
}

#Second: set names for transition functions
testo1 = paste (testo1,"c(")
for(i in 1:length(country_R)) {
  for(m in 1:length(country_R)) {
    if(m != i) {
      eq1 = paste(sprintf("'%s_%s_rates'", glue("{country_R[i]}"), glue("{country_R[m]}")))
      if (i == length(country_R) & m == (i-1)) {
        testo1 = paste(testo1,eq1)
      } else {
        testo1 = paste(testo1 ,eq1, ",")
      }
    }
    if (m == i & m == length(country_R)){
      testo1 = paste(testo1, "))")
    }
  }
}
eval(parse(text = testo1))

# ---------------------------------------------------------------------------
# Transition matrix for mortality and emigration
# ---------------------------------------------------------------------------

testo <- "absTransitions <- rbind("
for(i in 1:length(country_R)) {
  for(m in 1:length(sex)) {
    eq1 = paste(sprintf("c('%s/%s/dead', 'mortRates_%s_%s')", glue("{sex[m]}"), glue("{country_R[i]}") ,
                        glue("{sex[m]}"), glue("{country_R[i]}")),',',
                sprintf("c('%s/%s/rest', 'emigrRates_%s_%s')", glue("{sex[m]}"),glue("{country_R[i]}") ,
                        glue("{sex[m]}"), glue("{country_R[i]}")))
    if(i == length(country_R) & m == length(sex)) {
      testo = paste(testo, eq1, ")")
    } else {
      testo = paste(testo,eq1, ",")
    }
  }
}
eval(parse(text = testo))

# ---------------------------------------------------------------------------
# Combine all
# ---------------------------------------------------------------------------

allTransitions <- rbind(fertTrMatrix, intmigTrMatrix)
transitionMatrix <- buildTransitionMatrix(allTransitions=allTransitions,
                                          absTransitions=absTransitions, 
                                          stateSpace=stateSpace)

# ---------------------------------------------------------------------------
# Define transitions triggering a birth event
# ---------------------------------------------------------------------------

fertTr <- fertTrMatrix[,1]

## -----------------------------------------------------------------------------
pop <- micSim(initPop=initPop[1:500,], immigrPop=immigrPop[1:100,],
              transitionMatrix=transitionMatrix, absStates=absStates,
              varInitStates=varInitStates, initStatesProb=initStatesProb,
              fixInitStates=fixInitStates,
              maxAge=maxAge, simHorizon=simHorizon,fertTr=fertTr)
head(pop)


## -----------------------------------------------------------------------------
popLong <- convertToLongFormat(pop, migr=TRUE)
head(popLong)


## -----------------------------------------------------------------------------
popWide <- convertToWideFormat(pop)
head(popWide)

## -----------------------------------------------------------------------------

cores <- 3
seeds <- c(34,145,97)
pop <- micSimParallel(initPop=initPop, immigrPop=immigrPop,
                      transitionMatrix=transitionMatrix, absStates=absStates,
                      varInitStates=varInitStates, initStatesProb=initStatesProb,
                      fixInitStates=fixInitStates,
                      maxAge=maxAge, simHorizon=simHorizon,fertTr=fertTr,
                      cores=cores, seeds=seeds)
head(pop)