Pediatric Hepatitis C Patient Outcomes in a Tertiary Academic Medical Center Utilizing an Integrated Specialty Pharmacy Model: Cascade of Care

1 Study Details

1.1 Background

FDA approvals of hepatitis C medications (Harvoni, Epclusa, and Mavyret) for the pediatric population <18 years of age began in April of 2017. HSSP worked with the pediatric hepatology clinic to assist with starting treatment for these patients. HSSP pharmacists worked with patients on swallowing tablets and pellets, the insurance approval process, education of medication, on treatment monitoring, and ensuring labs completed for sustained virologic response (SVR); otherwise known as cure.

1.2 Methods

Descriptive reports of cascade of care among a pediatric population with Hepatitis C

1.2.1 Inclusion Criteria

  • Patients <18 years

  • Patients seen by a pediatric hepatologist at the Pediatric Hepatology Clinic at the Monroe Carrell Jr. Children’s Hospital at Vanderbilt

  • Clinic visits between January 2017 and September 2022.

1.3 Outcomes

Primary outcome: Rate of achievement of each individual stage of the cascade of care (CoC) following evaluation (clinic evaluation, specialty pharmacy referral, successful practice administration, medication approval, medication initiation, medication completion, SVR completed, SVR achieved)

Secondary outcomes:

  • Time from referral to the integrated specialty pharmacy to successful practice administration, to DAA approval, and to DAA initiation

  • Reasons for gaps or delays in between steps of the cascade of care

2 Analysis

2.1 Demographics

2.1.1 All patients

tab <- summaryM(age + sex + Race + height + weight + bmi +InsuranceType + referred_vsp + medication_initiation +time_to_initiation ~1,
                data=dat,continuous=5, na.include=T)
Hmisc::html(tab, 
      caption="Patient Characteristics",long=T,exclude1=F)
Patient Characteristics.
N
N=119
Age 119 3 5 8
Gender 119
    Female 0.47 (56)
    Male 0.53 (63)
Race 119
    Hispanic/Latino 0.01 ( 1)
    Black or African American 0.01 ( 1)
    White 0.73 (87)
    Unknown / Not Reported 0.24 (28)
    Other 0.01 ( 1)
    Arabic 0.01 ( 1)
Height (cm) 116 93.025 113.850 129.525
Weight (kilograms) 119 15.05 22.10 32.95
BMI 116 16.075 17.400 20.525
InsuranceType 119
    Medicaid 0.85 (101)
    Commercial 0.12 ( 14)
    None 0.03 ( 4)
Was patient referred to HSSP? 119
    Yes 0.61 (73)
    No 0.39 (46)
Medication Initiation 119
    Did not initate 0.46 (55)
    Initiated 0.54 (64)
Time from HSSP referral to daa b start
days
 64 22 43 129
a b c represent the lower quartile a, the median b, and the upper quartile c for continuous variables.   N is the number of non-missing values. Numbers after proportions are frequencies.

2.1.2 All patients grouped by FDA age groups

tab <- summaryM(age + sex + Race + height + weight + bmi +InsuranceType +  genotype+
                  referred_vsp + medication_initiation + Regimen + fibrosis_score  + viral_load_b +dosage_form + daa_lot +
                  treatment_naive  + ast_b + alt_b ~ age_group ,
                data=dat,continuous=5, na.include=T, overall=T)
Hmisc::html(tab, 
      caption="Patient Characteristics",long=T,exclude1=F, digits=3)
Patient Characteristics.
N
 0-5 years
N=61		 6-11 years
N=42		 12-17 years
N=16		 Combined
N=119
Age 119 1.0 3.0 5.0 6.0 7.0 9.0 14.0 15.0 16.2 3.0 5.0 8.0
Gender 119
    Female 0.30 (18) 0.60 (25) 0.81 (13) 0.47 (56)
    Male 0.70 (43) 0.40 (17) 0.19 ( 3) 0.53 (63)
Race 119
    Hispanic/Latino 0.00 ( 0) 0.00 ( 0) 0.06 ( 1) 0.01 ( 1)
    Black or African American 0.02 ( 1) 0.00 ( 0) 0.00 ( 0) 0.01 ( 1)
    White 0.70 (43) 0.76 (32) 0.75 (12) 0.73 (87)
    Unknown / Not Reported 0.28 (17) 0.24 (10) 0.06 ( 1) 0.24 (28)
    Other 0.00 ( 0) 0.00 ( 0) 0.06 ( 1) 0.01 ( 1)
    Arabic 0.00 ( 0) 0.00 ( 0) 0.06 ( 1) 0.01 ( 1)
