04-week04.Rmd

# Week 4 {#week4}

```{r, echo=FALSE, warning=FALSE, message=FALSE}
options(tigris_use_cache = TRUE)
pacman::p_load(demogR, demography, magrittr, knitr, kableExtra, readstata13, captioner,tigris, sf, tidyverse)

figure_nums <- captioner(prefix = "Figure")
table_nums <- captioner(prefix = "Table")

# path to this file name
if (!interactive()) {
    fnamepath <- current_input(dir = TRUE)
    fnamestr <- paste0(Sys.getenv("COMPUTERNAME"), ": ", fnamepath)
} else {
    fnamepath <- ""
}
```

<h2>Topics: Functions and sampling</h2>
This week's topics cover functions and sampling, the latter including a cursory treatment of loops and bootstrapping. We will revisit Ben's code for reading Human Mortality and Human Fertility databases to look at the use of `purrr::map()` as an alternative to `for()` looping in functions. Finally we will delve into `demogR` and `demography` packages for analysis of age-structured demographic models and forecasting mortality, fertility, migration and population data.

* [R environments](#renviron)
* [R functions](#rfunc)
* [Sampling in R](#rsampling)
* [Packages for running demographic analysis](#demogpackages)
    * [`demogR`](#demogR)
    * [`demography`](#demography)
    
Download [template.Rmd.txt](files/template.Rmd.txt) as a template for this lesson. Save in in your working folder for the course, renamed to `week_04.Rmd` (with no `.txt` extension). Make any necessary changes to the YAML header.


## Environments {#renviron}
Functions exist in `environments`, which are "frames" or collections containing objects including variables, functions, etc. There is a global environment (`.GlobalEnv`) that contains all of the data, functions, in the current R session. Those objects can be enumerated with the `ls()` function. 

The reason environments are mentioned at this time is because functions instantiate environments that exist while the function is running, but once the function is completed, the environment is removed from memory. More on this [below](#funcenviron).

## Functions {#rfunc}
Functions are sets of statements in R that are grouped together to perform a specific task or set of tasks. Functions are either built in, included in packages, or user-defined. Functions are used mainly to simplify running a series of individual commands or functions, for situations where the same process will need to be run multiple times on different inputs, or when control structures are needed (e.g., looping, logical branching).

### Function Components
The different parts of a function are:

1. (Usually) Name: This is the actual name of the function. It is stored in an R environment as an object with this name.
1. (Optional) Arguments: Arguments specify the inputs or options to the function. When a function is invoked, you pass a value to the argument. Arguments are optional; that is, a function may contain no arguments. Also arguments can have default values.
1. Body: The function body contains a collection of statements that defines what the function does.
1. Return value: The return value of a function is the last expression in the function body to be evaluated.

Another important concept for functions is environments.

#### Name
Most functions are created with code of the form

```
function_name <- function(argument(s)){
    statement(s)
}
```

For example, to square a vector of numerical values:

```{r}
f_square <- function(x){
    x^2
}
```

the function name is `f_square`.

Some functions are not named, and are referred to as "anonymous" functions. For example, functions can be used within the `apply` family of functions. Here is an oprationalized example.

```{r}
# create a list of three vectors of random numbers of different random lengths

# set.seed() makes the random process reproducible.
set.seed(10)
# vector lengths
v.len <- rnorm(n = 3, mean = 30, sd = 10) %>% round(0)

# make the random vectors
set.seed(5)
v1 <- rnorm(n = v.len[1])
set.seed(3)
v2 <- rnorm(n = v.len[2])
set.seed(6)
v3 <- rnorm(n = v.len[3])

# create the list
L <- list(v1, v2, v3)

# get the first value from each vector in the list
lapply(X = L, FUN = function(x) {x[1]})
```

The `lapply()` function is used to _apply_ a function to each element of a list. In this example, the last line of the code chunk is:

```
lapply(X = L, FUN = function(x) {x[1]})

```

in which the body of the function is `x[1]`, i.e., obtain the first element of a vector. In natural language, this translates to "for each element of the list _L_, obtain the first element of vector _x_". But the function itself is not named, and hence "anonymous."

#### Arguments
Most functions require arguments. Arguments are used to instantiate variables within the function's environment that can be used later in the body of the function. Each argument is named, and the name is used within the function as a local variable within the function's environment.

Following our example from above, `f_square` takes an argument named "x" that is a numeric vector.

Here, let's modify the function to demonstrate that within the environment of the function, `x` is a variable by using `print(x)`:

```{r}
f_square_2 <- function(x){
    message("input:")
    print(x)
    message("output:")
    x^2
}

f_square_2(c(1,2,3))
```

We can try running the original function using different (or no) arguments:

Here, using a vector of a single NA numeric

```{r}
f_square(as.numeric(NA))
```

... or a vector that contains a numeric NA (with the first two elements being numeric, the third element `NA` is automatically cast as a numeric):

```{r}
f_square(c(1, 2, NA))
```

... or a null:

```{r}
f_square(NULL)
```

... or a vector containing a null:

```{r}
f_square(c(1, 2, NULL))
```

... or with no argument at all:

```
f_square()
```

<font color="red">
```
## Error in f_square() : argument "x" is missing, with no default
```
</font>

Some functions do not require arguments, e.g., to get the current date or time:

```{r}
Sys.Date()
Sys.time()
```

... and if we try to use an argument we get an error:

```
Sys.Date(1)
```

<font color="red">
```
## Error in Sys.Date(1) : unused argument (1)
```
</font>

##### Default values for arguments
If you want an argument to have a default value, it is specified in the listed arguments in the form `argument = value`. 

Following our previous `f_square_2()` function, we can set the default value of `x` to `3`:

```{r}
f_square_3 <- function(x = 3){
    x^2
}
```

Because the default argument is set, the function can be run with no arguments, and the default will be substituted in:

```{r}
f_square_3()
```


A more meaningful example demonstrates stratification of counts into intensity bins using accelerometry data. We will be using accelerometry from one day's data from one subject in a study.

The cut points for accelerometry were identified at 0, 2690, 6167, and 9642 counts per minute, from Sasaki et al. (2011).

*Sasaki JE, John D, Freedson PS. Validation and comparison of ActiGraph activity monitors. J Sci Med Sport. 2011;14(5):411-416. doi:10.1016/j.jsams.2011.04.003"*

The variable `vm3` is the vector magnitude measured with the accelerometer for each minute. Data: [accelerometry.csv](files/accelerometry.csv).

```{r}
# read the accelerometry data
acc <- read.csv("files/accelerometry.csv")

# print first 6 lines
head(acc)
```


The following function codes intensity by the aforementioned cut points by default and using default labels:

```{r}
f_acc_intensity <- function(x,
                            cuts = c(
                                -Inf,
                                2690, 6167, 9642, Inf
                            ),
                            labels = c(
                                "sedentary/low",
                                "moderate", "vigorous", "very vigorous"
                            )) {
    cut(x = acc$vm3, breaks = cuts, labels = labels)
}
```

... and when run with the defaults to tabulate the minutes spent in different PA levels

```{r}
# recode
acc$intens_default <- f_acc_intensity(acc$vm3)

# tabulate
acc %>% 
    group_by(intens_default) %>% 
    summarise(n = n())
```

But we could run this with different thresholds and levels, where SPLA = "sedentary/low physical activity" and MVPA = "moderate-to-very vigorous physical activity):

```{r}
acc$intens_2lev <- f_acc_intensity(x = acc$vm3,
                                cuts = c(-Inf, 2690, Inf),
                                labels = c("SLPA", "MVVPA"))
acc %>% 
    group_by(intens_2lev) %>% 
    summarise(n = n())
```


##### The `...` argument
When functions do not have a known _a priori_ number or set of arguments, or when a large number of arguments is to be passed to another function the `...` argument is used. We will not cover this here, but you are encouraged to read more: [How to Use the Dots Argument in R](https://www.dummies.com/programming/r/how-to-use-the-dots-argument-in-r/); [The three-dots construct in R](https://www.r-bloggers.com/2013/01/the-three-dots-construct-in-r/).

#### Body
The function's body contains all of the code to perform the purpose of the function. Following our initial example, the body of the function is simply 

```
x^2
```

The body can be as simple or complicated as it needs to be in order to achieve the desired result.

### Logical testing for branching
Sometimes you want to vary what your code does based on some condition. If the condition is met, then execute a block of code. If the condition is not met, or some other condition is met, then execute other code. For a more complete tutorial, see [R if else elseif Statement](https://www.learnbyexample.org/r-if-else-elseif-statement/)


Following our previous `f_square_2()` function, we can modify to print the input based on the logical argument `verbose`. In this code, if the `verbose` object is set to `TRUE`, some text is printed in a message as well as the value that was supplied for `x`. In either case (`verbose = TRUE` or `verbose = FALSE`), the output value of `x^2` is returned.

```{r, collapse=TRUE}
f_square_4 <- function(x, verbose = FALSE){
    # only run the next lines if verbose is true
    if(verbose){
        message("input:")
        print(x)
        message("output:")
    }
    x^2
}
```

Here we run with the default option `verbose = FALSE`; only the final value `x^2` is in the output:

```{r}
f_square_4(x = c(1, 2, 3))
```

... and with `verbose = TRUE`, additional text is printed :


```{r}
f_square_4(x = c(1, 2, 3), verbose = TRUE)
```

There are additional ways to use `if()`, as nested statements, using:

```
if(condition1){
    if(condition2){
        statement
    }
}
```

Or with an alternative:

```
if(condition){
    statement1
    statement2
    ...
} else {
    statement3
    statement4
}
...
```

Or with an additional condition to check if the first condition is not met:


```
if(condition1){
    statement1
    statement2
    ...
} else if (condition2){
    statement3
    statement4
} else {
    statement5
    statement6
}
...
```

### Return value
The return value is either the last evaluated expression in the function or an object specified using the `return()` function. For functions that are intended to return only one value, by convention that is the last line in the function.

In our original `f_square()` function, the return value is `x^2` since no other `return()` value was specified, e.g., for a vector of one element:

```{r}
f_square <- function(x){
    x^2
}
f_square(3)
```

or a vector with multiple elements:

```{r}
f_square(c(1,2,3))
```

However, if it is possible that different outputs can be produced by a function based on some logical testing, one can explicitly use `return(object)` in the code; at that time the object will be output and the function will stop. A simple example of explicitly specifying return values is shown in this numerical comparison function:

```{r}
f_compare <- function(x, y){
    # either missing?
    if(nargs() != 2)
        return("invalid number of arguments")
    # numeric?
    if(!is.numeric(x) | !is.numeric(y)){
        return(sprintf("%s or %s is not numeric.", x, y))
    }
    # comparisons follow
    if(x > y){
        return(sprintf("%s is greater than %s", x, y))
    } 
    if(x < y) {
        return(sprintf("%s is less than %s", x, y))
    }
    if(x == y){
        return(sprintf("%s equals %s", x, y))
    }
}
```

Based on criteria such as the number of arguments or the relative value of the arguments `x` and `y`, different outputs are generated. Here are a few examples running the function:

```{r}
f_compare(1)

f_compare(1, 2)

f_compare(2, 1)
```

If you want to handle an expected error, you can print an informative message and use `return(invisible())`, which returns nothing at all (whereas `return()` results in a `NULL` object) e.g., here without `invisible()`:

```{r}
f_readfile <- function(fname){
    if(!file.exists(fname)){
        warning(paste(fname, "does not exist!"))
        return()
    } else {
        read.csv(fname)
    }
}

f_readfile("foobar.txt")
```

... and with `invisible()`:

```{r}
f_readfile <- function(fname){
    if(!file.exists(fname)){
        warning(paste(fname, "does not exist!"))
        return(invisible())
    } else {
        read.csv(fname)
    }
}

f_readfile("foobar.txt")
```

### Function environments {#funcenviron}
As mentioned before, functions instantiate environments that exist only while the function is being evaluated. This means that functions can include named variables that have the same name as a variable in a different environment. For example here is a function that only lists what objects are in the local and global environments:

```{r}
# declare a few variables
x <- 1
y <- "hello"

# a simple function
f <- function(x){
    # create a local variable
    y <- x + 2
    # another function inside this function
    g <- function(x){
        x * 3
    }
    # what variables are in this environment?
    print("----------")
    print("objects in this function's environment:")
    print(ls())
    # what is in the global env?
    print("----------")
    print("objects in the global environment:")
    print(ls(envir = .GlobalEnv))
    # return the output of the function
    print("----------")
    y
}

f(1)
```

Once the function completes, all objects in its local environment are purged. If you are running a complicated function that creates intermediate values that you want to examine for troubleshooting, you can create an environment in the `.GlobalEnv` (i.e., the overall environment in which R is running), and then assign a variable with a specific value to that environment created specifically for examining the intermediate products of the function:

```{r}
# a function to show how we can assign
g <- function(x){
    # code for a bunch of complicated operations
    # ...
    # create environment "foo"
    if(!exists("foo")){
        message("make foo")
        assign(x = "foo", value = new.env(), envir = .GlobalEnv)
    }    
    # generates an intermediate data frame named "bar"
    bar <- head(iris)
    # save to the foo env
    assign(x = "bar", value = bar, envir = foo)
    # more code to do more complicated stuff
    # ...
    foobar <- head(cars)
    # also assign 
    assign(x = "foobar", value = foobar, envir = foo)
    # yet more complicated stuff here
    # ...
}
```

When the function runs, the objects `bar` and `foobar` are placed in the `foo` environment. We can examine those:

```{r}
# run the function
g()

# what is in environment "foo"?
ls(envir = foo)
```

And we can view their values:

```{r}
print(foo$bar)
```

```{r}
print(foo$foobar)
```

But those objects will not appear in the `.GloalEnv`, even though the environment `foo` does; you could consider that the objects `bar` and `foobar` are "nested" in environment `foo` which is itself nested in `.GlobalEnv`.

```{r}
ls(envir = .GlobalEnv)
```

### Looping 
Loops are run when you want to perform a series of tasks over and over on a set of objects. 

#### Looping with `for()` loops
 A loop is constructed using the form ...

```
for (element in set){
    do something
}
```

where `element` represents the variable value of the current iteration and `set` is a group of objects, such as elements in a vector or list

A few simple examples follow.

Print each letter in the first 5 letters of the alphabet. In this example, the iterated element is a letter, and the variable `i` has the value of the letter in each iteration. The built in vector `letters` is used. So on the first iteration, `i = 1`, on the second iteration,m `i = 2` and so on.

```{r}
for (i in head(letters, 5)){
    print(i)
}
```


In this example, the object `i` takes on integer values 1 through 10. Within the loop, a message is printed that includes the term `i^2`, which evaluates to `1^2`, `2^2`, $\dots$, `10^2`.

```{r}
for(i in 1:10){
    message(paste(i, "squared equals", i^2))
}
```

Here is a similar approach, but using a numerical index referring to the position within a vector. The vector to be iterated over is `1:length(states_5)`, which evaluates to ``r 1:length(head(state.name, 5))``. In each iteration, we print the index of the iterator and the state name.

```{r}
# take the first 5 state names
states_5 <- head(state.name, 5)

# iterate over those
for (i in 1:length(states_5)){
    s <- states_5[i]
    message(paste0(i, ": ", s))
}
```

In this example, we will use an index that is the position of rows within a data frame.

```{r}
# initialize a value to hold a sum
mySum <- 0

# create a data frame of 5 rows
y <- iris %>% head(5)

# loop
for(x in 1:nrow(y)){
    message(paste("iteration:", x))
    # get the sepal length for this iteration
    sl <- y$Sepal.Length[x]
    # add to make a cumulative sum
    mySum <- mySum + sl
    # calculate the mean
    myMean <- mySum / x
    message(paste0("  sepal length = ", sl, 
                   "; cumulative sum = ", mySum, 
                   "; mean = ", myMean))
  
}
```

Here we will iterate over the elements of a list, also including the system run time by wrapping the loop in a `system.time()` call.

```{r}
# make it reproducible
set.seed(5)

# create a list
L <- list(
    v1 = rnorm(n = 10, mean = 1),
    v2 = rpois(n = 20, lambda = 5),
    v3 = runif(n = 25, min = 0, max = 10)
)

# run the loop
system.time(
    for(i in L){
        print(mean(i))
    }
)
```

#### Alternatives to loops: `apply()`. `lapply()`, `tapply()`, `sapply()`
There is a family of functions (`*apply()`) in R that are built to iterate over a matrix, or elements of vectors or lists.

The `apply` function performs a function over the rows or columns of a matrix. Here are some simple examples that get the sums of the rows and columns of a matrix:

```{r}
# a matrix
(m <- matrix(1:9, nrow = 3))


# row sums
message("row sums")
(rs <- apply(X = m, MARGIN = 1, FUN = sum))

# column sums
message("column sums")
(cs <- apply(X = m, MARGIN = 2, FUN = sum))

```

The `lapply()` function runs over elements of a list. Here we will run the same analysis as before (mean of the vectors in list `L`). The output of `lapply()` is a list with the same number of elements as the input; the output elements also have the original list element names.

```{r}
system.time(
    Lmean <- lapply(X = L, FUN = mean)
)

Lmean
```

Note the difference in `system.time()` when using `lapply()`. In general one should always use the `*apply()` functions rather than loops if possible. Loops are better for multi-step operations versus operations on single matrices, vectors, or lists.

You should look up the help documentation for `tapply()`, `lapply()`. There is also a nice tutorial, [Tutorial on the R Apply Family](https://www.datacamp.com/community/tutorials/r-tutorial-apply-family).

#### `purrr` 
`purrr` aims to simplify looping, particularly for multiple-step processes.

The main functions in `purrr` are `map_*()`, with specific variants to be run on data frames (`map_dfr()`), character vectors (`map_chr()`), logical (`map_lgl()`), integer (`map_int()`), and double precision numeric (`map_dbl()`).

Although a full treatment of `purrr` is beyond our scope, see a [great `purrr` tutorial with worked examples](https://jennybc.github.io/purrr-tutorial.)

##### Efficient download of census data using `purrr` and `tigris`

Here, following on the lessons from using the `tigirs` package for downloading US Census TIGER/Line data, we will download all counties for all states using `purrr::map_dfr()` and export to a GPKG database.

```{r}
# the tigris download function for counties
f_county <- function(state_name, year = 2019){
    # this downloads a single county
    counties(state = state_name, year)
}

# the map_dfr() functionality is wrapped in here
get_counties <- function(year = 2019) {
    # the output is a data frame. input (.x) is the state name, the function (.f) is the county download
    map_dfr(
        # input is the built in vector "state.name"
        .x = state.name,
        # each iteration runs the f_county() function over the iteration's state
        .f = function(x) {
            f_county(state_name = x, year = year)
        }
    )
}

# run the function
all_counties <- get_counties()

# export to GPKG
myTmpDir <- tempdir()
st_write(obj = all_counties, dsn = file.path(myTmpDir, "counties.gpkg"), layer = "us_counties_2019", delete_dsn = TRUE)
```

In `r figure_nums(name = "countymap", display = "cite")` we show the counties displayed in QGIS.

![](images/week04/2022-01-26 23_42_02-Window.png)
\    
_`r figure_nums(name = "countymap", caption = "Downloaded counties for the US")`_

This process could easily be applied to multiple years of data or different census units (block groups, blocks, tracts).

##### Revisiting Ben's HMDHFDplus example
Let's revisit Ben's [example using HMDHFDplus](#benhmdhfdplus)

```{r}
# Help function to list the available countries
# this generates a vector of all country abbreviations
countries <- HMDHFDplus::getHMDcountries()[1:2]

# Function to download a specified HMD data set item for a single county
# the country code is referenced as "CNTRY"
# the "item" is the base name of the link with ".txt" removed. For example,
# https://www.mortality.org/hmd/ISR/STATS/Mx_1x1.txt
#                                         Mx_1x1       <<- this is the item for 1 year x 1 year death rates
read_hmd_country <- function(CNTRY, item) {
  HMDHFDplus::readHMDweb(
    # the country from the function call
    CNTRY = CNTRY,
    # the item to download
    item = item,
    # the username from this key's record
    username = keyring::key_list("human-mortality-database")$username,
    # the password for this key's record
    password = keyring::key_get(
      service = "human-mortality-database",
      username = keyring::key_list("human-mortality-database")$username
    )
  )
}

# Help function to list the available countries
# this generates a vector of all country abbreviations
# for speed in class we will use only two countries
countries <- HMDHFDplus::getHMDcountries()[1:2]

# Download a data set iteratively for all countries using purrr::map()
# In this case, age-specific mortality in 1-year periods x 1-year age groups
# for all 1-year periods available
# output is a data frame named "mx_1x1"
mx_1x1 <- countries %>%
    # Returns a list of data.frames, adding a column for country code to each
    # the map() function performs a run of Ben's read_hmd_country() function for each listed country
    purrr::map_dfr(function(country) {
        # the item to read is 1 x 1 death rates
        read_hmd_country(CNTRY = country, item = "Mx_1x1") %>%
            # this adds the column "country" storing the country ISO code
            dplyr::mutate(country = country)
    }) %>%
    # Phil added this to make it a tibble
    tibble()
```

The key piece here is 

```
purrr::map_dfr(function(country) {
    # the item to read is 1 x 1 death rates
    read_hmd_country(country, "Mx_1x1")
```

which constructs a data frame by iterating over `countries` and applying the function `read_hmd_country()` to download the variable `Mx_1x1`. Here we plot mortality data from Australia and Austria over the time series (`r figure_nums(name = "aus_aut_mort", display = "cite")`).

```{r}
tmx <- mx_1x1 %>% 
    group_by(Year, country) %>% 
    dplyr::summarise(mortality = sum(Total, na.rm = TRUE), .groups = "keep")

tmx %>% ggplot(mapping = aes(x = Year, y = mortality)) +
    geom_line() +
    facet_grid(country ~ .) + 
    xlab("year")
```
\    
*`r figure_nums(name = "aus_aut_mort", caption = "Mortality data from Australia and Austria")`*

#### Bootstrapping 
Bootstrapping is used to estimate the characteristics of a population by repeating a large number of random samples (with replacement) and then calculating summary statistics on the set of samples. The reason replacement is used is that each sample represents a "snapshot" of the population. As the number of samples increases, the summary continues to approach the true population characteristics.

We will use a simple example from a hypothetical population. If we had no idea what the male/female split was, we could generate a large number of small-ish samples and the take the mean.

Back to the bootstrap: we will use 5,000 samples of size 100 from our hypothetical population to estimate the proportion of females.

```{r}
# create the population
# 1 indicates female and 0 indicates male
pop <- c(rep(1, 5e4 * 3 / 5), rep(0, 5e4 * 2 / 5))

# initialize a vector
F <- NULL

# run the bootstrap
for (i in seq(from = 1, to = 5000, by = 1)){
    # sample once
    s <- sample(x = pop, size = 100, replace = TRUE)
    # calculate percent female
    p <- sum(s) / length(s)
    # concatenate the result with the running result
    F <- c(F, p)
}

# mean and standard deviation of the bootstrap
mean(F)
sd(F)
```

Using the bootstrap results, we can estimate the 95% confidence interval around the mean using the `Rmisc::CI()` function, and make a density plot showing the mean and the confidence interval (`r figure_nums(name = "bootstrat_plot", display = "cite")).

```{r}
# 95% CI
ci_95 <- Rmisc::CI(x = F, ci = 0.95)

# plot with 95 % CI
plot(density(F), main = "")
abline(v = ci_95, col=c(2,1,2))
```
\    
*`r figure_nums(name = "bootstrat_plot", caption = "Density plot from bootstrapping")`*

If we already had an enumeration of the population, this method would be unnecessary. However, most survey-derived data are generated from samples. If the sample is representative of the underlying population,  the sample can be considered an acceptable proxy.


## Sampling {#rsampling}
The general topic of sampling is far greater than what we can cover in this course. However, a few examples may be helpful for students who have little experience in sampling.

The main R function for sampling is `sample()`, which draws a random sample from a vector or data frame. Options for `sample()` include the size of the sample, whether or not to replace sampled individuals to the pool, and a probability weight of the same length (or number of rows) as the population from which the sample is drawn.

### Sampling with replacement
Sampling with replacement is similar to rolling a die. Each roll of the die has a 1/6 probability of coming up with each value. For example we simulate rolling a die 100,00 times:

```{r}
# set.seed makes it so that random processes can be repeated with the same outputs
set.seed(5)

# create a vector to represent the faces of a die
d <- 1:6

# now make a sample of 100,000 rolls of the die, with replacement for independent rolls
s <- sample(x = d, size = 100000, replace = TRUE)

# tabulate the result
(tbs <- table(s))

# proportions
(prop.table(tbs))
```

The probability of rolling two of the same number in a row is 1/6 $\times$ 1/6, or $\approx$ 0.028

### Sampling without replacement
Sampling without replacement removes from the population the individual that was sampled. For example, the probability of selecting an ace from a whole deck of cards is 4/52, or $\approx$ 0.077. If we drew once from the deck and returned the card back to the deck to draw again, both samples would be independent. The probability of drawing two aces would be 4/52 $\times$ 4/52, or $\approx$ 0.0059. _That_ is an example of sampling _with_ replacement, similar to rolling a die.

Performing the same sample _without_ replacement, i.e., finding an ace, removing it, and then sampling again would have a probability of drawing two aces being 4/52 $\times$ 3/51, or $\approx$ `r round((4/52) * (3/51), 4)`.

If you were to design an study in which once a person was surveyed once they would not be eligible to be surveyed again, that would also be an example of samping *without* replacement.

For a population of sufficient size relative to the sample, sampling with or without replacement will lead to nearly identical results.

Let's use an example of a population of 50,000 persons, where 30,000 are female and 20,000 are male. If we were to sample two persons with replacement, the probability that they would both be female would be 30000/50000 $\times$ 30000/50000, or 0.36.

If we sample without replacement, the probability of selecting a female in the first sample of one would be 30000/50000, and the second sample of 1 would have a probability of selecting a female 29999/49999, with a joint probability of $\approx$ 0.359995.

## R packages for demographic analysis {#demogpackages}
There is always a growing list of available packages for analysis in R. Here we will briefly cover two packages for demographic analysis.

### `demogR` {#demogR}
Age-Structured Demographic Models

The [`demogR`](https://cran.r-project.org/web/packages/demogR/) package provides a number of functions for analysis of age-structured demographic models. In this exercise we will look at the `life.table()` function for creating life tables from enumerated deaths and mid-interval population estimates.

We will follow the example from the help page on `life.table()`. The example uses the `goodman` data set, which has aggregate data from Venezuela, Madagascar, and the USA from 1965, 1966, and 1967, respectively. The table contains these fields:

* age\
age classes
* ven.nKx\
mid-year population structure for Venezuela
* ven.nDx\
enumerated deaths for Venezuela
* ven.bx\
enumerated births for Venezuela
* mad.nKx\
mid-year population structure for Madagascar
* mad.nDx\
enumerated deaths for Madagascar
* mad.bx\
enumerated births for Madagascar
* usa.nKx\
mid-year population structure for the United States
* usa.nDx\
enumerated deaths for the United States
* usa.bx\
enumerated births for the United States

The first step is to create the life table itself, which should be run with at least three arguments, `x` (age at the beginning of the age classes of the life table), `nDx` (deaths), and `nKx` (population size). We will create the life table using the three countries' data and the default Keyfitz and Fleiger (`kf`) calculation type. Results are in `r table_nums(name = "kf", display = "cite")`.

```{r, warning=FALSE}
# load the Goodman data from the demogR package
data(goodman)

## default type="kf", Venezuela data
vlt <- with(goodman, life.table(x=age, nKx=ven.nKx, nDx=ven.nDx))

## US life table
ult <- with(goodman, life.table(x=age, nKx=usa.nKx, nDx=usa.nDx))

## Madagascar life table
mlt <- with(goodman, life.table(x=age, nKx=mad.nKx, nDx=mad.nDx))

# some values for the text
vlx35 <- vlt %>% filter(x == 35) %>% pull(lx) %>% round(2)
vlx80 <- vlt %>% filter(x == 80) %>% pull(lx) %>% round(2)
ulx35 <- ult %>% filter(x == 35) %>% pull(lx) %>% round(2)
ulx80 <- ult %>% filter(x == 80) %>% pull(lx) %>% round(2)
mlx35 <- mlt %>% filter(x == 35) %>% pull(lx) %>% round(2)
mlx80 <- mlt %>% filter(x == 80) %>% pull(lx) %>% round(2)

# combine these
lt <- bind_rows(vlt, ult, mlt)
```

*`r table_nums(name = "kf", caption = "Venezuela life table, Keyfitz and Fleiger method")`*

```{r, warning=FALSE}
lt %>% kable() %>% 
    kable_styling(bootstrap_options =
                      c("striped", "hover", "condensed", "responsive"), 
                  full_width = F, position = "left", font_size = 12, fixed_thead = T) %>% 
    pack_rows(group_label = "Venezuela", start_row = 1, end_row = 19) %>% 
    pack_rows(group_label = "USA", start_row = 20, end_row = 38) %>% 
    pack_rows(group_label = "Madagascar", start_row = 39, end_row = 57)
```

The results have nine output columns:

1. x\	
age at the beginning of the interval
1. nax\
person-years lived by those dying in the interval x to x+n
1. nMx\
period central death rate
1. nqx\
probability of death between ages x and x+n
1. lx\
probability of survival to exact age x
1. ndx\
proportion of deaths occurring between ages x and x+n
1. nLx\
person-years lived in the interval x to x+n
1. Tx\
person-years of life left in the cohort at age x
1. ex\
life expectancy at age x

For example, the life table estimates in Venezuela in 1965, a probability of reaching the age of 35 was `r vlx35`, whereas reaching 80 had a probability of `r vlx80`. In the US in 1967, the respective probabilities were `r ulx35` and `r ulx80`, and in Madagascar in 1966 the probabilities were `r mlx35` and `r mlx80`.

Another way to look at the data is to plot the life expectancy, shown in `r figure_nums(name = "demogrltplot", display = "cite")`. The difference between Venezuela and the US were less stark than the difference between Madagascar and the other two countries. In 1966, a woman from Madagascar had a less than 50% chance of reaching age 40! One wonders what these life expectancies look like with more recent data.

```{r}
lt %<>% mutate(
    country = c(rep("VEN", 19), rep("USA", 19), rep("MDG", 19))
)
ggplot(data = lt, mapping = aes(x = x, y = lx, col = country))+
    geom_line() +
    xlab("age") +
    ylab("probability of surviving to age on X axis")
```
\    
*`r figure_nums(name = "demogrltplot", caption = "Life expectancy, Madagascar (1966), Venezuela (1965), USA (1967)")`*

### `demography` {#demography}
The [`demography`](https://cran.r-project.org/web/packages/demography/index.html) also contains a large number of functions for demographic analysis. For other reference, see [An R intro to the demography package](https://rpubs.com/Timexpo/487053).

We will look at one example from the documentation.

```{md}
france.lt <- lifetable(fr.mort)
plot(france.lt)
lt1990 <- print(lifetable(fr.mort,year=1990))
```

The first step is creating the life table using the `lifetable()` function. The function requires data of a specific format, which can be provided in the package's `fr.mort` data set. What is in the `fr.mort` data set? Let's see:

```{r}
fr.mort
```

This tells us there were time series for females, males, and all persons, spanning years 1816 to 2006 and covering ages 0 to 110.

We can delve deeper:

```{r}
attributes(fr.mort)
```

Here we see there are specific objects within the data set. Let's look at some of them. `age`:

```{r}
fr.mort$age
```

`year`:

```{r}
fr.mort$year
```

The population data are stored in the `pop` object, which is composed of three matrices representing all persons, females, and males:

```{r}
fr.mort$pop %>% names()
```

Each of these matrices has the same number of rows and columns, corresponding to the ages (0-110, 111 rows) and years (1816-2006, columns). Here let's look at a sample of the total, the first 5 rows and columns:

```{r}
fr.mort$pop$total[1:5, 1:5]
```

This shows for example in 1816 there were estimated to be `r prettyNum(fr.mort$pop$total[1, 1])` persons less than one year of age.

Output from running the `lifetable()` function can be plotted. `r figure_nums(name = "france_ltplot", display = "cite") shows curves for each year using a rainbow palette. It is interesting to note that the earliest years in the time series had more positive life expectancies than some of the intermediate years particularly for persons below the age of 40.

```{r}
# min and max years
y0 <- min(fr.mort$year)
y1 <- max(fr.mort$year)

# make the life table and plot life expectancy graph
france.lt <- lifetable(fr.mort)
plot(france.lt)
```
*`r figure_nums(name = "france_ltplot", paste("Life expectancy in France from", y0, "to", y1))`*

Finally, let's look at the life table from 2005. 

```{r}
# create the life table
lt_france_2005 <- lifetable(fr.mort,year=2005)

# print using the native interface
print(lt_france_2005)
```
The output of `lifetable` prints in a tabular format but is not a data frame and so the various methods of making nice looking tables do not work out of the box. This also means standard exporting functions such as `write.csv()` will not work. However, with some work the specific objects returned from the function, we can construct a decent looking table.

```{r, message=FALSE, warning=FALSE}
# I looked at the individual vectors output from the lifetable function. These are the ones we want:
positions <- c(1,3:10)

# use lapply() across some of the list objects to
#   get the object (a vector), that is the anonymous "function(x) x"
# use "do.call" to do something to each of those vectors. that something is a "cbind"
# turn it into a data frame and then a tibble
lt_france_2005_df <- do.call(cbind, lapply(lt_france_2005[positions], function(x) x)) %>% 
    data.frame %>% tibble()

# set the column names from the names of the vectors
colnames(lt_france_2005_df) <- names(lt_france_2005)[positions]

# print a nice table
lt_france_2005_df %>% kable() %>% 
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = F, position = "left", font_size = 12, )
```

<hr>
Rendered at <tt>`r Sys.time()`</tt>

## Source code
File is at `r fnamestr`.

### R code used in this document
```{r ref.label=knitr::all_labels(), echo=TRUE, eval=FALSE}
```

### Complete Rmd code
```{r comment=''}
cat(readLines(fnamepath), sep = '\n')
```