Preliminaries

head(iris)
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#> 1          5.1         3.5          1.4         0.2  setosa
#> 2          4.9         3.0          1.4         0.2  setosa
#> 3          4.7         3.2          1.3         0.2  setosa
#> 4          4.6         3.1          1.5         0.2  setosa
#> 5          5.0         3.6          1.4         0.2  setosa
#> 6          5.4         3.9          1.7         0.4  setosa

library(surveytable)

class(namcs2019sv_df$AGER)
#> [1] "factor"
class(namcs2019sv_df$PAYNOCHG)
#> [1] "logical"
class(namcs2019sv_df$AGE)
#> [1] "numeric"

attr(namcs2019sv_df$AGE, "label")
#> [1] "Patient age in years (raw - use caution)"

mysurvey = survey::svydesign(ids = ~ CPSUM
  , strata = ~ CSTRATM
  , weights = ~ PATWT
  , data = namcs2019sv_df)

attr(mysurvey, "label") = "NAMCS 2019 PUF"

all.equal(namcs2019sv, mysurvey)
#> [1] TRUE

Begin analysis

First, specify the survey object that you’d like to analyze.

library(surveytable)
set_survey(namcs2019sv)

Check the survey label, survey design variables, and the number of observations to verify that it all looks correct.

For this example, we do want to turn on certain NCHS-specific options, such as identifying low-precision estimates. If you do not care about identifying low-precision estimates, you can skip this command. To turn on the NCHS-specific options:

set_opts(mode = "NCHS")
#> * Mode: NCHS.

var_list("age")

Tabulate categorical and logical variables

The main function of the surveytable package is tab(), which tabulates variables. It operates on categorical and logical variables, and presents both estimated counts, with their standard errors (SEs) and 95% confidence intervals (CIs), and percentages, with their SEs and CIs. For example, to tabulate AGER, type:

tab("AGER")

The table title shows the variable label (the long variable name) and the survey label.

For each level of the variable, the table shows:

• the estimated count, its standard error, and its 95% confidence interval; and
• the estimated percentage, its standard error, and its 95% confidence interval.

Low-precision estimates. Optionally, the tab() function, as well as the other tabulation functions that are discussed below, can automatically identify low-precision estimates using algorithms developed at NCHS. For counts, rates, and percentages, the functions flag estimates if, according to the algorithms, they should not be presented, should be reviewed by a clearance official, or should be presented with a footnote. If no estimates are flagged by the checks, the table has a footnote that indicates this. If the checks do identify an estimate, that is denoted in an additional column and in the table footnote.

Turn on this functionality using any of the following: set_opts(lpe = TRUE), set_opts(mode = "nchs"), set_survey(*, mode = "nchs"), or options(surveytable.find_lpe = TRUE).

As an example, let’s tabulate PAYNOCHG:

tab("PAYNOCHG")

This table tells us that the estimated number of visits in which there was no charge for the visit has low precision. Intuitively, we can see that the CI for this count estimate is very wide, indicating high uncertainty.

The CIs that are displayed are the ones that are used by the NCHS presentation standards. Specifically, for counts, the tables show the log Student’s t 95% CI, with adaptations for complex surveys; for percentages, they show the 95% Korn and Graubard CI.

Drop missing values. Some variables might contain missing values (NA). Consider the following variable, which is not part of the actual survey, but was constructed specifically for this example:

tab("SPECCAT.bad")

To calculate percentages based on the non-missing values only, use the drop_na argument:

tab("SPECCAT.bad", drop_na = TRUE)

The above table gives percentages based only on the knowns, that is, based only on non-NA values.

Multiple tables. Multiple tables can be created with a single command:

tab("MDDO", "SPECCAT", "MSA")

total()

tab_subset("AGER", "SEX")

tab_cross("AGER", "SEX")

Tabulate numeric variables

The tab() and tab_subset() functions also work with numeric variables, though with such variables, the output is different. To tabulate NUMMED (number of medications), a numeric variable, type:

tab("NUMMED")

As before, the table title shows the variable label (the long variable name) and the survey label.

The table shows the percentage of values that are not missing (not NA), the mean, the standard error of the mean (SEM), and the standard deviation (SD).

Subsetting works too:

tab_subset("NUMMED", "AGER")

Perform statistical hypothesis testing

The tab_subset() function makes it easy to perform hypothesis testing by using the test argument. When the argument is TRUE, a test of association is performed. In addition, t-tests for all pairs of levels are performed as well.

tab_subset("AGER", "SPECCAT", test = TRUE)

tab_subset("MRI", "SPECCAT", test = TRUE)

tab_subset("NUMMED", "AGER", test = TRUE)

tab_subset("NUMMED", "SPECCAT", test = TRUE)

tab("SPECCAT", test = TRUE)

Calculate rates

A rate is a ratio of count estimates based on the survey in question divided by population size, which is assumed to be known. For example, the number of physician visits per 100 people in the population is a rate: the number of physician visits is estimated from the namcs2019sv survey, while the number of people in the population comes from another source.

To calculate rates, in addition to the survey, we need a source of information on population size. You would typically use a function such as read.csv() to load the population figures and get them into the correct format. The surveytable package comes with an object called uspop2019 that contains several population figures for use in these examples.

Let’s examine uspop2019:

class(uspop2019)
#> [1] "list"
names(uspop2019)
#> [1] "total"       "MSA"         "AGER"        "Age group"   "SEX"        
#> [6] "AGER x SEX"  "Age group 5"

The overall population size for the country as a whole is:

uspop2019$total
#> [1] 323186697

Once we have the overall population size, the overall rate is:

total_rate(uspop2019$total)

