Introduction to PHEindicatormethods

Introduction

This vignette introduces the following functions from the PHEindicatormethods package and provides basic sample code to demonstrate their execution. The code included is based on the code provided within the ‘examples’ section of the function documentation. This vignette does not explain the methods applied in detail but these can (optionally) be output alongside the statistics or for a more detailed explanation, please see the references section of the function documentation.

The following packages must be installed and loaded if not already available

library(PHEindicatormethods)
library(dplyr)

Package functions

This vignette covers the following core functions available within PHEindicatormethods:

Function Type Description
phe_proportion Non-aggregate Performs a calculation on each row of data (unless data is grouped)
phe_rate Non-aggregate Performs a calculation on each row of data (unless data is grouped)
phe_mean Aggregate Performs a calculation on each grouping set
calculate_dsr Aggregate, standardised Performs a calculation on each grouping set and requires additional reference inputs
calculate_ISRatio Aggregate, standardised Performs a calculation on each grouping set and requires additional reference inputs
calculate_ISRate Aggregate, standardised Performs a calculation on each grouping set and requires additional reference inputs

Other functions are introduced in separate vignettes.

Non-aggregate functions

Create some test data for the non-aggregate functions

The following code chunk creates a data frame containing observed number of events and populations for 4 geographical areas over 2 time periods that is used later to demonstrate the PHEindicatormethods package functions:

df <- data.frame(
        area = rep(c("Area1","Area2","Area3","Area4"), 2),
        year = rep(2015:2016, each = 4),
        obs = sample(100, 2 * 4, replace = TRUE),
        pop = sample(100:200, 2 * 4, replace = TRUE))
df
#>    area year obs pop
#> 1 Area1 2015  97 172
#> 2 Area2 2015  50 167
#> 3 Area3 2015   3 191
#> 4 Area4 2015  90 130
#> 5 Area1 2016  49 135
#> 6 Area2 2016  91 105
#> 7 Area3 2016   4 104
#> 8 Area4 2016  37 144

Execute phe_proportion and phe_rate

INPUT: The phe_proportion and phe_rate functions take a single data frame as input with columns representing the numerators and denominators for the statistic. Any other columns present will be retained in the output.

OUTPUT: The functions output the original data frame with additional columns appended. By default the additional columns are the proportion or rate, the lower 95% confidence limit, the upper 95% confidence limit, the confidence level, the statistic name and the method.

OPTIONS: The functions also accept additional arguments to specify the level of confidence, the multiplier and a reduced level of detail to be output.

Here are some example code chunks to demonstrate these two functions and the arguments that can optionally be specified

# default proportion
phe_proportion(df, obs, pop)
#>    area year obs pop      value     lowercl    uppercl confidence
#> 1 Area1 2015  97 172 0.56395349 0.489247841 0.63586486        95%
#> 2 Area2 2015  50 167 0.29940120 0.235086550 0.37273698        95%
#> 3 Area3 2015   3 191 0.01570681 0.005355862 0.04515422        95%
#> 4 Area4 2015  90 130 0.69230769 0.608401892 0.76517443        95%
#> 5 Area1 2016  49 135 0.36296296 0.286680886 0.44682811        95%
#> 6 Area2 2016  91 105 0.86666667 0.788564932 0.91888608        95%
#> 7 Area3 2016   4 104 0.03846154 0.015056813 0.09474751        95%
#> 8 Area4 2016  37 144 0.25694444 0.192543683 0.33397614        95%
#>         statistic method
#> 1 proportion of 1 Wilson
#> 2 proportion of 1 Wilson
#> 3 proportion of 1 Wilson
#> 4 proportion of 1 Wilson
#> 5 proportion of 1 Wilson
#> 6 proportion of 1 Wilson
#> 7 proportion of 1 Wilson
#> 8 proportion of 1 Wilson

