Clinical outcome of a no. 2 suture (Dynacord): Supplementary analysis report

Author

Affiliation

Corey Scholes, PhD

EBM Analytics

1 Introduction

This analysis links to the manuscript of the Dynacord (Depuy-Mitek, USA) product as one of two companion publications assessing new-to-market hardware. The dataset is derived from the PRULO registry snapshot and live tables. A protocol has been previously prepared for the registry (Scholes et al. 2023).

1.1 Reporting

The study was reported according to the RECORD guidelines (Benchimol et al. 2015) and companion checklist.

The analysis was conducted in RStudio IDE (RStudio 2024.12.0+467 “Kousa Dogwood” Release) using Rbase, quarto and attached packages to perform the following;

• Data import and preparation

• Sample selection

• Describe and address missingness

• Data manipulation, modelling and visualisation of;

  • Patient characteristics

  • Pathology characteristics (diagnosis)

  • Management and surgical technique

  • Treatment and repair survival

  • Adverse events and complications

  • Patient reported outcomes

• Publish to posit connect for dissemination

1.2 Preparation

Packages were loaded initially with pacman package. Citations were applied to each library at first use in the text.

if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  "tidycmprsk",
  "stopwords",
  "wordcloud",
  "ggsurvfit",
  "tidytext",
  "dplyr",
  "flextable",
  "litedown",
  "grateful",
  "modelsummary",
  "quantreg",
  "readr",
  "knitr",
  "cardx",
  "forcats",
  "gargle",
  "googledrive",
  "googlesheets4",
  "openxlsx2",
  "tidyverse",
  "tidymodels",
  "lubridate",
  "gt",
  "consort",
  "gtsummary",
  "survival",
  "ggplot2",
  "ggdist",
  "ggfortify",
  "mice",
  "marginaleffects",
  "patchwork",
  "naniar",
  "quantreg",
  "broom",
  "epoxy",
  "broom.helpers",
  "broom.mixed",
  "lme4",
  "stringr"
  )

Table 1: Summary of package usage and citations

Package	Version	Citation
base	4.4.2	R Core Team (2024a)
broom	1.0.10	Robinson, Hayes, and Couch (2025)
broom.helpers	1.22.0	Larmarange and Sjoberg (2025)
broom.mixed	0.2.9.6	Bolker and Robinson (2024)
cardx	0.3.0	Sjoberg, Yogasekaram, and de la Rua (2025)
consort	1.2.2	Dayim (2024)
dplyr	1.1.4	Wickham et al. (2023)
epoxy	1.0.0	Aden-Buie (2023)
flextable	0.9.10	Gohel and Skintzos (2025)
forcats	1.0.1	Wickham (2025a)
gargle	1.6.0	Bryan, Citro, and Wickham (2025)
ggdist	3.3.3	Kay (2024); Kay (2025)
ggfortify	0.4.19	Tang, Horikoshi, and Li (2016); Horikoshi and Tang (2018)
ggplot2	4.0.0	Wickham (2016)
ggsurvfit	1.2.0	Sjoberg et al. (2025)
googledrive	2.1.2	D’Agostino McGowan and Bryan (2025)
googlesheets4	1.1.2	Bryan (2025)
grid	4.4.2	R Core Team (2024b)
gt	1.1.0	Iannone et al. (2025)
gtsummary	2.4.0	Sjoberg et al. (2021)
knitr	1.50	Xie (2014); Xie (2015); Xie (2025a)
litedown	0.8	Xie (2025b)
lme4	1.1.37	Bates et al. (2015)
lubridate	1.9.4	Grolemund and Wickham (2011)
marginaleffects	0.30.0	Arel-Bundock, Greifer, and Heiss (2024)
mice	3.18.0	van Buuren and Groothuis-Oudshoorn (2011)
modelsummary	2.5.0	Arel-Bundock (2022)
naniar	1.1.0	Tierney and Cook (2023)
openxlsx2	1.21	Barbone and Garbuszus (2025)
pacman	0.5.1	Rinker and Kurkiewicz (2018)
patchwork	1.3.2	Pedersen (2025)
quantreg	6.1	Koenker (2025)
readr	2.1.5	Wickham, Hester, and Bryan (2024)
rmarkdown	2.30	Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2025)
scales	1.4.0	Wickham, Pedersen, and Seidel (2025)
stopwords	2.3	Benoit, Muhr, and Watanabe (2021)
stringr	1.6.0	Wickham (2025b)
survival	3.8.3	Terry M. Therneau and Patricia M. Grambsch (2000); Therneau (2024)
tidycmprsk	1.1.0	Sjoberg and Fei (2024)
tidymodels	1.4.1	Kuhn and Wickham (2020)
tidyr	1.3.1	Wickham, Vaughan, and Girlich (2024)
tidytext	0.4.3	Silge and Robinson (2016)
tidyverse	2.0.0	Wickham et al. (2019)
wordcloud	2.6	Fellows (2018)

The packages drawn on to produce the following report are summarised in Table 1.

1.3 Authorisations

Access to PRULO datasets was pre-authorised using the gargle package and googledrive.

1.4 Functions for Processing

A function was generated to retrieve files using the googledrive package, to call on later in the analysis for processing data imports.

get_specific_snapshot <- function(folder_name, base_folder_id = base_folder_id1) {
  tryCatch({
    # Check if the folder exists in the base directory
    folder <- googledrive::drive_ls(as_id(base_folder_id), pattern = paste0("^", folder_name, "$"))
    
    if(nrow(folder) == 0) {
      stop(paste("Folder", folder_name, "not found"))
    }
    
    # Find the snapshot file in the specified folder
    snapshot_file <- googledrive::drive_ls(
      folder$id, 
      pattern = "Registry data snapshot\\.xlsx$"
    )
    
    if(nrow(snapshot_file) == 0) {
      stop("No snapshot file found in specified folder")
    }
    
    # Return both pieces of information as a list
    return(list(
      snapshot = snapshot_file,
      folder_name = folder$name
    ))
    
  }, error = function(e) {
    stop(paste("Error finding specified snapshot:", e$message))
  })
}

ReplaceTermFun <- function(String) {
  # Function must handle vectors - process each element
  map_chr(String, function(s) {
    # Find the matched term in TargetTerms2
    match <- filter(TargetTerms2, s == Term)
    
    # Check if a match is found
    if (nrow(match) == 1) {
      return(match$ReplaceTerm)
    } else {
      # Return the original string if no match is found
      return(s)
    }
  })
}

1.5 Analysis Aim

To describe the clinical and patient-reported outcomes, in patients presenting for surgical review of shoulder pathology and electing to undergo reconstruction or repair of soft-tissue structures with a biodegradable anchor (Healix Advance BR, Depuy-Mitek, USA), at a private, regional orthopaedic clinic between 2020 - 2024.

1.6 Analysis Hypotheses

It was hypothesised that i) a low incidence of adverse events would be observed and ii) that significant improvements in general function (QuickDASH) and pathology-specific (WORC) outcomes would be observed at up to 12months follow up.

2 Methods

2.0.1 RECORD [4] - Study Design

Subgroup analysis of a clinical registry embedded into private practice. Observational, cohort design.

2.1 Data Import and Preparation

Data was retrieved using googlesheets4 to retrieve live database tables. Source files were specified and stored as global variables to call on in further functions.

# Authenticate for sheets using the same token
gs4_auth(token = drive_token())



ComplicTable <- googlesheets4::read_sheet(
  ss = SheetIDs$DbSS,
  sheet = "Complications",
  col_names = TRUE,
  col_types = "TcccDccD"
  )

IntraComplic <- googlesheets4::range_read(
  ss = SheetIDs$DbSS,
  range = "A3:AQ",
  sheet = "IntraComplications", 
  col_names = TRUE, 
  col_types = "ccccliicicccccccccccccccciccccccccccccccccc"
  )

ImplantTable <- googlesheets4::range_read(
  ss = SheetIDs$ImplantSS,
  sheet = "Dynacord",
  range = "F1:I",
  col_names = TRUE,
  col_types = "nccc"
  )

# AcctData <- googlesheets4::range_read(
#   ss = SheetIDs$AcctSS,
#   sheet = "AcctType2015", 
#   range = "A1:G",
#   col_names = TRUE, 
#   col_types = "ccccDcc"
#   )

#To match to acctData
PatientTable <- googlesheets4::range_read(
  ss = SheetIDs$DbSS,
  sheet = "Patient",
  range = "A10:N",
  col_names = FALSE,
  col_types = "DccccDcccDcicc"
  )


Patient_Col <- c(
  "PatientCreationDate",
  "PatientID",
  "LastName",
  "FirstName",
  "AlternateID",
  "DateOfBirth",
  "Sex",
  "RegistryStatus",
  "RegistryStatusNotes",
  "DateRegistryStatus",
  "NotificationMethod",
  "NoTreatmentRecords",
  "Email",
  "Phone"
)

colnames(PatientTable) <- Patient_Col

A static registry snapshot was retrieved using the pre-specified function (see Functions for Processing) and formatted using openxlsx based on the fixed date of preparation of the snapshot (31-Mar-2024) and using tidyverse syntax and associated packages (dplyr, lubridate). Date columns were prepared for further analysis using lubridate.

# Authenticate for sheets using the same token
gs4_auth(token = drive_token())

# To get a snapshot from a specific folder (e.g., "20230415")
specific_snapshot <- get_specific_snapshot("20240331")

temp_file1 <- tempfile(fileext = ".xlsx")
drive_download(
  file = specific_snapshot$snapshot$id,
  path = temp_file1,
  overwrite = TRUE
)

# Correction to reset back to excel origin
DaysDiff <- as.numeric(as.duration(interval(ymd("1899-12-30"), ymd("1970-01-01"))),"days")

SnapshotGen <- openxlsx2::read_xlsx(
  temp_file1,
  sheet = "ShoulderGeneral",
  colNames = TRUE,
  detectDates = TRUE
  ) |> dplyr::mutate(
    Cohort = "General"
  )

SnapshotRC <- openxlsx2::read_xlsx(
  temp_file1,
  sheet = "RotatorCuff",
  colNames = TRUE,
  detectDates = TRUE
  ) |> dplyr::mutate(
    Cohort = "Rotator Cuff"
  )

SnapshotGH <- openxlsx2::read_xlsx(
  temp_file1,
  sheet = "GlenohumeralInstability",
  colNames = TRUE,
  detectDates = TRUE
  ) |> dplyr::mutate(
    Cohort = "Glenohumeral Instability"
  )


STROBEInput <- openxlsx2::read_xlsx(
  temp_file1,
  sheet = "Strobe_Input",
  colNames = TRUE,
  detectDates = TRUE
  )

Dataframes were combined into one for further analysis.

SnapshotComb <- SnapshotRC |> dplyr::bind_rows(
  dplyr::select(SnapshotGH, -ExternalStudyTag)
  ) |> bind_rows(
    dplyr::select(
      SnapshotGen, 
      -ExternalStudyTag)
    ) |> dplyr::mutate(
    PatientID = stringr::str_split_i(TreatmentUID,"\\.",1)
  ) |> relocate(
    PatientID, .before = TreatmentUID
  ) |> mutate(
    CombID = paste0(PatientID,".",AffectedSide)
  ) |> relocate(
    CombID, .after = TreatmentUID
  ) |> relocate(
    Cohort, .before = EligibleAtPreop
  ) |> group_by(CombID) |> mutate(
  CombIDn = row_number(DateTreatment),
  CombIDn = tidyr::replace_na(CombIDn,1)
) |> ungroup()

An additional export of account data was prepared and imported to the workspace using the readr package, as well as tidyverse syntax and stringr, to categorise text fields.

#Read in full text file
#
temp_file2 <- tempfile(fileext = ".txt")
drive_download(
  file = as_id(AcctNewFile),
  path = temp_file2,
  overwrite = TRUE
)

AcctDataNew <- readr::read_tsv(
  file = temp_file2,
  col_names = TRUE,
  #trim_ws = TRUE,
  col_types = list(
  Id = "c",
  AccountType = "c",
  Surname = "c",
  FirstName = "c",
  DOB = col_date(format = "%d/%m/%Y"),
  UsualProvider = "c",
  HealthFundName = "c",
  HccPensionNum = "c",
  DvaNum = "c"
),
col_select = c(
  AlternateID = Id,
  LastName = Surname,
  FirstName,
  DateOfBirth = DOB,
  HealthFundName
  )
) |> dplyr::mutate(
  DateOfBirth2 = as.numeric(DateOfBirth),
  HealthFund2 = str_to_lower(HealthFundName),
  LastName = stringr::str_to_title(LastName)
) |> unite(
  col = "CombID",
  sep = ".",
  c("FirstName","LastName","DateOfBirth2"),
  remove = FALSE
) |> mutate(
  CombID = stringr::str_squish(CombID),
    AccountType2 = case_when(
      stringr::str_detect(HealthFund2,"nil|uninsured") == TRUE ~ "Uninsured",
      stringr::str_detect(HealthFund2,"fund|pty|limit*|nib|hcf|bupa|medibank|hbf|ahm|health|hba|a.u|^unity$|cbhs|unity|^yes$")  == TRUE ~ "Private",
      stringr::str_detect(HealthFund2,"dva|vaff|vet|aff|def")  == TRUE ~ "DVA",
      stringr::str_detect(HealthFund2,"work*|wc|w//c*")  == TRUE ~ "WorkCover",
      stringr::str_detect(HealthFund2,"^tac$")  == TRUE ~ "TAC",
      stringr::str_detect(HealthFund2,"sf|self")  == TRUE ~ "SelfFund",
      )
    ) |> left_join(
  PatientTable |> dplyr::select(
  PatientID,
  AlternateID
  ),
  by = "AlternateID"
  ) |> dplyr::filter(
  PatientID %in% SnapshotComb$PatientID
)

2.1.1 RECORD [5] - Setting

The PRULO registry is based in a regional private practice for upper limb orthopaedics (Scholes et al. 2023).

The registry has 2681 treatment records with the first patient enrolled 13 October 2020 and the final treatment record created 26 March 2024.The registry snapshot was extracted on 31 March 2024. Patients are followed for up to 2 years after surgery to capture treatment outcomes and patient-reported outcome measures (PROMs).

2.2 Record [6] Participants

Record [6.1] Sample selection

Identify cases receiving the suture of interest. Cases were identified by SKUs identified from the SKU database maintained as part of implant tracking within the registry. Cases were not restricted by available follow up.

Inclusion criteria;

• Case involves anchor of interest

• Case is the index procedure within the registry (first use of suture)

• Patient has not withdrawn consent for inclusion of data in the registry

• Treatment record is eligible for surgery (it has occurred)

Data manipulation (add columns and filter tables based on column values) was performed with tidyverse and converted to display format using gt.

# ImplantTable2 |> gt(
#   rowname_col = "row"
#   )
#   

knitr::kable(
  ImplantTable
)

Table 2: Summary of SKUs (Reference) used to identify cases of interest from the PRULO registry

