Introduction to joinpointR

Joinpoint regression is commonly used in epidemiology to identify changes in temporal trends. The joinpointR package provides tools to fit joinpoint regression models using the grid-search method, estimate summary statistics, and present results.

Joinpoint regression models by group

The function model_jp_grid() allows to fit joinpoint regression models with a log-response by levels of up to two categorical variables, using the grid-search method and model selection based in the Bayesian Information Criteria (BIC). For the examples, we are going to use a dataset with the HIV rates by sex in Argentina by jurisdiction between 2010 and 2022:

# Load required packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(joinpointR)

# Load example data
data(hiv_data)

Model selection by one grouping variable

## Create a reduced dataset
data_sex <- hiv_data |>
  filter(admin == "ARG")

## Fit the joinpoint model
mods_sex <- model_jp_grid(
  data = data_sex,
  rate = hiv_rate,
  time = year,
  group = "sex"
)
#> Model selection based on the BIC.
#> Both.sexes | Joinpoints: None detected | BIC: -3.775
#> Female | Joinpoints: 2020 | BIC: -3.686
#> Male | Joinpoints: None detected | BIC: -3.85

Model selection by two grouping variables

## Create a reduced dataset
data_admin <- hiv_data |>
  filter(between(admin, "Buenos Aires", "Chaco"))

## Fit the joinpoint model
mods_admin <- model_jp_grid(
  data = data_admin,
  rate = hiv_rate,
  time = year,
  group = c("admin", "sex")
)
#> Model selection based on the BIC.
#> Buenos.Aires_Both.sexes | Joinpoints: None detected | BIC: -3.768
#> CABA_Both.sexes | Joinpoints: None detected | BIC: -3.621
#> Catamarca_Both.sexes | Joinpoints: 2015, 2018 | BIC: -1.338
#> Chaco_Both.sexes | Joinpoints: None detected | BIC: -0.526
#> Buenos.Aires_Female | Joinpoints: None detected | BIC: -3.537
#> CABA_Female | Joinpoints: 2020 | BIC: -3.104
#> Catamarca_Female | Joinpoints: 2015, 2018 | BIC: -0.599
#> Chaco_Female | Joinpoints: None detected | BIC: -0.518
#> Buenos.Aires_Male | Joinpoints: None detected | BIC: -3.815
#> CABA_Male | Joinpoints: None detected | BIC: -3.591
#> Catamarca_Male | Joinpoints: 2015, 2018 | BIC: -0.838
#> Chaco_Male | Joinpoints: None detected | BIC: -0.431

Change model selection method

By default model_jp_grid() uses the Bayesian Information Criteria (BIC) to select the best-fit model. Model selection method can be changed using the argument method.

# Fit the model using the weighted BIC
mods_wbic <- model_jp_grid(
  data = data_admin,
  rate = hiv_rate,
  time = year,
  group = c("sex", "admin"),
  method = "wbic"
)
#> Model selection based on the WBIC.
#> Both.sexes_Buenos.Aires | Joinpoints: None detected| WBIC: -3.768
#> Female_Buenos.Aires | Joinpoints: None detected| WBIC: -3.537
#> Male_Buenos.Aires | Joinpoints: None detected| WBIC: -3.815
#> Both.sexes_CABA | Joinpoints: None detected| WBIC: -3.621
#> Female_CABA | Joinpoints: None detected| WBIC: -3.032
#> Male_CABA | Joinpoints: None detected| WBIC: -3.591
#> Both.sexes_Catamarca | Joinpoints: 2015, 2018| WBIC: -0.997
#> Female_Catamarca | Joinpoints: 2015, 2018| WBIC: -0.311
#> Male_Catamarca | Joinpoints: 2015, 2018| WBIC: -0.514
#> Both.sexes_Chaco | Joinpoints: None detected| WBIC: -0.526
#> Female_Chaco | Joinpoints: None detected| WBIC: -0.518
#> Male_Chaco | Joinpoints: None detected| WBIC: -0.431

Summarise model results

Bayesian Information Criteria

Regardless of the method selected, all the BIC metrics can be accesed using the function bic_jp():

# BIC table for all the models
bic_jp(mods_admin)
#> # A tibble: 12 × 5
#>    model                       BIC   BIC3 Weight   WBIC
#>    <chr>                     <dbl>  <dbl>  <dbl>  <dbl>
#>  1 Buenos Aires: Both sexes -3.77  -3.77   0     -3.77 
#>  2 CABA: Both sexes         -3.62  -3.62   0     -3.62 
#>  3 Catamarca: Both sexes    -1.34  -0.943  0.863 -0.997
#>  4 Chaco: Both sexes        -0.526 -0.526  0     -0.526
#>  5 Buenos Aires: Female     -3.54  -3.54   0     -3.54 
#>  6 CABA: Female             -3.10  -2.91   0.373 -3.03 
#>  7 Catamarca: Female        -0.599 -0.205  0.730 -0.311
#>  8 Chaco: Female            -0.518 -0.518  0     -0.518
#>  9 Buenos Aires: Male       -3.82  -3.82   0     -3.82 
#> 10 CABA: Male               -3.59  -3.59   0     -3.59 
#> 11 Catamarca: Male          -0.838 -0.443  0.819 -0.514
#> 12 Chaco: Male              -0.431 -0.431  0     -0.431