# specify confidence level for proportion
phe_proportion(df, obs, pop, confidence = 99.8)
#>    area year obs pop      value     lowercl    uppercl confidence
#> 1 Area1 2015  97 172 0.56395349 0.446808065 0.67437098      99.8%
#> 2 Area2 2015  50 167 0.29940120 0.203183253 0.41731986      99.8%
#> 3 Area3 2015   3 191 0.01570681 0.003159054 0.07437558      99.8%
#> 4 Area4 2015  90 130 0.69230769 0.557697323 0.80059839      99.8%
#> 5 Area1 2016  49 135 0.36296296 0.248091378 0.49594101      99.8%
#> 6 Area2 2016  91 105 0.86666667 0.733299440 0.93889883      99.8%
#> 7 Area3 2016   4 104 0.03846154 0.009329569 0.14522444      99.8%
#> 8 Area4 2016  37 144 0.25694444 0.162049529 0.38207153      99.8%
#>         statistic method
#> 1 proportion of 1 Wilson
#> 2 proportion of 1 Wilson
#> 3 proportion of 1 Wilson
#> 4 proportion of 1 Wilson
#> 5 proportion of 1 Wilson
#> 6 proportion of 1 Wilson
#> 7 proportion of 1 Wilson
#> 8 proportion of 1 Wilson

# specify multiplier to output proportions as percentages
phe_proportion(df, obs, pop, multiplier = 100)
#>    area year obs pop     value    lowercl   uppercl confidence  statistic
#> 1 Area1 2015  97 172 56.395349 48.9247841 63.586486        95% percentage
#> 2 Area2 2015  50 167 29.940120 23.5086550 37.273698        95% percentage
#> 3 Area3 2015   3 191  1.570681  0.5355862  4.515422        95% percentage
#> 4 Area4 2015  90 130 69.230769 60.8401892 76.517443        95% percentage
#> 5 Area1 2016  49 135 36.296296 28.6680886 44.682811        95% percentage
#> 6 Area2 2016  91 105 86.666667 78.8564932 91.888608        95% percentage
#> 7 Area3 2016   4 104  3.846154  1.5056813  9.474751        95% percentage
#> 8 Area4 2016  37 144 25.694444 19.2543683 33.397614        95% percentage
#>   method
#> 1 Wilson
#> 2 Wilson
#> 3 Wilson
#> 4 Wilson
#> 5 Wilson
#> 6 Wilson
#> 7 Wilson
#> 8 Wilson

# specify multiplier for proportion, confidence level and remove metadata columns
phe_proportion(df, obs, pop, confidence = 99.8, multiplier = 100, type = "standard")
#>    area year obs pop     value    lowercl   uppercl
#> 1 Area1 2015  97 172 56.395349 44.6808065 67.437098
#> 2 Area2 2015  50 167 29.940120 20.3183253 41.731986
#> 3 Area3 2015   3 191  1.570681  0.3159054  7.437558
#> 4 Area4 2015  90 130 69.230769 55.7697323 80.059839
#> 5 Area1 2016  49 135 36.296296 24.8091378 49.594101
#> 6 Area2 2016  91 105 86.666667 73.3299440 93.889883
#> 7 Area3 2016   4 104  3.846154  0.9329569 14.522444
#> 8 Area4 2016  37 144 25.694444 16.2049529 38.207153

# default rate
phe_rate(df, obs, pop)
#>    area year obs pop     value    lowercl    uppercl confidence       statistic
#> 1 Area1 2015  97 172 56395.349 45731.7125  68798.215        95% rate per 100000
#> 2 Area2 2015  50 167 29940.120 22220.4194  39473.148        95% rate per 100000
#> 3 Area3 2015   3 191  1570.681   323.9121   4590.195        95% rate per 100000
#> 4 Area4 2015  90 130 69230.769 55668.1699  85097.279        95% rate per 100000
#> 5 Area1 2016  49 135 36296.296 26850.0523  47986.663        95% rate per 100000
#> 6 Area2 2016  91 105 86666.667 69776.6946 106408.585        95% rate per 100000
#> 7 Area3 2016   4 104  3846.154  1047.9475   9847.681        95% rate per 100000
#> 8 Area4 2016  37 144 25694.444 18088.8642  35417.428        95% rate per 100000
#>   method
#> 1  Byars
#> 2  Byars
#> 3  Exact
#> 4  Byars
#> 5  Byars
#> 6  Byars
#> 7  Exact
#> 8  Byars