Size (mm)	Description	Category	Reference
4.5	Healix Advance BR Dynacord (x2) with Needles	Anchor + Suture	10886705029440
4.5	Healix Advance BR Dynacord (x2)	Anchor + Suture	10886705029402
4.5	Healix Advance BR Dynacord (x3)	Anchor + Suture	10886705029396
5.5	Healix Advance BR Dynacord (x2)	Anchor + Suture	10886705029464
5.5	Healix Advance BR Dynacord (x3)	Anchor + Suture	10886705029457
5.5	Healix Advance BR Dynacord (x2) with Needles	Anchor + Suture	10886705029471
6.5	Healix Advance BR Dynacord (x2)	Anchor + Suture	10886705029525
6.5	Healix Advance BR Dynacord (x3)	Anchor + Suture	10886705029518
6.5	Healix Advance BR Dynacord (x2) with Needles	Anchor + Suture	10886705029532
4.5	Healix Advance PEEK Dynacord (x2) with Needles	Anchor + Suture	10886705029433
4.5	Healix Advance PEEK Dynacord (x2)	Anchor + Suture	10886705029426
4.5	Healix Advance PEEK Dynacord (x3)	Anchor + Suture	10886705029419
5.5	Healix Advance PEEK Dynacord (x2)	Anchor + Suture	10886705029495
5.5	Healix Advance PEEK Dynacord (x3)	Anchor + Suture	10886705029488
5.5	Healix Advance PEEK Dynacord (x2) with Needles	Anchor + Suture	10886705029501
6.5	Healix Advance PEEK Dynacord (x2)	Anchor + Suture	10886705029556
6.5	Healix Advance PEEK Dynacord (x3)	Anchor + Suture	10886705029549
6.5	Healix Advance PEEK Dynacord (x2) with Needles	Anchor + Suture	10886705029563
NA	Gryphon BR Dynacord BL	Anchor + Suture	10886705029877
NA	Gryphon BR Dynacord STR/BL	Anchor + Suture	10886705029884
NA	Gryphon P PEEK With Dynacord	Anchor + Suture	10886705029891
NA	Gryphon P PEEK DS Anchor with Dynacord	Anchor + Suture	10886705029907
NA	Dynacord #2 suture Pack Blue (with OS-6 needles)	Suture	222065
NA	Dynacord #2 suture Pack Blue (with MO-7 needles)	Suture	222066
NA	Dynacord #2 suture Pack Blue (without needles)	Suture	222067
NA	Dynacord #2 suture Pack Striped (without needles)	Suture	222068
NA	Dynacord #2 suture Pack Striped/Blue (without needles)	Suture	222069
NA	Dynacord #2 suture Pack Striped/Blue (with MO-7 needles)	Suture	222071
NA	Dynacord #2 suture Pack Striped/Blue (with OS-6 needles)	Suture	222073

A dataframe was prepared to generate a flow chart of record retrieval, screening and patient follow up within the sample of interest.

IncludeAny = paste(ImplantTable$Reference,sep = "|")

STROBEInput1 <- STROBEInput |> left_join(
  SnapshotComb |> dplyr::select(
    TreatmentUID,
    CombID,
    CombIDn,
    PatientID
  ),
  by = "TreatmentUID"
) |> dplyr::filter(
  stringr::str_detect(ImplantCodes,paste(IncludeAny,collapse = "|"))
) |> group_by(CombID) |> dplyr:: mutate(
  CombIDIndex = row_number(DateTreatment)
) |> ungroup()

#| label: Consort-Diagram
#| code-summary: "CONSORT|STROBE"

# Inclusion
# - Surgical treatment
# - With hardware of interest
# - Index procedure (first surgery of interest)
# After "induction"
# - Patient withdraws consent
# - Treatment fails before analysis date
# - Not eligible for 12m followup
# 
CurrentDate <- as.character("20240331")


STROBEFlow <- STROBEInput |> dplyr::filter(
  !is.na(TreatmentUID)
) |> dplyr::left_join(
  SnapshotComb |> dplyr::select(
    TreatmentUID,
    CombID,
    DateInitialExamination,
    EligibleAtPreop,
    EligibleAtx12months
  ),
  by = "TreatmentUID"
) |> left_join( #Bring in index procedure marker
  STROBEInput1 |> dplyr::select(
    TreatmentUID,
    CombIDIndex # This is the count of how many treatments with the CombID once dataset filtered for Dynacord SKUs
  ), 
  by = "TreatmentUID"
) |> dplyr::mutate(
  exclusion1 = case_when(
    is.na(SurgicalTreatment) ~ "Not a surgical treatment",
    SurgicalTreatment == "Surgical" ~ NA_character_,
    .default = "Not a surgical treatment"
  ),
  exclusion2 = case_when(
    is.na(exclusion1) & stringr::str_detect(RegistryStatus,"Opt-out") ~ "Patient Opt-Out",
    is.na(exclusion1) & (is.na(ImplantCodes) | stringr::str_detect(ImplantCodes,paste(IncludeAny,collapse = "|"), negate = TRUE))  ~ "No hardware of interest",
    is.na(exclusion1) & stringr::str_detect(ImplantCodes,paste(IncludeAny,collapse = "|")) & CombIDIndex != 1 ~ "Non-index procedure",
    is.na(exclusion1) & stringr::str_detect(ImplantCodes,paste(IncludeAny,collapse = "|")) & CombIDIndex == 1 ~ NA_character_
  ),
  followup = if_else(
    is.na(exclusion1) & is.na(exclusion2),
    TreatmentUID,
    NA_character_
  ),
  lost_followup = case_when(
    is.na(exclusion1) & is.na(exclusion2) & TreatmentStatus == "Failed" & (ymd(DateStatusChange) < ymd(CurrentDate)) ~ "Repair failure",
    is.na(exclusion1) & is.na(exclusion2) & TreatmentStatus == "No further followup" & (ymd(DateStatusChange) < ymd(CurrentDate)) ~ "Patient Opt-out",
    is.na(exclusion1) & is.na(exclusion2) & TreatmentStatus == "Ongoing" & is.na(EligibleAtx12months) ~ "Not eligible for followup"
  ),
  mitt = if_else(
    !is.na(followup) & is.na(lost_followup),
    TreatmentUID,
    NA_character_
)
) |> dplyr::rename(
  trialno = "TreatmentUID"
)

The combined snapshot dataframe was filtered using the results of the STROBE flowchart dataframe and the sample of interest retrieved.

Mastersheet <- SnapshotComb |> dplyr::select(
  Cohort,
  EligibleAtPreop:SignificantComorbidities_Preop,
  EligibleAtIntraop:OtherShoulderGirdle,
  starts_with("QDASH_TotalScore"),
  starts_with("EligibleAt"), 
  starts_with("WORC")
  ) |> filter(
    TreatmentUID %in% STROBEFlow$followup
  ) |> 
  relocate(
    Cohort, .before = EligibleAtPreop
  ) |> left_join(
    STROBEInput1 |> dplyr::select(
      TreatmentUID,
      CombIDn
    ),
    by = "TreatmentUID"
  ) |> relocate(
    CombIDn,
    .after = CombID
  )

Of the 255 records in the mastersheet, 0 treatment records had withdrawn consent for data inclusion and 3 had declined to participate in PROMs.

2.2.1 Record [6.2] Algorithm validation

Record selection code was cross-checked by manual record checking within the registry snapshot for a subset (N = 10) of cases.

2.2.2 Record [6.3] Data linkage

No data linkage was utilised for this analysis.

2.3 Record [7] Variables

Table 3: Summary of variables

Category	Variable	Comments	Citation
Patient Characteristics	Insurance Status	Recode from account data to insurance status
Pathology	Primary diagnosis	Free text coded using ICD-10 international
	CuffRetraction	Defined as per modified Patte grading	(Lädermann et al. 2016)
	CuffCondition	Fatty infiltration as assessed by Goutallier scale	(Fuchs et al. 1999)
	TearPattern	Shape the tear makes within the margins of the cuff as viewed in the transverse plane	(Lädermann et al. 2016)
	OtherShoulderPathology	Free-text coded as present [Yes] or not [No]
Management - Surgery	RepairAugment	Techniques used to augment the repair
	CuffTension	Surgeon perceived tension to restore anatomical footprint of repair
	RepairQuality	Surgeon subjective rating of the repair quality
Adverse Events	Modidifed sink grade	Modification of the Sink grading of complication severity	(Felsch et al. 2021)
Patient-Reported Outcomes	WORC Physical Q3	How much weakness do you experience in your shoulder?	(Kirkley, Alvarez, and Griffin 2003)

Key variables defined as part of this analysis are summarised Table 3.

2.4 Record [8] Data sources

Data was sourced directly from the PRULO clinical registry as described in (Scholes et al. 2023).Patient and treatment information were entered into the database through the registry interface and compiled into a data cube (snapshot) every quarter. Complications and adverse events captured into an online form (QuestionPro, USA) and linked using record identifier codes.

Adverse Events

The complications tables (intraop and postop) were also processed for further analysis.

# Join sufficient columns to complictable to perform calculations

ComplicTable1 <- ComplicTable |> rename_with(
  ~ (gsub(" ", "", .x))
  ) |> rename(
 TreatmentUID = "TreatmentID") |> left_join(
    SnapshotComb |> dplyr::select(
      CombID,
      TreatmentUID,
      DateTreatment),
    by = "TreatmentUID"
  ) |> mutate(
    ReopDelay = as.numeric(as.duration(interval(ymd(DateTreatment), ymd(ReoperationProcedureDate))),"weeks"),
    OccurDelay = as.numeric(as.duration(interval(ymd(DateTreatment), DateofOccurrence)),"weeks"),
    ComplicationNature = str_replace_all(str_to_lower(ComplicationNature),"[;.,-]",""),
    Intraop = case_when(
      ReopDelay == 0 ~ "Yes",
      str_detect(ComplicationNature,"intraop") ~ "Yes",
      .default = "No"
      ),
  ComplicationOccurrence2 = case_when(
    str_detect(str_to_lower(ComplicationOccurrence),"possible|yes") ~ "Yes",
    ComplicationOccurrence == "No" & !is.na(ComplicationNature) & str_detect(ComplicationNature,"nil", negate = TRUE) ~ "Yes",
    .default = ComplicationOccurrence
  )
) |> filter(
  stringr::str_detect(str_to_lower(ComplicationOccurrence),"yes|possible"),
  !is.na(CombID),
  OccurDelay <= 52
)

Intraoperative events were prepared to append to the Complications Table for further analysis.

IntraComp1 <- IntraComplic |> rename_with(
  ~ (gsub(" ", "", .x))
  ) |> filter(
    ICCount > 0
  ) |> dplyr::select(
    TreatmentID,
    ICNature1:ICIntervention1Technique # there is one case with two complications listed, skipped for now
  ) |> rename(
    TreatmentUID = "TreatmentID"
  ) |> mutate(
    ComplicationID = stringr::str_c(TreatmentUID,".","01"),
    OccurDelay = 0,
    ComplicationOccurrence2 = "Yes",
    Intraop = "Yes"
  ) |> left_join(
    SnapshotComb |> dplyr::select(
      TreatmentUID,
      CombID,
      DateTreatment
    ),
    by = "TreatmentUID"
  ) |> rename(
    DateofOccurrence = "DateTreatment",
    ComplicationNature = "ICNature1"
  ) |> filter(
    TreatmentUID %in% Mastersheet$TreatmentUID
  )

Reoperations were identified and subsetted to append to the Complications Table for survival analysis.

ComplicReop <- ComplicTable1 |> filter(
  ReoperationProcedure == "Yes"
) |> dplyr::select(
  !(c(
    OccurDelay,
    DateofOccurrence,
    ComplicationNature
  ))
) |> rename(
  DateofOccurrence = "ReoperationProcedureDate",
  OccurDelay = "ReopDelay"
) |> mutate(
  ComplicationNature = "Reoperation",
  ComplicationOccurrence2 = "Yes",
  Intraop = "No"
)   # need to add patientID to this table?

ComplicTable2 <- ComplicTable1 |> bind_rows(
  ComplicReop |> dplyr::select(
    CombID,
    ComplicationID,
    TreatmentUID,
    ComplicationOccurrence2,
    ComplicationNature,
    DateofOccurrence,
    OccurDelay,
    DateTreatment,
    Intraop
  ),
  IntraComp1 |> dplyr::select(
    CombID,
    ComplicationID,
    TreatmentUID,
    ComplicationOccurrence2,
    ComplicationNature,
    DateofOccurrence,
    OccurDelay,
    Intraop
  )
) |> dplyr::filter(
    TreatmentUID %in% Mastersheet$TreatmentUID & (DateofOccurrence < ymd(CurrentDate))
    )

Complication entries were written to an external file for co-author review.

googlesheets4::range_write(
  ss = SheetIDs$StudySS,
  data = ComplicTable2 |> dplyr::select(
    ComplicationID,
    TreatmentUID,
    CombID,
    DateTreatment,
    DateofOccurrence,
    ComplicationNature,
    ReoperationProcedure,
    ReoperationProcedureDate,
    ReopDelay,
    OccurDelay
    ),
  sheet = "Complications2",
  range = "A1:K",
  col_names = TRUE) 

The free-text describing the nature of the complication or adverse event was pre-processed using tidytext (v0.4.3) (Silge and Robinson 2016) to split into word tokens and remove stop words.

Stop <- tibble(get_stopwords())

ComplicTable3 <- tidytext::unnest_tokens(
  ComplicTable2,
  output = Term,
  input = ComplicationNature,
  token = "words",
  format = "text",
  to_lower = TRUE,
  drop = FALSE
) |> anti_join(
  Stop,
  by = c("Term" = "word")
) |> mutate(
  TermLength = stringr::str_length(Term)
)

Terms with less than five characters were extracted and reproduced in an external file for manual spelling of abbreviations. Terms with digits (e.g. L5) were removed.

# Retrieve terms less than 4 characters that likely need to be recast into full words

TargetTerms <- ComplicTable3 |> dplyr::select(
  ComplicationID,
  Term
  ) |> distinct(
  Term,
  .keep_all = TRUE
  ) |> mutate(
    TermLength = stringr::str_length(Term)
    ) |> filter(
     TermLength < 5 & stringr::str_detect(Term,"\\d",negate = TRUE)
  ) |> arrange(
    Term
  )

#Commented out after first use
googlesheets4::range_write(
ss= SheetIDs$StudySS,
data = TargetTerms |> dplyr::select(
  ComplicationID,
  Term
),
sheet = "ComplicTerm2",
range = paste0("A1:","B",nrow(TargetTerms)+1),
col_names = TRUE

)

The abbreviated terms with expanded definitions were read back into the workspace for replacement in the complication descriptions.

# read in new terms

TargetTerms2 <- googlesheets4::range_read(
  ss= SheetIDs$StudySS,
  sheet = "ComplicTerm2",
  range = "A1:C",
  col_names = TRUE,
  trim_ws = TRUE
) |> mutate(
  TargetTerm = paste0("\\b",Term,"\\b")
)

The terms were replaced and added to the dataframe containing complication data.

TargetTermsList <- str_c(TargetTerms2$TargetTerm, collapse = "|")

ComplicTable4 <- ComplicTable3 |> 
  mutate(
    Term1 = case_when(
      str_detect(Term, "\\d", negate = TRUE) ~ str_replace_all(Term, TargetTermsList, ReplaceTermFun),
      .default = Term
    )
  ) |> filter(
    str_detect(Term1, "\\d", negate = TRUE)
  ) |> tidytext::unnest_tokens(
  output = Term2,
  input = Term1,
  token = "words",
  format = "text",
  to_lower = TRUE,
  drop = FALSE
) |> anti_join(
  Stop,
  by = c("Term2" = "word")
) |> mutate(
  TermLength = stringr::str_length(Term2)
)

