Exploring the relationship between Air pollution and incident arrhythmia: a UKBAnalytica demo case study

This document demonstrates an end-to-end UKBAnalytica workflow using ukb_demo() to generate a RAP-style toy table at runtime. The case study evaluates associations between air pollution exposure and incident cardiac arrhythmia. The workflow covers raw field extraction, phenotype definition, baseline preprocessing, participant flow, baseline table generation, Cox regression, and sensitivity analyses.

The demo data are fully synthetic, generated from toy random distributions, and do not contain UK Biobank participant-level source records. Results are for workflow illustration only.

1 Load demo data

if (requireNamespace("pkgload", quietly = TRUE) && dir.exists("../R")) {
  pkgload::load_all("..", quiet = TRUE)
}
library(UKBAnalytica)
library(data.table)

demo_raw <- ukb_demo(seed = 20260618)
dim(demo_raw)
#> [1] 500 754
names(demo_raw)[1:12]
#>  [1] "eid"       "p31"       "p53_i0"    "p53_i1"    "p53_i2"    "p53_i3"   
#>  [7] "p21022"    "p21001_i0" "p21001_i1" "p21001_i2" "p21001_i3" "p21000_i0"

2 Specify variables from raw UKB field IDs

The analysis uses the following raw UKB fields:

  • NO2: p24016, p24018, and p24017; p24003 is used for the 2010-exposure sensitivity analysis.
  • NOx: p24004.
  • PM2.5: p24006.
  • PM10: p24019; p24005 is used for the 2010-exposure sensitivity analysis.
  • Covariates: sex (p31), baseline date (p53_i0), age (p21022), BMI (p21001_i0), smoking (p20116_i0), drinking (p20117_i0), Townsend deprivation index (p189), ethnicity (p21000_i0), and education (p6138_i0).

3 Define incident arrhythmia

The predefined Arrhythmia phenotype in UKBAnalytica captures broad cardiac arrhythmia using ICD-10 codes I44-I49 and procedure codes including K57.6, K59-K62, K64.1, and K72-K74 where available. Baseline history and prospective incident events are handled separately.

disease_catalog <- get_predefined_diseases()
arrhythmia_def <- disease_catalog["Arrhythmia"]
arrhythmia_def
#> $Arrhythmia
#> $Arrhythmia$name
#> [1] "Cardiac Arrhythmia"
#> 
#> $Arrhythmia$icd10_pattern
#> [1] "^(I44|I45|I46|I47|I48|I49)"
#> 
#> $Arrhythmia$icd9_pattern
#> NULL
#> 
#> $Arrhythmia$sr_codes
#> NULL
#> 
#> $Arrhythmia$death_icd10
#> [1] "^(I44|I45|I46|I47|I48|I49)"
#> 
#> $Arrhythmia$opcs4_pattern
#> [1] "^(K576|K59|K60|K61|K62|K641|K72|K73|K74)"
#> 
#> $Arrhythmia$first_occurrence_fields
#> [1] 131342 131344 131346 131348 131350 131352
#> 
#> $Arrhythmia$first_occurrence_source_fields
#> NULL
#> 
#> $Arrhythmia$cancer_icd10_pattern
#> NULL
#> 
#> $Arrhythmia$cancer_histology
#> NULL
#> 
#> $Arrhythmia$cancer_behaviour
#> NULL
#> 
#> $Arrhythmia$algo_date_field
#> NULL
#> 
#> $Arrhythmia$algo_source_field
#> NULL