Height (cm) 116 76.2 94.0 106.6 118.3 125.3 134.9 148.8 157.0 162.1 93.0 113.8 129.5
Weight (kilograms) 119 9.8 15.1 18.9 24.7 27.5 33.4 49.8 58.5 74.0 15.1 22.1 33.0
BMI 116 15.6 16.6 17.8 16.6 17.6 20.5 20.3 22.8 26.6 16.1 17.4 20.5
InsuranceType 119
    Medicaid 0.87 ( 53) 0.83 ( 35) 0.81 ( 13) 0.85 (101)
    Commercial 0.08 ( 5) 0.14 ( 6) 0.19 ( 3) 0.12 ( 14)
    None 0.05 ( 3) 0.02 ( 1) 0.00 ( 0) 0.03 ( 4)
Genotype 119
    Genotype 1 0.34 (21) 0.57 (24) 0.69 (11) 0.47 (56)
    Genotype 2 0.02 ( 1) 0.05 ( 2) 0.06 ( 1) 0.03 ( 4)
    Genotype 3 0.05 ( 3) 0.19 ( 8) 0.06 ( 1) 0.10 (12)
    Genotype 4 0.00 ( 0) 0.00 ( 0) 0.06 ( 1) 0.01 ( 1)
    NA 0.59 (36) 0.19 ( 8) 0.12 ( 2) 0.39 (46)
Was patient referred to HSSP? 119
    Yes 0.43 (26) 0.79 (33) 0.88 (14) 0.61 (73)
    No 0.57 (35) 0.21 ( 9) 0.12 ( 2) 0.39 (46)
Medication Initiation 119
    Did not initate 0.67 (41) 0.24 (10) 0.25 ( 4) 0.46 (55)
    Initiated 0.33 (20) 0.76 (32) 0.75 (12) 0.54 (64)
Regimen 119
    Epclusa 200/50mg Pellets x 12 weeks 0.02 ( 1) 0.00 ( 0) 0.00 ( 0) 0.01 ( 1)
    Epclusa 200/50mg Tablets x 12 weeks 0.05 ( 3) 0.17 ( 7) 0.00 ( 0) 0.08 (10)
    Epclusa 400/100mg Pellets x 12 weeks 0.00 ( 0) 0.05 ( 2) 0.00 ( 0) 0.02 ( 2)
    Epclusa 400/100mg Tablets x 12 weeks 0.00 ( 0) 0.10 ( 4) 0.00 ( 0) 0.03 ( 4)
    Harvoni 33.75/150mg Pellets x 12 weeks 0.05 ( 3) 0.00 ( 0) 0.00 ( 0) 0.03 ( 3)
    Harvoni 45/200mg Pellets x 12 weeks 0.05 ( 3) 0.02 ( 1) 0.00 ( 0) 0.03 ( 4)
    Harvoni 45/200mg Tablets x 12 weeks 0.18 (11) 0.36 (15) 0.00 ( 0) 0.22 (26)
    Harvoni 90/400mg Tablets x 12 weeks 0.00 ( 0) 0.07 ( 3) 0.50 ( 8) 0.09 (11)
    Mavyret 300/120mg Tablets x 8 weeks 0.00 ( 0) 0.00 ( 0) 0.25 ( 4) 0.03 ( 4)
    NA 0.66 (40) 0.24 (10) 0.25 ( 4) 0.45 (54)
Fibrosis Score 119
    F0 0.05 ( 3) 0.07 ( 3) 0.25 ( 4) 0.08 (10)
    F1 0.00 ( 0) 0.02 ( 1) 0.12 ( 2) 0.03 ( 3)
    F4 0.00 ( 0) 0.02 ( 1) 0.00 ( 0) 0.01 ( 1)
    F2-3 0.00 ( 0) 0.05 ( 2) 0.00 ( 0) 0.02 ( 2)
    no score available and no cirrhosis 0.34 (21) 0.64 (27) 0.50 ( 8) 0.47 (56)
    NA 0.61 (37) 0.19 ( 8) 0.12 ( 2) 0.39 (47)
Baseline Viral Load (IU/mL) 74 230500 651144 2993814 432619 1396356 3242928 285110 554147 1661144 329576 826275 3083763
Tablets or Pellets 119
    Tablets 0.23 (14) 0.69 (29) 0.75 (12) 0.46 (55)
    Pellets 0.11 ( 7) 0.07 ( 3) 0.00 ( 0) 0.08 (10)
    NA 0.66 (40) 0.24 (10) 0.25 ( 4) 0.45 (54)
Length of DAA Regimen 119
    8 weeks 0.00 ( 0) 0.00 ( 0) 0.25 ( 4) 0.03 ( 4)
    12 weeks 0.34 (21) 0.76 (32) 0.50 ( 8) 0.51 (61)
    NA 0.66 (40) 0.24 (10) 0.25 ( 4) 0.45 (54)
Treatment naïve 119
    Yes 0.43 (26) 0.81 (34) 0.88 (14) 0.62 (74)
    NA 0.57 (35) 0.19 ( 8) 0.12 ( 2) 0.38 (45)
AST (U/L) 74 45.0 50.5 62.8 34.0 44.0 57.2 28.8 37.0 47.5 35.2 46.5 61.2
ALT (U/L) 74 39.0 51.5 79.0 34.2 49.0 65.0 29.2 48.0 64.2 33.2 49.5 74.2
a b c represent the lower quartile a, the median b, and the upper quartile c for continuous variables.   N is the number of non-missing values. Numbers after proportions are frequencies.

2.2 Rate of achievment in each cascade of care

2.2.1 Completion of Cascade of Care

dat$clinic_visit <- ymd(dat$clinic_visit)
dat$vsp_referral_date <- ymd(dat$vsp_referral_date)
dat$date_successful_practice <- ymd(dat$date_successful_practice)
dat$date_medication_approval <- ifelse(dat$pa_outcome == "Approved", dat$pa_outcome_date,
                                ifelse(dat$pa_outcome %in% c("Denied", "Partial approval") & dat$appeal_outcome == "Approved", dat$appeal_outcome_date, NA))
dat$date_medication_approval[dat$study_id %in% c(8, 3475)] <- dat$appeal_outcome_date[dat$study_id %in% c(8, 3475)]
dat$date_medication_approval[dat$study_id %in% c(3438)] <- dat$pa_outcome_date[dat$study_id %in% c(3438)]

#id 8, 3438, and 3475 were given approval but did not have appeal approved
dat$date_medication_approval <-  ymd(dat$date_medication_approval)                                  
dat$daa_b_start_date <- ymd(dat$daa_b_start_date)
dat$end_date <- ymd(dat$end_date)
dat$svr_viral_load_date <- ymd(dat$svr_viral_load_date)
dat[dat$swallowing == "Yes" & !is.na(dat$swallowing), "date_successful_practice"] <- dat[dat$swallowing == "Yes"& !is.na(dat$swallowing), "vsp_referral_date"]
dat$`Complete Clinic Visit` <- ifelse(!is.na(dat$clinic_visit), "Yes" ,"No")
dat$`Complete VSP Referral` <- ifelse(!is.na(dat$vsp_referral_date), "Yes" ,"No")
dat$`Complete Successful Practice` <- ifelse(!is.na(dat$date_successful_practice), "Yes" ,"No")
dat$`Complete Medication Initiation` <- ifelse(!is.na(dat$daa_b_start_date), "Yes" ,"No")
dat$`Complete Medication End` <- ifelse(!is.na(dat$end_date), "Yes" ,"No")
dat$`Complete SVR Viral Load` <- ifelse(!is.na(dat$svr_viral_load_date), "Yes" ,"No")
dat$insurance_approval_date <- ifelse(dat$pa_outcome == "Approved", dat$pa_outcome_date, 
                               ifelse((dat$pa_outcome == "Partial approval" | 
                                         dat$pa_outcome == "Denied") &  
                                        (dat$appeal_outcome == "Approved"|
                                        dat$appeal_outcome == "Cancelled"), dat$appeal_outcome_date, NA))