A figure displaying term frequency was generated using ggplot2 (v4.0.0) (Wickham 2016) and formatted for reporting using (v1.50) (knitr?).

Figure1 <- ComplicTable4 |>
  count(Term2, sort = TRUE) |>
  filter(n > 5) |>
  mutate(Term2 = reorder(Term2, n)) |>
  ggplot(aes(n, Term2)) +
  geom_col() +
  labs(y = NULL)


knitr::knit_print(Figure1)

A wordcloud was generated using wordcloud (v2.6) (Fellows 2018) to express the most common terms in the complication description free text field.

Figure2 <- ComplicTable4 |>
  count(Term2) |>
  with(wordcloud::wordcloud(Term2, n, max.words = 50))

knitr::knit_print(Figure2)

NULL

Figure 1

ComplicTable5 <- left_join(
  ComplicTable2,
  ComplicTable4 |> group_by(
    ComplicationID
  ) |> summarise(
    Term3 = str_c(Term2, collapse = " ")
  ) |> ungroup(),
  by = "ComplicationID"
  ) |> mutate(
    #Ditch reoperation term from Term as it doubles up the cases when binding the tables together and conducting the regex to identify reoperation cases; also removed Healix as a term as it overlaps with "heal" to pick up fail to heal cases
  Term4 = case_when(
    ComplicationNature == "Reoperation" ~ "Reoperation",
    .default = str_replace_all(Term3,"reoperation|healix","")
  ),
Infection = case_when(
    stringr::str_detect(Term4,"sepsis|sinus|washout|infect*|antibiotic*|vac.*dress*|culture|cellulitis|staphylococcus|acnes") ~ "Infection", #added "culture" as term for Infection
    .default = "No",
  ),
  ImplantRemoval = case_when(
    stringr::str_detect(Term4,"explant|remov*") & OccurDelay > 0 ~ "ImplantRemoval",
    stringr::str_detect(Term4,"explant|remov*") & OccurDelay == 0 ~ "Explant",
    .default = "No",
  ),
RepairFailure = case_when(
    stringr::str_detect(Term4,"(repair fail*)|fail.*repair|recur*|rupture|retear|torn|retorn|heal") ~ "RepairFailure",
    .default = "No",
  ),
Loosening = case_when(
   stringr::str_detect(Term4,"loose*") ~ "Loosening",
    .default = "No"
    ),
Capsulitis = case_when(
    stringr::str_detect(Term4,"adhes*|capsul*|hydrodil*|froz*|stiff*|restrict*|range.*motion") ~ "Capsulitis",
    .default = "No"
  ),
Neurological = case_when(
    stringr::str_detect(Term4,"nerve|palsy|radiculo*") ~ "Neurological",
    .default = "No"
  ),
Hardware = case_when(
    stringr::str_detect(Term4,"brok*|snap*|break") ~ "Hardware",
    .default = "No"
  ),
HardwareIrritation = case_when(
    stringr::str_detect(Term4,"annoy*|irritat|promine*") ~ "Hardware Irritation",
    .default = "No"
  ),
PainOther = case_when(
  str_detect(Term4,"persistent|infarction|pain|complex") ~ "Pain - Other",
    .default = "No"
),
Explant = case_when(
  str_detect(Term4,"explant*|intraop*|regraft") ~ "Explant",
  .default = "No"
 ),
  Fracture = case_when(
    stringr::str_detect(Term4,"fract*") ~ "Fracture",
    .default = "No"
  ),
Wound = case_when(
    stringr::str_detect(Term4,"(wound|flap|absce*|stitch|ooze|skin|rash)") & Infection == "No" ~ "Wound-Skin",
    .default = "No"
  ),
Thrombosis = case_when(
    stringr::str_detect(Term4,"vein|thromb*|embo*|occlus*") ~ "Thrombosis",
    .default = "No"
  ),
Reoperation = case_when(
  ReoperationProcedure == "Yes" ~ "Reoperation",
  is.na(ReoperationProcedure) ~ "No",
  .default = "No"
),
  PatientID = stringr::str_split_i(TreatmentUID,"\\.",1)
)

ComplicCensor <- Mastersheet |> filter(
  !(TreatmentUID %in% ComplicTable5$TreatmentUID)
) |> dplyr::select(
  TreatmentUID
) |> mutate(
  Censor = "Censor",
  DateofOccurrence = ymd(CurrentDate)
)

ComplicMasterRev <- bind_rows(
  ComplicTable5 |> dplyr::filter(
    !(ComplicationNature == "Reoperation"),
    !(Explant == "Explant")
  ) |> dplyr::select(
    -Explant
  ),
  ComplicCensor
) |> group_by(
  TreatmentUID
) |> arrange(DateofOccurrence) |> mutate(
  RecordN = row_number()
) |> ungroup() |> arrange(TreatmentUID)

# Collapse complications to one row per TreatmentUID
ComplicMasterTotal <- ComplicMasterRev |>
  group_by(TreatmentUID) |>
  summarise(
    # Keep first occurrence for complication-specific columns
    Timestamp = first(Timestamp),
    ComplicationID = first(ComplicationID),
    ComplicationOccurrence = first(ComplicationOccurrence),
    DateofOccurrence = first(DateofOccurrence),
    ComplicationNature = first(ComplicationNature),
    ReoperationProcedure = first(ReoperationProcedure),
    ReoperationProcedureDate = first(ReoperationProcedureDate),
    CombID = first(CombID),
    DateTreatment = first(DateTreatment),
    ReopDelay = first(ReopDelay),
    OccurDelay = first(OccurDelay),
    Intraop = first(Intraop),
    ComplicationOccurrence2 = first(ComplicationOccurrence2),
    Term3 = first(Term3),
    
    # Concatenate Term4 column
    Term4 = paste(unique(na.omit(Term4)), collapse = "; "),
    
    # Any "Yes" approach for binary columns
    Infection = if_else(any(Infection == "Infection", na.rm = TRUE), "Infection", "No"),
    ImplantRemoval = if_else(any(ImplantRemoval == "ImplantRemoval", na.rm = TRUE), "Implant Removal", "No"),
    RepairFailure = if_else(any(RepairFailure == "RepairFailure", na.rm = TRUE), "Repair Failure", "No"),
    Loosening = if_else(any(Loosening == "Loosening", na.rm = TRUE), "Loosening", "No"),
    Capsulitis = if_else(any(Capsulitis == "Capsulitis", na.rm = TRUE), "Capsulitis", "No"),
    Neurological = if_else(any(Neurological == "Neurological", na.rm = TRUE), "Neurological", "No"),
    Hardware = if_else(any(Hardware == "Hardware", na.rm = TRUE), "Hardware", "No"),
    HardwareIrritation = if_else(any(HardwareIrritation == "Hardware Irritation", na.rm = TRUE), "Hardware Irritation", "No"),
    PainOther = if_else(any(PainOther == "Pain - Other", na.rm = TRUE), "Pain - Other", "No"),
    Fracture = if_else(any(Fracture == "Fracture", na.rm = TRUE), "Fracture", "No"),
    Wound = if_else(any(Wound == "Wound", na.rm = TRUE), "Wound", "No"),
    Thrombosis = if_else(any(Thrombosis == "Thrombosis", na.rm = TRUE), "Thrombosis", "No"),
    Reoperation = if_else(any(Reoperation == "Reoperation", na.rm = TRUE), "Reoperation", "No"),
    Censor = if_else(any(Censor == "Censor", na.rm = TRUE), "Censor", "No"),
    
    # Keep other administrative columns
    PatientID = first(PatientID),
    Censor = first(Censor),
    RecordN = first(RecordN)
  ) |>  ungroup() |> rows_insert( # add two records that seem to have fallen off the complictable 
    tibble(
      TreatmentUID = c("1615.1","781.1"),
      Censor = c("Censor","Censor"),
      DateofOccurrence = c(ymd(CurrentDate),ymd(CurrentDate))
    )
  ) |> mutate(across(Infection:Reoperation, ~replace_na(.x, "No")))

MastersheetFail <- Mastersheet |> filter(
  TreatmentStatus == "Failed"
) |> dplyr::select(
  -ComplicationOccurrence
) |> left_join(
  ComplicTable5 |> dplyr::select(
    TreatmentUID,
    ComplicationOccurrence2
  ),
  by = "TreatmentUID"
) |> distinct(
  TreatmentUID,
  .keep_all = TRUE
)

Adverse Events Grading

Complication events were graded in a separate table according to (Felsch et al. 2021) and retrieved into the workspace using googlesheets4.

# comment out after first run

# Authenticate for sheets using the same token
gs4_auth(token = drive_token())

 googlesheets4::range_write(
   ss = SheetIDs$StudySS,
   data = ComplicTable5 |> filter(
    Term4 != "Reoperation"
   ) |> dplyr::select(
     TreatmentUID,
     ComplicationID,
     ComplicationNature,
     OccurDelay,
     Reoperation,
     ReopDelay,
     Term4:Thrombosis
   ),
   sheet = "ComplicGrade2",
   range = "A1:T"
 )

# Authenticate for sheets using the same token
gs4_auth(token = drive_token())

ComplicGraded <- googlesheets4::range_read(
  ss = SheetIDs$StudySS,
  sheet = "ComplicGrade2",
  range = "U1:V",
  col_names = TRUE,
  col_types = "ic"
  )

#Ready for table presentation

The dataset was reshaped to a long format and the indicator columns for each adverse event type were combined into one column within the dataframe (Category).

The date of surgery for the index procedure was linked to each complication entry and subsequent quantitative variables such as the durations between;

• date of surgery (index procedure) and date of occurrence

• date of occurrence and date of reoperation

• date of surgery (index procedure) and date of reoperation

Whether the event was intraoperative or presented postoperatively was also flagged within the table.

ComplicTable6 <- ComplicTable5 |> dplyr::select(!(c(
  ReopDelay,
  OccurDelay))
  ) |> pivot_longer(
  cols = Infection:Reoperation,
  cols_vary = "slowest",
  names_to = "Category_Name",
  values_to = "Category_Value",
  values_drop_na = FALSE
) |> filter(
  !(Category_Value == "No"),
  DateofOccurrence <= ymd(CurrentDate),
  !((ComplicationID == "1539.2.3" & Category_Value == "ImplantRemoval") | (ComplicationID == "1539.2.2" & Category_Value == "Hardware") | (ComplicationID == "164.1.2" & Category_Value == "Capsulitis"))
) |> distinct( #duplicate entries for the same patient(Treatment)
  pick(CombID,
  DateofOccurrence,
  Category_Value
  ),
  .keep_all = TRUE
) |> mutate(
  PatientID = stringr::str_split_i(TreatmentUID,"\\.",1)
)

Records with no complication recorded, as well as the final period of right-censor for each record that did not undergo removal of surgery hardware at the end of the chart review period (censored) were generated and added to the complication table to enable reorganisation into a format appropriate for the analysis selected.

MasterEnd <- ComplicTable6 |> filter(
  Category_Value == "ImplantRemoval"
)

EndDiscrep <- MasterEnd |> filter(
  !(TreatmentUID %in% MastersheetFail$TreatmentUID)
)

MasterDiscrep <- Mastersheet |> filter(
  TreatmentUID == "439.2"
)

MastersheetCensor <- Mastersheet |> filter(
  !(TreatmentUID %in% ComplicTable6$TreatmentUID)
) |> dplyr::select(
  PatientID,
  TreatmentUID,
  CombID,
  DateTreatment
) |> mutate(
  ComplicationOccurrence2 = "No",
  ) |> bind_rows(
   ComplicTable6 |> dplyr::select(
    PatientID,
    TreatmentUID,
    CombID,
    DateTreatment,
    DateofOccurrence
  )
  ) |> distinct(
    CombID,
    .keep_all = TRUE
  ) |> filter(
    !(CombID %in% MasterEnd$CombID),
    !(CombID %in% MastersheetFail$CombID)
  ) |> mutate(
    DateofOccurrence = ymd(CurrentDate),
  Category_Value = "Censored"
)

Censored treatment records (with no complication recorded at all) were combined with records that were censored after one or more complication events to form the Censored component of the adverse events dataset.

The censored data records were integrated into the dataset, with the resultant new frame reorganised into a format appropriate for a multi-state model (see RECORD 12.5) of procedure survival after use of the suture of interest, as described in the survival package (v3.8.3) (survival?).

A duration variable was calculated to arrange the dataframe rows within each PatientID in descending order of occurrence to establish the transition patterns from one health state to the next. The start and stop times for certain events (mortality, amputation) were offset by one week to remove ties for recurrent events or different event types occurring on the same date for the same patient.The presence of each adverse event type were restricted to the first occurrence of each Category within a patient subsequent to an index procedure per date of occurrence.

ComplicMaster <- bind_rows(
  ComplicTable6,
  MastersheetCensor
) |> mutate(
  Category = fct(
    Category_Value,
    levels = c(
      "Censored",
      "Explant",
      "Capsulitis",
      "RepairFailure",
      "Reoperation",
      "ImplantRemoval",
      "Infection",
      "Thrombosis",
      "Loosening",
      "Neurological",
      "Pain - Other"
      )
  )
) |> group_by(
    CombID
  ) |> arrange(
    by_group = TRUE,
    DateofOccurrence
  ) |> ungroup() |> filter( # remove reoperations that occur with end states
    #Category_Value != "Reoperation",
    !(CombID %in% MasterEnd$CombID | (CombID %in% MastersheetFail$CombID))
  ) |> mutate(
    DateofOccurrence1 = case_when(
      Category_Value == "RepairFailure" | Category_Value == "ImplantRemoval" | Category_Value == "Reoperation" ~ DateofOccurrence + lubridate::days(3),
      .default = DateofOccurrence
    )
  ) |> relocate(
    DateofOccurrence1, .after = DateofOccurrence
  ) |> slice_min( # take first occurrence of each category per Patient
  by = c(
    CombID,
    Category_Value),
  DateofOccurrence,
  n = 1
  ) |> group_by(
    CombID
  ) |> mutate(
    Duration = case_when(
    Category != "Censored" ~ as.numeric(as.duration(interval(ymd(DateTreatment), ymd(DateofOccurrence1))),"weeks"),
    Category == "Censored" ~ as.numeric(as.duration(interval(ymd(DateTreatment), ymd(CurrentDate))),"weeks"),
    )
  ) |> arrange(
    CombID,
    Duration,
    .by_group = TRUE
  ) |> mutate(
    RowNum = row_number()
  )  |> mutate(
  DurationStart = case_when(
    RowNum > 1 ~ dplyr::lag(Duration),
    .default = 0
  )
  ) |> mutate(
    DurationDiff = as.numeric(Duration-DurationStart)
  ) |> rename(
    DurationStop = "Duration"
  ) |> relocate(
    DurationStop, .after = DurationDiff
  ) |> mutate(
    DurationStop1 = case_when(
      DurationDiff == 0 ~ DurationStop + 0.17,
      .default = DurationStop
      )
    ) |> mutate(
  DurationStart1 = case_when(
    RowNum > 1 ~ dplyr::lag(DurationStop1),
    .default = 0
    )
  ) |> relocate(
    DurationStop1, .after = DurationStart1
  ) |> dplyr::select(
    ComplicationID,
    CombID,
    DurationStart1,
    DurationStop1,
    Category,
    RowNum
  )

