Load Libarys

library(tidyverse)
library(plotly)
library(kableExtra)
library(janitor)

Load Data

base.reba=read.csv('data/trinity_reba_rroi_baseline.csv')
rf=read.csv('data/trinity_reba_rroi_rf.csv')
cfl=read.csv('data/trinity_reba_rroi_flp.csv')
rf.cfl=read.csv('data/trinity_reba_rroi_both.csv')

Experiment 1: Response function alterations

Step 1: Check for ROI differences

Experiment bounds step 1: filter for ast_de0 and treatment units common across all scenarios

base.reba=base.reba%>%
  filter(ast_de0>0)

rf=rf%>%
  filter(ast_de0>0)

cfl=cfl%>%
  filter(ast_de0>0)

rf.cfl=rf.cfl%>%
  filter(ast_de0>0)

Optional Experiment bounds step 2: filter for 5 treatments to test

base.reba=base.reba%>%
  filter(ini_rx %in% c('GFL', 'TA', 'TH', 'MASL', 'CC'))
rf=rf%>%
  filter(ini_rx %in% c('GFL', 'TA', 'TH', 'MASL', 'CC'))
cfl=cfl%>%
  filter(ini_rx %in% c('GFL', 'TA', 'TH', 'MASL', 'CC'))
rf.cfl=rf.cfl%>%
  filter(ini_rx %in% c('GFL', 'TA', 'TH', 'MASL', 'CC'))
fig1= plot_ly(alpha=.6)
fig1=fig1%>%
  add_histogram(data=base.reba, type='histogram', x=base.reba$ast_de0, name='base')
fig1=fig1%>%
    add_histogram(data=rf, type='histogram', x=rf$ast_de0, name='RF edits')
fig1 <- fig1 %>% layout(barmode = "overlay")
fig1

Step 2: Eval by treatment type assigned

base.trt=base.reba%>%
  count(ini_rx)%>%
  arrange(desc(n))

rf.trt=rf%>%
  count(ini_rx)%>%
  arrange(desc(n))
trt.select=left_join(base.trt, rf.trt, by="ini_rx")
trt.select=trt.select%>%
  mutate(delta=n.y-n.x)
colnames(trt.select)=c('ini_rx', 'base.reba', 'RF.Edits', 'delta')
trt.select.t=trt.select%>%
  adorn_totals("row")

tbl=kbl(trt.select.t, align='c', booktabs=TRUE)%>%
  column_spec(4, background=ifelse(trt.select$delta<0, 'coral', 'white'))
## Warning in ensure_len_html(background, nrows, "background"): The number of
## provided values in background does not equal to the number of rows.
tbl%>%
  kable_styling(full_width = TRUE)
ini_rx base.reba RF.Edits delta
TH 1297 1286 -11
GFL 749 746 -3
CC 476 470 -6
VDT 409 404 -5
SD 401 396 -5
MASL 338 336 -2
BR 123 124 1
AI 73 71 -2
URC 73 73 0
GMP 29 30 1
TFF 2 3 1
Total 3970 3939 -31

Delta Treatment Selection:

fig2= plot_ly(alpha=.8)
fig2=fig2%>%
  add_bars(data=base.trt, type='bar', x=base.trt$ini_rx, y=base.trt$n, name='base')
fig2=fig2%>%
    add_bars(data=rf.trt, type='bar', x=rf.trt$ini_rx, y=rf.trt$n, name='RF edits')
fig2

Treatment Selection Change due to Response Function Edits

fig3= plot_ly(alpha=.6)
fig3=fig3%>%
  add_bars(data=trt.select, x=trt.select$ini_rx, y=trt.select$delta, name='base')

fig3

Takeaways:

  1. In this experiment we had 31 fewer units with ast_de0 values
  2. Our RF edits most affected CC , TH & SD by deselecting those treatments in the REBA optimization, but this shift didn’t have a large effect in our overall treatment selection method.

Experiment 2: Conditional Flame Length Evaluation

Step 1: Visual check for repeat ROI values:

cfl=cfl%>%
  filter(ast_de0>0)
fig4=plot_ly(alpha=.8)
fig4=fig4%>%
  add_trace(cfl, type='scatter', y=cfl$ast_de0, x=~as.factor(cfl$unit_id),  name=~cfl$ini_rx, color=~cfl$ini_rx)
fig4

Step 2: Evaluate number of unique values between methods

cfl.u=length(unique(cfl$ast_de0))
cfl.u
## [1] 1603
base.u=length(unique(base.reba$ast_de0))
base.u
## [1] 1651
delta.u=cfl.u-base.u
delta.u
## [1] -48

Evaluate a single treatment: ‘TH’

cfl.TH=cfl%>%
  filter(ini_rx=="TH")
base.TH=base.reba%>%
  filter(ini_rx=="TH")

Evaluate distribution of ast_de0 values between methods:

fig1= plot_ly(alpha=.6)
fig1=fig1%>%
  add_trace(data=base.TH, type='violin', y=base.TH$ast_de0, name='base',     box = list(
      visible = T
    ),
    meanline = list(
      visible = T
    ))
fig1=fig1%>%
    add_trace(data=cfl.TH, type='violin', y=cfl.TH$ast_de0, name='CFL',     box = list(
      visible = T
    ),
    meanline = list(
      visible = T
    ))