dat$insurance_approval_date[dat$study_id == 3438] <- "2020-2-21"
dat$insurance_approval_date[dat$study_id == 4006] <- "2020-05-12"
#dat %>% filter(study_id %in%  c(4006,8, 3438, 3475)) %>% select(1,pa_outcome, pa_required, pa_outcome_date, appeal_outcome, appeal_outcome_date, insurance_approval_date)
dat$"Complete Insurance Approval" <- ifelse(!is.na(dat$insurance_approval_date), "Yes", "No")
dat$"Complete Achieved SVR" <- ifelse(dat$svr_achieved == "Yes", "Yes", "No")
dat$`Complete Achieved SVR`[is.na(dat$`Complete Achieved SVR`)] <- "No"
tab <- summaryM(`Complete Clinic Visit`+
                 `Complete VSP Referral` +`Complete Successful Practice`+ `Complete Insurance Approval`+
                  `Complete Medication Initiation` +`Complete Medication End`+ 
                  `Complete SVR Viral Load`+`Complete Achieved SVR`~1,
                data=dat,continuous=5, na.include=T)
Hmisc::html(tab, 
      caption="Completion rates of CoC",long=T,exclude1=F)
Completion rates of CoC.

N=119
Complete Clinic Visit
    Yes 1 (119)
    NO 0 ( 0)
Complete VSP Referral
    No 0.39 (46)
    Yes 0.61 (73)
Complete Successful Practice
    No 0.42 (50)
    Yes 0.58 (69)
Complete Insurance Approval
    No 0.44 (52)
    Yes 0.56 (67)
Complete Medication Initiation
    No 0.46 (55)
    Yes 0.54 (64)
Complete Medication End
    No 0.46 (55)
    Yes 0.54 (64)
Complete SVR Viral Load
    No 0.49 (58)
    Yes 0.51 (61)
Complete Achieved SVR
    No 0.5 (60)
    Yes 0.5 (59)
Numbers after proportions are frequencies.

2.2.2 Modeling time from referral to medication initiation

  • The model outcome is time to initiation from VSP referral in days

  • Only the 64 patients with a medication initiation are included in the model

  • An ordinal logistic regression model was fit due to the strong skew of the outcome

2.2.2.1 Distribution of the outcome

ggplot(dat=dat, aes(x=time_to_initiation)) +
  geom_histogram(color="Steelblue4", fill="steelblue3") +
  xlab("Time to initiation from VSP referral (Days)") + ylab("Count") +
  ggthemes::theme_calc()

2.2.2.2 Model Summary

2.2.2.2.1 Without vsp availability
mod_dat <- dat[!is.na(dat$daa_b_start_date),]
mod_dat$time_to_initiation <- as.numeric(ymd(mod_dat$daa_b_start_date)  - ymd(mod_dat$vsp_referral_date))
mod_dat$pa_outcome <- droplevels(mod_dat$pa_outcome)
dd <- datadist(mod_dat %>% select(time_to_initiation,swallowing,vsp_available,pa_outcome))
options(datadist='dd')
fit <- orm(time_to_initiation~swallowing+pa_outcome, data=mod_dat)


res <- data.frame(Variable = grep('y>',names(coef(fit)),inv=T, v = T), 
             lor = coef(fit)[grep('y>',names(coef(fit)),inv=T, v = T)], 
             rse =sqrt(diag(vcov(fit)))[grep('y>',names(coef(fit)),inv=T, v = T)]) %>% 
    mutate(lower = (lor - 1.96*rse), upper = (lor + 1.96*rse), 
           `p-value` = format.pval( round(pchisq((lor/rse)^2,1, lower.tail = FALSE),4), digits = 3, eps = 0.001)) %>% 
    mutate(`Log Odds Ratio` = lor, `Odds Ratio` = exp(lor), Lower = exp(lower), Upper = exp(upper)) %>% 
    select(Variable, `Log Odds Ratio`,`Odds Ratio`, Lower, Upper, `p-value`)
res$Variable <- factor(res$Variable,
  levels=c("swallowing=No",
           "vsp_available=No",
           "pa_outcome=Denied",
           "pa_outcome=Partial approval"),
  labels=c("No",
           "Unavailable",
           "Denied",
           "Partially approved")
)
rownames(res) <- NULL
res %>%
  arrange(Variable) %>%
  kable(digits=2) %>% 
  kable_paper("hover", full_width = F)%>%
  pack_rows("Swallowing (Reference: Yes)", 1, 1) %>%
  pack_rows("PA Outcome (Reference: Approved)", 2, 3) %>%
  add_header_above(c(" " = 3, "95% Confidence interval" = 2, " " = 1)) %>%
  kable_classic(full_width = F, html_font = "Cambria")