2.5 Record [9] Bias

For a discussion of biases in the context of the clinical registry utilised for this analysis, refer to (Scholes et al. 2023). Specific to this analysis, the following considerations were noted;

Table 4: Biases in analysis of observational cohort of a clinical registry

Bias	Definition	Source	Mitigation
Misclassification	Treatment record labelled into incorrect cohort. PROMs package not aligned to clinical presentation	(Benchimol et al. 2015)	Clinical notes reviewed by experienced reviewer and matched to ICD10 code by definition.
Confounder	An variable of interest and a target outcome simultaneously influenced by a third variable	(Tennant et al. 2020)	PROMs analysis incorporated adjustment for age and sex
Missing data	The absence of a data value where a treatment record is eligible to have a data value collected	(Carroll, Morris, and Keogh 2020)	Multiple imputation utilised
Prevalent user	Follow-up starts after eligible individuals have started the treatment. The follow-up time is left-truncated	(Nguyen et al. 2021)	Eligibility and enrollment is performed prior to treatment offering for any patient or new presentation. Index procedures identified for analysis are followed prior to surgery occurring.
Selection	Treatments are selected based on post-treatment criteria	(Nguyen et al. 2021)	Unable to be mitigated fully - records are identified by presence of hardware code associated with suture of interest
Immortal time	Individuals need to meet eligibility criteria that can only be assessed after follow-up has started	(Nguyen et al. 2021)	Patients enrolled at time of diagnosis
Pseudoreplication	Analyse data while ignoring dependency between observations. Inadequate model specification.	(Davies and Gray 2015; Lazic 2010)	Cluster for patient in survival (all-cause failure and retear). Utilise mixed effects linear model (lme4::lmer) for PROMs analysis with treatment identifier as random effect

2.6 Record [10] Sample size

Sample size was derived from the available records of the Registry at the time of analysis.

2.7 Record [11] Quantitative variables

The anterior-posterior (AP) and mediolateral (ML) dimensions of the cuff tear were reported and multiplied to calculate tear area (mm^2). The tear was also classified according to (Rashid et al. 2017).

• Small tears were defined as full-thickness defects in the supraspinatus tendon under 1 cm in the anterior–posterior (AP) dimension.

• Medium tears were defined as full-thickness defects in the supraspinatus tendon only, greater than 1 cm and less than 3 cm in the AP dimension.

• Large tears involved full-thickness defects of both the supraspinatus and infraspinatus tendons, greater than 3 cm, and less than 5 cm in the AP dimension.

• Massive tears involved all 3 tendons (supraspinatus, infraspinatus, and subscapularis) and were greater than 5 cm in the AP dimension.

Partial tears were left labelled as partial. Ultimately recoded tear classification based on AP tear length, as the involvement of other tendons for tears of small length was not adequately defined in the original paper.

## Slice inputs for columns and rows

##Sort out variable types; Calculate wait time

Mastersheet1 <- Mastersheet |> rename_with(
   ~sub("_TotalScore_","_",.), 
   contains("_TotalScore_")
   ) |> rename_with(
     ~sub("EligibleAtx","Eligible_",.), 
     starts_with("EligibleAtx")
     ) |> rename_with(
       ~sub("EligibleAt","Eligible_",.), 
       starts_with("EligibleAt")
       ) |> rename_with(
         ~sub("WORC_","WORC",.), 
         starts_with("WORC_")
         ) |> mutate(
    CuffTearSizeML = as.numeric(CuffTearSizeML),
    CuffTearSizeAP = as.numeric(CuffTearSizeAP),
    WaitTime = as.numeric(as.duration(interval(ymd(DateInitialExamination), ymd(DateTreatment))),"weeks"),
    across(starts_with("WORC"),as.numeric),
    across(starts_with("QDASH"),as.numeric),
    across(where(is.character) & !contains("ID", ignore.case = TRUE), as.factor),
    Sex2 = case_when(
      Sex == "F" ~ "Female",
      Sex == "M" ~ "Male"),
    Surgeon2 = case_when(
      Surgeon == "KE" ~ "A",
      Surgeon == "RP" ~ "B",
      Surgeon == "GB" ~ "C",
      .default = NA_character_),
    TearClass = case_when(
      CuffStatus == "Partial Tear" ~ "Partial",
      CuffStatus == "Full Tear" & CuffTearSizeAP <=10 ~ "Small",
      CuffStatus == "Full Tear" & between(CuffTearSizeAP,11,30) ~ "Medium",
      CuffStatus == "Full Tear" & between(CuffTearSizeAP,31,50) ~ "Large",
      CuffStatus == "Full Tear" & CuffTearSizeAP > 50 ~ "Massive",
      .default = NA_character_),
    TendonsInvolved = case_when(
      CuffTendonsInvolved == "No other tendon involved" ~ "Supraspinatus (isolated)",
      .default = CuffTendonsInvolved
    ),
    TearArea = CuffTearSizeAP * CuffTearSizeML,
    OtherShoulderPathology = case_when(
      stringr::str_detect(str_to_lower(OtherShoulderGirdle),"no",negate = TRUE) & !is.na(OtherShoulderGirdle) ~ "Yes",
      .default = OtherShoulderGirdle
      ),
    RepairAugment2 = case_when(
      stringr::str_detect(str_to_lower(AdjunctProcedures), "scr|superior") ~ "Superior Capsular",
      RepairAugmentation == "Nil" ~ "None",
      .default = RepairAugmentation
      ),
    TreatStatus = case_when(
    TreatmentStatus == "Failed" ~ 1,
    .default = 0
  ),
  TreatEndDate = coalesce(DateStatusChange,ymd(CurrentDate)),
  TreatDuration = as.numeric(as.duration(interval(ymd(DateTreatment), ymd(TreatEndDate))),"weeks"),
  TreatmentStatus2 = case_when(
    stringr::str_detect(str_to_lower(TreatmentStatus),"pend|further") ~ "Ongoing",
      .default = TreatmentStatus
    )
  ) |> relocate(
      TearClass, 
      .after = CuffCondition) |> relocate(
        TearArea, 
        .after = CuffCondition
        )

Data was read in from database table to determine account type.

PatientTable1 <- PatientTable |> filter(
  !is.na(PatientCreationDate) & !is.na(DateOfBirth)
) |> mutate(
  DateOfBirth2 = as.numeric(DateOfBirth),
  LastName = stringr::str_to_title(LastName)
) |> unite(
  col = "CombID",
  sep = ".",
  c("FirstName","LastName","DateOfBirth2"),
  remove = FALSE
)

Mastersheet1b <- Mastersheet1 |> left_join(
  PatientTable1 |> dplyr::select(
    FirstName,
    LastName,
    PatientID
  ),
  by = "PatientID"
)  |> mutate(
  DateOfBirth2 = as.numeric(as_date(DateOfBirth))
  ) |> unite(
    col = "CombID",
    sep = ".",
    c("FirstName","LastName","DateOfBirth2"),
    remove = FALSE
) |> left_join(
  AcctDataNew |> dplyr::select(
    AccountType2,
    HealthFund2,
    PatientID
    ),
  by = "PatientID"
  )

Procedure details were extracted from the master table and processed to enable presentations in summary tables.

OpData <- Mastersheet1b |> dplyr::select(
  TreatmentUID,
  SurgicalTreatment,
  PatientPosition,
  ArthroscopicApproach:OtherShoulderGirdle
  ) |> dplyr::mutate(
  AdjunctProcedures = as.character(AdjunctProcedures)
) |> separate_longer_delim(
    AdjunctProcedures,
    stringr::regex("[;,]")
    ) |> dplyr::mutate(
      AdjunctProcedures2 = str_squish(AdjunctProcedures)
      ) |> dplyr::mutate(
      AdjunctProcedures3 = case_when(
      AdjunctProcedures2 == "None performed" ~ "None",
      str_detect(str_to_lower(AdjunctProcedures2), "resection|osteot") ~ "Clavicle Resection",
str_detect(str_to_lower(AdjunctProcedures2), "release|debulk") ~ "Ligament|Capsule Release",
str_detect(str_to_lower(AdjunctProcedures2), "transfer") ~ "Tendon Transfer",
str_detect(str_to_lower(AdjunctProcedures2), "debride") ~ "Interval|Acromion|Capsule Debridement",
str_detect(str_to_lower(AdjunctProcedures2), "burs|excis") ~ "Bursectomy",
str_detect(str_to_lower(AdjunctProcedures2), "scr|superior|calcific|incorporat|supraspin|cable") ~ "None",
      .default = AdjunctProcedures2
    )
  )

OpData2 <- OpData |>
  group_by(TreatmentUID) |>
    summarise(AdjunctProcedures4 = paste(AdjunctProcedures3,collapse = "; ")
              ) |> mutate(
                AdjunctProcedures5 = case_when(
                  stringr::str_detect(AdjunctProcedures4,"None") & str_length(AdjunctProcedures4) > 4 ~ str_replace_all(AdjunctProcedures4,"(None|None[;\\s]*|[;\\s]*None|[;\\s]+)",""),
                  .default = AdjunctProcedures4
                )
              ) |> mutate(
                AdjunctProcedures6 = str_squish(AdjunctProcedures5)
              ) |> dplyr::select(-AdjunctProcedures4,-AdjunctProcedures5)

Mastersheet2 <- left_join(
  Mastersheet1b,
  OpData2,
  by = "TreatmentUID"
)

Tables were rearranged with tidyverse to prepare patient-reported outcomes (PROMs) for analysis in the long format. Separate dataframes were created for the QuickDASH and the WORC, as the QuickDASH was collected at 3months and the WORC was not.

MasterPROM <- Mastersheet2 |> dplyr::select(
  TreatmentUID,
  starts_with("QDASH"),
  starts_with("WORC"),
  starts_with("Eligible")
) |> pivot_longer(
  cols = !TreatmentUID,
  names_to = c(".value","TimePoint"),
  names_sep = "_",
  values_drop_na = TRUE
) |> mutate(
  TimePoint = factor(TimePoint, levels = c("Preop","3months","6months","12months"), ordered = TRUE, exclude = NA),
) |> filter(!is.na(TimePoint))

MasterPROMWORC <- MasterPROM |> filter(TimePoint != "3months")

Tables were modified to track anchor usage.

AnchorUsage <- Mastersheet2 |> dplyr::select(
  TreatmentUID,
  ImplantCodes
) |> separate_longer_delim(cols = ImplantCodes,", ") |> left_join(
  ImplantTable |> dplyr::select(
    Reference,
    Category
  ),
  by = join_by("ImplantCodes" == "Reference")
)

2.8 Record [12] Statistical methods

A number of analytical techniques were employed to i) clean the data inputs as well as ii) evaluate missingness in the dataset and iii) complete the descriptive analysis of;

• Patient characteristics

• Pathology details

• Patient, implant and adverse event time to event

• Patient-reported outcomes

2.8.1 Record [12.1] Access to population

The registry system represents all cases presenting to the rooms of a surgical group within Geelong, Australia using the implant of interest from the inception of the clinical registry to the analysis date. All reviewed charts from the operating surgeons practice records (electronic medical record) were entered into database and the present analysis draws data from a regular compilation of the registry records (snapshot) produced quarterly by the registry administration team.

2.8.2 Record [12.2] Data cleaning methods

Complication descriptions were pre-processed to remove relational terms (stopwords) and expand abbreviations to improve clarity.

Dates of events (preceding and subsequent surgical records; adverse events including mortality) relative to index surgery date were assessed using coded checks to flag anomalies and were resolved by further manual review to resolve inconsistencies or discrepancies with the chart review input data stored in the registry database.

Diagnosis and complication description free text fields were pre-processed to remove relational terms (stopwords) and expand abbreviations to improve clarity.

The dataset used as input for the survival analysis of adverse outcomes was assessed survival analysis, with visual assessment of the transitions table to ensure procedure endstates (mortality, implant removal) did not have subsequent states and that the numbers of events and unique identifiers matched the numbers in the dataframe.

SurvCheck <- survival::survcheck(Surv(DurationStart1, DurationStop1, Category) ~1, ComplicMaster, id=CombID) 

knitr::knit_print(SurvCheck)

Call:
survival::survcheck(formula = Surv(DurationStart1, DurationStop1, 
    Category) ~ 1, data = ComplicMaster, id = CombID)

Unique identifiers       Observations        Transitions 
               234                292                 58 
4 observations removed due to missing 

Transitions table:
               to
from            Capsulitis RepairFailure Reoperation Infection Neurological
  (s0)                  31            12           0         1            2
  Capsulitis             0             3           1         0            0
  RepairFailure          0             0           1         0            0
  Reoperation            0             0           0         0            0
  Infection              0             0           1         0            0
  Neurological           0             0           0         0            0
  Pain - Other           0             0           0         0            0
               to
from            Pain - Other (censored)
  (s0)                     4        184
  Capsulitis               1         26
  RepairFailure            1         13
  Reoperation              0          3
  Infection                0          0
  Neurological             0          2
  Pain - Other             0          6

Number of subjects with 0, 1, ... transitions to each state:
               count
state             0  1 2
  Capsulitis    203 31 0
  RepairFailure 219 15 0
  Reoperation   231  3 0
  Infection     233  1 0
  Neurological  232  2 0
  Pain - Other  228  6 0
  (any)         184 42 8

2.8.3 Record [12.3] Data linkage

Not applicable

2.8.4 Record [12.4] Missingness

Evaluation

Missingness was assessed with visualisation and table functions in the naniar package and compiled into figures using patchwork.

# Assessing missingness from Mastersheet2

#Patient characteristics

MissFig1 <- Mastersheet2 |>
 dplyr::select(
    AgeAtTreatment, 
    Sex2, 
    IndexSide,
    Surgeon2,
    BMI,
    BilateralStatus,
    WaitTime,
    AccountType2) |> rename(
      `Age at Surgery` = "AgeAtTreatment", 
    Sex = "Sex2", 
    `Affected Side` = "IndexSide",
    Surgeon = "Surgeon2",
    `Body Mass Index` = "BMI",
    `Bilateral Status` = "BilateralStatus",
    `Initial Consult to Surgery` = "WaitTime",
    `Insurance Status` = "AccountType2"
    ) |> gg_miss_var(show_pct = TRUE) + labs(title = "Patient Characteristics")

#Pathology characteristics

MissFig2 <- Mastersheet2 |>
  dplyr::select(
    TreatmentType,
    CuffStatus,
    CuffCondition,
    CuffTearRetraction,
    CuffTendonDelaminated,
    TendonsInvolved,
    CuffTearSizeAP, 
    CuffTearSizeML,
    TearArea, 
    TearClass, 
    CuffTearPattern,
    OtherShoulderPathology
    ) |> rename(
   `Primary Presentation` = "TreatmentType",
    `Full Tear` = "CuffStatus",
    `Cuff Condition` = "CuffCondition",
    `Tendon Retraction` = "CuffTearRetraction",
    `Tendon Delamination` = "CuffTendonDelaminated",
    `Tendons Involved (+Ssp)` = "TendonsInvolved",
    `Tear Size AP (mm)` = "CuffTearSizeAP", 
    `Tear Size ML (mm)` = "CuffTearSizeML",
    `Tear Area (mm^2)` = "TearArea", 
    `Tear Classification` = "TearClass", 
     `Tear Pattern` = "CuffTearPattern",
    `Other Pathology` = "OtherShoulderPathology" 
    ) |> gg_miss_var(show_pct = TRUE) + labs(title = "Pathology Characteristics")

