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"

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

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

identical(namcs2019sv, mysurvey)
#> [1] TRUE

Begin analysis

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

set_survey(namcs2019sv)
#> * To adjust how counts are rounded, see ?set_count_int
#>                        _                                                                    
#> Survey name            NAMCS 2019 PUF                                                       
#> Number of variables    33                                                                   
#> Number of observations 8250                                                                 
#> Info1                  Stratified 1 - level Cluster Sampling design (with replacement)      
#> Info2                  With (398) clusters.                                                 
#> Info3                  survey::svydesign(ids = ~CPSUM, strata = ~CSTRATM, weights = ~PATWT, 
#> Info4                      data = namcs2019sv_df)

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

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.

NCHS presentation standards. The tab() function also applies the National Center for Health Statistics (NCHS) presentation standards for counts and percentages, and flags estimates if, according to the standards, they should be suppressed, footnoted, or reviewed by an analyst. 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.

One does not need to do anything extra to perform presentation standards checking – it is performed automatically. For example, let’s tabulate PAYNOCHG:

tab("PAYNOCHG")

This table tells us that, according to the NCHS presentation standards, the estimated number of visits in which there was no charge for the visit should be suppressed due to low precision. However, the lack of a percentage flag indicates that the estimated percentage of such visits can be shown.

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, 4, 14, 64, Inf)
        , c("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.

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.

The variables data frame. Recall that survey objects have an element called variables, which is a data frame that contains the survey variables.

class(namcs2019sv$variables)
#> [1] "data.frame"

More advanced users can create or modify variables in the variables data frame directly. After they modify these variables, they must call set_survey() again. For example:

namcs2019sv$variables$`Medicare and Medicaid` = ( 
  namcs2019sv$variables$PAYMCARE & namcs2019sv$variables$PAYMCAID)
set_survey(namcs2019sv)
#> * To adjust how counts are rounded, see ?set_count_int
#>                        _                                                                    
#> Survey name            NAMCS 2019 PUF                                                       
#> Number of variables    34                                                                   
#> Number of observations 8250                                                                 
#> Info1                  Stratified 1 - level Cluster Sampling design (with replacement)      
#> Info2                  With (398) clusters.                                                 
#> Info3                  survey::svydesign(ids = ~CPSUM, strata = ~CSTRATM, weights = ~PATWT, 
#> Info4                      data = namcs2019sv_df)

tab("Medicare and Medicaid")

Note, however, that the var_*() functions do not modify the survey object specified in set_survey() directly. Rather, they modify variables inside the following data frame: surveytable:::env$survey$variables. If you use the var_*() functions and then need access to the modified / created variables, that’s where you should look. For example:

var_cross("newvar", "MAJOR", "AGER")
# This should give NULL. The new variable does not exist here:
namcs2019sv$variables$newvar
#> NULL
# Rather, the new variable is here:
str(surveytable:::env$survey$variables$newvar)
#>  Factor w/ 42 levels "Blank : Under 15 years",..: 17 31 24 38 38 24 30 31 24 38 ...
#>  - attr(*, "label")= chr "(Major reason for this visit) x (Patient age recode)"

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_output() 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:

tmp_file = tempfile(fileext = ".csv")
suppressMessages( set_output(csv = tmp_file) )

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

set_output(csv = "")
#> * Turning off CSV output.
#> * ?set_output for other options.

If the tabulation functions are called from within an R Markdown notebook, they produce HTML tables. 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.