95% Confidence interval
Variable Log Odds Ratio Odds Ratio Lower Upper p-value
Swallowing (Reference: Yes)
No 1.36 3.89 1.58 9.61 0.003
PA Outcome (Reference: Approved)
Denied 0.88 2.42 0.77 7.64 0.132
Partially approved -0.32 0.72 0.22 2.35 0.591 
options(prType="html")
print(fit)
Logistic (Proportional Odds) Ordinal Regression Model 
orm(formula = time_to_initiation ~ swallowing + pa_outcome, data = mod_dat)
Model Likelihood
Ratio Test		 Discrimination
Indexes		 Rank Discrim.
Indexes
Obs 64 LR χ2 11.78 R2 0.168 ρ 0.422
Distinct Y 50 d.f. 3 R23,64 0.128
Y0.5 43 Pr(>χ2) 0.0082 R23,63.9 0.128
max |∂log L/∂β| 6×10-5 Score χ2 11.92 |Pr(Y ≥ median)-½| 0.166
Pr(>χ2) 0.0077
β S.E. Wald Z Pr(>|Z|)
swallowing=No   1.3588  0.4615 2.94 0.0032
pa_outcome=Denied   0.8841  0.5862 1.51 0.1315
pa_outcome=Partial approval  -0.3220  0.5995 -0.54 0.5912 
print(anova(fit))
Wald Statistics for time_to_initiation
χ2 d.f. P
swallowing 8.67 1 0.0032
pa_outcome 2.98 2 0.2256
TOTAL 11.15 3 0.0109
2.2.2.2.2 With vsp availability
mod_dat <- dat[!is.na(dat$daa_b_start_date),]
mod_dat$time_to_initiation <- as.numeric(ymd(mod_dat$daa_b_start_date)  - ymd(mod_dat$vsp_referral_date))
mod_dat$pa_outcome <- droplevels(mod_dat$pa_outcome)
dd <- datadist(mod_dat %>% select(time_to_initiation,swallowing,vsp_available,pa_outcome))
options(datadist='dd')
fit <- orm(time_to_initiation~swallowing+vsp_available+pa_outcome, data=mod_dat)


res <- data.frame(Variable = grep('y>',names(coef(fit)),inv=T, v = T), 
             lor = coef(fit)[grep('y>',names(coef(fit)),inv=T, v = T)], 
             rse =sqrt(diag(vcov(fit)))[grep('y>',names(coef(fit)),inv=T, v = T)]) %>% 
    mutate(lower = (lor - 1.96*rse), upper = (lor + 1.96*rse), 
           `p-value` = format.pval( round(pchisq((lor/rse)^2,1, lower.tail = FALSE),4), digits = 3, eps = 0.001)) %>% 
    mutate(`Log Odds Ratio` = lor, `Odds Ratio` = exp(lor), Lower = exp(lower), Upper = exp(upper)) %>% 
    select(Variable, `Log Odds Ratio`,`Odds Ratio`, Lower, Upper, `p-value`)
res$Variable <- factor(res$Variable,
  levels=c("swallowing=No",
           "vsp_available=No",
           "pa_outcome=Denied",
           "pa_outcome=Partial approval"),
  labels=c("No",
           "Unavailable",
           "Denied",
           "Partially approved")
)
rownames(res) <- NULL
res %>%
  arrange(Variable) %>%
  kable(digits=2) %>% 
  kable_paper("hover", full_width = F)%>%
  pack_rows("Swallowing (Reference: Yes)", 1, 1) %>%
  pack_rows("VSP Availability (Reference: Available)", 2, 2) %>%
  pack_rows("PA Outcome (Reference: Approved)", 3, 4) %>%
  add_header_above(c(" " = 3, "95% Confidence interval" = 2, " " = 1)) %>%
  kable_classic(full_width = F, html_font = "Cambria")
95% Confidence interval
Variable Log Odds Ratio Odds Ratio Lower Upper p-value
Swallowing (Reference: Yes)
No 1.37 3.94 1.56 9.98 0.004
VSP Availability (Reference: Available)
Unavailable 3.73 41.47 9.51 180.87 <0.001
PA Outcome (Reference: Approved)
Denied 0.99 2.68 0.80 9.03 0.112
Partially approved 0.06 1.06 0.33 3.42 0.918 
options(prType="html")
print(fit)
Logistic (Proportional Odds) Ordinal Regression Model 
orm(formula = time_to_initiation ~ swallowing + vsp_available + 
    pa_outcome, data = mod_dat)
Model Likelihood
Ratio Test		 Discrimination
Indexes		 Rank Discrim.
Indexes
Obs 64 LR χ2 42.33 R2 0.484 ρ 0.604
Distinct Y 50 d.f. 4 R24,64 0.451
Y0.5 43 Pr(>χ2) <0.0001 R24,63.9 0.451
max |∂log L/∂β| 0.0005 Score χ2 48.67 |Pr(Y ≥ median)-½| 0.228
Pr(>χ2) <0.0001
β S.E. Wald Z Pr(>|Z|)
swallowing=No  1.3721  0.4739 2.90 0.0038
vsp_available=No  3.7251  0.7514 4.96 <0.0001
pa_outcome=Denied  0.9862  0.6197 1.59 0.1115
pa_outcome=Partial approval  0.0614  0.5958 0.10 0.9179 
print(anova(fit))
Wald Statistics for time_to_initiation
χ2 d.f. P
swallowing 8.38 1 0.0038
vsp_available 24.58 1 <0.0001
pa_outcome 2.57 2 0.2760
TOTAL 33.42 4 <0.0001

2.2.3 Median time between steps in cascade of care

require(gtsummary)
dat$time_referral_admin <- dat$date_successful_practice - dat$vsp_referral_date 
#View(dat %>% select(time_referall_admin, vsp_referral_date, date_successful_practice))
dat$insurance_approval_date <- ymd(dat$insurance_approval_date)
dat$time_admin_insurance <- dat$insurance_approval_date - dat$date_successful_practice
dat$time_insurance_initiation <- dat$daa_b_start_date - dat$insurance_approval_date
dat$time_referral_initiation <- dat$daa_b_start_date - dat$vsp_referral_date

#dat %>% select(1,date_successful_practice, insurance_approval_date,time_admin_insurance)
dat %>% 
  select(time_referral_admin, time_admin_insurance, time_insurance_initiation, time_referral_initiation) %>% 
  tbl_summary(type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{N_nonmiss}",
      "{median} ({p25}, {p75})",
      "{min}, {max}"
    ))
Characteristic N = 119
time_referral_admin
    N 69
    Median (IQR) 0 (0, 51)
    Range 0, 676
    Unknown 50
time_admin_insurance
    N 67
    Median (IQR) 15 (5, 52)
    Range -14, 283
    Unknown 52
time_insurance_initiation
    N 64
    Median (IQR) 11 (7, 20)
    Range -42, 59
    Unknown 55
time_referral_initiation
    N 64
    Median (IQR) 43 (22, 129)
    Range 10, 710
    Unknown 55

2.2.4 Visualization of Cascade of care

label(dat$clinic_visit) <- NULL
label(dat$vsp_referral_date) <- NULL
label(dat$date_successful_practice) <- NULL
label(dat$daa_b_start_date) <- NULL
label(dat$end_date) <- NULL
label(dat$svr_viral_load_date) <- NULL
stages <- c("Complete Clinic Visit",
                 "Complete VSP Referral","Complete Successful Practice", "Complete Insurance Approval",
                  "Complete Medication Initiation" ,"Complete Medication End",
                  "Complete SVR Viral Load", "Complete Achieved SVR")
plot_dat <- dat[,stages]
plot_dat <- sapply(plot_dat, function(x){ifelse(x == "Yes", 1, 0)})
stage_sum <- colSums(plot_dat)
plot_dat <- data.frame(Stage=names(stage_sum),
                       Count = unname(stage_sum))
plot_dat$Stage <- factor(plot_dat$Stage,
                         levels=c("Complete Clinic Visit",
                                  "Complete VSP Referral",
                                  "Complete Successful Practice",
                                  "Complete Insurance Approval",
                                  "Complete Medication Initiation" ,
                                  "Complete Medication End",
                                  "Complete SVR Viral Load",
                                  "Complete Achieved SVR"),
                         labels=c("Clinic Visit",
                                  "VSP Referral",
                                  "Successful Practice",
                                  "Insurance Approval",
                                  "Medication Initiation" ,
                                  "Medication End",
                                  "SVR Viral Load",
                                  "Achieved SVR"))
my.labels <- c(                   "Clinic \nEvaluation",
                                  "Referral to HSSP",
                                  "Successful \nAdministration",
                                  "Insurance \nApproval",
                                  "Treatment \nInitiation" ,
                                  "Treatment \nCompletion",
                                  "SVR Lab \nCompletion",
                                  "Achieved SVR")
plot_dat$Percentage <- plot_dat$Count/max(plot_dat$Count) 
ggplot(plot_dat) +
  geom_bar(aes(y=Percentage, x=Stage), stat= "identity", alpha=1, color="black", fill="#0C3D77") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_discrete(labels= my.labels)+
  xlab("Stage of Care") + 
  geom_text(aes(label=paste0("n= ",Count),x=Stage, y=Percentage+0.025), size=4) +
  #ylim(c(0,1.10)) +
  ggthemes::theme_calc()

ggsave("COC_overall.png", dpi=300)

#ix <- which(dat$`Complete Medication Initiation` == "Yes" & dat$`Complete Insurance Approval` == "No") 
#id <- dat$study_id[ix]
#dat %>% filter(study_id %in% id) %>% select(study_id ,pa_outcome, appeal_outcome, pa_outcome_date, appeal_outcome_date, insurance_approval_date, daa_b_start_date) %>% View()
stages <- c("Complete Clinic Visit",
                                  "Complete VSP Referral",
                                  "Complete Successful Practice",
                                  "Complete Insurance Approval",
                                  "Complete Medication Initiation" ,
                                  "Complete Medication End",
                                  "Complete SVR Viral Load",
                                  "Complete Achieved SVR")


####################### 0-5

plot_dat <- dat[dat$age_group == "0-5 years",stages]
plot_dat <- sapply(plot_dat, function(x){ifelse(x == "Yes", 1, 0)})
stage_sum <- colSums(plot_dat)
plot_dat <- data.frame(Stage=names(stage_sum),
                       Count = unname(stage_sum))
plot_dat$Percentage <- plot_dat$Count/max(plot_dat$Count) 

plot_dat$Stage <- factor(plot_dat$Stage,
                         levels=c("Complete Clinic Visit",
                                  "Complete VSP Referral",
                                  "Complete Successful Practice",
                                  "Complete Insurance Approval",
                                  "Complete Medication Initiation" ,
                                  "Complete Medication End",
                                  "Complete SVR Viral Load",
                                  "Complete Achieved SVR"),
                         labels=c("Clinic Visit",
                                  "VSP Referral",
                                  "Successful Practice",
                                  "Insurance Approval",
                                  "Medication Initiation" ,
                                  "Medication End",
                                  "SVR Viral Load",
                                  "Achieved SVR"))
plot_dat1 <- plot_dat
#######################
####################### 6-11

plot_dat <- dat[dat$age_group == "6-11 years",stages]
plot_dat <- sapply(plot_dat, function(x){ifelse(x == "Yes", 1, 0)})
stage_sum <- colSums(plot_dat)
plot_dat <- data.frame(Stage=names(stage_sum),
                       Count = unname(stage_sum))
plot_dat$Percentage <- plot_dat$Count/max(plot_dat$Count) 

plot_dat$Stage <- factor(plot_dat$Stage,
                         levels=c("Complete Clinic Visit",
                                  "Complete VSP Referral",
                                  "Complete Successful Practice",
                                  "Complete Insurance Approval",
                                  "Complete Medication Initiation" ,
                                  "Complete Medication End",
                                  "Complete SVR Viral Load",
                                  "Complete Achieved SVR"),
                         labels=c("Clinic Visit",
                                  "VSP Referral",
                                  "Successful Practice",
                                  "Insurance Approval",
                                  "Medication Initiation" ,
                                  "Medication End",
                                  "SVR Viral Load",
                                  "Achieved SVR"))
plot_dat2 <- plot_dat
#######################
####################### 12-17

plot_dat <- dat[dat$age_group == "12-17 years",stages]
plot_dat <- sapply(plot_dat, function(x){ifelse(x == "Yes", 1, 0)})
stage_sum <- colSums(plot_dat)
plot_dat <- data.frame(Stage=names(stage_sum),
                       Count = unname(stage_sum))
plot_dat$Percentage <- plot_dat$Count/max(plot_dat$Count) 

plot_dat$Stage <- factor(plot_dat$Stage,
                         levels=c("Complete Clinic Visit",
                                  "Complete VSP Referral",
                                  "Complete Successful Practice",
                                  "Complete Insurance Approval",
                                  "Complete Medication Initiation" ,
                                  "Complete Medication End",
                                  "Complete SVR Viral Load",
                                  "Complete Achieved SVR"),
                         labels=c("Clinic Visit",
                                  "VSP Referral",
                                  "Successful Practice",
                                  "Insurance Approval",
                                  "Medication Initiation" ,
                                  "Medication End",
                                  "SVR Viral Load",
                                  "Achieved SVR"))
plot_dat3 <- plot_dat
#######################
my.labels <- c(                   "Clinic \nEvaluation",
                                  "Referral \nto \nHSSP",
                                  "Successful \nAdministration",
                                  "Insurance \nApproval",
                                  "Treatment \nInitiation" ,
                                  "Treatment \nCompletion",
                                  "SVR \nLab \nCompletion",
                                  "Achieved \nSVR")


g1 <- ggplot(plot_dat1) +
  geom_bar(aes(y=Percentage, x=Stage), stat= "identity", 
           alpha=1, color="black",  fill="#0C3D77") +
  scale_x_discrete(labels= my.labels)+
  xlab("Stage of Care") +
    ggtitle("0-5 years")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +

  ggthemes::theme_calc()+
 theme(text = element_text(size=20))
  
g2 <- ggplot(plot_dat2) +
  geom_bar(aes(y=Percentage, x=Stage), stat= "identity", 
           alpha=1, color="black",  fill="#0C3D77") +
  scale_x_discrete(labels= my.labels)+
  xlab("Stage of Care") +
    ggtitle("6-11 years")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +

  ggthemes::theme_calc()+
  theme(text = element_text(size=20))
  
g3 <- ggplot(plot_dat3) +
  geom_bar(aes(y=Percentage, x=Stage), stat= "identity", 
           alpha=1, color="black",  fill="#0C3D77") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_discrete(labels= my.labels)+
  xlab("Stage of Care") +
  ggtitle("12-17 years")+
    #ylim(c(0,1))+

  ggthemes::theme_calc()+
  theme(text = element_text(size=20))
        #axis.text.x = element_text(angle=35, hjust=0.5)) 
cowplot::plot_grid(g1,g2,g3, cols=3)

d <- 2.6
ggsave("COC_age.png", dpi=100, width=12*d, height=4*d, limitsize=F)

2.2.5 Referral information

tab <- summaryM(referred_vsp + vsp_not_referred + start +
                  reason_for_not_start~1,
                data=dat,continuous=5, na.include=T)
Hmisc::html(tab, 
      caption="Referral Information",long=T,exclude1=F)
Referral Information.

N=119
Was patient referred to HSSP?
    Yes 0.61 (73)
    No 0.39 (46)
Why was patient not referred to VSP?
    Lost to follow-up in clinic 0.18 (22)
    Will be followed by local provider 0.02 ( 2)
    Enrolled in HCV study 0.03 ( 3)
    HCV RNA negative x2 readings 0.07 ( 8)
    HCV Ab- after 18 months 0.05 ( 6)
    < 3 yo, will continue to follow until DAA is appropriate if applicable 0.02 ( 2)
    < 18 months, will continue to follow to see if HCV RNA+ after 18 months of age 0.02 ( 2)
    Prescriber monitoring swallowing (not yet sent to VIP and not able to swallow) 0.01 ( 1)
    NA 0.61 (73)
Did Patient Start Therapy at VUMC?
    Yes 0.54 (64)
    No 0.08 ( 9)
    NA 0.39 (46)
Reason Not Started
    Patient unable to administer prescribed dosage form 0.03 ( 4)
    Patient lost to follow-up 0.03 ( 3)
    Other 0.01 ( 1)
    Started treatment with outside provider and will be followed at OSF 0.01 ( 1)
    NA 0.92 (110)
Numbers after proportions are frequencies.