# specify multiplier for rate and confidence level
phe_rate(df, obs, pop, confidence = 99.8, multiplier = 100)
#>    area year obs pop     value    lowercl   uppercl confidence    statistic
#> 1 Area1 2015  97 172 56.395349 40.3314147  76.44153      99.8% rate per 100
#> 2 Area2 2015  50 167 29.940120 18.5237775  45.49060      99.8% rate per 100
#> 3 Area3 2015   3 191  1.570681  0.0997557   6.83887      99.8% rate per 100
#> 4 Area4 2015  90 130 69.230769 48.8361086  94.89841      99.8% rate per 100
#> 5 Area1 2016  49 135 36.296296 22.3360064  55.37209      99.8% rate per 100
#> 6 Area2 2016  91 105 86.666667 61.2617114 118.59930      99.8% rate per 100
#> 7 Area3 2016   4 104  3.846154  0.4120696  14.22514      99.8% rate per 100
#> 8 Area4 2016  37 144 25.694444 14.5674857  41.63684      99.8% rate per 100
#>   method
#> 1  Byars
#> 2  Byars
#> 3  Exact
#> 4  Byars
#> 5  Byars
#> 6  Byars
#> 7  Exact
#> 8  Byars

# specify multiplier for rate, confidence level and remove metadata columns
phe_rate(df, obs, pop, type = "standard", confidence = 99.8, multiplier = 100)
#>    area year obs pop     value    lowercl   uppercl
#> 1 Area1 2015  97 172 56.395349 40.3314147  76.44153
#> 2 Area2 2015  50 167 29.940120 18.5237775  45.49060
#> 3 Area3 2015   3 191  1.570681  0.0997557   6.83887
#> 4 Area4 2015  90 130 69.230769 48.8361086  94.89841
#> 5 Area1 2016  49 135 36.296296 22.3360064  55.37209
#> 6 Area2 2016  91 105 86.666667 61.2617114 118.59930
#> 7 Area3 2016   4 104  3.846154  0.4120696  14.22514
#> 8 Area4 2016  37 144 25.694444 14.5674857  41.63684

These functions can also return aggregate data if the input dataframes are grouped:

# default proportion - grouped
df %>%
  group_by(year) %>%
  phe_proportion(obs, pop)
#> # A tibble: 2 × 9
#> # Groups:   year [2]
#>    year   obs   pop value lowercl uppercl confidence statistic       method
#>   <int> <int> <int> <dbl>   <dbl>   <dbl> <chr>      <chr>           <chr> 
#> 1  2015   240   660 0.364   0.328   0.401 95%        proportion of 1 Wilson
#> 2  2016   181   488 0.371   0.329   0.415 95%        proportion of 1 Wilson

# default rate - grouped
df %>%
  group_by(year) %>%
  phe_rate(obs, pop)
#> # A tibble: 2 × 9
#> # Groups:   year [2]
#>    year   obs   pop  value lowercl uppercl confidence statistic       method
#>   <int> <int> <int>  <dbl>   <dbl>   <dbl> <chr>      <chr>           <chr> 
#> 1  2015   240   660 36364.  31908.  41267. 95%        rate per 100000 Byars 
#> 2  2016   181   488 37090.  31883.  42905. 95%        rate per 100000 Byars



Aggregate functions

The remaining functions aggregate the rows in the input data frame to produce a single statistic. It is also possible to calculate multiple statistics in a single execution of these functions if the input data frame is grouped - for example by indicator ID, geographic area or time period (or all three). The output contains only the grouping variables and the values calculated by the function - any additional unused columns provided in the input data frame will not be retained in the output.

The df test data generated earlier can be used to demonstrate phe_mean:

Execute phe_mean

INPUT: The phe_mean function take a single data frame as input with a column representing the numbers to be averaged.

OUTPUT: By default, the function outputs one row per grouping set containing the grouping variable values (if applicable), the mean, the lower 95% confidence limit, the upper 95% confidence limit, the confidence level, the statistic name and the method.

