DemoDecomp (version 1.0.1)

horiuchi: Numeric Approximation of Continuous Decomposition

Description

This is an exact R implementation of the decomposition code in Matlab offered by the authors in the supplementary material given here: <http://www.demog.berkeley.edu/~jrw/Papers/decomp.suppl.pdf>. The difference between DecompContinuous() and this function is that DecompContinuousOrig takes rates1 and rates2 as single vectors, rather than as matrices, and output is also returned as a vector. This difference makes the function more flexible, but may add a step when writing the function to be decomposed. See examples.

Usage

horiuchi(func, pars1, pars2, N, ...)

Value

returns effectmat, a matrix of the variable effects that is organized in the same way as pars1 and pars2.

Arguments

func

A function specified by the user. This must be able to take the vectors rates1 or rates2 as its argument, and to return the value of the function, y, when evaluated for these rates. It may also have additional arguments, not to be decomposed.

pars1

vector of covariates to be passed on as arguments to func(). Covariates can be in any order, as long as func() knows what to do with them. pars1 is for time 1 (or population 1).

pars2

is the same as pars2 but for time/population 2.

N

The number of intervals to integrate over.

...

optional parameters to pass on to func(). These are not decomposed.

Details

The decomposition works by assuming a linear change in all covariates between time 1 and time 2 (or population 1 and population 2). At each small time step approaching time 2 (the size of which is the inverse of N) each covariate is moved forward along its linear trajectory. One at a time, each covariate (of which there are ages*variables of) is switched out twice, once for its value at 1/(2N) forward and once for its value at 1/(2N) backward in time. The difference between func() evaluated with these two rate matrices is the change in yattributable to that particular covariate and that particular time step. Summing over all N time steps, we get the contribution to the difference of each covariate, effectmat. The sum of effectmat should come very close to func(rates2)-func(rates1). The error decreases with larger N, but there is not much point in having an N higher than 100, and 20 is usually sufficient. This ought to be able to handle a very wide variety of functions.

If pars1 are observations from 2005 and pars2 are observations from 2006 an N of 20 would imply a delta of 1/20 of a year for each integration step. Higher N provides finer results (a lower total residual), but takes longer to compute. In general, there are decreasing returns to higher N. sum(effectmat) ought to approximate func(rates2)-func(rates1).

References

andreev2002algorithmDemoDecomp andreev2002algorithmDemoDecomp

Examples

data(rates1)
data(rates2)

# we need rates1 and rates2 as vectors
rates1 <- c(rates1)
rates2 <- c(rates2)
# look at the function:
R0vec
# 2 things to point out:
# 1) it has an argument pfem, proportion female of births (1/(1+SRB)), 
#    that must be specified, but that we don't care about decomposing
# 2) x is a single vector. The the inside of the function needs to 
#    either refer to parts of it by indexing, as done here, or else 
#    re-assign x to various objects. In this case x[1:l] is Lx and 
#    x[(l+1):(2*l)] is Fx...
A <- horiuchi(func = R0vec,
              pars1 = rates1,
              pars2 = rates2,
              N = 10,
              pfem = .4886)
# the output, A, is also a single vector. Each element corresponds 
# to the effect of changes in that particular covariate toward the 
# overall change in the function value. sum(A) should be close to
# original difference
(check1 <- R0vec(rates2) - R0vec(rates1)) 
(check2 <- sum(A))

# \dontshow{
# de facto unit test. In this case the residual is very tiny,
# but if differences are very large and there are very many components,
# then the residual can be larger albeit trivial. In that case
# increase N and go make a coffee.
stopifnot(abs(check1 - check2) < 1e-6)
# }

# This package does not supply default plotting functions, but one 
# strategy might be the following:

# reorder A into a matrix (sideways):
A <- t(matrix(A,ncol=2))
# call barplot() twice, once for positive values and again for
# negative values
Apos <- A * .5 * (sign(A) + 1)      
Aneg <- A * .5 * abs(sign(A) - 1)   
if (FALSE) {
barplot(Apos, 
        width = rep(1, length(A) / 2),
        space = 0, 
        ylim = range(A), 
        main = "A fake decomposition of R0",
        col=c("yellow","green"),
        axisnames = FALSE,
        xlim=c(0, 90), 
        ylab = "contrib to change in R0",
        cex.axis = .7)
barplot(Aneg, 
        width = rep(1, length(A) / 2),
        add = TRUE, 
        space = 0,
        col = c("yellow", "green"),
        axes = FALSE, axisnames = FALSE)
segments(seq(from=0,to=90,by=10),0,seq(from=0,to=90,by=10),-.027,lty=2,col="grey")
text(seq(from=0,to=90,by=10),-.027,seq(from=0,to=90,by=10),pos=1,xpd=T)
legend("bottomright",fill=c("yellow","green"),legend=c("contrib from change in Lx",
"contrib from change in Fx"),title="age specific contrib of changes in Fx and Lx",bg="white") 
}