knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The sample data on migration between Spain, Sweden and the Netherlands were prepared by Claudio Bosco (European Commission), Daniela Ghio (European Commission), Maurizio Teobaldelli (European Commission) and Sabine Zinn (German Socio-Economic Panel, Humboldt University Berlin). EUROSTAT data were used as the data source. The sample was intentionally kept small. MicSim can also handle larger numbers of cases, but to be efficient in terms of run times, this requires a bit more computing power with more CPUs.
# 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")
# 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)
Beware: Rates have to be given at least for age [0,maxAge) and for all years within the simulation horizon. At this the exact value of maxAge is excluded, i.e. here 100.00 but not e.g. age=99.9999.
# 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]
For illustration purpose, the subsequent run is limited to the first 500 people and to 100 migrants. However, just remove the restriction to run the whole sample, i.e. using initPop instead of initPop[1:500,] and immigrPop instead of immigrPop[1:100,].
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)
Try this to speed up your simulation run. The more cores you can use the faster the simulation will be executed. This example uses three cores. Give as many seeds as you use cores. That way you can replicate your results.
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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.