OPTIONS: The function also accepts additional arguments to specify the level of confidence and a reduced level of detail to be output.

Here are some example code chunks to demonstrate the phe_mean function and the arguments that can optionally be specified

# default mean
phe_mean(df,obs)
#>   value_sum value_count    stdev  value  lowercl  uppercl confidence statistic
#> 1       421           8 37.66369 52.625 21.13737 84.11263        95%      mean
#>                     method
#> 1 Student's t-distribution

# multiple means in a single execution with 99.8% confidence
df %>%
    group_by(year) %>%
        phe_mean(obs, confidence = 0.998)
#> # A tibble: 2 × 10
#> # Groups:   year [2]
#>    year value_sum value_count stdev value lowercl uppercl confidence statistic
#>   <int>     <int>       <int> <dbl> <dbl>   <dbl>   <dbl> <chr>      <chr>    
#> 1  2015       240           4  43.3  60     -161.    281. 99.8%      mean     
#> 2  2016       181           4  35.9  45.2   -138.    229. 99.8%      mean     
#> # ℹ 1 more variable: method <chr>

# multiple means in a single execution with 99.8% confidence and data-only output
df %>%
    group_by(year) %>%
        phe_mean(obs, type = "standard", confidence = 0.998)
#> # A tibble: 2 × 7
#> # Groups:   year [2]
#>    year value_sum value_count stdev value lowercl uppercl
#>   <int>     <int>       <int> <dbl> <dbl>   <dbl>   <dbl>
#> 1  2015       240           4  43.3  60     -161.    281.
#> 2  2016       181           4  35.9  45.2   -138.    229.

Standardised Aggregate functions

Create some test data for the standardised aggregate functions

The following code chunk creates a data frame containing observed number of events and populations by age band for 4 areas, 5 time periods and 2 sexes:

df_std <- data.frame(
            area = rep(c("Area1", "Area2", "Area3", "Area4"), each = 19 * 2 * 5),
            year = rep(2006:2010, each = 19 * 2),
            sex = rep(rep(c("Male", "Female"), each = 19), 5),
            ageband = rep(c(0, 5,10,15,20,25,30,35,40,45,
                           50,55,60,65,70,75,80,85,90), times = 10),
            obs = sample(200, 19 * 2 * 5 * 4, replace = TRUE),
            pop = sample(10000:20000, 19 * 2 * 5 * 4, replace = TRUE))
head(df_std)
#>    area year  sex ageband obs   pop
#> 1 Area1 2006 Male       0 180 10251
#> 2 Area1 2006 Male       5 188 17675
#> 3 Area1 2006 Male      10 168 10268
#> 4 Area1 2006 Male      15  39 14612
#> 5 Area1 2006 Male      20  73 15771
#> 6 Area1 2006 Male      25 148 13248

Execute calculate_dsr

INPUT: The minimum input requirement for the calculate_dsr function is a single data frame with columns representing the numerators and denominators and standard populations for each standardisation category. The standard populations must be appended to the input data frame by the user prior to execution of the function. The 2013 European Standard Population is provided within the package in vector form (esp2013), which you can join to your dataset. Alternative standard populations can be used but must be provided by the user.

OUTPUT: By default, the function outputs one row per grouping set containing the grouping variable values, the total count, the total population, the dsr, the lower 95% confidence limit, the upper 95% confidence limit, the confidence level, the statistic name and the method.

OPTIONS: The function also accepts additional arguments to specify the level of confidence, the multiplier and a reduced level of detail to be output. It is also possible to calculate CIs when we can’t assume events are independent - further details can be found in the DSR vignette.

Here are some example code chunks to demonstrate the calculate_dsr function and the arguments that can optionally be specified


# Append the standard populations to the data frame
# calculate separate dsrs for each area, year and sex
df_std %>%
    mutate(refpop = rep(esp2013, 40)) %>%
    group_by(area, year, sex) %>%
    calculate_dsr(obs,pop, stdpop = refpop)