# BIC table for a single model
bic_jp(mods_admin[1])
#> # A tibble: 1 × 5
#>   model                      BIC  BIC3 Weight  WBIC
#>   <chr>                    <dbl> <dbl>  <dbl> <dbl>
#> 1 Buenos Aires: Both sexes -3.77 -3.77      0 -3.77

Summary tables

The functions get_summary() provides a summary table containing: - model: Model name. - jp: The number of fitted joinpoints. - period: The time segments defined by the observed joinpoints. Defaults to the whole period when no joinpoints were detected. - apc,apc_lower, apc_upper, apc_sig: The Annual Percent Change (APC) with its confidence interval and significance stars. - aapc,aapc_lower, aapc_upper, apc_sig: The Average Annual Percent Change (APC) with its confidence interval and significance stars.

get_summary(mods_sex)
#> # A tibble: 4 × 12
#>   group_var   n_jp segment period      apc apc_lower apc_upper apc_sig  aapc
#>   <chr>      <dbl>   <int> <chr>     <dbl>     <dbl>     <dbl> <chr>   <dbl>
#> 1 Both sexes     0       1 2010-2022 -4.65     -6.73     -2.52 "*"     -4.65
#> 2 Female         1       1 2010-2020 -7.35     -9.66     -4.98 "*"     -4.28
#> 3 Female         1       2 2020-2022 12.7      -3.07     30.9  ""      -4.28
#> 4 Male           0       1 2010-2022 -4.26     -6.27     -2.20 "*"     -4.26
#> # ℹ 3 more variables: aapc_lower <dbl>, aapc_upper <dbl>, aapc_sig <chr>

By default, get_summary() displays both the APC and AAPC, to show only the APC use:

get_summary(mods_sex, stats = "apc")
#> # A tibble: 4 × 8
#>   group_var   n_jp segment period      apc apc_lower apc_upper apc_sig
#>   <chr>      <dbl>   <int> <chr>     <dbl>     <dbl>     <dbl> <chr>  
#> 1 Both sexes     0       1 2010-2022 -4.65     -6.73     -2.52 "*"    
#> 2 Female         1       1 2010-2020 -7.35     -9.66     -4.98 "*"    
#> 3 Female         1       2 2020-2022 12.7      -3.07     30.9  ""     
#> 4 Male           0       1 2010-2022 -4.26     -6.27     -2.20 "*"

get_apc(mods_sex)
#> # A tibble: 4 × 8
#>   group_var   n_jp segment period      apc apc_lower apc_upper apc_sig
#>   <chr>      <dbl>   <int> <chr>     <dbl>     <dbl>     <dbl> <chr>  
#> 1 Both sexes     0       1 2010-2022 -4.65     -6.73     -2.52 "*"    
#> 2 Female         1       1 2010-2020 -7.35     -9.66     -4.98 "*"    
#> 3 Female         1       2 2020-2022 12.7      -3.07     30.9  ""     
#> 4 Male           0       1 2010-2022 -4.26     -6.27     -2.20 "*"

To display only the AAPC use:

get_summary(mods_sex, stats = "aapc")
#> # A tibble: 4 × 8
#>   group_var   n_jp segment period     aapc aapc_lower aapc_upper aapc_sig
#>   <chr>      <dbl>   <int> <chr>     <dbl>      <dbl>      <dbl> <chr>   
#> 1 Both sexes     0       1 2010-2022 -4.65      -6.73      -2.52 *       
#> 2 Female         1       1 2010-2020 -4.28      -6.50      -2.01 *       
#> 3 Female         1       2 2020-2022 -4.28      -6.50      -2.01 *       
#> 4 Male           0       1 2010-2022 -4.26      -6.27      -2.20 *

get_aapc(mods_sex)
#> # A tibble: 4 × 8
#>   group_var   n_jp segment period     aapc aapc_lower aapc_upper aapc_sig
#>   <chr>      <dbl>   <int> <chr>     <dbl>      <dbl>      <dbl> <chr>   
#> 1 Both sexes     0       1 2010-2022 -4.65      -6.73      -2.52 *       
#> 2 Female         1       1 2010-2020 -4.28      -6.50      -2.01 *       
#> 3 Female         1       2 2020-2022 -4.28      -6.50      -2.01 *       
#> 4 Male           0       1 2010-2022 -4.26      -6.27      -2.20 *

The confidence interval or significance stars can be hided using the argument hide:

#|id: get_summary-4
# Hide the confidence interval
get_summary(mods_sex, hide = "ci")
#> # A tibble: 4 × 8
#>   group_var   n_jp segment period      apc apc_sig  aapc aapc_sig
#>   <chr>      <dbl>   <int> <chr>     <dbl> <chr>   <dbl> <chr>   
#> 1 Both sexes     0       1 2010-2022 -4.65 "*"     -4.65 *       
#> 2 Female         1       1 2010-2020 -7.35 "*"     -4.28 *       
#> 3 Female         1       2 2020-2022 12.7  ""      -4.28 *       
#> 4 Male           0       1 2010-2022 -4.26 "*"     -4.26 *

# Hide the significance stars
get_summary(mods_sex, hide = "sig")
#> # A tibble: 4 × 10
#>   group_var   n_jp segment period      apc apc_lower apc_upper  aapc aapc_lower
#>   <chr>      <dbl>   <int> <chr>     <dbl>     <dbl>     <dbl> <dbl>      <dbl>
#> 1 Both sexes     0       1 2010-2022 -4.65     -6.73     -2.52 -4.65      -6.73
#> 2 Female         1       1 2010-2020 -7.35     -9.66     -4.98 -4.28      -6.50
#> 3 Female         1       2 2020-2022 12.7      -3.07     30.9  -4.28      -6.50
#> 4 Male           0       1 2010-2022 -4.26     -6.27     -2.20 -4.26      -6.27
#> # ℹ 1 more variable: aapc_upper <dbl>

The confidence level defaults to 95% and can be changed using the argument level.ci:

get_summary(mods_sex, level.ci = .9)
#> # A tibble: 4 × 12
#>   group_var   n_jp segment period      apc apc_lower apc_upper apc_sig  aapc
#>   <chr>      <dbl>   <int> <chr>     <dbl>     <dbl>     <dbl> <chr>   <dbl>
#> 1 Both sexes     0       1 2010-2022 -4.65    -6.35      -2.92 "*"     -4.65
#> 2 Female         1       1 2010-2020 -7.35    -9.24      -5.42 "*"     -4.28
#> 3 Female         1       2 2020-2022 12.7     -0.312     27.3  ""      -4.28
#> 4 Male           0       1 2010-2022 -4.26    -5.90      -2.59 "*"     -4.26
#> # ℹ 3 more variables: aapc_lower <dbl>, aapc_upper <dbl>, aapc_sig <chr>

Results can also be displayed as a flextable:

get_summary(mods_sex, as.ft = TRUE)

group

n_jp

segment

period

apc

apc_lower

apc_upper

apc_sig

aapc

aapc_lower

aapc_upper

aapc_sig

Both sexes

0

1

2010-2022

-4.65

-6.73

-2.52

*

-4.65

-6.73

-2.52

*

Female

1

1

2010-2020

-7.35

-9.66

-4.98

*

-4.28

-6.50

-2.01

*

2

2020-2022

12.66

-3.07

30.95

Male

0

1

2010-2022

-4.26

-6.27

-2.20

*

-4.26

-6.27

-2.20

*

Regression plots

Default plot

One grouping variable

gg_jpoint(mods_sex)

Two grouping variables

gg_jpoint(mods_admin)

Change the facets layout

gg_jpoint(mods_admin, facets = "grid")

Reverse the facet grid

gg_jpoint(mods_admin, facets = "grid2")

Change the colors by APC trend

gg_jpoint(mods_admin, facets = "grid2", color.by = "trend")

Change the colors by time segment

gg_jpoint(mods_admin, facets = "grid2", color.by = "segment")

Change the colors by time period

gg_jpoint(mods_admin, facets = "grid2", color.by = "period")

Show Average Annual Percent Change (AAPC)

gg_jpoint(mods_admin, facets = "grid2", aapc = TRUE)

Hide joinpoint positions

gg_jpoint(mods_admin, facets = "grid2", jp = FALSE)
#> Position of joinpoints will not be displayed.

Plot fitted lines

The observed log-rates can be hidden using geom = "line"

gg_jpoint(mods_admin, facets = "grid2", geom = "line")

Alternatively, the same results can be achieved using gg_jpoint_line()

gg_jpoint_line(mods_admin, facets = "grid2")

Plot fitted lines with background

gg_jpoint(mods_sex, facets = "grid2", geom = "area")

The shortcut gg_jpoint_area()produces the same result

gg_jpoint_area(mods_sex, facets = "grid2")

Change the default colors

By default, gg_jpoint() uses the default colorblind-friendly palette "viridis" for coloring the plots. Colorblind-friendly palettes can be changed using the functions scale_cbpal_color() and scale_cbpal_fill().

gg_jpoint_line(mods_sex) +
  scale_cbpal_color(palette = "algae")
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.


gg_jpoint_area(mods_sex, color.by = "trend") +
  scale_cbpal_fill(palette = "blue_fluoride")
#> Scale for fill is already present.
#> Adding another scale for fill, which will replace the existing scale.

The available colorblind-friendly palettes can be checked using plot_cbpal():

plot_cbpal(type = "div")


plot_cbpal(series = "scico")