Skip to content

Instantly share code, notes, and snippets.

@mhermans
Last active January 2, 2016 07:38
Show Gist options
  • Select an option

  • Save mhermans/8270906 to your computer and use it in GitHub Desktop.

Select an option

Save mhermans/8270906 to your computer and use it in GitHub Desktop.
Reproducable example of projection using mstls
library(devtools)
dev_mode(on=TRUE)
install_github(username='mhermans', repo='mstls')
library(mstls)
data(YU61female) # included dataset in mstls: 2-region female population (Yugoslavia, 1961)
# calculate observed rates
obs_rates <- with(
YU61female,
cbind(
to_SL/population,
to_RYU/population,
deaths/population,
births/population
)
)
obs_rates <- cbind(
rep(seq(0, 85, 5), 2),
rep(seq(4, 89, 5), 2),
rep(c(1,2), each=18),
obs_rates)
colnames(obs_rates) <- c('x', 'xn', 'state', 'Mx_i1', 'Mx_i2', 'Mdx', 'Fx')
Mx.ex <- transfer_matrix(obs_rates[,-7], multiple=TRUE) # drop Fx
Px.ex <- transfer_prob(Mx.ex, multiple=TRUE)
lx.ex <- expected_survivors(Px.ex, 100000)
Lx.ex <- years_lived(lx.ex, Mx.ex)
Sx.ex <- survivor_prop(years_lived=Lx.ex)
Bx.ex <- birth_prop(birth_rates=obs_rates[,7], transition_probs=Px.ex, survivor_prop=Sx.ex)
G.ex <- projection_matrix(birth_prop=Bx.ex, survivor_prop=Sx.ex)
dim(G.ex)
# add labels to G
a <- rep(seq(0,80,5), each=2)
b <- rep(1:2, 17)
rlbls <- c(c('Bx_i1', 'Bx_i2'), paste('S', a, '_i', b, sep=''))
a <- rep(seq(0,85,5), each=2)
b <- rep(1:2, 18)
clbls <- paste(a, '_', b, 'j', sep='')
dimnames(G.ex) <- list(rlbls, clbls)
rm(a, b, clbls, rlbls)
# project 8 steps
N.projected <- project(YU61female$population, n_states=2, pmat=G.ex, n_steps=8)
rownames(N.projected) <- rep(seq(0, 85, 5), 2)
N.projected
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment