Counts and percentages

This table shows the overall estimated count as well as the counts and percentages by type of doctor, physician specialty, and metropolitan statistical area.

The variables that are necessary for creating this table are already in the survey, making the commands very straightforward.

total()

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

Rates

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

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

Here is the overall population estimate:

uspop2019$total
## [1] 323186697

Once we have the overall population estimate, 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 variable called Level that matches the levels of the variable in the survey, and a variable called Population that gives the population size (which is assumed to be a constant rather than a random variable).

For MSA, we can see the levels of the variables just by using the tab() command, just as we did above. Thus, to calculate rates, we need a data frame as follows:

uspop2019$MSA
##                                 Level Population
## 1 MSA (Metropolitan Statistical Area)  277229518
## 2                             Non-MSA   45957179

Now that we have the appropriate population estimates, the rate is:

tab_rate("MSA", uspop2019$MSA)

We can also calculate rates of a specific variable based on the entire population:

tab_rate("MDDO", uspop2019$total)
## * Rate based on the entire population.

tab_rate("SPECCAT", uspop2019$total)
## * Rate based on the entire population.

Counts and percentages

This table presents estimates for each age group, as well as for each age group by sex.

var_list("age")

The survey has a couple of relevant age-related variables. AGE is the patient age in years. AGER is a categorical variable based on AGE. However, for this table, in addition to AGER, we need another age group variable, with different age categories. We create it using the var_cut function.

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

Now that we’ve created the Age group variable, we can create the tables:

tab("AGER", "Age group", "SEX")

tab_cross("AGER", "SEX")

Rates

tab_rate("AGER", uspop2019$AGER)

tab_rate("Age group", uspop2019$`Age group`)
## * Population for some levels not defined: 15-64

tab_rate("SEX", uspop2019$SEX)

To calculate the rates for one variable (AGER) by another variable (SEX), we need population estimates 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

Once we have these population estimates, the rates are:

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

More advanced coding

In addition, for each age-sex category, the published table shows the percentage of preventive care visits made to primary care physicians.

To calculate these percentages, a slightly more involved for loop is needed. Below is the code, followed by an explanation:

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

for (vr in c("AGER", "Age group", "SEX", "Age x Sex")) {
    var_cross("tmp", "MAJOR", vr)
    for (lvl in levels(surveytable:::env$survey$variables[,vr])) {
        tab_subset("SPECCAT", "tmp", paste0("Preventive care: ", lvl))
    }
}
## Warning in var_cross("tmp", "MAJOR", vr): tmp: overwriting a variable that
## already exists.
## Warning in var_cross("tmp", "MAJOR", vr): tmp: overwriting a variable that
## already exists.
## Warning in var_cross("tmp", "MAJOR", vr): tmp: overwriting a variable that
## already exists.
set_output(csv = "")
## * Turning off CSV output.
## * ?set_output for other options.

• Since tab_subset() is called from within a for loop, if we wanted to print to the screen, we would need to use print( tab_subset(*) ). Since we don’t want to print to the screen, a call to print() is omitted.
• Since so many tables are being produced, the output is being sent to a CSV file.
• As before, the loop goes through the age, sex, and age / sex interaction variables, calling each of these variables vr.
• MAJOR and vr are crossed, with the result stored in a variable called tmp.
• Next, the inner loop goes through all levels of vr, calling each of these levels lvl.
• The code tabulates SPECCAT (Type of specialty – Primary, Medical, Surgical) on a subset in which tmp (which is MAJOR crossed with vr) is restricted to "Preventive care: " followed by lvl, which is some level of vr, such as “Under 15 years” for AGER.
• Finally, CSV output is turned off.

If you run this code, all of the tables should be stored in the CSV file. To give you an idea of what the tables should look like, here is just one of the tables:

vr = "AGER"
var_cross("tmp", "MAJOR", vr)
## Warning in var_cross("tmp", "MAJOR", vr): tmp: overwriting a variable that
## already exists.
lvl = levels(surveytable:::env$survey$variables[,vr])[1]
tab_subset("SPECCAT", "tmp", paste0("Preventive care: ", lvl))

To match the percentage in the published table, see the “Primary care specialty” row. Be sure to check the presentation standards flags.