#fig1 <- fig1 %>% layout(barmode = "overlay")
fig1
fig1= plot_ly(alpha=.6)
fig1=fig1%>%
  add_trace(data=base.TH, type='violin', y=base.TH$crb_de0, name='base',     box = list(
      visible = T
    ),
    meanline = list(
      visible = T
    ))
fig1=fig1%>%
    add_trace(data=cfl.TH, type='violin', y=cfl.TH$crb_de0, name='CFL',     box = list(
      visible = T
    ),
    meanline = list(
      visible = T
    ))
#fig1 <- fig1 %>% layout(barmode = "overlay")
fig1

Evaluate a single treatment: ‘GFL’

cfl.GFL=cfl%>%
  filter(ini_rx=="GFL")
base.GFL=base.reba%>%
  filter(ini_rx=="GFL")

For a single treatment type, check that each Unit_ID has a single ast_de0

fig5=plot_ly(alpha=.8)
fig5=fig5%>%
  add_trace(cfl.GFL, type='scatter', y=cfl.GFL$ast_de0, x=as.character(cfl.GFL$unit_id), color=cfl.GFL$ast_de0)
fig5=fig5%>%
  layout(xaxis=list(range=list(0, 20)))
fig5
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Takeaways

  1. No statistically significant difference in ast_de0 between the base and CFL versions of the REBA.
  2. Instead of adding unique values, we end up with 48 fewer ast_de0 values.

Are there differences in ast_de0 across treatment types?

i.e. Are we ‘breaking ties’?

fig6= plot_ly(alpha=.4)
fig6=fig6%>%
  add_histogram(data=cfl, type='histogram', x=cfl$ast_de0, name=~cfl$ini_rx)
fig6 <- fig6 %>% layout(barmode = "overlay")
fig6

Evaluate iterative treatment selection change by each edit

cfl.trt=cfl%>%
  #group_by(ini_mthd)%>%
  count(ini_rx)%>%
  arrange(desc(n))

rf.cfl.trt=rf.cfl%>%
  #group_by(ini_mthd)%>%
  count(ini_rx)%>%
  arrange(desc(n))
trt.select=left_join(base.trt, rf.trt, by="ini_rx")
trt.select=left_join(trt.select, cfl.trt,  by='ini_rx')
trt.select=left_join(trt.select, rf.cfl.trt,  by='ini_rx')
colnames(trt.select)=c('ini_rx', 'base.reba', 'RF.Edits', 'CFL.Edits', 'RF.CFL.Edits')
trt.select=trt.select%>% 
  mutate(dRF=RF.Edits-base.reba)%>%
  mutate(dCFL=CFL.Edits-base.reba)%>%
  mutate(dRF.CFL=RF.CFL.Edits-base.reba)
trt.summ=trt.select%>%
  adorn_totals("row")%>%
  mutate(prc.RF=(RF.Edits/tail(RF.Edits, n=1))*100)%>%
  mutate(prc.CFL=(CFL.Edits/tail(RF.Edits, n=1))*100, nm.rm=TRUE)%>%
  mutate(prc.RF.CFL=(RF.CFL.Edits/tail(RF.Edits, n=1))*100, nm.rm=TRUE)%>%
  select(ini_rx, dRF, prc.RF, dCFL, prc.CFL, dRF.CFL, prc.RF.CFL)
trt.summ=trt.summ%>%
  mutate(across(where(is.numeric), round, 0))
tbl=kbl(trt.summ, align='c', booktabs=TRUE)

tbl%>%
  kable_styling(full_width = TRUE)
ini_rx dRF prc.RF dCFL prc.CFL dRF.CFL prc.RF.CFL
TH -11 33 -14 33 -12 33
GFL -3 19 19 19 24 20
CC -6 12 -7 12 -7 12
VDT -5 10 -2 10 -2 10
SD -5 10 0 10 0 10
MASL -2 9 1 9 5 9
BR 1 3 -114 0 -114 0
AI -2 2 0 2 0 2
URC 0 2 1 2 1 2
GMP 1 1 2 1 3 1
TFF 1 0 NA NA NA NA
Total -31 100 -114 98 -102 98 
trt.summ=trt.summ%>%
  slice(1:11)
dFig=plot_ly(alpha = .6)
dFig=dFig%>%
  add_bars(data=trt.summ, x=trt.summ$ini_rx, y=trt.summ$prc.RF, name='RF')%>%
  add_bars(data=trt.summ, x=trt.summ$ini_rx, y=trt.summ$prc.CFL, name='CFL')%>%
  add_bars(data=trt.summ, x=trt.summ$ini_rx, y=trt.summ$prc.RF.CFL, name='RF-CFL')
dFig=dFig%>%
  layout(title=list(text='% Treatment Allocated by Method'))
  
dFig
## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations
fig8= plot_ly(alpha=.6)
fig8=fig8%>%
  add_bars(data=trt.select, x=trt.select$ini_rx, y=trt.select$dRF, name='dRF')
fig8=fig8%>%
  add_bars(data=trt.select, x=trt.select$ini_rx, y=trt.select$dCFL, name='dCFL')
fig8=fig8%>%
  add_bars(data=trt.select, x=trt.select$ini_rx, y=trt.select$dRF.CFL, name='dRF-CFL')



fig8
## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations