The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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")

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.