#> # A tibble: 40 × 11
#>    area   year sex    total_count total_pop value lowercl uppercl confidence
#>    <chr> <int> <chr>        <int>     <int> <dbl>   <dbl>   <dbl> <chr>     
#>  1 Area1  2006 Female        2017    266445  775.    739.    812. 95%       
#>  2 Area1  2006 Male          2482    291977  859.    823.    896. 95%       
#>  3 Area1  2007 Female        2110    283861  822.    785.    861. 95%       
#>  4 Area1  2007 Male          1761    303140  611.    580.    644. 95%       
#>  5 Area1  2008 Female        1607    298530  607.    576.    640. 95%       
#>  6 Area1  2008 Male          2203    278454  783.    748.    819. 95%       
#>  7 Area1  2009 Female        2085    290196  754.    720.    788. 95%       
#>  8 Area1  2009 Male          2131    285076  752.    717.    788. 95%       
#>  9 Area1  2010 Female        1573    259135  627.    595.    661. 95%       
#> 10 Area1  2010 Male          1988    271954  746.    711.    783. 95%       
#> # ℹ 30 more rows
#> # ℹ 2 more variables: statistic <chr>, method <chr>


# Append the standard populations to the data frame
# calculate separate dsrs for each area, year and sex and drop metadata fields from output
df_std %>%
    mutate(refpop = rep(esp2013, 40)) %>%
    group_by(area, year, sex) %>%
    calculate_dsr(obs,pop, stdpop = refpop, type = "standard")
#> # A tibble: 40 × 8
#>    area   year sex    total_count total_pop value lowercl uppercl
#>    <chr> <int> <chr>        <int>     <int> <dbl>   <dbl>   <dbl>
#>  1 Area1  2006 Female        2017    266445  775.    739.    812.
#>  2 Area1  2006 Male          2482    291977  859.    823.    896.
#>  3 Area1  2007 Female        2110    283861  822.    785.    861.
#>  4 Area1  2007 Male          1761    303140  611.    580.    644.
#>  5 Area1  2008 Female        1607    298530  607.    576.    640.
#>  6 Area1  2008 Male          2203    278454  783.    748.    819.
#>  7 Area1  2009 Female        2085    290196  754.    720.    788.
#>  8 Area1  2009 Male          2131    285076  752.    717.    788.
#>  9 Area1  2010 Female        1573    259135  627.    595.    661.
#> 10 Area1  2010 Male          1988    271954  746.    711.    783.
#> # ℹ 30 more rows

# calculate for under 75s by filtering out records for 75+ from input data frame and standard population
df_std %>%
  filter(ageband <= 70) %>%
  mutate(refpop = rep(esp2013[1:15], 40)) %>%
  group_by(area, year, sex) %>%
  calculate_dsr(obs, pop, stdpop = refpop)
#> # A tibble: 40 × 11
#>    area   year sex    total_count total_pop value lowercl uppercl confidence
#>    <chr> <int> <chr>        <int>     <int> <dbl>   <dbl>   <dbl> <chr>     
#>  1 Area1  2006 Female        1570    222259  755.    717.    795. 95%       
#>  2 Area1  2006 Male          1842    226309  839.    801.    879. 95%       
#>  3 Area1  2007 Female        1707    215448  838.    797.    879. 95%       
#>  4 Area1  2007 Male          1389    229806  619.    586.    654. 95%       
#>  5 Area1  2008 Female        1309    228223  613.    579.    648. 95%       
#>  6 Area1  2008 Male          1665    213212  775.    737.    814. 95%       
#>  7 Area1  2009 Female        1855    230066  789.    753.    827. 95%       
#>  8 Area1  2009 Male          1459    225044  712.    676.    751. 95%       
#>  9 Area1  2010 Female        1195    201975  590.    557.    626. 95%       
#> 10 Area1  2010 Male          1443    213375  724.    686.    762. 95%       
#> # ℹ 30 more rows
#> # ℹ 2 more variables: statistic <chr>, method <chr>

Execute calculate_ISRatio and calculate_ISRate

INPUT: These functions take a single data frame as input, with columns representing the numerators and denominators for each standardisation category, plus reference numerators and denominators for each standardisation category.