#Surgical Technique

MissFig3 <- Mastersheet2 |>
  dplyr::select(
    ArthroscopicApproach,
    PatientPosition,
    CuffTendonsTreated,
    RepairType,
    AnchorFixation,
    RepairAugment2,
    CuffRepairTension,
    CuffRepairQuality
    ) |> rename(
          Arthroscopy = "ArthroscopicApproach",
         `Beachchair Position` = "PatientPosition",
         `Supraspinatus Isolated Repair` = "CuffTendonsTreated",
         `Double Row Repair` = "RepairType",
         `Knotted Anchor Fixation` = "AnchorFixation",
         `Superior Capsular Augment` = "RepairAugment2",
         `Low Repair Tension` = "CuffRepairTension",
         `Anatomic Repair` = "CuffRepairQuality"  
    ) |> gg_miss_var(show_pct = TRUE) + labs(title = "Surgical Details")

#PROMs
MissFig4 <- MasterPROM |>
  dplyr::select(QDASH,
                TimePoint) |>  gg_miss_var(show_pct = TRUE,
                                        facet = TimePoint) + labs(title = "Patient Reported\n Outcomes - 1")


MissFig5 <- MasterPROMWORC |>
  dplyr::select(WORCNorm,
                TimePoint) |>  gg_miss_var(show_pct = TRUE,
                                        facet = TimePoint) + labs(title = "Patient Reported\n Outcomes - 2")

MissFig1 + MissFig2 + MissFig3 + MissFig4 + plot_layout(ncol = 2)

MissFig5

Figure 2: Missingness rates of patient, pathology, management and patient-reported outcomes

Figure 3: Missingness rates of patient, pathology, management and patient-reported outcomes

PROMsmiss <- MasterPROM |> dplyr::select(-TreatmentUID) |>
  group_by(TimePoint) |>
  miss_var_summary() |>
  pivot_wider(names_from = TimePoint, 
              values_from = c(pct_miss))

The compliance for the QuickDASH at baseline was 63.5% and for the WORCNorm it was 50.2%.

The compliance for the QuickDASH was 46.3% and for the WORCNorm it was 41.3% at 12months.

Management

The data tables were reduced to the required columns (PROMs and adjunct columns) in preparation for multiple imputation using chained equations (White, Royston, and Wood 2010) with the mice package. One patient with bilateral records in the sample had one field (EducationLevel_Preop) mirrored from one side record to the other, where it was missing. Character fields were converted to factors and the dataset was filtered to those cases that were eligible for 12months followup.

QDASHPROMInput <- MasterPROM |> left_join(
  dplyr::select(
    Mastersheet2,
    TreatmentUID,
    WaitTime,
    IndexSide,
    BilateralStatus,
    EducationLevel_Preop,
    AgeAtTreatment,
    Sex2,
    ComplicationOccurrence
    ),
  by = "TreatmentUID"
) |> dplyr::mutate(
  Complication2 = case_when(
    ComplicationOccurrence != "No" ~ "Yes",
    .default = "No"
  ),
  across(where(is.character) & !contains("ID", ignore.case = TRUE), as.factor),
  TimePoint = fct_relevel(
      as.factor(TimePoint), 
      c("Preop", "3months", "6months", "12months")
    ),
    Sex2 = fct_relevel(
      as.factor(Sex2), 
      c("Female", "Male")
    ),
    BilateralStatus = fct_relevel(
      as.factor(BilateralStatus), 
      c("No", "Yes")
    ),
    EducationLevel_Preop = fct_relevel(
      as.factor(EducationLevel_Preop), 
      c(
        "Up to Secondary Year 10", 
        "Secondary - Year 12",
        "Post-secondary trade certificate or diploma",
        "Undergraduate degree",
        "Postgraduate degree")
      ),
    TreatmentInt = as.integer(str_replace(TreatmentUID, "\\.", ""))
) |> dplyr::select(
  !c(
    ComplicationOccurrence,
    Eligible,
    WORCPhysical:WORCEmotions,
    WORCPhysicalQ3,
    WORCNorm
    )
  ) |> dplyr::filter(
    TreatmentUID %in% STROBEFlow$mitt
  ) |> mutate(
    
  ) |> dplyr::rows_update(
    tibble(
    TreatmentUID = "253.2",
    EducationLevel_Preop = "Undergraduate degree"
    ),
    by = "TreatmentUID"
  )

A row was inserted for one case that did not return an entry for the 3month timepoint. The dataframe was reordered to create a visitsequence for the multiple imputation function.

# 
QDASHPROMInput2 <- dplyr::rows_append(
  QDASHPROMInput,
  tibble(TreatmentUID = "1553.1",
  TimePoint = "3months",
  BilateralStatus = "No",
  AgeAtTreatment = 74,
  Sex2 = "Female",
  TreatmentInt = 22,
  IndexSide = "Dominant",
  WaitTime = 24.8571429,
  Complication2 = "No",
  EducationLevel_Preop = NA_character_
  )
) |> dplyr::select(
    -TreatmentUID
  ) |> dplyr::relocate( # reorder columns for visitsequence
    all_of(
      c(
        "TreatmentInt",
        "TimePoint",
        "Sex2",
        "AgeAtTreatment",
        "BilateralStatus",
        "WaitTime",
        "IndexSide",
        "EducationLevel_Preop",
        "Complication2",
        "QDASH"
      )
    )
    
  )

A predictor matrix was generated to specify the combination of variables to be drawn on for the imputation of each column in the dataset. In addition, a method matrix was created to specify varying univariate imputations to account for the multilevel nature of the dataset (van Buuren and Groothuis-Oudshoorn 2011). Patient level variables (education, sex, bilateralstatus) were imputed as level-2 variables and the PROMs columns treated as level-1 variables.

# TreatmentInt - no missing
# TimePoint - no missing
# Sex2 - no missing
# AgeAtTreatment - no missing
# BilateralStatus - no missing
# WaitTime - no missing
# IndexSide - Level 2 missing
# EducationLevel_Preop - level 2 missing
# Complication2 - not missing
# QDASH - level 2 missing


# Make predictor matrix
# 
QDASHPred <- make.predictorMatrix(
  QDASHPROMInput2
)

QDASHPred["TreatmentInt", ] <- 0
QDASHPred["TimePoint", ] <- 0
QDASHPred["Sex2", ] <- 0
QDASHPred["AgeAtTreatment", ] <- 0
QDASHPred["BilateralStatus", ] <- 0
QDASHPred["WaitTime", ] <- 0
QDASHPred["Complication2", ] <- 0
QDASHPred["IndexSide", ] <- c(-2,0,1,1,1,1,0,1,1,1)
QDASHPred["EducationLevel_Preop", ] <- c(-2,0,1,1,1,1,1,0,1,1)
QDASHPred["QDASH", ] <- c(-2,1,1,1,1,1,1,1,1,0)

# Specify Method
# Method specification
QDASHMeth <- make.method(QDASHPROMInput2)

# For time-invariant variables (level-2)
QDASHMeth["EducationLevel_Preop"] <- "2lonly.pmm"  # Or "2lonly.polyreg" for categorical
QDASHMeth["IndexSide"] <- "2lonly.bin"  # Or "2lonly.polyreg" for categorical


# For time-varying variables (level-1)
QDASHMeth["QDASH"] <- "pmm"  # Multilevel imputation allowing for within-person variation
# Visit sequence
# 
QDASHseq <- mice::make.visitSequence(QDASHPROMInput2)

QDASHImputed <- mice::mice(
  data = QDASHPROMInput2,
  m = 10,
  predictorMatrix = QDASHPred,
  visitSequence = QDASHseq,
  method = QDASHMeth,
  maxit= 15,
  printFlag = FALSE
)

plot(QDASHImputed)

Figure 4: Stability of imputed variables over iterations for QuickDASH dataset

WORCPROMInput <- MasterPROM |> left_join(
  dplyr::select(
    Mastersheet2,
    TreatmentUID,
    WaitTime,
    IndexSide,
    BilateralStatus,
    EducationLevel_Preop,
    AgeAtTreatment,
    Sex2,
    ComplicationOccurrence
    ),
  by = "TreatmentUID"
) |> dplyr::mutate(
  Complication2 = case_when(
    ComplicationOccurrence != "No" ~ "Yes",
    .default = "No"
  ),
  across(where(is.character) & !contains("ID", ignore.case = TRUE), as.factor),
  TimePoint = fct_relevel(
      as.factor(TimePoint), 
      c("Preop", "3months", "6months", "12months")
    ),
    Sex2 = fct_relevel(
      as.factor(Sex2), 
      c("Female", "Male")
    ),
    BilateralStatus = fct_relevel(
      as.factor(BilateralStatus), 
      c("No", "Yes")
    ),
    EducationLevel_Preop = fct_relevel(
      as.factor(EducationLevel_Preop), 
      c(
        "Up to Secondary Year 10", 
        "Secondary - Year 12",
        "Post-secondary trade certificate or diploma",
        "Undergraduate degree",
        "Postgraduate degree")
      ),
    TreatmentInt = as.integer(str_replace(TreatmentUID, "\\.", ""))
) |> dplyr::select(
  !c(
    ComplicationOccurrence,
    Eligible,
    WORCPhysical:WORCEmotions
    )
  ) |> dplyr::filter(
    TreatmentUID %in% STROBEFlow$mitt,
    TimePoint != "3months"
  ) |> mutate(
    
  ) |> dplyr::rows_update(
    tibble(
    TreatmentUID = "253.2",
    EducationLevel_Preop = "Undergraduate degree"
    ),
    by = "TreatmentUID"
  )

# Fix up missing row (General switched to RotatorCuff)
# 
WORCPROMInput2 <- WORCPROMInput |> dplyr::select(
    -TreatmentUID
  ) |> dplyr::relocate( # reorder columns for visitsequence
    all_of(
      c(
        "TreatmentInt",
        "TimePoint",
        "Sex2",
        "AgeAtTreatment",
        "BilateralStatus",
        "WaitTime",
        "IndexSide",
        "EducationLevel_Preop",
        "Complication2",
        "QDASH",
        "WORCPhysicalQ3",
        "WORCNorm"
      )
    )
    
  )

# TreatmentInt - no missing
# TimePoint - no missing
# Sex2 - no missing
# AgeAtTreatment - no missing
# BilateralStatus - no missing
# WaitTime - no missing
# IndexSide - Level 2 missing
# EducationLevel_Preop - level 2 missing
# Complication2 - not missing
# QDASH - level 2 missing


# Make predictor matrix
# 
WORCPred <- make.predictorMatrix(
  WORCPROMInput2
)

WORCPred["TreatmentInt", ] <- 0
WORCPred["TimePoint", ] <- 0
WORCPred["Sex2", ] <- 0
WORCPred["AgeAtTreatment", ] <- 0
WORCPred["BilateralStatus", ] <- 0
WORCPred["WaitTime", ] <- 0
WORCPred["Complication2", ] <- 0
WORCPred["IndexSide", ] <- c(-2,0,1,1,1,1,0,1,1,1,1,1)
WORCPred["EducationLevel_Preop", ] <- c(-2,0,1,1,1,1,1,0,1,1,1,1)
WORCPred["QDASH", ] <- c(-2,1,1,1,1,1,1,1,1,0,1,1)
WORCPred["WORCPhysicalQ3", ] <- c(-2,1,1,1,1,1,1,1,1,1,0,1)
WORCPred["WORCNorm", ] <- c(-2,1,1,1,1,1,1,1,1,1,1,0)

# Specify Method
# Method specification
WORCMeth <- make.method(WORCPROMInput2)

# For time-invariant variables (level-2)
WORCMeth["EducationLevel_Preop"] <- "2lonly.pmm"  # Or "2lonly.polyreg" for categorical
WORCMeth["IndexSide"] <- "2lonly.bin"  # Or "2lonly.polyreg" for categorical


# For time-varying variables (level-1)
WORCMeth["QDASH"] <- "pmm"  # Multilevel imputation allowing for within-person variation
WORCMeth["WORCPhysicalQ3"] <- "pmm"  # Multilevel imputation allowing for within-person variation
WORCMeth["WORCNorm"] <- "pmm"  # Multilevel imputation allowing for within-person variation

# Visit sequence
# 
WORCseq <- mice::make.visitSequence(WORCPROMInput2)

WORCImputed <- mice::mice(
  data = WORCPROMInput2,
  m = 10,
  predictorMatrix = WORCPred,
  visitSequence = WORCseq,
  method = WORCMeth,
  maxit = 15,
  printFlag = FALSE
)

plot(WORCImputed)

Figure 5: Stability of imputed variables over iterations for WORC dataset

Figure 6: Stability of imputed variables over iterations for WORC dataset

2.9 Record [12.5] Analysis

A flow chart was created with the consort package (v1.2.2) (Dayim 2024) to describe the inclusion and exclusion of records into the sample pool for the present analysis to be drawn from. Patient demographics, pathology characteristics and surgical details were summarised using gtsummary (v2.4.0) (Sjoberg et al. 2021). Alpha was set for all significance tests at 5%, with confidence intervals of 95% used to bound point estimates for central tendency and model coefficients.

2.9.1 Adverse Events

The analysis of adverse events and treatment/patient survival after arthroplasty remains a challenging endeavour, made more so by the complexities of ortho-oncology. Attempts have been made to standardize reporting of adverse outcomes after rotator cuff surgery [citations]. However, a key challenge of reporting incidence rates of these outcomes in a given sample is the variability in follow up from one patient to another in the same analytical sample. With variation in follow up, the uni-dimensional estimate of incidence (number with condition/total available sample) leads to considerable underestimation of the true rate, since some cases have not yet reached sufficient follow up to experience the event of interest. For this reason time-to-event (survival) analysis provides superior incidence estimates - however, there are additional aspects of the present analysis that preclude the use of traditional Kaplan-Meier analysis.

The first is that each patient can experience multiple adverse events after the index procedure (recurring events) which adds a element of dependency to the structure of the adverse event data (thenmozhi2019?), which is not accounted for in a KM curve. The second is that certain events (e.g. implant removal) preclude the appearance of subsequent adverse events. When these records are subsequently censored (removed from the pool available records) it can bias estimates of other events of interest upward to impossible values (Coemans2022?). In the present dataset, where these elements exist simultaneously, traditional (simplistic) methods can lead to analytical decisions that remove a considerable amount of information from the dataset (e.g. analysis of first occurrence of any type) or biased estimates.

To address these issues within the analysis, the survival and tidycmprsk (v1.1.0) (Sjoberg and Fei 2024) packages were utilised to deploy a multi-state survival model (see Table 2) to estimate time-varying incidences of competing events such as;

• Implant removal (competing)

• Tendon retear | Hardware breakage

• Infection

• Adhesive capsulitis

• Dislocation - Instability

• Other events

The survival model was expressed in the form;

CRModelRCR <- survfit2(Surv(DurationStart1, DurationStop1, Category) ~ 1, 
                         data = ComplicMaster, 
                         id = CombID
                         )

SummaryRCR <- summary(CRModelRCR,
                  times = c(
                    4,
                    14,
                    38,
                    52,
                    104,
                    156
                  )
                  )

