Introduction to PHEindicatormethods

Georgina Anderson

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  23 135
#> 2 Area2 2015 100 132
#> 3 Area3 2015  43 170
#> 4 Area4 2015  35 177
#> 5 Area1 2016  34 196
#> 6 Area2 2016  35 137
#> 7 Area3 2016   1 141
#> 8 Area4 2016 100 188

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  23 135 0.170370370 0.116293355 0.24268774        95%
#> 2 Area2 2015 100 132 0.757575758 0.677858179 0.82272537        95%
#> 3 Area3 2015  43 170 0.252941176 0.193551591 0.32324952        95%
#> 4 Area4 2015  35 177 0.197740113 0.145756520 0.26256500        95%
#> 5 Area1 2016  34 196 0.173469388 0.126873824 0.23261844        95%
#> 6 Area2 2016  35 137 0.255474453 0.189808603 0.33447920        95%
#> 7 Area3 2016   1 141 0.007092199 0.001253046 0.03907697        95%
#> 8 Area4 2016 100 188 0.531914894 0.460663888 0.60188776        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  23 135 0.170370370 0.0930914703 0.29120264      99.8%
#> 2 Area2 2015 100 132 0.757575758 0.6275393887 0.85285781      99.8%
#> 3 Area3 2015  43 170 0.252941176 0.1649733715 0.36718916      99.8%
#> 4 Area4 2015  35 177 0.197740113 0.1217788226 0.30464698      99.8%
#> 5 Area1 2016  34 196 0.173469388 0.1056259311 0.27165313      99.8%
#> 6 Area2 2016  35 137 0.255474453 0.1589438166 0.38387289      99.8%
#> 7 Area3 2016   1 141 0.007092199 0.0006190634 0.07609679      99.8%
#> 8 Area4 2016 100 188 0.531914894 0.4206537088 0.64009055      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  23 135 17.0370370 11.6293355 24.268774        95% percentage
#> 2 Area2 2015 100 132 75.7575758 67.7858179 82.272537        95% percentage
#> 3 Area3 2015  43 170 25.2941176 19.3551591 32.324952        95% percentage
#> 4 Area4 2015  35 177 19.7740113 14.5756520 26.256500        95% percentage
#> 5 Area1 2016  34 196 17.3469388 12.6873824 23.261844        95% percentage
#> 6 Area2 2016  35 137 25.5474453 18.9808603 33.447920        95% percentage
#> 7 Area3 2016   1 141  0.7092199  0.1253046  3.907697        95% percentage
#> 8 Area4 2016 100 188 53.1914894 46.0663888 60.188776        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  23 135 17.0370370  9.30914703 29.120264
#> 2 Area2 2015 100 132 75.7575758 62.75393887 85.285781
#> 3 Area3 2015  43 170 25.2941176 16.49733715 36.718916
#> 4 Area4 2015  35 177 19.7740113 12.17788226 30.464698
#> 5 Area1 2016  34 196 17.3469388 10.56259311 27.165313
#> 6 Area2 2016  35 137 25.5474453 15.89438166 38.387289
#> 7 Area3 2016   1 141  0.7092199  0.06190634  7.609679
#> 8 Area4 2016 100 188 53.1914894 42.06537088 64.009055