To calculate the rates for a particular variable, we need to provide a data frame with a column called Level that matches the levels of the variable in the survey, and a column called Population that gives the size of the population for that level.

For example, for AGER, this data frame as follows:

uspop2019$AGER
#>               Level Population
#> 1    Under 15 years   60526656
#> 2       15-24 years   41718700
#> 3       25-44 years   85599410
#> 4       45-64 years   82562049
#> 5       65-74 years   31260202
#> 6 75 years and over   21519680

Now that we have the appropriate population figures, the rates table is obtained by typing:

tab_rate("AGER", uspop2019$AGER)

To calculate the rates for one variable (AGER) by another variable (SEX), we need population figures in the following format:

uspop2019$`AGER x SEX`
#>                Level Subset Population
#> 1     Under 15 years Female   29604762
#> 2        15-24 years Female   20730118
#> 3        25-44 years Female   43192143
#> 4        45-64 years Female   42508901
#> 5        65-74 years Female   16673240
#> 6  75 years and over Female   12421444
#> 7     Under 15 years   Male   30921894
#> 8        15-24 years   Male   20988582
#> 9        25-44 years   Male   42407267
#> 10       45-64 years   Male   40053148
#> 11       65-74 years   Male   14586962
#> 12 75 years and over   Male    9098236

With this data frame, the rates table is obtained by typing:

tab_subset_rate("AGER", "SEX", uspop2019$`AGER x SEX`)

Create or modify variables

In some situations, it might be necessary to modify survey variables, or to create new ones. This section describes how to do this.

Convert factor to logical. The variable MAJOR (major reason for this visit) has several levels.

tab("MAJOR")

Notice that one of the levels is called "Preventive care". Suppose an analyst is only interested in whether or not a visit is a preventive care visit – they are not interested in the other visit types. They can create a new variable called Preventive care visits that is TRUE for preventive care visits and FALSE for all other types of visits, as follows:

var_case("Preventive care visits", "MAJOR", "Preventive care")
tab("Preventive care visits")

This creates a logical variable that is TRUE for preventive care visits and then tabulates it. When using the var_case() function, specify the name of the new logical variable to be created, an existing factor variable, and one or more levels of the factor variable that should be set to TRUE in the logical variable.

Thus, if an analyst is interested in surgery-related visits, which are indicated by two different levels of MAJOR, they could type:

var_case("Surgery-related visits"
  , "MAJOR"
  , c("Pre-surgery", "Post-surgery"))
tab("Surgery-related visits")

Collapse levels. The variable PRIMCARE (whether the physician is this patient’s primary care provider) has levels Unknown and Blank, among others.

tab("PRIMCARE")

To collapse Unknown and Blank into a single level, type:

var_collapse("PRIMCARE", "Unknown if PCP", c("Unknown", "Blank"))
tab("PRIMCARE")

Convert numeric to factor. The variable AGE is numeric.

tab("AGE")

To create a new variable of age categories based on AGE, type:

var_cut("Age group"
   , "AGE"
   , c(-Inf, -0.1, 0, 4, 14, 64, Inf)
   , c(NA, "Under 1", "1-4", "5-14", "15-64", "65 and over"))
tab("Age group")

In the var_cut() command, specify the following information:

• name of the new categorical variable;
• name of the existing numeric variable;
• cut points – note that the intervals are inclusive on the right; and
  
• category labels.

Be cognizant of any “special values” that the numeric variable might have. In some data systems, negative values indicate unknowns, which should be coded as NA. That’s what we do here – any value between -Inf and -0.1 gets coded as missing (NA). Though in this particular data, there are no unknowns and no “special values”.

Check whether any variable is true. For a series of logical variables, you can check whether any of them are TRUE using the var_any() command.

A physician visit is considered to be an “imaging services” visit if it had any of a number of imaging services ordered or provided. Imaging services are indicated using logical variables, such as MRI and XRAY. To create the Imaging services variable, type:

var_any("Imaging services"
  , c("ANYIMAGE", "BONEDENS", "CATSCAN", "ECHOCARD", "OTHULTRA"
  , "MAMMO", "MRI", "XRAY", "OTHIMAGE"))
tab("Imaging services")

Interact variables. The tab_cross() function creates a table of an interaction of two variables, but it does not save the interacted variable. To create the interacted variable, use the var_cross() command:

var_cross("Age x Sex", "AGER", "SEX")

Specify the name of the new variable as well as names of the two variables to interact.

Copy a variable. Create a new variable that is a copy of another variable using var_copy(). You can modify the copy, while the original remains unchanged. For example:

var_copy("Age group", "AGER")
#> Warning in var_copy("Age group", "AGER"): Age group: overwriting a variable
#> that already exists.
var_collapse("Age group", "65+", c("65-74 years", "75 years and over"))
var_collapse("Age group", "25-64", c("25-44 years", "45-64 years"))
tab("AGER", "Age group")

Here, the AGER variable remains unchanged, while the Age group variable has fewer categories.

Save the output

The tab* and total* functions have an argument called csv that specifies the name of a comma-separated values (CSV) file to save the output to. Alternatively, you can name the default CSV output file using the set_opts() function. For example, the following directs surveytable to send all future output to a CSV file, create some tables, and then turn off sending output to the file:

set_opts(csv = "output.csv")

tab("MDDO")

set_opts(csv = "")
#> * Turning off CSV output.

If the tabulation functions are called from within an R Markdown notebook or a Quarto document, they produce HTML or LaTeX tables, as appropriate. This makes it easy to incorporate the output of the surveytable package directly into documents, presentations, “shiny” web apps, and other output types.

Finally, the tabulation functions return the tables that they produce. More advanced analysts can use this functionality to integrate surveytable into other programming tasks.