# Subset relevant information from Table1
SubsetRCR <- SummaryRCR[c("time", "pstate", "lower", "upper","states")]

SubsetRCRTrans <- lapply(SubsetRCR, t)

SummaryProbRCR <- as.data.frame(SubsetRCRTrans$pstate) %>% rename(
  T4Pr = "V1",
  T14Pr = "V2",
  T38Pr = "V3",
  T52Pr = "V4",
  T104Pr = "V5",
  T156Pr = "V6",
)

SummaryLCLRCR <- as.data.frame(SubsetRCRTrans$lower) %>% rename(
  T4LCL = "V1",
  T14LCL = "V2",
  T38LCL = "V3",
  T52LCL = "V4",
  T104LCL = "V5",
  T156LCL = "V6",
)

SummaryUCLRCR <- as.data.frame(SubsetRCRTrans$upper) %>% rename(
  T4UCL = "V1",
  T14UCL = "V2",
  T38UCL = "V3",
  T52UCL = "V4",
  T104UCL = "V5",
  T156UCL = "V6",
)

3 Analysis Results

3.1 Record [13] Participants

The initial export from the registry returned 2681 records of all types.

3.1.1 Record [13.1] Treatment selection

A flow chart of individual treatment episodes (treatments) was generated using the consort package and prepared for display with the knitr package.

The diagram below summarises recruitment and categorisation of patients into the PRULO registry.

STROBEPlot <- consort::consort_plot(
  data = STROBEFlow,
  orders = c(
    trialno = "Population",
    exclusion1 = "Excluded",
    trialno = "Received Surgery",
    exclusion2 = "Excluded from Sample",
    followup = "Study Sample at Baseline",
    lost_followup = "Unavailable 12m data",
    mitt = "Final Analysis"),
side_box = c(
  "exclusion1", 
  "exclusion2", 
  "lost_followup"
  ),
cex = 0.9
)

knitr::knit_print(STROBEPlot)

Figure 7: Flowchart of extraction and followup of sample from the Registry¹

The table below summarises patient diagnoses in the PRULO registry.

TableData <- Mastersheet2 |>
  mutate(
    ICD10 = stringr::str_extract(DiagnosisPrimary,"^[A-Za-z]+[0-9.]+")
  ) |> count(ICD10, sort = TRUE) |>
  slice_head(n = 5)

knitr::kable(TableData)

Table 5: Summary of diagnoses by ICD-10 code

ICD10	n
M75.1	151
S46.0	75
S43.43	13
S43.0	6
M24.4	2

3.2 Record [14] Patient and record characteristics

TablePatientChar <- Mastersheet2 |>
 dplyr::select(
    AgeAtTreatment, 
    Sex2, 
    IndexSide,
    #Surgeon2,
    #BMI,
    BilateralStatus,
    WaitTime,
    AccountType2
    ) |>
  tbl_summary(
    label = list(
    AgeAtTreatment ~ "Age at Surgery",
    #Surgeon2 ~ "Surgeon",
    #BMI ~ "Body Mass Index",
    BilateralStatus ~ "Bilateral",
    Sex2 ~ "Female",
    WaitTime ~ "Exam to surgery delay (weeks)",
    IndexSide ~ "Non-dominant",
    AccountType2 ~ "Insurance Type"
    ),
    type = list(
    Sex2 ~ "dichotomous",
    IndexSide ~ "dichotomous"),
    value = list(
      Sex2 ~ "Female",
      IndexSide ~ "Non-dominant"),
    statistic = list(
      AgeAtTreatment ~ "{mean} ({sd})",
      #BMI ~ "{mean} ({sd})",
      WaitTime ~ "{mean} ({sd})",
      all_categorical() ~ "{p}% ({n})"),
    missing = "no") |>
  add_n() |>
  add_ci(statistic = list(all_categorical() ~ "{conf.low} - {conf.high}",
                          all_continuous() ~ "{conf.low} - {conf.high}")) |>
   add_stat_label(
    location = "row"
  ) |> modify_table_styling(
      columns = label,
      rows = label == "DVA",
      footnote = "DVA = Department of Veterans Affairs"
    ) |> modify_table_styling(
      columns = label,
      rows = label == "TAC",
      footnote = "TAC = Transport Accident Commission"
    )

gtsummary::as_flex_table(TablePatientChar)

Table 6: Summary of patient characteristics

Characteristic	N	N = 255	95% CI
Age at Surgery, Mean (SD)	255	59 (11)	57 - 60
Female, % (n)	255	28% (71)	23 - 34
Non-dominant, % (n)	189	35% (67)	29 - 43
Bilateral, % (n)	255	9.0% (23)	5.9 - 13
Exam to surgery delay (weeks), Mean (SD)	255	19 (41)	14 - 24
Insurance Type, % (n)	255
DVA1		1.2% (3)	0.30 - 3.7
Private		71% (182)	65 - 77
TAC2		2.4% (6)	0.96 - 5.3
Uninsured		5.9% (15)	3.4 - 9.7
WorkCover		19% (49)	15 - 25
1DVA = Department of Veterans Affairs
2TAC = Transport Accident Commission
Abbreviation: CI = Confidence Interval

Patient characteristics for cases receiving the anchor of interest are summarised in Table 6.

3.2.1 Record [14.1] Pathology characteristics

TablePatientPath <- Mastersheet2 |>
  dplyr::select(
    TreatmentType,
    CuffStatus,
    CuffCondition,
    CuffTearRetraction,
    CuffTendonDelaminated,
    TendonsInvolved,
    CuffTearSizeAP, 
    CuffTearSizeML,
    TearArea, 
    TearClass, 
    CuffTearPattern,
    OtherShoulderPathology
    ) |>
  tbl_summary(
    label = list(
     TreatmentType ~ "Primary Presentation",
     CuffStatus ~ "Full Tear",
     CuffCondition ~ "Fatty Infiltration",
     CuffTearRetraction ~ "Tendon Retraction",
     CuffTendonDelaminated ~ "Tendon Delamination",
     TendonsInvolved ~ "Tendons Involved (+Supraspinatus)",
     CuffTearSizeAP ~ "Tear Size AP (mm)", 
     CuffTearSizeML ~ "Tear Size ML (mm)",
     TearArea  ~ "Tear Area (mm^2)", 
     TearClass ~ "Tear Classification", 
     CuffTearPattern  ~ "Tear Pattern",
     OtherShoulderPathology  ~ "Other Pathology"
    ),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p} ({n})"
    ),
    type = list(
     TreatmentType ~ "dichotomous",
     CuffStatus ~ "dichotomous",
     OtherShoulderPathology ~ "dichotomous"
     ),
    value = list(
      TreatmentType ~ "Primary",
      CuffStatus ~ "Full Tear",
      OtherShoulderPathology ~ "Yes"
    ),
    missing = "no"
  ) |>
  add_n() |>
  add_ci(statistic = list(
        all_categorical() ~ "{conf.low} - {conf.high}", 
        all_continuous() ~ "{conf.low} - {conf.high}")
  ) |> modify_header(
  label = "Characteristic",
  n = "Available\n Sample",
  stat_0 = "Summary\n Statistic",
  ci_stat_0 = "95% CI"
) |> add_stat_label(
    location = "row"
  ) |> modify_table_styling(
      columns = label,
      rows = var_label == "Fatty Infiltration",
      footnote = "Fuchs et al 1999"
    ) |> modify_table_styling(
      columns = label,
      rows = var_label == "Tendon Retraction",
      footnote = "Modified Patte Grading (Lädermann et al., 2016)"
    ) |> modify_table_styling(
      columns = label,
      rows = var_label == "Tear Classification",
      footnote = "(Rashid et al., 2017)"
    )

gtsummary::as_flex_table(TablePatientPath)

Table 7: Summary of presenting pathology characteristics

Characteristic	Available Sample	Summary Statistic	95% CI
Primary Presentation, % (n)	255	97 (247)	94 - 99
Full Tear, % (n)	221	92 (204)	88 - 95
Fatty Infiltration, % (n)1	221
01		56 (123)	49 - 62
11		30 (67)	24 - 37
21		12 (27)	8.3 - 17
31		1.8 (4)	0.58 - 4.9
Tendon Retraction, % (n)2	221
I2		39 (86)	33 - 46
II2		38 (84)	32 - 45
III2		15 (34)	11 - 21
IV2		6.8 (15)	4.0 - 11
No retraction2		0.9 (2)	0.16 - 3.6
Tendon Delamination, % (n)	221	56 (123)	49 - 62
Tendons Involved (+Supraspinatus), % (n)	221
Infraspinatus		14 (31)	9.9 - 19
Infraspinatus; Subscapularis		4.1 (9)	2.0 - 7.8
Infraspinatus; Teres Minor; Subscapularis		0.5 (1)	0.02 - 2.9
Subscapularis		16 (36)	12 - 22
Subscapularis (isolated)		9.0 (20)	5.8 - 14
Supraspinatus (isolated)		56 (124)	49 - 63
Tear Size AP (mm), Mean (SD)	219	23 (11)	22 - 25
Tear Size ML (mm), Mean (SD)	220	21 (9)	19 - 22
Tear Area (mm^2), Mean (SD)	219	547 (482)	483 - 611
Tear Classification, % (n)3	219
Large3		15 (33)	11 - 21
Massive3		1.4 (3)	0.35 - 4.3
Medium3		68 (148)	61 - 74
Partial3		7.3 (16)	4.4 - 12
Small3		8.7 (19)	5.4 - 13
Tear Pattern, % (n)	220
Crescent		46 (101)	39 - 53
L		21 (46)	16 - 27
Partial articular side		3.2 (7)	1.4 - 6.7
Partial bursal side		1.8 (4)	0.58 - 4.9
Reverse L		12 (27)	8.4 - 18
U		15 (33)	11 - 21
V		0.9 (2)	0.16 - 3.6
Other Pathology, % (n)	174	30 (52)	23 - 37
1Fuchs et al 1999
2Modified Patte Grading (Lädermann et al., 2016)
3(Rashid et al., 2017)
Abbreviation: CI = Confidence Interval

Pathology characteristics for cases receiving the anchor of interest are summarised in Table 7.

3.2.2 Record [14.2] Management summary

TableManage <- Mastersheet2 |>
  dplyr::select(
    ArthroscopicApproach,
    PatientPosition,
    CuffTendonsTreated,
    RepairType,
    AnchorFixation,
    RepairAugment2,
    CuffRepairTension,
    CuffRepairQuality
    ) |>
  tbl_summary(
    label = list(
     ArthroscopicApproach ~ "Arthroscopy",
     PatientPosition ~ "Beachchair Position",
     CuffTendonsTreated ~ "Supraspinatus (isolated) Repair",
     RepairType ~ "Double Row Repair",
     AnchorFixation ~ "Knotted Anchor Fixation",
     RepairAugment2 ~ "Superior Capsular Augment",
     CuffRepairTension ~ "Low Repair Tension",
     CuffRepairQuality ~ "Anatomic Repair"
    ),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p} ({n})"
    ),
    type = list(
     ArthroscopicApproach ~ "dichotomous",
     PatientPosition ~ "dichotomous",
     CuffTendonsTreated ~ "dichotomous",
     RepairType ~ "dichotomous",
     AnchorFixation ~ "dichotomous",
     RepairAugment2 ~ "dichotomous",
     CuffRepairTension ~ "dichotomous",
     CuffRepairQuality ~ "dichotomous"
     ),
    value = list(
      ArthroscopicApproach ~ "Yes",
      PatientPosition ~ "Beachchair",
      CuffTendonsTreated ~ "None",
      RepairType ~ "Double",
      AnchorFixation ~ "Knot",
      RepairAugment2 ~ "Other",
      CuffRepairTension ~ "Low",
      CuffRepairQuality ~ "Anatomic"
    ),
    missing = "no"
  ) |>
  add_n() |>
  add_ci(statistic = list(
        all_categorical() ~ "{conf.low} - {conf.high}", 
        all_continuous() ~ "{conf.low} - {conf.high}")
  ) |> 
  add_stat_label(location = "row") |>
  modify_header(
  label = "**Characteristic**",
  n = "**Available\n Sample**",
  stat_0 = "**Summary\n Statistic**",
  ci_stat_0 = "**95% CI**"
)

gtsummary::as_flex_table(TableManage)

Table 8: Summary of management and surgical details

Characteristic	**Available Sample**	**Summary Statistic**	95% CI
Arthroscopy, % (n)	221	87 (193)	82 - 91
Beachchair Position, % (n)	229	86 (197)	81 - 90
Supraspinatus (isolated) Repair, % (n)	217	59 (128)	52 - 66
Double Row Repair, % (n)	216	92 (198)	87 - 95
Knotted Anchor Fixation, % (n)	216	63 (136)	56 - 69
Superior Capsular Augment, % (n)	226	4.0 (9)	2.0 - 7.7
Low Repair Tension, % (n)	217	80 (174)	74 - 85
Anatomic Repair, % (n)	217	90 (196)	85 - 94
Abbreviation: CI = Confidence Interval

Details of surgical management are summarised in Table 8.

3.2.3 Record [14.3] Follow up

TableFollowup <- Mastersheet2 |>
  dplyr::select(
    TreatDuration,
    TreatmentStatus
    ) |>
  tbl_summary(
    by = TreatmentStatus,
    statistic = list(all_continuous() ~ "{mean} ({sd})")
  ) |> add_overall()

gtsummary::as_flex_table(TableFollowup)

Table 9: Summary of case follow up (weeks) at the time of analysis

Characteristic	Overall N = 2551	Failed N = 181	No further followup N = 21	Ongoing N = 2351
TreatDuration	94 (50)	25 (14)	43 (32)	100 (47)
1Mean (SD)

The followup of the cohort is summarised in Table 9.

FollowupSum <- Mastersheet2 |>
ggplot(aes(y = TreatmentStatus2, x = TreatDuration)) +
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
  scale_fill_brewer(palette = "Set2") +
  ggtitle('Follow up duration by treatment status')

knitr::knit_print(FollowupSum)

Figure 8: Summary of follow up duration for the included sample.

The follow up varied by the type of adverse event observed - as shown below in the Figure.

FigureFU <- ComplicMaster %>%
ggplot(aes(y = Category, x = DurationStop1)) +
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
  scale_fill_brewer(palette = "Set2") +
  xlab("Follow up (Months)")

knitr::knit_print(FigureFU)

Figure 9

3.3 Record [15] Outcomes

3.3.1 Record [15.3] Adverse events and complications

# Split out Hardware       Other          Infection      Dislocation    Mortality      Wound          Neurological   TumourLocal   Fracture       Loosening      Amputation     TumourSystemic Thrombosis     ImplantRemoval 

# Get factor levels (excluding "Censored")
# CategoryNoCensor <- ComplicTable6 |>
#   filter(Category_Value != "Censored") |>
#   pull(Category_Value) |>
#   unique()