demo_surv <- build_survival_dataset(
  dt = demo_raw,
  disease_definitions = arrhythmia_def,
  prevalent_sources = c("ICD10", "OPCS4"),
  outcome_sources = c("ICD10"),
  baseline_col = "p53_i0",
  primary_disease = "Arrhythmia",
  censor_date = as.Date("2023-10-31"),
  show_flow = TRUE
)
#>                                     step n_before n_after
#>                                   <char>    <int>   <int>
#> 1:                            Raw cohort      500     500
#> 2:          After build_survival_dataset      500     500
#> 3:       Keep non-missing outcome_status      500     478
#> 4: Exclude baseline prevalent Arrhythmia      478     478
#> 5:          Keep valid outcome_surv_time      478     464
#>                                             exclusion_rule excluded
#>                                                     <char>    <int>
#> 1:                                                    None        0
#> 2:                         Function output row count check        0
#> 3:                            Exclude outcome_status is NA       22
#> 4:                            Keep Arrhythmia_history == 0        0
#> 5: Keep !is.na(outcome_surv_time) & outcome_surv_time >= 0       14
#>    retained_from_prev retained_from_raw
#>                <char>            <char>
#> 1:            100.00%           100.00%
#> 2:            100.00%           100.00%
#> 3:             95.60%            95.60%
#> 4:            100.00%            95.60%
#> 5:             97.07%            92.80%

table(demo_surv$outcome_status, useNA = "ifany")
#> 
#>    0    1 <NA> 
#>  441   37   22
participant_flow <- attr(demo_surv, "participant_flow")
participant_flow
#>                                     step n_before n_after
#>                                   <char>    <int>   <int>
#> 1:                            Raw cohort      500     500
#> 2:          After build_survival_dataset      500     500
#> 3:       Keep non-missing outcome_status      500     478
#> 4: Exclude baseline prevalent Arrhythmia      478     478
#> 5:          Keep valid outcome_surv_time      478     464
#>                                             exclusion_rule excluded
#>                                                     <char>    <int>
#> 1:                                                    None        0
#> 2:                         Function output row count check        0
#> 3:                            Exclude outcome_status is NA       22
#> 4:                            Keep Arrhythmia_history == 0        0
#> 5: Keep !is.na(outcome_surv_time) & outcome_surv_time >= 0       14
#>    retained_from_prev retained_from_raw
#>                 <num>             <num>
#> 1:          1.0000000             1.000
#> 2:          1.0000000             1.000
#> 3:          0.9560000             0.956
#> 4:          1.0000000             0.956
#> 5:          0.9707113             0.928

outcome_status = 1 indicates incident arrhythmia, outcome_status = 0 indicates censoring, and NA indicates baseline prevalent arrhythmia.

4 Preprocess exposures and covariates

preprocess_baseline() maps raw UKB field names to standardized analysis variables. Townsend deprivation index is added through custom_mapping because it is not part of the default baseline variable set.

exposure_vars <- c(
  "no2_2005", "no2_2006", "no2_2007", "no2_2010",
  "nox", "pm25", "pm10_2007", "pm10_2010"
)

covariates <- c("age", "sex", "tdi", "ethnicity", "education", "smoking", "drinking")

custom_mapping <- list(
  tdi = list(
    ukb_col = "p189",
    description = "Townsend deprivation index at recruitment"
  )
)

demo_prepared <- preprocess_baseline(
  df = demo_surv,
  variables = c(exposure_vars, covariates),
  custom_mapping = custom_mapping,
  missing_action = "keep"
)

demo_prepared[, c(exposure_vars, covariates), with = FALSE][1:6]
#>    no2_2005 no2_2006 no2_2007 no2_2010   nox  pm25 pm10_2007 pm10_2010   age
#>       <num>    <num>    <num>    <num> <num> <num>     <num>     <num> <num>
#> 1:    12.81    13.20    14.16    15.47 41.40  9.39     14.16     13.52  67.6
#> 2:    31.55    33.77    34.00    31.29 52.86 10.10     16.64     16.66  54.4
#> 3:    24.51    21.24    22.58    24.58 37.58  9.37     25.85     24.10  47.0
#> 4:    26.26    26.13    23.33    22.41 55.11 10.64     17.47     19.95  64.5
#> 5:    24.87    24.60    26.70    23.37 35.48  8.34     16.43     17.08  71.0
#> 6:    27.45    30.21    27.86    29.32 53.44  6.98      7.68      7.76  45.4
#>       sex   tdi ethnicity education  smoking drinking
#>    <fctr> <num>    <fctr>    <fctr>   <fctr>   <fctr>
#> 1: Female  1.64     White    Medium    Never Previous
#> 2:   Male -0.78     White      High    Never  Current
#> 3:   Male  4.36     White    Medium    Never Previous
#> 4:   Male  2.13     White      High    Never  Current
#> 5:   Male  2.50     White      High    Never Previous
#> 6: Female -1.87     White    Medium Previous Previous

NO2 is summarized as the mean of 2005, 2006, and 2007 estimates for the main analysis. The main PM10 exposure uses 2007 estimates. The 2010 NO2 and PM10 fields are retained for sensitivity analysis.

demo_prepared[, no2_main := rowMeans(
  cbind(no2_2005, no2_2006, no2_2007),
  na.rm = TRUE
)]
demo_prepared[is.nan(no2_main), no2_main := NA_real_]

main_exposures <- c("no2_main", "nox", "pm25", "pm10_2007")
sensitivity_exposures_2010 <- c("no2_2010", "nox", "pm25", "pm10_2010")

main_exposures_z <- paste0(main_exposures, "_z")
sensitivity_exposures_2010_z <- paste0(sensitivity_exposures_2010, "_z")

for (var in unique(c(main_exposures, sensitivity_exposures_2010))) {
  demo_prepared[, (paste0(var, "_z")) := as.numeric(scale(get(var)))]
}

demo_prepared[, c(main_exposures, main_exposures_z), with = FALSE][1:6]
#>    no2_main   nox  pm25 pm10_2007  no2_main_z       nox_z     pm25_z
#>       <num> <num> <num>     <num>       <num>       <num>      <num>
#> 1: 13.39000 41.40  9.39     14.16 -1.98107965 -0.04713709 -0.2325075
#> 2: 33.10667 52.86 10.10     16.64  1.33380183  0.96928388  0.1078715
#> 3: 22.77667 37.58  9.37     25.85 -0.40293828 -0.38594408 -0.2420956
#> 4: 25.24000 55.11 10.64     17.47  0.01121175  1.16884297  0.3667513
#> 5: 25.39000 35.48  8.34     16.43  0.03643063 -0.57219923 -0.7358849
#> 6: 28.50667 53.44  6.98      7.68  0.56042289  1.02072577 -1.3878785
#>    pm10_2007_z
#>          <num>
#> 1:  -0.9767831
#> 2:  -0.4112271
#> 3:   1.6890837
#> 4:  -0.2219483
#> 5:  -0.4591169
#> 6:  -2.4545262

5 Assemble the analytic cohort

The prospective cohort excludes participants with baseline prevalent arrhythmia, missing outcome status, invalid follow-up time, or missing values in the main exposures and adjustment covariates. The stepwise complete-case flow is retained for reporting.

analysis_data <- demo_prepared[
  !is.na(outcome_status) &
    !is.na(outcome_surv_time) &
    outcome_surv_time > 0
]

complete_case_vars <- c(main_exposures_z, covariates)

analysis_data <- sensitivity_exclude_missing_covariates(
  data = analysis_data,
  covariates = complete_case_vars,
  stepwise = TRUE,
  verbose = TRUE
)
#>      step    variable n_before n_missing_in_remaining n_after n_removed
#>     <int>      <char>    <int>                  <int>   <int>     <int>
#>  1:     1  no2_main_z      464                      0     464         0
#>  2:     2       nox_z      464                     15     449        15
#>  3:     3      pm25_z      449                     13     436        13
#>  4:     4 pm10_2007_z      436                     12     424        12
#>  5:     5         age      424                      0     424         0
#>  6:     6         sex      424                      0     424         0
#>  7:     7         tdi      424                     17     407        17
#>  8:     8   ethnicity      407                     13     394        13
#>  9:     9   education      394                     23     371        23
#> 10:    10     smoking      371                     20     351        20
#> 11:    11    drinking      351                     15     336        15
#>     retained_from_previous retained_from_input
#>                      <num>               <num>
#>  1:              1.0000000           1.0000000
#>  2:              0.9676724           0.9676724
#>  3:              0.9710468           0.9396552
#>  4:              0.9724771           0.9137931
#>  5:              1.0000000           0.9137931
#>  6:              1.0000000           0.9137931
#>  7:              0.9599057           0.8771552
#>  8:              0.9680590           0.8491379
#>  9:              0.9416244           0.7995690
#> 10:              0.9460916           0.7564655
#> 11:              0.9572650           0.7241379

complete_case_flow <- attr(analysis_data, "complete_case_flow")
complete_case_flow
#>      step    variable n_before n_missing_in_remaining n_after n_removed
#>     <int>      <char>    <int>                  <int>   <int>     <int>
#>  1:     1  no2_main_z      464                      0     464         0
#>  2:     2       nox_z      464                     15     449        15
#>  3:     3      pm25_z      449                     13     436        13
#>  4:     4 pm10_2007_z      436                     12     424        12
#>  5:     5         age      424                      0     424         0
#>  6:     6         sex      424                      0     424         0
#>  7:     7         tdi      424                     17     407        17
#>  8:     8   ethnicity      407                     13     394        13
#>  9:     9   education      394                     23     371        23
#> 10:    10     smoking      371                     20     351        20
#> 11:    11    drinking      351                     15     336        15
#>     retained_from_previous retained_from_input
#>                      <num>               <num>
#>  1:              1.0000000           1.0000000
#>  2:              0.9676724           0.9676724
#>  3:              0.9710468           0.9396552
#>  4:              0.9724771           0.9137931
#>  5:              1.0000000           0.9137931
#>  6:              1.0000000           0.9137931
#>  7:              0.9599057           0.8771552
#>  8:              0.9680590           0.8491379
#>  9:              0.9416244           0.7995690
#> 10:              0.9460916           0.7564655
#> 11:              0.9572650           0.7241379

table(analysis_data$outcome_status)
#> 
#>   0   1 
#> 309  27

6 Baseline characteristics

Baseline characteristics are summarized by incident arrhythmia status. Air pollutants are shown on their original scales, while Cox models use standardized exposure variables.

table1 <- create_baseline_table(
  data = analysis_data,
  case_col = "outcome_status",
  factor_cols = c("sex", "ethnicity", "education", "smoking", "drinking"),
  continuous_cols = c("age", "tdi", main_exposures),
  test = TRUE
)

print(table1, quote = FALSE, noSpaces = TRUE)
#>                         Stratified by outcome_status
#>                          0             1             p     test
#>   n                      309           27                      
#>   sex = Male (%)         165 (53.4)    10 (37.0)     0.152     
#>   ethnicity = Others (%) 43 (13.9)     2 (7.4)       0.511     
#>   education = High (%)   147 (47.6)    17 (63.0)     0.182     
#>   smoking (%)                                        0.183     
#>      Never               173 (56.0)    12 (44.4)               
#>      Previous            95 (30.7)     8 (29.6)                
#>      Current             41 (13.3)     7 (25.9)                
#>   drinking (%)                                       0.629     
#>      Never               28 (9.1)      1 (3.7)                 
#>      Previous            209 (67.6)    19 (70.4)               
#>      Current             72 (23.3)     7 (25.9)                
#>   age (mean (SD))        57.33 (8.56)  57.50 (10.29) 0.925     
#>   tdi (mean (SD))        0.09 (3.07)   0.85 (3.51)   0.221     
#>   no2_main (mean (SD))   25.00 (6.07)  27.78 (5.65)  0.023     
#>   nox (mean (SD))        42.13 (10.91) 35.87 (13.70) 0.005     
#>   pm25 (mean (SD))       9.93 (2.09)   9.96 (2.48)   0.948     
#>   pm10_2007 (mean (SD))  18.59 (4.19)  18.03 (4.93)  0.518

7 Cox regression

Each pollutant is evaluated in a separate Cox proportional hazards model, adjusted for age, sex, Townsend deprivation index, ethnicity, education, smoking, and drinking. Hazard ratios are interpreted per one-standard-deviation increase in exposure.

pollutant_labels <- data.frame(
  variable = main_exposures_z,
  pollutant = c("NO2", "NOx", "PM2.5", "PM10"),
  stringsAsFactors = FALSE
)

cox_results <- runmulti_cox(
  data = analysis_data,
  main_var = main_exposures_z,
  covariates = covariates,
  endpoint = c("outcome_surv_time", "outcome_status")
)

cox_results <- merge(cox_results, pollutant_labels, by = "variable", all.x = TRUE)
cox_results$p_bh <- p.adjust(cox_results$pvalue, method = "BH")
cox_results$p_bonferroni <- p.adjust(cox_results$pvalue, method = "bonferroni")
cox_results[order(cox_results$pvalue), ]
#>      variable    contrast         coef        se           z        HR
#> 2       nox_z       nox_z -0.571962396 0.1946940 -2.93775115 0.5644167
#> 1  no2_main_z  no2_main_z  0.409119575 0.2001311  2.04425836 1.5054917
#> 3 pm10_2007_z pm10_2007_z -0.105484213 0.1969988 -0.53545613 0.8998887
#> 4      pm25_z      pm25_z -0.005581028 0.1906796 -0.02926914 0.9944345
#>     lower95   upper95      pvalue   n n_event pollutant       p_bh p_bonferroni
#> 2 0.3853684 0.8266538 0.003306023 336      27       NOx 0.01322409   0.01322409
#> 1 1.0170130 2.2285903 0.040928030 336      27       NO2 0.08185606   0.16371212
#> 3 0.6116502 1.3239588 0.592334486 336      27      PM10 0.78977931   1.00000000
#> 4 0.6843361 1.4450502 0.976649936 336      27     PM2.5 0.97664994   1.00000000

8 Sensitivity analysis 1: exclude early events

To reduce potential reverse causation, early incident arrhythmia events are excluded at 1-, 2-, and 4-year lag windows, and the adjusted Cox models are refitted.

fit_pollution_cox <- function(data, exposures, labels, lag_years = 0) {
  out <- runmulti_cox(
    data = data,
    main_var = exposures,
    covariates = covariates,
    endpoint = c("outcome_surv_time", "outcome_status")
  )
  out$lag_years <- lag_years
  merge(out, labels, by = "variable", all.x = TRUE)
}

lag_results <- lapply(c(0, 1, 2, 4), function(lag) {
  lag_data <- if (lag == 0) {
    analysis_data
  } else {
    sensitivity_exclude_early_events(
      data = analysis_data,
      endpoint = c("outcome_surv_time", "outcome_status"),
      n_years = lag,
      verbose = TRUE
    )
  }
  fit_pollution_cox(lag_data, main_exposures_z, pollutant_labels, lag_years = lag)
})

lag_results <- data.table::rbindlist(lag_results, use.names = TRUE)
lag_results[, p_bh := p.adjust(pvalue, method = "BH"), by = lag_years]
as.data.frame(lag_results[order(pollutant, lag_years)])
#>       variable    contrast         coef        se           z        HR
#> 1   no2_main_z  no2_main_z  0.409119575 0.2001311  2.04425836 1.5054917
#> 2   no2_main_z  no2_main_z  0.331617070 0.2052272  1.61585358 1.3932192
#> 3   no2_main_z  no2_main_z  0.345371242 0.2144938  1.61016880 1.4125142
#> 4   no2_main_z  no2_main_z  0.070924639 0.2635680  0.26909425 1.0735003
#> 5        nox_z       nox_z -0.571962396 0.1946940 -2.93775115 0.5644167
#> 6        nox_z       nox_z -0.519883339 0.2024629 -2.56779585 0.5945899
#> 7        nox_z       nox_z -0.540419935 0.2084588 -2.59245423 0.5825036
#> 8        nox_z       nox_z -0.487907336 0.2571520 -1.89734964 0.6139098
#> 9  pm10_2007_z pm10_2007_z -0.105484213 0.1969988 -0.53545613 0.8998887
#> 10 pm10_2007_z pm10_2007_z -0.133536912 0.2073504 -0.64401582 0.8749952
#> 11 pm10_2007_z pm10_2007_z -0.165176210 0.2212433 -0.74658172 0.8477443
#> 12 pm10_2007_z pm10_2007_z -0.221675303 0.2906503 -0.76268723 0.8011755
#> 13      pm25_z      pm25_z -0.005581028 0.1906796 -0.02926914 0.9944345
#> 14      pm25_z      pm25_z  0.050700808 0.2019440  0.25106369 1.0520081
#> 15      pm25_z      pm25_z  0.173807521 0.2189226  0.79392217 1.1898265
#> 16      pm25_z      pm25_z  0.483974513 0.3007344  1.60930856 1.6225103
#>      lower95   upper95      pvalue   n n_event lag_years pollutant       p_bh
#> 1  1.0170130 2.2285903 0.040928030 336      27         0       NO2 0.08185606
#> 2  0.9318152 2.0830954 0.106125974 333      24         1       NO2 0.21225195
#> 3  0.9277167 2.1506527 0.107361011 331      22         2       NO2 0.21472202
#> 4  0.6404020 1.7994993 0.787857153 323      14         4       NO2 0.78785715
#> 5  0.3853684 0.8266538 0.003306023 336      27         0       NOx 0.01322409
#> 6  0.3998350 0.8842076 0.010234741 333      24         1       NOx 0.04093896
#> 7  0.3871312 0.8764741 0.009529386 331      22         2       NOx 0.03811754
#> 8  0.3708654 1.0162317 0.057781808 323      14         4       NOx 0.21509777
#> 9  0.6116502 1.3239588 0.592334486 336      27         0      PM10 0.78977931
#> 10 0.5827855 1.3137194 0.519565176 333      24         1      PM10 0.69275357
#> 11 0.5494678 1.3079391 0.455316088 331      22         2      PM10 0.45531609
#> 12 0.4532373 1.4162165 0.445649948 323      14         4      PM10 0.59419993
#> 13 0.6843361 1.4450502 0.976649936 336      27         0     PM2.5 0.97664994
#> 14 0.7081480 1.5628386 0.801764866 333      24         1     PM2.5 0.80176487
#> 15 0.7747050 1.8273886 0.427240743 331      22         2     PM2.5 0.45531609
#> 16 0.8999157 2.9253180 0.107548887 323      14         4     PM2.5 0.21509777

9 Sensitivity analysis 2: 2010 exposure definition

NO2 and PM10 are redefined using 2010 exposure fields (p24003 and p24005), while NOx and PM2.5 are kept unchanged.

pollutant_labels_2010 <- data.frame(
  variable = sensitivity_exposures_2010_z,
  pollutant = c("NO2_2010", "NOx", "PM2.5", "PM10_2010"),
  stringsAsFactors = FALSE
)

analysis_data_2010 <- sensitivity_exclude_missing_covariates(
  data = demo_prepared[
    !is.na(outcome_status) &
      !is.na(outcome_surv_time) &
      outcome_surv_time > 0
  ],
  covariates = c(sensitivity_exposures_2010_z, covariates),
  stepwise = FALSE,
  verbose = TRUE
)

cox_results_2010 <- fit_pollution_cox(
  data = analysis_data_2010,
  exposures = sensitivity_exposures_2010_z,
  labels = pollutant_labels_2010,
  lag_years = 0
)
cox_results_2010$p_bh <- p.adjust(cox_results_2010$pvalue, method = "BH")
cox_results_2010
#>      variable    contrast        coef        se          z        HR   lower95
#> 1  no2_2010_z  no2_2010_z  0.26568811 0.2022776  1.3134826 1.3043282 0.8774207
#> 2       nox_z       nox_z -0.59030583 0.2039304 -2.8946431 0.5541578 0.3715759
#> 3 pm10_2010_z pm10_2010_z -0.15366181 0.2089349 -0.7354532 0.8575620 0.5694031
#> 4      pm25_z      pm25_z -0.02003608 0.1980889 -0.1011469 0.9801633 0.6647905
#>     upper95      pvalue   n n_event lag_years pollutant       p_bh
#> 1 1.9389468 0.189020396 326      25         0  NO2_2010 0.37804079
#> 2 0.8264552 0.003795899 326      25         0       NOx 0.01518359
#> 3 1.2915499 0.462063509 326      25         0 PM10_2010 0.61608468
#> 4 1.4451471 0.919433867 326      25         0     PM2.5 0.91943387