The event and population data are at the core of the BYM-based models
used in the RSTr package. They work alongside the adjacency information
to generate smoothed estimates. In this vignette, we’ll discuss
requirements for event and population data and walk through an example
with a data.frame.
Data must be a list object with names Y
and n for the event counts and for the population counts,
respectively;
Y and n are intended to be
entire-population data. While it is possible to use RSTr to analyze
survey data or datasets that don’t include all members of a population
of interest, RSTr does not currently allow for the inclusion of survey
weights and thus assumes that each Y / n is an unbiased
estimate of the underlying event rate;
Y and n must contain real numbers.
Negative and infinite counts are not allowed, but suppressed data
containing NA’s is acceptable for the Y
values. Note, however, that n must have all population
counts;
For a given CAR model, Y must have at least one
total event. The CAR model will not be able to smooth information if
there is no event data present. For a CAR model, this includes any set
of all regions; for an MCAR model, this is any set of all regions and
groups; and for the MSTCAR model, this is the entire dataset;
Y and n can be up to a
three-dimensional array: the first margin (rows) specifies the region,
the second margin (columns) specifies the groups of interest, and the
third margin (matrix slices) specifies the time period;
Time periods, regions, and groups must be consistent. If your data contains counts for all regions in a specified set of groups for 1979 and 1981, for example, it must also include counts for all regions and all groups for 1980 as well, even if those counts have zero events;
Groups of many types are allowed as long as your sociodemographic groups are combined in the appropriate margin. For example, your groups may include just age groups, a mixture of age-sex groups, or even a mix of age-race-sex groups;
Finally, Y and n must have dimension
names associated with them. This makes for easy identification of
counties, groups, and time periods, and is necessary should you want to
age-standardize data using RSTr’s additional functionality.
To walk through the data setup from a data.frame to the
final array list, we will use data generated
by CDC WONDER’s Underlying Cause of Death Compressed Mortality, ICD-9
database, found at https://wonder.cdc.gov/cmf-icd9.html:
head(maexample)
#> Notes Year Year.Code County County.Code Sex Sex.Code Deaths
#> 1 1979 1979 Barnstable County, MA 25001 Female F 15
#> 2 1979 1979 Barnstable County, MA 25001 Male M 57
#> 3 1979 1979 Berkshire County, MA 25003 Female F 11
#> 4 1979 1979 Berkshire County, MA 25003 Male M 63
#> 5 1979 1979 Bristol County, MA 25005 Female F 52
#> 6 1979 1979 Bristol County, MA 25005 Male M 191
#> Population Crude.Rate
#> 1 25239 59.4 (Unreliable)
#> 2 21261 268.1
#> 3 24884 44.2 (Unreliable)
#> 4 22465 280.4
#> 5 80171 64.9
#> 6 71943 265.5Our example dataset contains acute myocardial infarction (ICD-9: 410)
mortality and population data in all counties of Massachusetts for men
and women aged 35-64 from 1979 to 1981. This dataset also includes some
notes in the bottom rows describing the dataset. maexample
contains several variables:
Notes: Provides general information about the
dataset, starting at row 85;
Year and Year.Code specify the
year;
County and County.Code specify the
county name and associated FIPS code;
Sex and Sex.Code specify the sex
group;
Deaths contains our mortality counts of
interest;
Population contains our population counts of
interest;
Crude.Rate shows the rates per 100,000 in each
year-county-sex group. For our purposes, this column can be
ignored.
The first thing we want to do with our dataset is remove the notes
from the bottom rows - while they are useful for getting acquainted with
the dataset, they will ultimately mess up our population arrays. Since
Year does not have information in rows with notes, we can
use that to filter our data:
The above code searches for values in maexample$Year
that aren’t NA and creates a new dataset containing only
those rows. Before we start generating our arrays, let’s take stock of
how our data is listed out:
head(ma_mort)
#> Notes Year Year.Code County County.Code Sex Sex.Code Deaths
#> 1 1979 1979 Barnstable County, MA 25001 Female F 15
#> 2 1979 1979 Barnstable County, MA 25001 Male M 57
#> 3 1979 1979 Berkshire County, MA 25003 Female F 11
#> 4 1979 1979 Berkshire County, MA 25003 Male M 63
#> 5 1979 1979 Bristol County, MA 25005 Female F 52
#> 6 1979 1979 Bristol County, MA 25005 Male M 191
#> Population Crude.Rate
#> 1 25239 59.4 (Unreliable)
#> 2 21261 268.1
#> 3 24884 44.2 (Unreliable)
#> 4 22465 280.4
#> 5 80171 64.9
#> 6 71943 265.5RSTr offers a long_to_list_matrix() function which can
transform this dataset into mortality and population arrays with
properly oriented margins:
If you want to manually set up the data, you can create
Y and n arrays using the xtabs()
function and consolidate them into a list to be used with
the model:
Y <- xtabs(Deaths ~ County.Code + Sex.Code + Year.Code, data = ma_mort)
n <- xtabs(Population ~ County.Code + Sex.Code + Year.Code, data = ma_mort)
ma_data <- list(Y = Y, n = n)Note that you must specify the names of each array element as above,
as creating a list with just the objects will not name each element, and
the names Y and n are necessary for RSTr to
know how to use the data.
If you have multiple types of groups, such as race and sex, it can take a little finessing to set up your group data, such as creating a combined race-sex group variable, but data setup will follow the same principles as above.
The above dataset is prepared specifically for an MSTCAR model. But what if we want to prepare data for an MCAR or even a CAR model? We can filter the original dataset and follow a similar procedure to prepare our data for the MCAR model:
ma_mort_mcar <- ma_mort[ma_mort$Year == 1979, ] # filter dataset to only show 1979 data
ma_data_mcar <- long_to_list_matrix(ma_mort_mcar, Deaths, Population, County.Code, Sex.Code)Note that xtabs() works by aggregating data along the
specified variables in the expression argument. In the case of the MCAR
model, we filter down to the year we want because otherwise, it would
give us the mortality and population counts for all years in our dataset
instead of just for 1979.
For the CAR model, setup is similar:
In this vignette, we used data generated from CDC WONDER to construct
our event and population counts, remove unnecessary rows using
filter(), and construct our list using
long_to_list_matrix(). Setting up the data for RSTr can
seem daunting at first, but with a few quick tricks in R, it can be easy
to have your data organized for analysis.