#Create summary table
ComplicIncTable <- ComplicMasterTotal |>
  select(
    Infection,
    RepairFailure,
    ImplantRemoval,
    Capsulitis,
    Loosening,
    Neurological,
    Thrombosis,
    Reoperation,
    PainOther
  ) |> tbl_summary(
    statistic = all_categorical() ~ "{n} ({p})",
    type = all_categorical() ~ "dichotomous",
    value = list(
    Infection ~ "Infection",
    RepairFailure ~ "Repair Failure",
    ImplantRemoval ~ "Implant Removal",
    Capsulitis ~ "Capsulitis",
    Loosening ~ "Loosening",
    Neurological ~ "Neurological",
    Thrombosis ~ "Thrombosis",
    Reoperation ~ "Reoperation",
    PainOther ~ "Pain - Other"
    )
  ) |>
  add_ci(
    statistic = all_categorical() ~ "{conf.low} - {conf.high}"
  ) |>
  add_stat_label(
    label = all_categorical() ~ "n (%)",
    location = "row"
  )

as_flex_table(ComplicIncTable)

Table 10

Characteristic	N = 255	95% CI
Infection, n (%)	6 (2.4)	0.96 - 5.3
RepairFailure, n (%)	30 (12)	8.2 - 17
ImplantRemoval, n (%)	6 (2.4)	0.96 - 5.3
Capsulitis, n (%)	34 (13)	9.5 - 18
Loosening, n (%)	1 (0.4)	0.02 - 2.5
Neurological, n (%)	2 (0.8)	0.14 - 3.1
Thrombosis, n (%)	1 (0.4)	0.02 - 2.5
Reoperation, n (%)	15 (5.9)	3.4 - 9.7
PainOther, n (%)	6 (2.4)	0.96 - 5.3
Abbreviation: CI = Confidence Interval

The cohort displayed retear or failure to heal of N = 30 (12)% with (95%CI 8.2 - 17)%, as well as infection (N = 6 (2.4)%, 95%CI 0.96 - 5.3%), implant removal (N = 6 (2.4)%, 95%CI 0.96 - 5.3%) and capsulitis (N = 34 (13)%, 95%CI 9.5 - 18%).

Table 9:

ComplicGraded2 <- ComplicGraded |> left_join(
  ComplicTable6 |> dplyr::select(
    ComplicationID,
    Category_Value,
    TreatmentUID
  ) |> filter(
    ComplicationID %in% ComplicGraded$ComplicationID
  ),
  by = "ComplicationID"
) |> filter(
  !(Category_Value == "Reoperation")
)

ComplicGradeTable <- ComplicGraded2 %>% dplyr::select(
  Category_Value,
  GradeFelsch
  ) %>% tbl_summary(
    by = GradeFelsch,
    statistic = all_categorical() ~ "{n} ({p})"
    ) %>% add_stat_label(
    # update default statistic label for continuous variables
    label = all_categorical() ~ "n (%)",
    location = "row"
  ) %>% add_overall()

as_flex_table(ComplicGradeTable)

Table 11

Characteristic	Overall N = 981	1 N = 39	2 N = 29	3 N = 26	4 N = 4
Category_Value, n (%)
Capsulitis	38 (39)	29 (74)	7 (24)	2 (7.7)	0 (0)
Explant	3 (3.1)	3 (7.7)	0 (0)	0 (0)	0 (0)
ImplantRemoval	6 (6.1)	0 (0)	0 (0)	6 (23)	0 (0)
Infection	10 (10)	0 (0)	0 (0)	10 (38)	0 (0)
Loosening	1 (1.0)	0 (0)	0 (0)	1 (3.8)	0 (0)
Neurological	2 (2.0)	0 (0)	1 (3.4)	0 (0)	1 (25)
Pain - Other	7 (7.1)	5 (13)	1 (3.4)	1 (3.8)	0 (0)
RepairFailure	30 (31)	2 (5.1)	20 (69)	5 (19)	3 (75)
Thrombosis	1 (1.0)	0 (0)	0 (0)	1 (3.8)	0 (0)
1n (%)

Incidence rates were altered when viewed within the context of the multistate survival model. The cumulative incidences, when expressed at set follow up times, showed early peak incidence (<12months of surgery) for infection (Table 10). Cumulative tendon retear also peaked at 20.6% by the 3 year followup.

TableRCR <- dplyr::bind_cols(
  SummaryProbRCR,
  SummaryLCLRCR,
  SummaryUCLRCR
) |>
  # Apply gt formatting
  gt::gt(
    rownames_to_stub = TRUE
  ) |>
  gt::fmt_number(
    scale_by = 100,
    decimals = 1
  ) |>
  tab_spanner(
    label = "W4",
    columns = starts_with("T4")
  ) |>
  tab_spanner(
    label = "Wk14",
    columns = starts_with("T14")
  ) |>
  tab_spanner(
    label = "Wk38",
    columns = starts_with("T38")
  ) |>
   tab_spanner(
    label = "Wk52",
    columns = starts_with("T52")
  ) |>
   tab_spanner(
    label = "Wk104",
    columns = starts_with("T104")
  ) |>
   tab_spanner(
    label = "Wk156",
    columns = starts_with("T156")
  ) |>
  cols_label(
    contains("Pr") ~ "CumIncid",
    contains("LCL") ~ "CI Lower",
    contains("UCL") ~ "CI Upper"
  ) 

knitr::knit_print(TableRCR)

Table 12: Summary of cumulative incidences of adverse events after rotator cuff repair with suture of interest

	W4			Wk14			Wk38			Wk52			Wk104			Wk156
	CumIncid	CI Lower	CI Upper	CumIncid	CI Lower	CI Upper	CumIncid	CI Lower	CI Upper	CumIncid	CI Lower	CI Upper	CumIncid	CI Lower	CI Upper	CumIncid	CI Lower	CI Upper
(s0)	98.7	97.3	100.0	88.6	84.6	92.8	79.1	74.0	84.6	77.6	72.3	83.3	77.6	72.3	83.3	77.6	72.3	83.3
Explant	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA
Capsulitis	0.4	0.1	3.0	8.7	5.7	13.3	11.4	7.9	16.4	11.4	7.9	16.4	11.4	7.9	16.4	11.4	7.9	16.4
RepairFailure	0.0	NA	NA	0.5	0.1	3.3	5.4	3.1	9.4	6.0	3.5	10.2	6.0	3.5	10.2	6.0	3.5	10.2
Reoperation	0.0	NA	NA	0.9	0.2	3.5	0.9	0.2	3.5	1.4	0.4	4.2	1.4	0.4	4.2	1.4	0.4	4.2
ImplantRemoval	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA
Infection	0.4	0.1	3.0	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA
Thrombosis	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA
Loosening	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA	0.0	NA	NA
Neurological	0.0	NA	NA	0.4	0.1	3.1	0.9	0.2	3.7	0.9	0.2	3.7	0.9	0.2	3.7	0.9	0.2	3.7
Pain - Other	0.4	0.1	3.0	0.9	0.2	3.5	2.2	0.9	5.3	2.7	1.2	6.0	2.7	1.2	6.0	2.7	1.2	6.0

The following figure illustrates the different incidence trajectories for adverse events within each cohort, when taking into account retear and implant removal as competing risks.

Figure7 <- ggcuminc(CRModelRCR, outcome = c("Capsulitis","Infection","RepairFailure","Reoperation")) +
  #add_confidence_interval() +
  add_risktable(
    times = c(0,26,52,104,156),
    risktable_stats = "{n.risk} ({cum.event})"
  ) +
  #theme_ggsurvfit_KMunicate() +
    labs(
      x = "Follow up (weeks)",
       y = "Cumulative Percentage Incidence"
      ) +
    scale_color_discrete(labels = c("Capsulitis", "Reoperation", "Infection", "RepairFailure")) +
  scale_linetype_discrete(labels = c("Capsulitis", "Reoperation", "Infection", "RepairFailure")) +
  scale_x_continuous(
    limits = c(0, 158)) +
  scale_y_continuous(
    #limits = c(0, 0.35),
    labels = scales::percent, 
    expand = c(0.01, 0)
  ) + theme_ggdist()

knitr::knit_print(Figure7)

Figure 10

3.3.2 Record [15.4] Patient-reported outcome measures

The QuickDASH total score and WORC Normalised Index, as well as Question 3 of the Physical sub-scale of the WORC were visualised using the ggdist and ggplot2 packages. Plots were arranged using the patchwork package.

CCADASH <- ggplot(data = MasterPROM, mapping = aes(y = QDASH, x = TimePoint)) +
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "QuickDASH",
    y = "QuickDASH Total Score"
  )

CCAWORC1 <- ggplot(data = MasterPROMWORC, mapping = aes(y = WORCNorm, x = TimePoint)) +
 stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
 stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
 scale_fill_brewer(palette = "Set2") +
 ggtitle('WORC Normalised')

CCAWORC2 <- ggplot(data = MasterPROMWORC, mapping = aes(y = WORCPhysicalQ3, x = TimePoint)) +
 stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
 stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
 scale_fill_brewer(palette = "Set2") +
 ggtitle('WORC Physical Q3')

CCADASH + CCAWORC1 / CCAWORC2

Figure 11: Patient reported outcomes (QDASH, WORC) trajectories by timepoint

3.4 Record [16] Main results

The imputed datasets for QDASH and WORC were modeled with a linear mixed effects model in lme4 and summarised with broom.mixed. Up to a 38.7 point improvement in QuickDASH total score was observed (Table 11), as well as 47.1 and 54 point improvements in WORC Index Normalised and WORC Physical Question3 respectively (Table 12). Distributions of model-predicted results illustrated variability in recovery trajectories within all PROMs measures (Figure 7).

QDASHfitimp <- with(
  QDASHImputed,
  exp = lme4::lmer(
    QDASH ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2 + (1|TreatmentInt)
  )
                 )

TableQDASHfit <- gtsummary::tbl_regression(
  QDASHfitimp, 
  tidy_fun = pool_and_tidy_mice, 
  show_single_row = "Sex2",
  label = list(
    Sex2 ~ "Male vs Female",
    AgeAtTreatment ~ "Age at Surgery",
    `factor(TimePoint, ordered = FALSE)` ~ "TimePoint"),
  estimate_fun = function(x) style_number(x, digits = 2),
  pvalue_fun = function(x) style_pvalue(x, digits = 3))

gtsummary::as_flex_table(TableQDASHfit)

Table 13: Pooled linear mixed effects model for QuickDASH

Characteristic	Beta	95% CI	p-value
TimePoint
Preop	—	—
3months	-9.14	-13.96, -4.32	<0.001
6months	-24.80	-29.26, -20.34	<0.001
12months	-30.41	-34.45, -26.36	<0.001
Age at Surgery	0.03	-0.22, 0.29	0.785
Male vs Female	-3.48	-8.57, 1.61	0.176
Abbreviation: CI = Confidence Interval

TableQDASHCCA <- tbl_summary(
  QDASHPROMInput2 |> dplyr::select(
    TimePoint,
    QDASH
  ),
  by = TimePoint,
  label = list(
    QDASH ~ "QuickDASH"
  ),
  statistic = list(all_continuous() ~ "{median} ({p25} - {p75})"),
  missing = "no"
    )

gtsummary::as_flex_table(TableQDASHCCA)

Table 14: Summary of QuickDASH complete case analysis

Characteristic	Preop N = 1971	3months N = 1951	6months N = 1971	12months N = 1971
QuickDASH	45 (36 - 57)	34 (23 - 55)	18 (9 - 30)	9 (2 - 20)
1Median (Q1 - Q3)

# Fit the lm model(s)
WORCNormfit <- with(WORCImputed,exp = lmer(WORCNorm ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2 + (1|TreatmentInt)))

WORCNormfitsum <- tbl_regression(WORCNormfit, tidy_fun = pool_and_tidy_mice, show_single_row = "Sex2",
               label = list(Sex2 ~ "Male vs Female", AgeAtTreatment ~ "Age at Surgery", `factor(TimePoint, ordered = FALSE)` ~ "TimePoint"),
               estimate_fun = function(x) style_number(x, digits = 2), pvalue_fun = function(x) style_pvalue(x, digits = 3)) 

WORCQ3fit <- with(WORCImputed,exp = lmer(WORCPhysicalQ3 ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2 + (1|TreatmentInt)))

WORCQ3fitsum <- tbl_regression(WORCQ3fit, tidy_fun = pool_and_tidy_mice, show_single_row = "Sex2",
               label = list(Sex2 ~ "Male vs Female", AgeAtTreatment ~ "Age at Surgery", `factor(TimePoint, ordered = FALSE)` ~ "TimePoint"),
               estimate_fun = function(x) style_number(x, digits = 2), pvalue_fun = function(x) style_pvalue(x, digits = 3))

TableWORCfit <- tbl_merge(tbls = list(WORCNormfitsum, WORCQ3fitsum),
                        tab_spanner = c("Normalised Index", "Physical Q3")
)

gtsummary::as_flex_table(TableWORCfit)

Table 15: Pooled linear mixed effects model for WORC normalised total score and WORC Physical sub-scale Question 3

	Normalised Index			Physical Q3
Characteristic	Beta	95% CI	p-value	Beta	95% CI	p-value
TimePoint
Preop	—	—		—	—
6months	32.79	28.35, 37.23	<0.001	-35.72	-42.12, -29.32	<0.001
12months	38.38	32.61, 44.15	<0.001	-43.43	-50.14, -36.72	<0.001
Age at Surgery	0.13	-0.08, 0.35	0.222	-0.23	-0.49, 0.03	0.080
Male vs Female	3.32	-2.81, 9.45	0.281	-1.61	-9.03, 5.81	0.662
Abbreviation: CI = Confidence Interval

QDASHPredict <-  marginaleffects::predictions(QDASHfitimp)
WORCNormPredict <- marginaleffects::predictions(WORCNormfit)
WORCQ3Predict <- marginaleffects::predictions(WORCQ3fit)

TableQDASHImp <- tbl_summary(
  QDASHPredict |> dplyr::select(
    TimePoint,
    QDASH
  ),
  by = TimePoint,
  label = list(
    QDASH ~ "QuickDASH"
  ),
  statistic = list(all_continuous() ~ "{median} ({p25} - {p75})")
    )

gtsummary::as_flex_table(TableQDASHImp)

Table 16: Summary of model-predicted QuickDASH by TimePoint

Characteristic	Preop N = 1971	3months N = 1951	6months N = 1971	12months N = 1971
QuickDASH	48 (36 - 59)	32 (20 - 48)	16 (7 - 30)	11 (5 - 25)
1Median (Q1 - Q3)

TableWORCsum <- tbl_summary(
  WORCNormPredict |> dplyr::select(
    TimePoint,
    WORCNorm
  ),
  by = TimePoint,
  label = list(
    WORCNorm ~ "WORC Normalised"
  ),
  statistic = list(all_continuous() ~ "{median} ({p25} - {p75})")
    )


TableWORCQ3sum <- tbl_summary(
  WORCQ3Predict |> dplyr::select(
    TimePoint,
    WORCPhysicalQ3
  ),
  by = TimePoint,
  label = list(
    WORCPhysicalQ3 ~ "WORC Physical Q3"
  ),
  statistic = list(all_continuous() ~ "{median} ({p25} - {p75})")
    )

WORCfitsum <- tbl_stack(tbls = list(TableWORCsum, TableWORCQ3sum)
                        )

gtsummary::as_flex_table(WORCfitsum)

Table 17: Summary of model-predicted WORC Normalised Total Score and WORC Physical sub-scale Question 3 by TimePoint

Characteristic	Preop N = 1971	6months N = 1971	12months N = 1971
WORC Normalised	37 (24 - 52)	75 (58 - 87)	85 (70 - 91)
WORC Physical Q3	72 (51 - 87)	30 (16 - 54)	18 (6 - 33)
1Median (Q1 - Q3)

QDASHprPlot <- ggplot(data = QDASHPredict, mapping = aes(y = QDASH, x = TimePoint)) + stat_slabinterval(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
  scale_fill_brewer(palette = "Set2") +
  ggtitle('QuickDASH')
  

  
WORCNormprPlot <- ggplot(data =  WORCNormPredict, mapping = aes(y = WORCNorm, x = TimePoint)) +
 stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
 stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
 scale_fill_brewer(palette = "Set2") +
 ggtitle('WORC Normalised')



WORCQ3prPlot <- ggplot(data = WORCQ3Predict, mapping = aes(y = WORCPhysicalQ3, x = TimePoint)) +
 stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
 stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
 scale_fill_brewer(palette = "Set2") +
 ggtitle('WORC Physical Q3')


(QDASHprPlot) +
  (WORCNormprPlot / WORCQ3prPlot) +
  plot_layout(ncol = 2)

Figure 12: Model predicted PROMs (QuickDASH, WORC) trajectories across time points.

3.5 Record [17] Sensitivity analyses

Based on the distribution changes in QDASH and WORC over time, a sensitivity analysis was performed on the model structure using the complete case dataset. A comparison was made between quantile regression using the quantreg package and an ordinary least squares linear model from stats and a linear mixed effects model with lme4. Results were tabulated using the modelsummary package as rq models are not supported in gtsummary.

QDASHfit <- list(
  "RQ" = rq(
    QDASH ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2,
    tau = 0.5,
    method = "fn",
    data = QDASHPROMInput2
                 ),
  "LM" = lm(
    QDASH ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2,
    data = QDASHPROMInput2
  ),
  "ME" = lme4::lmer(
    QDASH ~ factor(TimePoint, ordered = FALSE) + AgeAtTreatment + Sex2 + (1|TreatmentInt),
    data = QDASHPROMInput2
    )
)

SensTablesum <- modelsummary::modelsummary(
  QDASHfit,
  statistic = c(
    "se = {std.error}",
    "conf.int"),
  fmt = fmt_decimal(1),
  output = "flextable",
  coef_rename = coef_rename
             )

knitr::knit_print(SensTablesum)

Table 18: Comparison of linear model types to assess QuickDASH by Timepoints.

	RQ	LM	ME
(Intercept)	44.3	44.9	44.5
	se = 6.8	se = 5.9	se = 7.7
	[30.9, 57.7]	[33.3, 56.6]	[29.3, 59.6]
3months	-9.6	-8.5	-9.7
	se = 3.3	se = 2.3	se = 1.8
	[-16.0, -3.2]	[-13.1, -3.9]	[-13.2, -6.3]
6months	-27.5	-24.9	-26.0
	se = 2.9	se = 2.3	se = 1.7
	[-33.1, -21.8]	[-29.4, -20.4]	[-29.4, -22.6]
12months	-36.1	-31.4	-30.5
	se = 2.4	se = 2.4	se = 1.8
	[-40.8, -31.3]	[-36.2, -26.6]	[-34.1, -26.9]
AgeAtTreatment	0.1	0.1	0.1
	se = 0.1	se = 0.1	se = 0.1
	[-0.1, 0.3]	[-0.1, 0.2]	[-0.1, 0.3]
Sex2Male	-4.8	-3.5	-3.7
	se = 2.3	se = 2.1	se = 2.9
	[-9.3, -0.2]	[-7.6, 0.6]	[-9.4, 2.0]
SD (Intercept TreatmentInt)			13.0
SD (Observations)			12.5
Num.Obs.	437	437	437
R2	0.318	0.336
R2 Adj.		0.329
R2 Marg.			0.323
R2 Cond.			0.674
AIC	3771.2	3763.5	3665.3
BIC	3795.7	3792.1	3697.9
ICC			0.5
Log.Lik.		-1874.771
RMSE	17.90	17.66	10.58

The comparison between models revealed an underestimate of the difference in 12month score to preoperative baseline of 5.6 points for the QuickDASH (15.5%) in the mixed effects linear model, compared to the quantile regression (50th percentile).

4 References

Aden-Buie, Garrick. 2023. epoxy: String Interpolation for Documents, Reports and Apps. https://CRAN.R-project.org/package=epoxy.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2025. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Arel-Bundock, Vincent. 2022. “modelsummary: Data and Model Summaries in R.” Journal of Statistical Software 103 (1): 1–23. https://doi.org/10.18637/jss.v103.i01.

Arel-Bundock, Vincent, Noah Greifer, and Andrew Heiss. 2024. “How to Interpret Statistical Models Using marginaleffects for R and Python.” Journal of Statistical Software 111 (9): 1–32. https://doi.org/10.18637/jss.v111.i09.

Barbone, Jordan Mark, and Jan Marvin Garbuszus. 2025. Openxlsx2: Read, Write and Edit “xlsx” Files. https://janmarvin.github.io/openxlsx2/.

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Benchimol, Eric I., Liam Smeeth, Astrid Guttmann, Katie Harron, David Moher, Irene Petersen, Henrik T. Sørensen, Erik von Elm, and Sinéad M. Langan. 2015. “The REporting of Studies Conducted Using Observational Routinely-Collected Health Data (RECORD) Statement.” PLOS Medicine 12 (10): e1001885. https://doi.org/10.1371/journal.pmed.1001885.

Benoit, Kenneth, David Muhr, and Kohei Watanabe. 2021. stopwords: Multilingual Stopword Lists. https://CRAN.R-project.org/package=stopwords.

Bolker, Ben, and David Robinson. 2024. broom.mixed: Tidying Methods for Mixed Models. https://CRAN.R-project.org/package=broom.mixed.

Bryan, Jennifer. 2025. Googlesheets4: Access Google Sheets Using the Sheets API V4. https://CRAN.R-project.org/package=googlesheets4.

Bryan, Jennifer, Craig Citro, and Hadley Wickham. 2025. gargle: Utilities for Working with Google APIs. https://CRAN.R-project.org/package=gargle.

Carroll, Orlagh U., Tim P. Morris, and Ruth H. Keogh. 2020. “How Are Missing Data in Covariates Handled in Observational Time-to-Event Studies in Oncology? A Systematic Review.” BMC Medical Research Methodology 20 (1). https://doi.org/10.1186/s12874-020-01018-7.

D’Agostino McGowan, Lucy, and Jennifer Bryan. 2025. googledrive: An Interface to Google Drive. https://CRAN.R-project.org/package=googledrive.

Davies, G. Matt, and Alan Gray. 2015. “Don’t Let Spurious Accusations of Pseudoreplication Limit Our Ability to Learn from Natural Experiments (and Other Messy Kinds of Ecological Monitoring).” Ecology and Evolution 5 (22): 5295–5304. https://doi.org/10.1002/ece3.1782.

Dayim, Alim. 2024. consort: Create Consort Diagram. https://CRAN.R-project.org/package=consort.

Fellows, Ian. 2018. wordcloud: Word Clouds. https://CRAN.R-project.org/package=wordcloud.

Felsch, Quinten, Victoria Mai, Holger Durchholz, Matthias Flury, Maximilian Lenz, Carl Capellen, and Laurent Audigé. 2021. “Complications Within 6 Months After Arthroscopic Rotator Cuff Repair: Registry-Based Evaluation According to a Core Event Set and Severity Grading.” Arthroscopy: The Journal of Arthroscopic & Related Surgery 37 (1): 50–58. https://doi.org/10.1016/j.arthro.2020.08.010.

Fuchs, Bruno, Dominik Weishaupt, Marco Zanetti, Juerg Hodler, and Christian Gerber. 1999. “Fatty Degeneration of the Muscles of the Rotator Cuff: Assessment by Computed Tomography Versus Magnetic Resonance Imaging.” Journal of Shoulder and Elbow Surgery 8 (6): 599–605. https://doi.org/10.1016/s1058-2746(99)90097-6.

Gohel, David, and Panagiotis Skintzos. 2025. flextable: Functions for Tabular Reporting. https://CRAN.R-project.org/package=flextable.

Grolemund, Garrett, and Hadley Wickham. 2011. “Dates and Times Made Easy with lubridate.” Journal of Statistical Software 40 (3): 1–25. https://www.jstatsoft.org/v40/i03/.

Horikoshi, Masaaki, and Yuan Tang. 2018. ggfortify: Data Visualization Tools for Statistical Analysis Results. https://CRAN.R-project.org/package=ggfortify.

Iannone, Richard, Joe Cheng, Barret Schloerke, Ellis Hughes, Alexandra Lauer, JooYoung Seo, Ken Brevoort, and Olivier Roy. 2025. gt: Easily Create Presentation-Ready Display Tables. https://CRAN.R-project.org/package=gt.

Kay, Matthew. 2024. “ggdist: Visualizations of Distributions and Uncertainty in the Grammar of Graphics.” IEEE Transactions on Visualization and Computer Graphics 30 (1): 414–24. https://doi.org/10.1109/TVCG.2023.3327195.

———. 2025. ggdist: Visualizations of Distributions and Uncertainty. https://doi.org/10.5281/zenodo.3879620.

Kirkley, Alexandra, Christine Alvarez, and Sharon Griffin. 2003. “The Development and Evaluation of a Disease-Specific Quality-of-Life Questionnaire for Disorders of the Rotator Cuff: The Western Ontario Rotator Cuff Index.” Clinical Journal of Sport Medicine 13 (2): 84–92. https://doi.org/10.1097/00042752-200303000-00004.

Koenker, Roger. 2025. quantreg: Quantile Regression. https://CRAN.R-project.org/package=quantreg.

Kuhn, Max, and Hadley Wickham. 2020. Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles. https://www.tidymodels.org.

Lädermann, Alexandre, Stephen S. Burkhart, Pierre Hoffmeyer, Lionel Neyton, Philippe Collin, Evan Yates, and Patrick J. Denard. 2016. “Classification of Full-Thickness Rotator Cuff Lesions: A Review.” EFORT Open Reviews 1 (12): 420–30. https://doi.org/10.1302/2058-5241.1.160005.

Larmarange, Joseph, and Daniel D. Sjoberg. 2025. broom.helpers: Helpers for Model Coefficients Tibbles. https://CRAN.R-project.org/package=broom.helpers.

Lazic, Stanley E. 2010. “The Problem of Pseudoreplication in Neuroscientific Studies: Is It Affecting Your Analysis?” BMC Neuroscience 11 (1). https://doi.org/10.1186/1471-2202-11-5.

Nguyen, Van Thu, Mishelle Engleton, Mauricia Davison, Philippe Ravaud, Raphael Porcher, and Isabelle Boutron. 2021. “Risk of Bias in Observational Studies Using Routinely Collected Data of Comparative Effectiveness Research: A Meta-Research Study.” BMC Medicine 19 (1). https://doi.org/10.1186/s12916-021-02151-w.

Pedersen, Thomas Lin. 2025. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.

R Core Team. 2024a. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

———. 2024b. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rashid, Mustafa S, Cushla Cooper, Jonathan Cook, David Cooper, Stephanie G Dakin, Sarah Snelling, and Andrew J Carr. 2017. “Increasing Age and Tear Size Reduce Rotator Cuff Repair Healing Rate at 1 Year.” Acta Orthopaedica 88 (6): 606–11. https://doi.org/10.1080/17453674.2017.1370844.

Rinker, Tyler W., and Dason Kurkiewicz. 2018. pacman: Package Management for R. Buffalo, New York. http://github.com/trinker/pacman.

Robinson, David, Alex Hayes, and Simon Couch. 2025. broom: Convert Statistical Objects into Tidy Tibbles. https://CRAN.R-project.org/package=broom.

Scholes, Corey, Kevin Eng, Meredith Harrison-Brown, Milad Ebrahimi, Graeme Brown, Stephen Gill, and Richard Page. 2023. “Patient Registry of Upper Limb Outcomes (PRULO) a Protocol for an Orthopaedic Clinical Quality Registry to Monitor Treatment Outcomes.” http://dx.doi.org/10.1101/2023.02.01.23284494.

Silge, Julia, and David Robinson. 2016. “tidytext: Text Mining and Analysis Using Tidy Data Principles in r.” JOSS 1 (3). https://doi.org/10.21105/joss.00037.

Sjoberg, Daniel D., Mark Baillie, Charlotta Fruechtenicht, Steven Haesendonckx, and Tim Treis. 2025. ggsurvfit: Flexible Time-to-Event Figures. https://CRAN.R-project.org/package=ggsurvfit.

Sjoberg, Daniel D., and Teng Fei. 2024. tidycmprsk: Competing Risks Estimation. https://CRAN.R-project.org/package=tidycmprsk.

Sjoberg, Daniel D., Karissa Whiting, Michael Curry, Jessica A. Lavery, and Joseph Larmarange. 2021. “Reproducible Summary Tables with the Gtsummary Package.” The R Journal 13: 570–80. https://doi.org/10.32614/RJ-2021-053.

Sjoberg, Daniel D., Abinaya Yogasekaram, and Emily de la Rua. 2025. cardx: Extra Analysis Results Data Utilities. https://CRAN.R-project.org/package=cardx.

Tang, Yuan, Masaaki Horikoshi, and Wenxuan Li. 2016. “ggfortify: Unified Interface to Visualize Statistical Result of Popular r Packages.” The R Journal 8 (2): 474–85. https://doi.org/10.32614/RJ-2016-060.

Tennant, Peter W G, Eleanor J Murray, Kellyn F Arnold, Laurie Berrie, Matthew P Fox, Sarah C Gadd, Wendy J Harrison, et al. 2020. “Use of Directed Acyclic Graphs (DAGs) to Identify Confounders in Applied Health Research: Review and Recommendations.” International Journal of Epidemiology 50 (2): 620–32. https://doi.org/10.1093/ije/dyaa213.

Terry M. Therneau, and Patricia M. Grambsch. 2000. Modeling Survival Data: Extending the Cox Model. New York: Springer.

Therneau, Terry M. 2024. A Package for Survival Analysis in r. https://CRAN.R-project.org/package=survival.

Tierney, Nicholas, and Dianne Cook. 2023. “Expanding Tidy Data Principles to Facilitate Missing Data Exploration, Visualization and Assessment of Imputations.” Journal of Statistical Software 105 (7): 1–31. https://doi.org/10.18637/jss.v105.i07.

van Buuren, Stef, and Karin Groothuis-Oudshoorn. 2011. “mice: Multivariate Imputation by Chained Equations in r.” Journal of Statistical Software 45 (3): 1–67. https://doi.org/10.18637/jss.v045.i03.

White, Ian R., Patrick Royston, and Angela M. Wood. 2010. “Multiple Imputation Using Chained Equations: Issues and Guidance for Practice.” Statistics in Medicine 30 (4): 377–99. https://doi.org/10.1002/sim.4067.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

———. 2025a. forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.

———. 2025b. stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Romain François, Lionel Henry, Kirill Müller, and Davis Vaughan. 2023. dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, Jim Hester, and Jennifer Bryan. 2024. readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2025. scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Wickham, Hadley, Davis Vaughan, and Maximilian Girlich. 2024. tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2025a. knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.

———. 2025b. litedown: A Lightweight Version of r Markdown. https://CRAN.R-project.org/package=litedown.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Footnotes

1. Non-index procedure refers to a procedure that has a preceding treatment using the hardware of interest