# default rate
phe_rate(df, obs, pop)
#>    area year obs pop      value     lowercl  uppercl confidence       statistic
#> 1 Area1 2015  23 135 17037.0370 10796.57087 25565.04        95% rate per 100000
#> 2 Area2 2015 100 132 75757.5758 61637.95870 92142.38        95% rate per 100000
#> 3 Area3 2015  43 170 25294.1176 18303.62989 34071.85        95% rate per 100000
#> 4 Area4 2015  35 177 19774.0113 13771.30535 27501.70        95% rate per 100000
#> 5 Area1 2016  34 196 17346.9388 12011.40359 24241.37        95% rate per 100000
#> 6 Area2 2016  35 137 25547.4453 17792.12443 35531.39        95% rate per 100000
#> 7 Area3 2016   1 141   709.2199    17.95589  3951.52        95% rate per 100000
#> 8 Area4 2016 100 188 53191.4894 43277.71569 64695.71        95% rate per 100000
#>   method
#> 1  Byars
#> 2  Byars
#> 3  Byars
#> 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  23 135 17.0370370 8.096878e+00  31.15473      99.8% rate per 100
#> 2 Area2 2015 100 132 75.7575758 5.447249e+01 102.23002      99.8% rate per 100
#> 3 Area3 2015  43 170 25.2941176 1.500964e+01  39.64842      99.8% rate per 100
#> 4 Area4 2015  35 177 19.7740113 1.101123e+01  32.45781      99.8% rate per 100
#> 5 Area1 2016  34 196 17.3469388 9.567179e+00  28.66924      99.8% rate per 100
#> 6 Area2 2016  35 137 25.5474453 1.422619e+01  41.93455      99.8% rate per 100
#> 7 Area3 2016   1 141  0.7092199 7.095747e-04   6.54852      99.8% rate per 100
#> 8 Area4 2016 100 188 53.1914894 3.824664e+01  71.77852      99.8% rate per 100
#>   method
#> 1  Byars
#> 2  Byars
#> 3  Byars
#> 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  23 135 17.0370370 8.096878e+00  31.15473
#> 2 Area2 2015 100 132 75.7575758 5.447249e+01 102.23002
#> 3 Area3 2015  43 170 25.2941176 1.500964e+01  39.64842
#> 4 Area4 2015  35 177 19.7740113 1.101123e+01  32.45781
#> 5 Area1 2016  34 196 17.3469388 9.567179e+00  28.66924
#> 6 Area2 2016  35 137 25.5474453 1.422619e+01  41.93455
#> 7 Area3 2016   1 141  0.7092199 7.095747e-04   6.54852
#> 8 Area4 2016 100 188 53.1914894 3.824664e+01  71.77852

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   201   614 0.327   0.291   0.365 95%        proportion of 1 Wilson
#> 2  2016   170   662 0.257   0.225   0.291 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   201   614 32736.  28367.  37588. 95%        rate per 100000 Byars 
#> 2  2016   170   662 25680.  21964.  29844. 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       371           8 35.41564 46.375 16.76678 75.98322        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       201           4  34.2  50.2   -124.    225. 99.8%      mean     
#> 2  2016       170           4  41.5  42.5   -169.    254. 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       201           4  34.2  50.2   -124.    225.
#> 2  2016       170           4  41.5  42.5   -169.    254.

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  43 12994
#> 2 Area1 2006 Male       5 118 16129
#> 3 Area1 2006 Male      10  78 14863
#> 4 Area1 2006 Male      15 117 10399
#> 5 Area1 2006 Male      20  75 10560
#> 6 Area1 2006 Male      25  62 11231

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        1932    282476  646.    615.    678. 95%       
#>  2 Area1  2006 Male          1665    269331  677.    644.    712. 95%       
#>  3 Area1  2007 Female        1868    270858  683.    650.    716. 95%       
#>  4 Area1  2007 Male          1875    284623  714.    678.    752. 95%       
#>  5 Area1  2008 Female        1772    273428  665.    632.    699. 95%       
#>  6 Area1  2008 Male          2064    306517  671.    640.    702. 95%       
#>  7 Area1  2009 Female        2151    259869  803.    767.    841. 95%       
#>  8 Area1  2009 Male          1882    276191  751.    716.    788. 95%       
#>  9 Area1  2010 Female        1975    290984  695.    662.    729. 95%       
#> 10 Area1  2010 Male          2047    280129  828.    791.    866. 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        1932    282476  646.    615.    678.
#>  2 Area1  2006 Male          1665    269331  677.    644.    712.
#>  3 Area1  2007 Female        1868    270858  683.    650.    716.
#>  4 Area1  2007 Male          1875    284623  714.    678.    752.
#>  5 Area1  2008 Female        1772    273428  665.    632.    699.
#>  6 Area1  2008 Male          2064    306517  671.    640.    702.
#>  7 Area1  2009 Female        2151    259869  803.    767.    841.
#>  8 Area1  2009 Male          1882    276191  751.    716.    788.
#>  9 Area1  2010 Female        1975    290984  695.    662.    729.
#> 10 Area1  2010 Male          2047    280129  828.    791.    866.
#> # ℹ 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        1389    228672  609.    576.    643. 95%       
#>  2 Area1  2006 Male          1510    212452  716.    680.    754. 95%       
#>  3 Area1  2007 Female        1462    208576  682.    647.    719. 95%       
#>  4 Area1  2007 Male          1328    216793  716.    677.    756. 95%       
#>  5 Area1  2008 Female        1442    217785  676.    640.    712. 95%       
#>  6 Area1  2008 Male          1534    246925  650.    617.    683. 95%       
#>  7 Area1  2009 Female        1547    212625  748.    711.    787. 95%       
#>  8 Area1  2009 Male          1700    214864  801.    762.    841. 95%       
#>  9 Area1  2010 Female        1492    227411  680.    645.    717. 95%       
#> 10 Area1  2010 Male          1672    225088  823.    784.    865. 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:

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   738 117578
#> 2       5  1116 111392
#> 3      10   644 112054
#> 4      15   980 117658
#> 5      20   889 125921
#> 6      25   829 105114

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     1932    1913. 1.01    0.965   1.06  95%       
#>  2 Area1  2006 Male       1665    1804. 0.923   0.879   0.968 95%       
#>  3 Area1  2007 Female     1868    1819. 1.03    0.981   1.07  95%       
#>  4 Area1  2007 Male       1875    1911. 0.981   0.937   1.03  95%       
#>  5 Area1  2008 Female     1772    1870. 0.948   0.904   0.993 95%       
#>  6 Area1  2008 Male       2064    2106. 0.980   0.938   1.02  95%       
#>  7 Area1  2009 Female     2151    1780. 1.21    1.16    1.26  95%       
#>  8 Area1  2009 Male       1882    1870. 1.01    0.961   1.05  95%       
#>  9 Area1  2010 Female     1975    2008. 0.983   0.940   1.03  95%       
#> 10 Area1  2010 Male       2047    1883. 1.09    1.04    1.14  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     1932    1913. 1.01    0.965   1.06 
#>  2 Area1  2006 Male       1665    1804. 0.923   0.879   0.968
#>  3 Area1  2007 Female     1868    1819. 1.03    0.981   1.07 
#>  4 Area1  2007 Male       1875    1911. 0.981   0.937   1.03 
#>  5 Area1  2008 Female     1772    1870. 0.948   0.904   0.993
#>  6 Area1  2008 Male       2064    2106. 0.980   0.938   1.02 
#>  7 Area1  2009 Female     2151    1780. 1.21    1.16    1.26 
#>  8 Area1  2009 Male       1882    1870. 1.01    0.961   1.05 
#>  9 Area1  2010 Female     1975    2008. 0.983   0.940   1.03 
#> 10 Area1  2010 Male       2047    1883. 1.09    1.04    1.14 
#> # ℹ 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…     1932    1913.     673.  679.    650.    710. 95%       
#>  2 Area1  2006 Male      1665    1804.     673.  621.    592.    652. 95%       
#>  3 Area1  2007 Fema…     1868    1819.     673.  691.    660.    723. 95%       
#>  4 Area1  2007 Male      1875    1911.     673.  660.    631.    691. 95%       
#>  5 Area1  2008 Fema…     1772    1870.     673.  638.    608.    668. 95%       
#>  6 Area1  2008 Male      2064    2106.     673.  660.    631.    689. 95%       
#>  7 Area1  2009 Fema…     2151    1780.     673.  813.    779.    848. 95%       
#>  8 Area1  2009 Male      1882    1870.     673.  677.    647.    708. 95%       
#>  9 Area1  2010 Fema…     1975    2008.     673.  662.    633.    691. 95%       
#> 10 Area1  2010 Male      2047    1883.     673.  731.    700.    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     1932    1913.     673.  679.    650.    710.
#>  2 Area1  2006 Male       1665    1804.     673.  621.    592.    652.
#>  3 Area1  2007 Female     1868    1819.     673.  691.    660.    723.
#>  4 Area1  2007 Male       1875    1911.     673.  660.    631.    691.
#>  5 Area1  2008 Female     1772    1870.     673.  638.    608.    668.
#>  6 Area1  2008 Male       2064    2106.     673.  660.    631.    689.
#>  7 Area1  2009 Female     2151    1780.     673.  813.    779.    848.
#>  8 Area1  2009 Male       1882    1870.     673.  677.    647.    708.
#>  9 Area1  2010 Female     1975    2008.     673.  662.    633.    691.
#> 10 Area1  2010 Male       2047    1883.     673.  731.    700.    764.
#> # ℹ 30 more rows