The reference data can either be provided in a separate data frame/vectors or as columns within the input data frame:

  • reference data provided as a data frame or as vectors - the data frame/vectors and the input data frame must both contain rows for the same standardisation categories, and both must be sorted, within each grouping set, by these standardisation categories in the same order.

  • reference data provided as columns within the input data frame - the reference numerators and denominators can be appended to the input data frame prior to execution of the function - if the data is grouped to generate multiple indirectly standardised rates or ratios then the reference data will need to be repeated and appended to the data rows for every grouping set.

OUTPUT: By default, the functions output one row per grouping set containing the grouping variable values, the observed and expected counts, the reference rate (ISRate only), the indirectly standardised rate or ratio, the lower 95% confidence limit, and the upper 95% confidence limit, the confidence level, the statistic name and the method.

OPTIONS: If reference data are being provided as columns within the input data frame then the user must specify this as the function expects vectors by default. The function also accepts additional arguments to specify the level of confidence, the multiplier and a reduced level of detail to be output.

The following code chunk creates a data frame containing the reference data - this example uses the all area data for persons in the baseline year:

df_ref <- df_std %>%
    filter(year == 2006) %>%
    group_by(ageband) %>%
    summarise(obs = sum(obs),
              pop = sum(pop),
              .groups = "drop_last")
    
head(df_ref)
#> # A tibble: 6 × 3
#>   ageband   obs    pop
#>     <dbl> <int>  <int>
#> 1       0  1189 112241
#> 2       5   771 123596
#> 3      10   780 112929
#> 4      15   691 135211
#> 5      20   603 123704
#> 6      25   826 115815

Here are some example code chunks to demonstrate the calculate_ISRatio function and the arguments that can optionally be specified

# calculate separate smrs for each area, year and sex
# standardised against the all-year, all-sex, all-area reference data
df_std %>%
    group_by(area, year, sex) %>%
    calculate_ISRatio(obs, pop, df_ref$obs, df_ref$pop)
#> # A tibble: 40 × 11
#> # Groups:   area, year, sex [40]
#>    area   year sex    observed expected value lowercl uppercl confidence
#>    <chr> <int> <chr>     <int>    <dbl> <dbl>   <dbl>   <dbl> <chr>     
#>  1 Area1  2006 Female     2017    1828. 1.10    1.06    1.15  95%       
#>  2 Area1  2006 Male       2482    2017. 1.23    1.18    1.28  95%       
#>  3 Area1  2007 Female     2110    1986. 1.06    1.02    1.11  95%       
#>  4 Area1  2007 Male       1761    2137. 0.824   0.786   0.863 95%       
#>  5 Area1  2008 Female     1607    2053. 0.783   0.745   0.822 95%       
#>  6 Area1  2008 Male       2203    1928. 1.14    1.10    1.19  95%       
#>  7 Area1  2009 Female     2085    2051. 1.02    0.973   1.06  95%       
#>  8 Area1  2009 Male       2131    1993. 1.07    1.02    1.12  95%       
#>  9 Area1  2010 Female     1573    1815. 0.867   0.824   0.911 95%       
#> 10 Area1  2010 Male       1988    1881. 1.06    1.01    1.10  95%       
#> # ℹ 30 more rows
#> # ℹ 2 more variables: statistic <chr>, method <chr>

# calculate the same smrs by appending the reference data to the data frame
# and drop metadata columns from output
df_std %>%
    mutate(refobs = rep(df_ref$obs,40),
           refpop = rep(df_ref$pop,40)) %>%
    group_by(area, year, sex) %>%
    calculate_ISRatio(obs, pop, refobs, refpop, refpoptype = "field",
                      type = "standard")
#> # A tibble: 40 × 8
#> # Groups:   area, year, sex [40]
#>    area   year sex    observed expected value lowercl uppercl
#>    <chr> <int> <chr>     <int>    <dbl> <dbl>   <dbl>   <dbl>
#>  1 Area1  2006 Female     2017    1828. 1.10    1.06    1.15 
#>  2 Area1  2006 Male       2482    2017. 1.23    1.18    1.28 
#>  3 Area1  2007 Female     2110    1986. 1.06    1.02    1.11 
#>  4 Area1  2007 Male       1761    2137. 0.824   0.786   0.863
#>  5 Area1  2008 Female     1607    2053. 0.783   0.745   0.822
#>  6 Area1  2008 Male       2203    1928. 1.14    1.10    1.19 
#>  7 Area1  2009 Female     2085    2051. 1.02    0.973   1.06 
#>  8 Area1  2009 Male       2131    1993. 1.07    1.02    1.12 
#>  9 Area1  2010 Female     1573    1815. 0.867   0.824   0.911
#> 10 Area1  2010 Male       1988    1881. 1.06    1.01    1.10 
#> # ℹ 30 more rows

The calculate_ISRate function works exactly the same way but instead of expressing the result as a ratio of the observed and expected rates the result is expressed as a rate and the reference rate is also provided. Here are some examples:

# calculate separate indirectly standardised rates for each area, year and sex
# standardised against the all-year, all-sex, all-area reference data
df_std %>%
    group_by(area, year, sex) %>%
    calculate_ISRate(obs, pop, df_ref$obs, df_ref$pop)
#> # A tibble: 40 × 12
#> # Groups:   area, year, sex [40]
#>    area   year sex   observed expected ref_rate value lowercl uppercl confidence
#>    <chr> <int> <chr>    <int>    <dbl>    <dbl> <dbl>   <dbl>   <dbl> <chr>     
#>  1 Area1  2006 Fema…     2017    1828.     692.  764.    731.    798. 95%       
#>  2 Area1  2006 Male      2482    2017.     692.  852.    818.    886. 95%       
#>  3 Area1  2007 Fema…     2110    1986.     692.  735.    704.    767. 95%       
#>  4 Area1  2007 Male      1761    2137.     692.  570.    544.    598. 95%       
#>  5 Area1  2008 Fema…     1607    2053.     692.  542.    515.    569. 95%       
#>  6 Area1  2008 Male      2203    1928.     692.  791.    758.    825. 95%       
#>  7 Area1  2009 Fema…     2085    2051.     692.  703.    674.    734. 95%       
#>  8 Area1  2009 Male      2131    1993.     692.  740.    709.    772. 95%       
#>  9 Area1  2010 Fema…     1573    1815.     692.  600.    570.    630. 95%       
#> 10 Area1  2010 Male      1988    1881.     692.  731.    699.    764. 95%       
#> # ℹ 30 more rows
#> # ℹ 2 more variables: statistic <chr>, method <chr>

# calculate the same indirectly standardised rates by appending the reference data to the data frame
# and drop metadata columns from output
df_std %>%
    mutate(refobs = rep(df_ref$obs,40),
           refpop = rep(df_ref$pop,40)) %>%
    group_by(area, year, sex) %>%
    calculate_ISRate(obs, pop, refobs, refpop, refpoptype = "field",
                     type = "standard")
#> # A tibble: 40 × 9
#> # Groups:   area, year, sex [40]
#>    area   year sex    observed expected ref_rate value lowercl uppercl
#>    <chr> <int> <chr>     <int>    <dbl>    <dbl> <dbl>   <dbl>   <dbl>
#>  1 Area1  2006 Female     2017    1828.     692.  764.    731.    798.
#>  2 Area1  2006 Male       2482    2017.     692.  852.    818.    886.
#>  3 Area1  2007 Female     2110    1986.     692.  735.    704.    767.
#>  4 Area1  2007 Male       1761    2137.     692.  570.    544.    598.
#>  5 Area1  2008 Female     1607    2053.     692.  542.    515.    569.
#>  6 Area1  2008 Male       2203    1928.     692.  791.    758.    825.
#>  7 Area1  2009 Female     2085    2051.     692.  703.    674.    734.
#>  8 Area1  2009 Male       2131    1993.     692.  740.    709.    772.
#>  9 Area1  2010 Female     1573    1815.     692.  600.    570.    630.
#> 10 Area1  2010 Male       1988    1881.     692.  731.    699.    764.
#> # ℹ 30 more rows