Reproducible supplement

Load wrangled dataset and define key functions

Show code

df <- read.xlsx("dfs/20220331_df.xlsx",sheet=1)

# Plot calibration curves for ordinal models with 2 intercepts
fit_check_ord2 <- function(fit,predictor){
  par(mfrow=c(1,2))
  plot(rms::calibrate(fit, kint = 1,B=200,group=df$crs_g_r))
  plot(rms::calibrate(fit, kint = 2,B=200,group=df$crs_g_r))
}

# Plot calibration curves for ordinal models with 3 intercepts
fit_check_ord3 <- function(fit,predictor){
  par(mfrow=c(1,3))
  plot(rms::calibrate(fit, kint = 1,B=200,group=df$nt_g_r))
  plot(rms::calibrate(fit, kint = 2,B=200,group=df$nt_g_r))
  plot(rms::calibrate(fit, kint = 3,B=200,group=df$nt_g_r))
}

# Sparser function to extract coefficients and p values from rms fits
get_model_stats = function(x) {
  cap = capture.output(print(x))
  
  #model stats
  stats = c()
  stats$R2.adj = str_match(cap, "R2 adj\\s+ (\\d\\.\\d+)") %>% na.omit() %>% .[, 2] %>% as.numeric()
  
  #coef stats lines
  coef_lines = cap[which(str_detect(cap, "Coef\\s+S\\.E\\.")):(length(cap) - 1)]
  
  #parse
  coef_lines_table = suppressWarnings(readr::read_table(coef_lines %>% stringr::str_c(collapse = "\n")))
  colnames(coef_lines_table)[1] = "Predictor"
  
  list(
    stats = stats,
    coefs = coef_lines_table
  )
}

# Sparse regression coefficients, CI and P values from univariate and multivariable fits (minimal adjustment set of covariates) for response
resp_table_sparser <- function(fituni,fitmulti){
#Parse ORs and 95% CIs from univariate model
d <- summary(fituni)
d <- data.frame(d)
resp_unadj <- rbind(
  cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-")),
  cbind(round(d[4,4],2),paste(round(d[4,6],2),round(d[4,7],2),sep="-"))
)
colnames(resp_unadj) <- c("unadjusted OR","95% CI")
rownames(resp_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")
#Parse p values from univariate model
unadjstats_d <- get_model_stats(fituni)
#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fitmulti))
resp_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep="-"))
)
colnames(resp_adj) <- c("adjusted OR","95% CI")
rownames(resp_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")
#Parse p values from multivariable model
adjstats_d <- get_model_stats(fitmulti)
#Create table object to merge later
resp_t <- cbind(resp_unadj,unadjstats_d$coefs[2:3,5],resp_adj,adjstats_d$coefs[2:3,5])
resp_t
}


# Impute missing data
set.seed(12345)
df$preleuka_LDH_l <- log10(df$preleuka_LDH)
df$preleuka_ALC_l <- log10(df$preleuka_ALC)
df$preLD_LDH_l <- log10(df$preLD_LDH)
df$preLD_ALC_l <- log10(df$preLD_ALC)

df <- df %>% 
  mutate(product=fct_relevel(product,"axicel","tisacel","jcar014"))
dd <- datadist(df)
options(datadist='dd')
mi <- aregImpute(formula = ~ product+preleuka_LDH_l+age+hct_score+preleuka_bulk_s+preleuka_ALC_l,
                 data=df, n.impute = 10, burnin=3, nk=3, type='pmm', pmmtype=1)

Iteration 1 
Iteration 2 
Iteration 3 
Iteration 4 
Iteration 5 
Iteration 6 
Iteration 7 
Iteration 8 
Iteration 9 
Iteration 10 
Iteration 11 
Iteration 12 
Iteration 13

Show code

# Impute missing data excluding JCAR014 patient
set.seed(12346)
df2 <- df[df$product!="jcar014",]
df2$product <- as.factor(as.character(df2$product))
df2 <- df2 %>% 
  select(
    product,
    preleuka_LDH_l,
    age,
    hct_score,
    preleuka_bulk_s,
    preleuka_ALC_l,
    crs_g_r,
    nt_g_r,
    bestCR_b,
    bestCR_CT_b,
    bestCR_PET_b
  )
dd2 <- datadist(df2)
options(datadist='dd2')
mi2 <- aregImpute(formula = ~ product+preleuka_LDH_l+age+hct_score+preleuka_bulk_s+preleuka_ALC_l,
                 data=df2, n.impute = 10, burnin=3, nk=3, type='pmm', pmmtype=1)

Iteration 1 
Iteration 2 
Iteration 3 
Iteration 4 
Iteration 5 
Iteration 6 
Iteration 7 
Iteration 8 
Iteration 9 
Iteration 10 
Iteration 11 
Iteration 12 
Iteration 13

Table 1. Patient and disease characteristics

Show code

df %>%
  select(`CAR T-cell product`=product,
           Age=age,
           Sex=sex,
         `HCT-CI`=hct_score,
         `Lymphoma histology`=lymph_type,
         `Number of prior therapies`=prior_tx,
         `Relapse pattern`=relapse_pattern,
         `LDH (U/L)`=preleuka_LDH,
         `Largest lesion diameter (cm)`=preleuka_bulk_s,
         `Ann Arbor Stage`=preleuka_annarbor,
         `Extranodal disease`=preleuka_extranodal,
           `ALC (G/L)`=preleuka_ALC,
         `Bridging therapy post-leukapheresis`=bridge,
         ) %>%
  mutate(`CAR T-cell product`=fct_relevel(`CAR T-cell product`,"axicel","tisacel","jcar014")) %>%
  mutate(`CAR T-cell product`=recode_factor(`CAR T-cell product`,`axicel`="axicel",`tisacel`="tisacel",`jcar014`="JCAR014")) %>%
  tbl_summary(
    by = `CAR T-cell product`,
    type = list(all_continuous() ~ "continuous2",
                `HCT-CI` ~ 'continuous2',
                `Number of prior therapies` ~ 'continuous2'),
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})", 
                                     "{min}, {max}"),
    missing_text = "Missing data",
    sort = all_categorical() ~ "frequency") %>%
  add_p(list(all_continuous() ~ "kruskal.test",
             all_categorical() ~ "fisher.test"),
             simulate.p.value=TRUE) %>% 
  bold_labels()

Characteristic	axicel, N = 68¹	tisacel, N = 31¹	JCAR014, N = 30¹	p-value²
Age				0.077
Median (IQR)	62 (50, 66)	64 (56, 72)	60 (53, 64)
Range	25, 79	23, 81	27, 71
Sex				0.9
M	47 (69%)	21 (68%)	19 (63%)
F	21 (31%)	10 (32%)	11 (37%)
HCT-CI				0.5
Median (IQR)	1.00 (0.00, 3.00)	1.00 (0.00, 3.00)	1.00 (0.00, 2.75)
Range	0.00, 9.00	0.00, 6.00	0.00, 6.00
Lymphoma histology				0.004
Diffuse large B-cell lymphoma	50 (74%)	18 (58%)	13 (43%)
Transformed	14 (21%)	12 (39%)	8 (27%)
Other large B-cell histologies	4 (5.9%)	1 (3.2%)	7 (23%)
Burkitt lymphoma	0 (0%)	0 (0%)	2 (6.7%)
Number of prior therapies				0.004
Median (IQR)	3.00 (2.00, 4.00)	3.00 (2.00, 4.00)	4.00 (3.25, 5.75)
Range	2.00, 9.00	2.00, 9.00	2.00, 9.00
Relapse pattern				0.4
Relapsed	33 (49%)	11 (35%)	11 (37%)
Secondary refractory	22 (32%)	16 (52%)	14 (47%)
Primary refractory	13 (19%)	4 (13%)	5 (17%)
LDH (U/L)				>0.9
Median (IQR)	231 (168, 379)	221 (176, 372)	227 (184, 307)
Range	105, 2,812	145, 549	110, 641
Largest lesion diameter (cm)				0.5
Median (IQR)	4.60 (2.50, 7.30)	4.25 (2.50, 6.70)	4.05 (2.58, 5.12)
Range	1.10, 13.10	1.50, 24.80	0.80, 10.90
Missing data	9	11	4
Ann Arbor Stage				0.2
4	38 (58%)	14 (54%)	18 (60%)
3	16 (25%)	5 (19%)	8 (27%)
2	11 (17%)	4 (15%)	4 (13%)
1	0 (0%)	3 (12%)	0 (0%)
Missing data	3	5	0
Extranodal disease	40 (62%)	21 (75%)	18 (60%)	0.4
Missing data	3	3	0
ALC (G/L)				0.8
Median (IQR)	0.72 (0.42, 1.10)	0.64 (0.45, 0.99)	0.77 (0.53, 1.03)
Range	0.11, 3.82	0.28, 2.69	0.27, 2.51
Bridging therapy post-leukapheresis	40 (59%)	22 (71%)	6 (20%)	<0.001
¹ n (%)
² Kruskal-Wallis rank sum test; Fisher's exact test

Table 2. CRS, ICANS and response after CD19 CAR T-cell therapy

Show code

t2a <- df %>%
  mutate(crs_g3 = case_when(crs_g3 == "≥3" ~ 1,
                            crs_g3 == "<3" ~ 0)) %>%
  transmute(
    `CAR T-cell product` = product,
    `Grade ≥1 CRS` = crs_g1_d,
    `Grade ≥2 CRS` = crs_g2_d,
    `Grade ≥3 CRS` = crs_g3,
    `CRS duration (days)` = CRS_duration,
    `Grade ≥1 ICANS` = nt_g1_d,
    `Grade ≥2 ICANS` = nt_g2_d,
    `Grade ≥3 ICANS` = nt_g3_d,
    `ICANS duration (days)` = ICANS_duration,
    `ICANS duration (days)` = ICANS_duration,
    `Tocilizumab` = ifelse(toci_nb==0,"No","Yes"),
    `Corticosteroids` = ifelse(steroids_cum_dex_d==0,"No","Yes")
  ) %>%
  mutate(`CAR T-cell product` = fct_relevel(`CAR T-cell product`, "axicel", "tisacel", "jcar014")) %>%
  mutate(
    `CAR T-cell product` = recode_factor(
      `CAR T-cell product`,
      `axicel` = "axicel",
      `tisacel` = "tisacel",
      `jcar014` = "JCAR014"
    )
  ) %>%
  tbl_summary(
    by = `CAR T-cell product`,
    type = list(`CRS duration (days)` ~ 'continuous2',
                `ICANS duration (days)` ~ 'continuous2'
                # `Tocilizumab (number of administrations)` ~ 'continuous2',
                # `Cumulative corticosteroids dose (mg of dexamethasone)` ~ 'continuous2'
                ),
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})",
                                     "{min}, {max}"),
    missing = "no"
  ) %>%
  add_p(test = list(
        `Grade ≥1 CRS` ~ "fisher.test",
        `Grade ≥2 CRS` ~ "fisher.test",
        `Grade ≥3 CRS` ~ "fisher.test",
        `Grade ≥1 ICANS` ~ "fisher.test",
        `Grade ≥2 ICANS` ~ "fisher.test",
        `Grade ≥3 ICANS` ~ "fisher.test",
        `CRS duration (days)` ~ "kruskal.test",
        `ICANS duration (days)` ~ "kruskal.test",
        `Tocilizumab` ~ "fisher.test",
        `Corticosteroids` ~ "fisher.test"
        ),
        simulate.p.value = TRUE) %>%
  bold_labels()

#### Table 2b response ####
t2b <- df %>%
  filter(!is.na(bestresp_b)) %>% # to use the right denominator for percentage calculation (pts available for response)
  transmute(
    `CAR T-cell product` = fct_relevel(product, "axicel", "tisacel", "jcar014"),
    `Overall response rate (best response)` = as.numeric(bestresp_b),
    `Overall response rate (best response - CT-only criteria)` = as.numeric(bestresp_CT_b),
    `Overall response rate (best response - PET-only criteria)` = as.numeric(bestresp_PET_b),
    `Complete response rate (best response)` = bestCR_b,
    `Complete response rate (best response - CT-only criteria)` = as.numeric(bestCR_CT_b),
    `Complete response rate (best response - PET-only criteria)` = as.numeric(bestCR_PET_b)
  ) %>%
  mutate(
    `CAR T-cell product` = recode_factor(
      `CAR T-cell product`,
      `axicel` = "axicel",
      `tisacel` = "tisacel",
      `jcar014` = "JCAR014"
    )
  ) %>%
  tbl_summary(
    by = `CAR T-cell product`,
    # type = list(`CRS grade` ~ 'continuous2',
    #             `ICANS grade` ~ 'continuous2'),
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})",
                                     "{min}, {max}"),
    missing = "no"
  ) %>%
  add_p(test = everything() ~ "fisher.test",
        simulate.p.value = TRUE) %>%
  bold_labels() 

#### Stack table 2a and 2b ####
tbl_stack(
  tbls = list(t2a, t2b),group_header = c("Toxicity","Response"))

Characteristic	axicel, N = 68¹	tisacel, N = 31¹	JCAR014, N = 30¹	p-value²
Toxicity
Grade ≥1 CRS	59 (87%)	22 (71%)	14 (47%)	<0.001
Grade ≥2 CRS	33 (49%)	9 (29%)	7 (23%)	0.030
Grade ≥3 CRS	5 (7.4%)	0 (0%)	0 (0%)	0.2
CRS duration (days)				0.024
Median (IQR)	6.0 (4.0, 9.0)	5.0 (3.0, 6.0)	4.5 (3.0, 5.0)
Range	2.0, 22.0	1.0, 14.0	2.0, 6.0
Grade ≥1 ICANS	42 (62%)	7 (23%)	6 (20%)	<0.001
Grade ≥2 ICANS	35 (51%)	7 (23%)	5 (17%)	<0.001
Grade ≥3 ICANS	20 (29%)	4 (13%)	3 (10%)	0.045
ICANS duration (days)				0.2
Median (IQR)	6 (3, 14)	5 (2, 7)	3 (2, 4)
Range	1, 101	1, 27	1, 5
Tocilizumab	33 (49%)	13 (42%)	2 (6.7%)	<0.001
Corticosteroids	41 (60%)	5 (16%)	3 (10%)	<0.001
Response
Overall response rate (best response)	49 (75%)	18 (58%)	17 (57%)	0.10
Overall response rate (best response - CT-only criteria)	41 (63%)	16 (52%)	14 (50%)	0.4
Overall response rate (best response - PET-only criteria)	47 (80%)	13 (68%)	17 (65%)	0.3
Complete response rate (best response)	36 (55%)	10 (32%)	12 (40%)	0.077
Complete response rate (best response - CT-only criteria)	16 (25%)	7 (23%)	3 (11%)	0.3
Complete response rate (best response - PET-only criteria)	34 (58%)	9 (47%)	12 (46%)	0.5
¹ n (%)
² Fisher's exact test; Kruskal-Wallis rank sum test

Table S1. Details of aggressive lymphoma types treated with JCAR014

Show code

df %>% 
  filter(product=="jcar014") %>% 
  select(`Lymphoma histologies`=lymph_type2) %>% 
  tbl_summary(
    # type = list(`CRS grade` ~ 'continuous2',
    #             `ICANS grade` ~ 'continuous2'),
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})",
                                     "{min}, {max}"),missing_text = "Missing data",
    sort = all_categorical() ~ "frequency")

Characteristic	N = 30¹
Lymphoma histologies
Diffuse large B-cell lymphoma, NOS	20 (67%)
B-cell lymphoma, unclassifiable	2 (6.7%)
Burkitt lymphoma	2 (6.7%)
High-grade B-cell lymphoma, with MYC and BCL2 and/or BCL6 rearrangements	2 (6.7%)
EBV+ diffuse large B-cell lymphoma, NOS	1 (3.3%)
PMBCL	1 (3.3%)
Primary cutaneous DLBCL, leg type	1 (3.3%)
T-cell/histiocyte-rich large B-cell lymphoma	1 (3.3%)
¹ n (%)

Table S2. Time from leukapheresis to CAR T-cell infusion

Show code

df %>% 
  transmute(
    `CAR T-cell product` = product,
    `Time from leukapheresis to CAR T-cell infusion` = as.numeric(veintovein)
  ) %>% 
  mutate(`CAR T-cell product` = fct_relevel(`CAR T-cell product`, "axicel", "tisacel", "jcar014")) %>%
  mutate(
    `CAR T-cell product` = recode_factor(
      `CAR T-cell product`,
      `axicel` = "axicel",
      `tisacel` = "tisacel",
      `jcar014` = "JCAR014"
    )
  ) %>%
  tbl_summary(
    by = `CAR T-cell product`,
    type = `Time from leukapheresis to CAR T-cell infusion` ~ 'continuous2',
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})",
                                     "{min}, {max}"),
    missing_text = "Missing data"
  ) %>%
  add_p(test = everything() ~ "kruskal.test",
        simulate.p.value = TRUE) %>%
  bold_labels()

Characteristic	axicel, N = 68	tisacel, N = 31	JCAR014, N = 30	p-value¹
Time from leukapheresis to CAR T-cell infusion				<0.001
Median (IQR)	27 (26, 31)	40 (36, 48)	17 (17, 19)
Range	19, 50	27, 238	16, 437
¹ Kruskal-Wallis rank sum test

Table S3. Rates of severe infections and cytopenia after CD19 CAR T-cell therapy

Show code

df %>% 
  select(
    `CAR T-cell product` = product,
    `Grade ≥3 infection` = day30_sinf,
    `Grade ≥3 neutropenia` = day30G34neutropenia,
    `Grade ≥3 thrombocytopenia` = day30G34thrombocytopenia,
    `Grade ≥3 anemia` = day30G34anemia
  ) %>% 
  mutate(`CAR T-cell product` = fct_relevel(`CAR T-cell product`, "axicel", "tisacel", "jcar014")) %>%
  mutate(
    `CAR T-cell product` = recode_factor(
      `CAR T-cell product`,
      `axicel` = "axicel",
      `tisacel` = "tisacel",
      `jcar014` = "JCAR014"
    )
  ) %>%
  tbl_summary(
    by = `CAR T-cell product`,
    # type = list(`CRS duration (days)` ~ 'continuous2',
    #             `ICANS duration (days)` ~ 'continuous2'
    #             # `Tocilizumab (number of administrations)` ~ 'continuous2',
    #             # `Cumulative corticosteroids dose (mg of dexamethasone)` ~ 'continuous2'
    #             ),
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})",
                                     "{min}, {max}"),
    missing_text = "Missing data"
  ) %>%
  add_p(test = everything() ~ "fisher.test",
         simulate.p.value = TRUE) %>%
  bold_labels()

Characteristic	axicel, N = 68¹	tisacel, N = 31¹	JCAR014, N = 30¹	p-value²
Grade ≥3 infection	6 (8.8%)	1 (3.2%)	0 (0%)	0.2
Grade ≥3 neutropenia	31 (48%)	11 (35%)	6 (21%)	0.052
Missing data	3	0	2
Grade ≥3 thrombocytopenia	31 (47%)	12 (39%)	6 (21%)	0.061
Missing data	2	0	2
Grade ≥3 anemia	8 (12%)	5 (16%)	0 (0%)	0.087
Missing data	2	0	2
¹ n (%)
² Fisher's exact test

Figure 2. CRS and ICANS rates across products

Show code

df <- df %>% 
  mutate(
    product_m = case_when(
      product== "axicel" ~ "axicel",
      product=="tisacel" ~ "tisacel/JCAR014",
      product=="jcar014" ~ "tisacel/JCAR014"
    )
  )
theme_set(theme_bw())
#### CRS rates by product ####
CRS_p <- df %>% 
  mutate(crs_g=as.character(crs_g),
         product=fct_relevel(product,"axicel","tisacel","jcar014")) %>% 
  mutate(product=recode(product,"jcar014"="JCAR014")) %>% 
  group_by(product) %>%
  count(crs_g) %>% 
  mutate(perc = round((n / sum(n)) * 100, 1), labs = paste0(perc, "%")) %>% 
  ggplot(aes(product, perc, group = product)) +
  geom_col(aes(fill = crs_g))+
  geom_text(
  aes(label = labs),
  position = position_stack(vjust = 0.5),
  col = "white",
  size = 4
) +
  scale_fill_manual(
    name = "CRS grade",
    values = c("#fc9272", "#fb6a4a", "#ef3b2c", "#cb181d", "#99000d")
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_text(size = 12),
    legend.position = "right",
    legend.text = element_text(size = 10),
    axis.text.x = element_text(size = 18, color = "black")
  ) +
  scale_y_continuous(name = "Percentage")


#### ICANS rates by product ####
ICANS_p <- df %>% 
  mutate(nt_g=as.character(nt_g),
         product=fct_relevel(product,"axicel","tisacel","jcar014")) %>% 
  mutate(product=recode(product,"jcar014"="JCAR014")) %>% 
  group_by(product) %>%
  count(nt_g) %>%
  mutate(perc = round((n / sum(n)) * 100, 1), labs = paste0(perc, "%")) %>%
  ggplot(aes(product, perc, group = product)) +
  geom_col(aes(fill = nt_g)) +
  geom_text(
    aes(label = labs),
    position = position_stack(vjust = 0.5),
    col = "white",
    size = 4
  ) +
  scale_fill_manual(
    name = "ICANS grade",
    values = c("#66c2a4", "#41ae76", "#238b45", "#006d2c", "#00441b")
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_text(size = 12),
    legend.position = "right",
    legend.text = element_text(size = 10),
    axis.text.x = element_text(size = 18, color = "black")
  ) +
  scale_y_continuous(name = "Percentage")

fig1 <- CRS_p + ICANS_p
fig1 + plot_annotation(
  title = 'Figure 2.',
  tag_levels = "A",
  theme = theme(plot.title = element_text(face='bold'))
) &
  theme(plot.tag = element_text(face = 'bold'))

Show code

# Save to eps
# cairo_ps("Manuscript/Blood_resub2/Fig2.eps", width = 12, height = 7)
# fig1 + plot_annotation(
#   title = 'Figure 2.',
#   tag_levels = "A",
#   theme = theme(plot.title = element_text(size = 24,face='bold'))
# ) &
#   theme(plot.tag = element_text(size=22,face = 'bold'))
# dev.off()

Univariate and multivariable ordinal logistic regression model for CRS grade

Show code

options(datadist='dd')
#Univariate
fit_crs <- lrm(crs_g_r ~ product,x=T,y=T,data=df)
options(prType="html")
print(
  fit_crs,
  labels = c("tisacel", "JCAR014"),
  digits = 2,
  coefs = T,
  title = "Univariate ordinal proportional odds logistic regression model for CRS grade"
)

Univariate ordinal proportional odds logistic regression model for CRS grade

 lrm(formula = crs_g_r ~ product, data = df, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 15.44	R² 0.127	C 0.645
0 34	d.f. 2	g 0.709	D_xy 0.290
1 46	Pr(>χ²) 0.0004	g_r 2.033	γ 0.454
2-5 49		g_p 0.140	τ_a 0.193
max \|∂log L/∂β\| 1×10^-10		Brier 0.169

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	1.68	0.28	5.93	<0.0001
y≥2-5	0.00	0.23	0.00	0.9980
product=tisacel	-0.84	0.40	-2.08	0.0375
product=jcar014	-1.64	0.44	-3.71	0.0002

Show code

#Multivariable
fit_mi <-
  fit.mult.impute(
    crs_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x = T,
    y = T
  )

Variance Inflation Factors Due to Imputation:

       y>=1          y>=2-5 product=tisacel product=jcar014 
       1.00            1.00            1.00            1.00

preleuka_LDH_l preleuka_bulk_s age hct_score 1.01 1.13 1.01 1.00 preleuka_ALC_l 1.00

Rate of Missing Information:

       y>=1          y>=2-5 product=tisacel product=jcar014 
       0.00            0.00            0.00            0.00

preleuka_LDH_l preleuka_bulk_s age hct_score 0.01 0.12 0.01 0.00 preleuka_ALC_l 0.00

d.f. for t-distribution for Tests of Single Coefficients:

       y>=1          y>=2-5 product=tisacel product=jcar014 
 1474346.17      1483149.27      1884670.51      1221633.48

preleuka_LDH_l preleuka_bulk_s age hct_score 76673.05 644.82 105773.52 3428713.54 preleuka_ALC_l 5552666.82

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

print(
  fit_mi,
  labels = c("tisacel", "JCAR014"),
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = crs_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 21.04	R² 0.170	C 0.684
0 34	d.f. 7	g 0.920	D_xy 0.369
1 46	Pr(>χ²) 0.0037	g_r 2.510	γ 0.369
2-5 49		g_p 0.165	τ_a 0.245
max \|∂log L/∂β\| 3×10^-9		Brier 0.161

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-0.25	2.15	-0.12	0.9063
y≥2-5	-1.98	2.15	-0.92	0.3565
product=tisacel	-0.75	0.41	-1.81	0.0698
product=jcar014	-1.66	0.45	-3.67	0.0002
preleuka_LDH_l	1.26	0.81	1.56	0.1196
preleuka_bulk_s	-0.03	0.06	-0.48	0.6293
age	-0.02	0.01	-1.08	0.2802
hct_score	0.02	0.10	0.22	0.8251
preleuka_ALC_l	0.57	0.59	0.96	0.3392

Show code

forestplot(summary(fit_mi),digits=2,title="OR for CRS grade")

Show code

print(summary(fit_mi,age=c(50,60)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `crs_g_r`
preleuka_LDH_l	2.20	2.600	0.32	0.410	0.26	-0.110	0.920
Odds Ratio	2.20	2.600	0.32	1.500		0.900	2.500
preleuka_bulk_s	2.50	6.700	4.20	-0.120	0.25	-0.600	0.360
Odds Ratio	2.50	6.700	4.20	0.890		0.550	1.400
age	50.00	60.000	10.00	-0.150	0.14	-0.430	0.120
Odds Ratio	50.00	60.000	10.00	0.860		0.650	1.100
hct_score	0.00	3.000	3.00	0.068	0.31	-0.540	0.670
Odds Ratio	0.00	3.000	3.00	1.100		0.590	2.000
preleuka_ALC_l	-0.37	0.017	0.38	0.220	0.23	-0.230	0.660
Odds Ratio	-0.37	0.017	0.38	1.200		0.800	1.900
product --- tisacel:axicel	1.00	2.000		-0.750	0.41	-1.600	0.061
Odds Ratio	1.00	2.000		0.470		0.210	1.100
product --- jcar014:axicel	1.00	3.000		-1.700	0.45	-2.600	-0.770
Odds Ratio	1.00	3.000		0.190		0.078	0.460

Show code

print(anova(fit_mi),digits = 2)

	χ²	d.f.	P
Wald Statistics for `crs_g_r`
product	14.00	2	0.0009
preleuka_LDH_l	2.42	1	0.1196
preleuka_bulk_s	0.23	1	0.6293
age	1.17	1	0.2802
hct_score	0.05	1	0.8251
preleuka_ALC_l	0.91	1	0.3392
TOTAL	18.81	7	0.0088

Show code

plot(anova(fit_mi))

Show code

fit_check_ord2(fit_mi)

n=129 Mean absolute error=0.019 Mean squared error=0.00055 0.9 Quantile of absolute error=0.035

n=129 Mean absolute error=0.057 Mean squared error=0.0054 0.9 Quantile of absolute error=0.105

Show code

options(prType=NULL)

#Sparse ORs and 95% CIs from univariate model
crs_d <- data.frame(summary(fit_crs))
crs_unadj <- rbind(
  cbind(round(crs_d[2,4],2),paste(round(crs_d[2,6],2),round(crs_d[2,7],2),sep="-")),
  cbind(round(crs_d[4,4],2),paste(round(crs_d[4,6],2),round(crs_d[4,7],2),sep="-"))
)
colnames(crs_unadj) <- c("unadjusted OR","95% CI")
rownames(crs_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
crs_unadjstats_d <- get_model_stats(fit_crs)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
crs_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep = "-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep = "-"))
)
colnames(crs_adj) <- c("adjusted OR","95% CI")
rownames(crs_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from multivariable model
crs_adjstats_d <- get_model_stats(fit_mi)

### Create object to merge tables later
CRS_t <-
  cbind(crs_unadj,
        crs_unadjstats_d$coefs[3:4, 5],
        crs_adj,
        crs_adjstats_d$coefs[3:4, 5])

Sensitivity analysis: univariate and multivariable ordinal logistic regression model for CRS grade excluding JCAR014 patients

Show code

options(datadist='dd2')
#Univariate
fit_crs_sub <- lrm(crs_g_r ~ product,x=T,y=T,data=df2)
options(prType="html")
print(
  fit_crs,
  labels = c("tisacel", "JCAR014"),
  digits = 2,
  coefs = T,
  title = "Univariate ordinal proportional odds logistic regression model for CRS grade"
)

Univariate ordinal proportional odds logistic regression model for CRS grade

 lrm(formula = crs_g_r ~ product, data = df, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 15.44	R² 0.127	C 0.645
0 34	d.f. 2	g 0.709	D_xy 0.290
1 46	Pr(>χ²) 0.0004	g_r 2.033	γ 0.454
2-5 49		g_p 0.140	τ_a 0.193
max \|∂log L/∂β\| 1×10^-10		Brier 0.169

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	1.68	0.28	5.93	<0.0001
y≥2-5	0.00	0.23	0.00	0.9980
product=tisacel	-0.84	0.40	-2.08	0.0375
product=jcar014	-1.64	0.44	-3.71	0.0002

Show code

#Multivariate
fit_mi_sub <-
  fit.mult.impute(
    crs_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df2,
    xtrans = mi2,
    x = T,
    y = T
  )

Variance Inflation Factors Due to Imputation:

       y>=1          y>=2-5 product=tisacel  preleuka_LDH_l 
       1.00            1.00            1.00            1.01

preleuka_bulk_s age hct_score preleuka_ALC_l 1.12 1.00 1.00 1.00

Rate of Missing Information:

       y>=1          y>=2-5 product=tisacel  preleuka_LDH_l 
       0.00            0.00            0.00            0.01

preleuka_bulk_s age hct_score preleuka_ALC_l 0.11 0.00 0.00 0.00

d.f. for t-distribution for Tests of Single Coefficients:

       y>=1          y>=2-5 product=tisacel  preleuka_LDH_l 
 8286063.36      6704358.53      5789491.73       120171.67

preleuka_bulk_s age hct_score preleuka_ALC_l 738.18 400335.41 25144833.66 7334750.50

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

print(
  fit_mi_sub,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = crs_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi2, data = df2, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 99	LR χ² 8.28	R² 0.092	C 0.642
0 18	d.f. 6	g 0.648	D_xy 0.285
1 39	Pr(>χ²) 0.2200	g_r 1.913	γ 0.285
2-5 42		g_p 0.091	τ_a 0.182
max \|∂log L/∂β\| 9×10^-11		Brier 0.139

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-0.57	2.33	-0.24	0.8082
y≥2-5	-2.50	2.34	-1.07	0.2866
product=tisacel	-0.86	0.42	-2.02	0.0433
preleuka_LDH_l	1.37	0.88	1.56	0.1197
preleuka_bulk_s	-0.05	0.06	-0.75	0.4520
age	-0.01	0.02	-0.57	0.5677
hct_score	-0.04	0.11	-0.36	0.7225
preleuka_ALC_l	0.36	0.64	0.56	0.5765

Show code

print(summary(fit_mi_sub,age=c(50,60)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `crs_g_r`
preleuka_LDH_l	2.20	2.600	0.34	0.470	0.30	-0.12	1.100
Odds Ratio	2.20	2.600	0.34	1.600		0.89	2.900
preleuka_bulk_s	2.50	7.200	4.70	-0.220	0.29	-0.79	0.350
Odds Ratio	2.50	7.200	4.70	0.800		0.45	1.400
age	50.00	60.000	10.00	-0.088	0.15	-0.39	0.210
Odds Ratio	50.00	60.000	10.00	0.920		0.68	1.200
hct_score	0.00	3.000	3.00	-0.120	0.34	-0.78	0.540
Odds Ratio	0.00	3.000	3.00	0.890		0.46	1.700
preleuka_ALC_l	-0.38	0.021	0.40	0.140	0.26	-0.36	0.650
Odds Ratio	-0.38	0.021	0.40	1.200		0.70	1.900
product --- tisacel:axicel	1.00	2.000		-0.860	0.42	-1.70	-0.026
Odds Ratio	1.00	2.000		0.430		0.19	0.970

Show code

print(anova(fit_mi_sub),digits = 2)

	χ²	d.f.	P
Wald Statistics for `crs_g_r`
product	4.08	1	0.0433
preleuka_LDH_l	2.42	1	0.1197
preleuka_bulk_s	0.57	1	0.4520
age	0.33	1	0.5677
hct_score	0.13	1	0.7225
preleuka_ALC_l	0.31	1	0.5765
TOTAL	7.53	6	0.2750

Show code

plot(anova(fit_mi_sub))

Show code

# fit_check_ord2(fit_mi_sub)
options(prType=NULL)

#Sparse ORs and 95% CIs from univariate model
crs_d <- data.frame(summary(fit_crs_sub))
crs_unadj <- cbind(round(crs_d[2,4],2),paste(round(crs_d[2,6],2),round(crs_d[2,7],2),sep="-"))
colnames(crs_unadj) <- c("unadjusted OR","95% CI")
rownames(crs_unadj) <- c("Tisacel versus axicel")

#Parse p values from univariate model
crs_unadjstats_d <- get_model_stats(fit_crs_sub)

#Parse p values from multivariable model
crs_adjstats_d <- get_model_stats(fit_mi_sub)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi_sub))
crs_adj <- cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep = "-"))
colnames(crs_adj) <- c("adjusted OR","95% CI")
rownames(crs_adj) <- c("Tisacel versus axicel")

### Create object to merge tables later
CRS_sub_t <-
  cbind(crs_unadj,
        crs_unadjstats_d$coefs[3, 5],
        crs_adj,
        crs_adjstats_d$coefs[3, 5])

Sensitivity analysis: multivariable ordinal logistic regression model for CRS including bridging

Show code

options(datadist='dd')
fit_mi <-
  fit.mult.impute(
    crs_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x = T,
    y = T
  )


Variance Inflation Factors Due to Imputation:

           y>=1          y>=2-5 product=tisacel product=jcar014 
           1.00            1.00            1.00            1.00 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
           1.01            1.13            1.01            1.00 
 preleuka_ALC_l      bridge=yes 
           1.00            1.00 

Rate of Missing Information:

           y>=1          y>=2-5 product=tisacel product=jcar014 
           0.00            0.00            0.00            0.00 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
           0.01            0.11            0.01            0.00 
 preleuka_ALC_l      bridge=yes 
           0.00            0.00 

d.f. for t-distribution for Tests of Single Coefficients:

           y>=1          y>=2-5 product=tisacel product=jcar014 
     1395803.87      1400672.09      2549622.19      1405355.31 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
       97417.66          688.84       112741.92      4037078.01 
 preleuka_ALC_l      bridge=yes 
     5595900.70     34488058.60 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  labels = c("tisacel", "JCAR014"),
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for CRS grade"
)

Multivariable ordinal proportional odds logistic regression model for CRS grade

 fit.mult.impute(formula = crs_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge, 
     fitter = lrm, xtrans = mi, data = df, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 21.10	R² 0.170	C 0.683
0 34	d.f. 8	g 0.923	D_xy 0.366
1 46	Pr(>χ²) 0.0069	g_r 2.517	γ 0.366
2-5 49		g_p 0.166	τ_a 0.243
max \|∂log L/∂β\| 3×10^-9		Brier 0.161

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-0.16	2.19	-0.07	0.9435
y≥2-5	-1.89	2.19	-0.86	0.3891
product=tisacel	-0.77	0.42	-1.82	0.0681
product=jcar014	-1.63	0.47	-3.46	0.0005
preleuka_LDH_l	1.19	0.86	1.38	0.1685
preleuka_bulk_s	-0.03	0.06	-0.48	0.6310
age	-0.01	0.01	-1.06	0.2904
hct_score	0.02	0.10	0.24	0.8112
preleuka_ALC_l	0.57	0.59	0.96	0.3369
bridge=yes	0.09	0.39	0.23	0.8172

Show code

options(prType=NULL)

#Sparse ORs and 95% CIs from univariate model
crs_d <- data.frame(summary(fit_crs))
crs_unadj <- rbind(
  cbind(round(crs_d[2,4],2),paste(round(crs_d[2,6],2),round(crs_d[2,7],2),sep="-")),
  cbind(round(crs_d[4,4],2),paste(round(crs_d[4,6],2),round(crs_d[4,7],2),sep="-"))
)
colnames(crs_unadj) <- c("unadjusted OR","95% CI")
rownames(crs_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
crs_unadjstats_d <- get_model_stats(fit_crs)

#Parse p values from multivariable model
crs_adjstats_d <- get_model_stats(fit_mi)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
crs_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep = "-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep = "-"))
)
colnames(crs_adj) <- c("adjusted OR","95% CI")
rownames(crs_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

### Create object to merge tables later
CRS_bridge_t <-
  cbind(crs_unadj,
        crs_unadjstats_d$coefs[3:4, 5],
        crs_adj,
        crs_adjstats_d$coefs[3:4, 5])

Sensitivity analysis: multivariable ordinal logistic regression model for CRS including preLD ALC instead of preleukapheresis ALC

Show code

options(datadist='dd')
fit_mi <-
  fit.mult.impute(
    crs_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x = T,
    y = T
  )


Variance Inflation Factors Due to Imputation:

           y>=1          y>=2-5 product=tisacel product=jcar014 
           1.00            1.00            1.00            1.00 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
           1.01            1.13            1.01            1.00 
    preLD_ALC_l      bridge=yes 
           1.00            1.00 

Rate of Missing Information:

           y>=1          y>=2-5 product=tisacel product=jcar014 
           0.00            0.00            0.00            0.00 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
           0.01            0.11            0.01            0.00 
    preLD_ALC_l      bridge=yes 
           0.00            0.00 

d.f. for t-distribution for Tests of Single Coefficients:

           y>=1          y>=2-5 product=tisacel product=jcar014 
     1520503.26      1513649.57      2025017.03      1584498.18 
 preleuka_LDH_l preleuka_bulk_s             age       hct_score 
       91551.83          718.21       178660.14      2516157.74 
    preLD_ALC_l      bridge=yes 
     1910762.40     25291395.73 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  labels = c("tisacel", "JCAR014"),
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for CRS grade"
)

Multivariable ordinal proportional odds logistic regression model for CRS grade

 fit.mult.impute(formula = crs_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge, 
     fitter = lrm, xtrans = mi, data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

         crs_g_r         product  preleuka_LDH_l preleuka_bulk_s 
               0               0               0               0 
             age       hct_score     preLD_ALC_l          bridge 
               0               0               1               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 128	LR χ² 20.43	R² 0.166	C 0.678
0 34	d.f. 8	g 0.900	D_xy 0.355
1 45	Pr(>χ²) 0.0089	g_r 2.460	γ 0.355
2-5 49		g_p 0.164	τ_a 0.236
max \|∂log L/∂β\| 3×10^-9		Brier 0.163

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-0.33	2.18	-0.15	0.8777
y≥2-5	-2.04	2.18	-0.93	0.3501
product=tisacel	-0.75	0.42	-1.78	0.0755
product=jcar014	-1.59	0.47	-3.39	0.0007
preleuka_LDH_l	1.21	0.86	1.40	0.1608
preleuka_bulk_s	-0.03	0.06	-0.50	0.6158
age	-0.01	0.01	-0.96	0.3386
hct_score	0.04	0.10	0.34	0.7338
preLD_ALC_l	0.24	0.35	0.67	0.5019
bridge=yes	0.08	0.39	0.21	0.8331

Show code

options(prType=NULL)

#Sparse ORs and 95% CIs from univariate model
crs_d <- data.frame(summary(fit_crs))
crs_unadj <- rbind(
  cbind(round(crs_d[2,4],2),paste(round(crs_d[2,6],2),round(crs_d[2,7],2),sep="-")),
  cbind(round(crs_d[4,4],2),paste(round(crs_d[4,6],2),round(crs_d[4,7],2),sep="-"))
)
colnames(crs_unadj) <- c("unadjusted OR","95% CI")
rownames(crs_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
crs_unadjstats_d <- get_model_stats(fit_crs)

#Parse p values from multivariable model
crs_adjstats_d <- get_model_stats(fit_mi)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
crs_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep = "-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep = "-"))
)
colnames(crs_adj) <- c("adjusted OR","95% CI")
rownames(crs_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

### Create object to merge tables later
CRS_preLDALC_t <-
  cbind(crs_unadj,
        crs_unadjstats_d$coefs[3:4, 5],
        crs_adj,
        crs_adjstats_d$coefs[3:4, 5])

Univariate and multivariable proportional odds logistic regression model for ICANS grade

Show code

options(datadist='dd')
#Univariate
fit_nt <- lrm(nt_g_r ~ product,x=T,y=T,data=df)
options(prType="html")
print(
  fit_nt,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 lrm(formula = nt_g_r ~ product, data = df, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  74   8  20  27

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 19.79	R² 0.160	C 0.666
max \|∂log L/∂β\| 1×10^-9	d.f. 2	g 0.839	D_xy 0.332
	Pr(>χ²) <0.0001	g_r 2.314	γ 0.549
		g_p 0.176	τ_a 0.200
		Brier 0.206

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	0.38	0.24	1.60	0.1092
y≥2	0.07	0.23	0.31	0.7537
y≥3-5	-0.78	0.25	-3.13	0.0017
product=tisacel	-1.52	0.49	-3.14	0.0017
product=jcar014	-1.73	0.51	-3.40	0.0007

Show code

#Multivariable
fit_mi <-
  fit.mult.impute(
    nt_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

       y>=1            y>=2          y>=3-5 product=tisacel 
       1.00            1.00            1.00            1.01

product=jcar014 preleuka_LDH_l preleuka_bulk_s age 1.00 1.02 1.23 1.01 hct_score preleuka_ALC_l 1.00 1.00

Rate of Missing Information:

       y>=1            y>=2          y>=3-5 product=tisacel 
       0.00            0.00            0.00            0.01

product=jcar014 preleuka_LDH_l preleuka_bulk_s age 0.00 0.02 0.19 0.01 hct_score preleuka_ALC_l 0.00 0.00

d.f. for t-distribution for Tests of Single Coefficients:

       y>=1            y>=2          y>=3-5 product=tisacel 
15504028.07     15488201.77     15930170.81        49514.81

product=jcar014 preleuka_LDH_l preleuka_bulk_s age 21973632.32 37932.93 256.98 307653.67 hct_score preleuka_ALC_l 4792663.68 14041663.68

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

print(
  fit_mi, 
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for ICANS grade"
)

Multivariable ordinal proportional odds logistic regression model for ICANS grade

 fit.mult.impute(formula = nt_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  74   8  20  27

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 25.97	R² 0.205	C 0.716
max \|∂log L/∂β\| 4×10^-7	d.f. 7	g 1.121	D_xy 0.432
	Pr(>χ²) 0.0005	g_r 3.067	γ 0.432
		g_p 0.228	τ_a 0.261
		Brier 0.193

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-3.21	2.16	-1.49	0.1370
y≥2	-3.53	2.17	-1.63	0.1026
y≥3-5	-4.44	2.18	-2.03	0.0420
product=tisacel	-1.75	0.52	-3.35	0.0008
product=jcar014	-1.78	0.52	-3.42	0.0006
preleuka_LDH_l	1.01	0.74	1.36	0.1730
preleuka_bulk_s	0.02	0.06	0.28	0.7810
age	0.02	0.02	1.43	0.1527
hct_score	-0.02	0.11	-0.19	0.8465
preleuka_ALC_l	0.87	0.60	1.45	0.1476

Show code

forestplot(summary(fit_mi),digits=2,title="OR for ICANS grade")

Show code

print(summary(fit_mi,age=c(50,60)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `nt_g_r`
preleuka_LDH_l	2.20	2.600	0.32	0.320	0.24	-0.140	0.79
Odds Ratio	2.20	2.600	0.32	1.400		0.870	2.20
preleuka_bulk_s	2.50	6.700	4.20	0.074	0.27	-0.450	0.60
Odds Ratio	2.50	6.700	4.20	1.100		0.640	1.80
age	50.00	60.000	10.00	0.230	0.16	-0.084	0.54
Odds Ratio	50.00	60.000	10.00	1.300		0.920	1.70
hct_score	0.00	3.000	3.00	-0.061	0.32	-0.680	0.56
Odds Ratio	0.00	3.000	3.00	0.940		0.510	1.70
preleuka_ALC_l	-0.37	0.017	0.38	0.330	0.23	-0.120	0.78
Odds Ratio	-0.37	0.017	0.38	1.400		0.890	2.20
product --- tisacel:axicel	1.00	2.000		-1.700	0.52	-2.800	-0.73
Odds Ratio	1.00	2.000		0.170		0.063	0.48
product --- jcar014:axicel	1.00	3.000		-1.800	0.52	-2.800	-0.76
Odds Ratio	1.00	3.000		0.170		0.061	0.47

Show code

print(anova(fit_mi),digits = 2)

	χ²	d.f.	P
Wald Statistics for `nt_g_r`
product	19.11	2	<0.0001
preleuka_LDH_l	1.86	1	0.1730
preleuka_bulk_s	0.08	1	0.7810
age	2.04	1	0.1527
hct_score	0.04	1	0.8465
preleuka_ALC_l	2.10	1	0.1476
TOTAL	22.20	7	0.0023

Show code

plot(anova(fit_mi))

Show code

fit_check_ord3(fit_mi)

n=129 Mean absolute error=0.026 Mean squared error=0.00106 0.9 Quantile of absolute error=0.052

n=129 Mean absolute error=0.035 Mean squared error=0.0019 0.9 Quantile of absolute error=0.069

n=129 Mean absolute error=0.039 Mean squared error=0.00235 0.9 Quantile of absolute error=0.074

Show code

options(prType=NULL)

#Parse ORs and 95%CI from univariate model
d <- summary(fit_nt)
d <- data.frame(d)
nt_unadj <- rbind(
  cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-")),
  cbind(round(d[4,4],2),paste(round(d[4,6],2),round(d[4,7],2),sep="-"))
)
colnames(nt_unadj) <- c("unadjusted OR","95% CI")
rownames(nt_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_nt)

#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
nt_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep="-"))
)
colnames(nt_adj) <- c("adjusted OR","95% CI")
rownames(nt_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Create object to merge later
NT_t <-
  cbind(nt_unadj,
        round(unadjstats_d$coefs[4:5, 5], 4),
        nt_adj,
        round(adjstats_d$coefs[4:5, 5], 4))

Sensitivity analysis: univariate and multivariable proportional odds logistic regression model for ICANS grade excluding JCAR014 patients

Show code

options(datadist='dd2')
#Univariate
fit_nt_sub <- lrm(nt_g_r ~ product,x=T,y=T,data=df2)
options(prType="html")
print(
  fit_nt_sub,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 lrm(formula = nt_g_r ~ product, data = df2, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  50   7  18  24

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 99	LR χ² 11.50	R² 0.121	C 0.622
max \|∂log L/∂β\| 2×10^-10	d.f. 1	g 0.669	D_xy 0.245
	Pr(>χ²) 0.0007	g_r 1.952	γ 0.572
		g_p 0.144	τ_a 0.160
		Brier 0.227

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	0.40	0.24	1.67	0.0942
y≥2	0.08	0.23	0.34	0.7313
y≥3-5	-0.81	0.25	-3.18	0.0015
product=tisacel	-1.54	0.49	-3.16	0.0016

Show code

#Multivariate
fit_mi_sub <-
  fit.mult.impute(
    nt_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df2,
    xtrans = mi2,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

       y>=1            y>=2          y>=3-5 product=tisacel 
       1.00            1.00            1.00            1.01

preleuka_LDH_l preleuka_bulk_s age hct_score 1.02 1.23 1.00 1.00 preleuka_ALC_l 1.00

Rate of Missing Information:

       y>=1            y>=2          y>=3-5 product=tisacel 
       0.00            0.00            0.00            0.01

preleuka_LDH_l preleuka_bulk_s age hct_score 0.02 0.19 0.00 0.00 preleuka_ALC_l 0.00

d.f. for t-distribution for Tests of Single Coefficients:

       y>=1            y>=2          y>=3-5 product=tisacel 
66097322.73     71061564.11    110606919.99       112886.93

preleuka_LDH_l preleuka_bulk_s age hct_score 33707.95 253.50 535731.68 5881558.60 preleuka_ALC_l 275035062.15

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

print(
  fit_mi_sub,
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for ICANS grade"
)

Multivariable ordinal proportional odds logistic regression model for ICANS grade

 fit.mult.impute(formula = nt_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi2, data = df2, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  50   7  18  24

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 99	LR χ² 19.00	R² 0.193	C 0.705
max \|∂log L/∂β\| 3×10^-7	d.f. 6	g 1.102	D_xy 0.411
	Pr(>χ²) 0.0042	g_r 3.010	γ 0.411
		g_p 0.230	τ_a 0.269
		Brier 0.208

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-3.50	2.32	-1.51	0.1317
y≥2	-3.84	2.33	-1.65	0.0986
y≥3-5	-4.81	2.35	-2.05	0.0406
product=tisacel	-1.84	0.53	-3.48	0.0005
preleuka_LDH_l	0.89	0.78	1.15	0.2515
preleuka_bulk_s	0.03	0.07	0.46	0.6468
age	0.03	0.02	1.81	0.0702
hct_score	0.01	0.12	0.05	0.9615
preleuka_ALC_l	0.94	0.63	1.50	0.1345

Show code

print(summary(fit_mi_sub,age=c(50,60)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `nt_g_r`
preleuka_LDH_l	2.20	2.600	0.34	0.300	0.27	-0.220	0.83
Odds Ratio	2.20	2.600	0.34	1.400		0.810	2.30
preleuka_bulk_s	2.50	7.200	4.70	0.140	0.31	-0.460	0.74
Odds Ratio	2.50	7.200	4.70	1.200		0.630	2.10
age	50.00	60.000	10.00	0.320	0.17	-0.026	0.66
Odds Ratio	50.00	60.000	10.00	1.400		0.970	1.90
hct_score	0.00	3.000	3.00	0.017	0.35	-0.660	0.70
Odds Ratio	0.00	3.000	3.00	1.000		0.520	2.00
preleuka_ALC_l	-0.38	0.021	0.40	0.380	0.25	-0.120	0.88
Odds Ratio	-0.38	0.021	0.40	1.500		0.890	2.40
product --- tisacel:axicel	1.00	2.000		-1.800	0.53	-2.900	-0.80
Odds Ratio	1.00	2.000		0.160		0.057	0.45

Show code

print(anova(fit_mi_sub),digits = 2)

	χ²	d.f.	P
Wald Statistics for `nt_g_r`
product	12.10	1	0.0005
preleuka_LDH_l	1.32	1	0.2515
preleuka_bulk_s	0.21	1	0.6468
age	3.28	1	0.0702
hct_score	0.00	1	0.9615
preleuka_ALC_l	2.24	1	0.1345
TOTAL	15.60	6	0.0160

Show code

plot(anova(fit_mi_sub))

Show code

options(prType=NULL)

#Parse ORs and 95%CI from univariate model
d <- summary(fit_nt_sub)
d <- data.frame(d)
nt_unadj <- cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-"))
colnames(nt_unadj) <- c("unadjusted OR","95% CI")
rownames(nt_unadj) <- c("Tisacel versus axicel")

#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_nt_sub)

#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi_sub)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi_sub))
nt_adj <- cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-"))
colnames(nt_adj) <- c("adjusted OR","95% CI")
rownames(nt_adj) <- c("Tisacel versus axicel")

#Create object to merge later
NT_sub_t <-
  cbind(nt_unadj,
        round(unadjstats_d$coefs[4, 5], 4),
        nt_adj,
        round(adjstats_d$coefs[4, 5], 4))

Sensitivity analysis: fit multivariable model for ICANS including bridging therapy

Show code

options(datadist='dd')
fit_mi <-
  fit.mult.impute(
    nt_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

           y>=1            y>=2          y>=3-5 product=tisacel 
           1.00            1.00            1.00            1.01 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
           1.00            1.02            1.24            1.01 
      hct_score  preleuka_ALC_l      bridge=yes 
           1.00            1.00            1.00 

Rate of Missing Information:

           y>=1            y>=2          y>=3-5 product=tisacel 
           0.00            0.00            0.00            0.01 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
           0.00            0.02            0.20            0.01 
      hct_score  preleuka_ALC_l      bridge=yes 
           0.00            0.00            0.00 

d.f. for t-distribution for Tests of Single Coefficients:

           y>=1            y>=2          y>=3-5 product=tisacel 
     9608649.10      9571788.16      9654165.99        49822.01 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
     9854844.26        35883.70          234.70       205141.92 
      hct_score  preleuka_ALC_l      bridge=yes 
     2198939.87     17189259.07     46816900.34 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = nt_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge, 
     fitter = lrm, xtrans = mi, data = df, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  74   8  20  27

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 26.88	R² 0.211	C 0.719
max \|∂log L/∂β\| 4×10^-7	d.f. 8	g 1.151	D_xy 0.438
	Pr(>χ²) 0.0007	g_r 3.162	γ 0.438
		g_p 0.233	τ_a 0.264
		Brier 0.192

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-2.93	2.17	-1.35	0.1761
y≥2	-3.25	2.17	-1.50	0.1339
y≥3-5	-4.16	2.19	-1.90	0.0569
product=tisacel	-1.82	0.53	-3.44	0.0006
product=jcar014	-1.65	0.54	-3.06	0.0022
preleuka_LDH_l	0.75	0.78	0.97	0.3341
preleuka_bulk_s	0.02	0.07	0.31	0.7589
age	0.02	0.02	1.52	0.1278
hct_score	-0.02	0.11	-0.18	0.8574
preleuka_ALC_l	0.90	0.60	1.50	0.1342
bridge=yes	0.40	0.42	0.95	0.3418

Show code

options(prType=NULL)
d <- summary(fit_nt)
d <- data.frame(d)

#Parse ORs and 95%CI from univariate model
nt_unadj <- rbind(
  cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-")),
  cbind(round(d[4,4],2),paste(round(d[4,6],2),round(d[4,7],2),sep="-"))
)
colnames(nt_unadj) <- c("unadjusted OR","95% CI")
rownames(nt_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_nt)

#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
nt_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep="-"))
)
colnames(nt_adj) <- c("adjusted OR","95% CI")
rownames(nt_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Create object to merge later
NT_bridge_t <-
  cbind(nt_unadj,
        round(unadjstats_d$coefs[4:5, 5], 4),
        nt_adj,
        round(adjstats_d$coefs[4:5, 5], 4))

Sensitivity analysis: fit multivariable model for ICANS including preLD ALC instead of preleukapheresis ALC

Show code

options(datadist='dd')
fit_mi <-
  fit.mult.impute(
    nt_g_r ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

           y>=1            y>=2          y>=3-5 product=tisacel 
           1.00            1.00            1.00            1.01 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
           1.00            1.01            1.20            1.01 
      hct_score     preLD_ALC_l      bridge=yes 
           1.00            1.01            1.00 

Rate of Missing Information:

           y>=1            y>=2          y>=3-5 product=tisacel 
           0.00            0.00            0.00            0.01 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
           0.00            0.01            0.17            0.01 
      hct_score     preLD_ALC_l      bridge=yes 
           0.00            0.01            0.00 

d.f. for t-distribution for Tests of Single Coefficients:

           y>=1            y>=2          y>=3-5 product=tisacel 
     9932112.01      9978272.10     10073552.28        47307.73 
product=jcar014  preleuka_LDH_l preleuka_bulk_s             age 
     9501228.25        54442.96          327.34       237368.57 
      hct_score     preLD_ALC_l      bridge=yes 
     1828888.37       306247.05      8120345.25 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = nt_g_r ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge, 
     fitter = lrm, xtrans = mi, data = df, x = T, y = T)

Frequencies of Responses

   0   1   2 3-5 
  73   8  20  27

Frequencies of Missing Values Due to Each Variable

          nt_g_r         product  preleuka_LDH_l preleuka_bulk_s 
               0               0               0               0 
             age       hct_score     preLD_ALC_l          bridge 
               0               0               1               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 128	LR χ² 25.65	R² 0.204	C 0.712
max \|∂log L/∂β\| 4×10^-7	d.f. 8	g 1.120	D_xy 0.424
	Pr(>χ²) 0.0012	g_r 3.065	γ 0.424
		g_p 0.228	τ_a 0.257
		Brier 0.195

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-3.35	2.17	-1.55	0.1223
y≥2	-3.67	2.17	-1.69	0.0911
y≥3-5	-4.57	2.19	-2.09	0.0367
product=tisacel	-1.74	0.53	-3.29	0.0010
product=jcar014	-1.60	0.54	-2.99	0.0028
preleuka_LDH_l	0.78	0.78	1.01	0.3140
preleuka_bulk_s	0.03	0.06	0.42	0.6771
age	0.03	0.02	1.80	0.0715
hct_score	0.00	0.11	-0.04	0.9677
preLD_ALC_l	0.41	0.37	1.12	0.2646
bridge=yes	0.40	0.42	0.96	0.3376

Show code

options(prType=NULL)
d <- summary(fit_nt)
d <- data.frame(d)

#Parse ORs and 95%CI from univariate model
nt_unadj <- rbind(
  cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-")),
  cbind(round(d[4,4],2),paste(round(d[4,6],2),round(d[4,7],2),sep="-"))
)
colnames(nt_unadj) <- c("unadjusted OR","95% CI")
rownames(nt_unadj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_nt)

#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi)

#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi))
nt_adj <- rbind(
  cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-")),
  cbind(round(d[14,4],2),paste(round(d[14,6],2),round(d[14,7],2),sep="-"))
)
colnames(nt_adj) <- c("adjusted OR","95% CI")
rownames(nt_adj) <- c("Tisacel versus axicel","JCAR014 versus axicel")

#Create object to merge later
NT_preLDALC_t <-
  cbind(nt_unadj,
        round(unadjstats_d$coefs[4:5, 5], 4),
        nt_adj,
        round(adjstats_d$coefs[4:5, 5], 4))

Univariate and multivariable logistic regression model for CR (best response by PET or CT)

Show code

options(datadist='dd')
#Univariate
fit_uni <- lrm(bestCR_b ~ product,x=T,y=T,data=df)
options(prType="html")
print(
  fit_uni,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 lrm(formula = bestCR_b ~ product, data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

 bestCR_b  product 
        3        0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 126	LR χ² 5.16	R² 0.054	C 0.606
0 68	d.f. 2	g 0.439	D_xy 0.212
1 58	Pr(>χ²) 0.0758	g_r 1.551	γ 0.339
max \|∂log L/∂β\| 2×10^-8		g_p 0.106	τ_a 0.106
		Brier 0.238

	β	S.E.	Wald Z	Pr(>\|Z\|)
Intercept	0.22	0.25	0.87	0.3862
product=tisacel	-0.96	0.46	-2.09	0.0365
product=jcar014	-0.62	0.45	-1.39	0.1657

Show code

fit_uniCT <- lrm(bestCR_CT_b ~ product,x=T,y=T,data=df)
fit_uniPET <- lrm(bestCR_PET_b ~ product,x=T,y=T,data=df)

#Multivariate
fit_mi <-
  fit.mult.impute(
    bestCR_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       1.04            1.20            1.11            1.02

preleuka_bulk_s age hct_score preleuka_ALC_l 1.84 1.06 1.01 1.09

Rate of Missing Information:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       0.04            0.17            0.10            0.02

preleuka_bulk_s age hct_score preleuka_ALC_l 0.46 0.06 0.01 0.08

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
    5968.60          312.62          875.13        17132.55

preleuka_bulk_s age hct_score preleuka_ALC_l 43.00 2536.71 51613.38 1424.91

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression for CR (best response)"
)

Multivariable ordinal proportional odds logistic regression for CR (best response)

 fit.mult.impute(formula = bestCR_b ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

        bestCR_b         product  preleuka_LDH_l preleuka_bulk_s 
               3               0               0               0 
             age       hct_score  preleuka_ALC_l 
               0               0               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 126	LR χ² 54.27	R² 0.467	C 0.848
0 68	d.f. 7	g 2.230	D_xy 0.696
1 58	Pr(>χ²) <0.0001	g_r 9.494	γ 0.696
max \|∂log L/∂β\| 7×10^-6		g_p 0.351	τ_a 0.349
		Brier 0.159

	β	S.E.	Wald Z	Pr(>\|Z\|)
Intercept	13.25	3.57	3.71	0.0002
product=tisacel	-1.49	0.64	-2.34	0.0191
product=jcar014	-1.55	0.63	-2.46	0.0138
preleuka_LDH_l	-4.63	1.31	-3.54	0.0004
preleuka_bulk_s	-0.30	0.14	-2.20	0.0275
age	0.00	0.02	0.02	0.9815
hct_score	0.10	0.14	0.71	0.4768
preleuka_ALC_l	2.29	0.91	2.52	0.0116

Show code

print(summary(fit_mi),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `bestCR_b`
preleuka_LDH_l	2.20	2.600	0.32	-1.5000	0.42	-2.300	-0.67
Odds Ratio	2.20	2.600	0.32	0.2200		0.098	0.51
preleuka_bulk_s	2.50	6.700	4.20	-1.3000	0.57	-2.400	-0.14
Odds Ratio	2.50	6.700	4.20	0.2900		0.094	0.87
age	53.00	67.000	14.00	0.0065	0.28	-0.540	0.55
Odds Ratio	53.00	67.000	14.00	1.0000		0.580	1.70
hct_score	0.00	3.000	3.00	0.3000	0.42	-0.530	1.10
Odds Ratio	0.00	3.000	3.00	1.4000		0.590	3.10
preleuka_ALC_l	-0.37	0.017	0.38	0.8800	0.35	0.200	1.60
Odds Ratio	-0.37	0.017	0.38	2.4000		1.200	4.80
product --- tisacel:axicel	1.00	2.000		-1.5000	0.64	-2.700	-0.24
Odds Ratio	1.00	2.000		0.2300		0.065	0.78
product --- jcar014:axicel	1.00	3.000		-1.5000	0.63	-2.800	-0.32
Odds Ratio	1.00	3.000		0.2100		0.062	0.73

Show code

print(anova(fit_mi),digits = 2)

	χ²	d.f.	P
Wald Statistics for `bestCR_b`
product	8.00	2	0.0184
preleuka_LDH_l	12.52	1	0.0004
preleuka_bulk_s	4.86	1	0.0275
age	0.00	1	0.9815
hct_score	0.51	1	0.4768
preleuka_ALC_l	6.37	1	0.0116
TOTAL	23.17	7	0.0016

Show code

plot(anova(fit_mi))

Show code

plot(calibrate(fit_mi))

n=126 Mean absolute error=0.058 Mean squared error=0.00436 0.9 Quantile of absolute error=0.102

Show code

options(prType=NULL)

#Response by CT only criteria
fit_mi_CT <-
  fit.mult.impute(
    bestCR_CT_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       1.02            1.08            1.03            1.08

preleuka_bulk_s age hct_score preleuka_ALC_l 2.28 1.02 1.06 1.08

Rate of Missing Information:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       0.02            0.08            0.03            0.08

preleuka_bulk_s age hct_score preleuka_ALC_l 0.56 0.02 0.06 0.07

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
   15691.91         1500.05         8307.83         1519.90

preleuka_bulk_s age hct_score preleuka_ALC_l 28.48 15825.25 2855.76 1831.22

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

#Response by PET only criteria
fit_mi_PET <-
  fit.mult.impute(
    bestCR_PET_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       1.01            1.05            1.05            1.04

preleuka_bulk_s age hct_score preleuka_ALC_l 1.42 1.02 1.02 1.02

Rate of Missing Information:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       0.01            0.05            0.05            0.04

preleuka_bulk_s age hct_score preleuka_ALC_l 0.30 0.02 0.02 0.02

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel product=jcar014  preleuka_LDH_l 
  107641.60         3721.71         3619.76         7228.56

preleuka_bulk_s age hct_score preleuka_ALC_l 103.34 28962.56 17309.33 15553.39

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

CR_t <- resp_table_sparser(fit_uni,fit_mi)
CR_CT_t <- resp_table_sparser(fit_uniCT,fit_mi_CT)
CR_PET_t <- resp_table_sparser(fit_uniPET,fit_mi_PET)

Figure 3. Forest plot of predictors of complete response after CD19 CAR T-cell therapy

Show code

newnames <- c("Preleukapheresis LDH", "Preleukapheresis largest lesion diameter", "Age", "HCT-CI", "Preleukapheresis ALC","tisacel versus axicel","JCAR014 versus axicel")
newhead <- c("Odds Ratio", "Lower CI", "Upper CI" )
forestplot(summary(fit_mi),digits=2,title="OR for complete response (best response by PET or CT)",row.names.y = newnames,header = newhead)

Show code

#Save as .eps
# cairo_ps("manuscript/Blood_resub2/Fig3.eps", width = 12, height = 4)
# newnames <- c("Preleukapheresis LDH", "Preleukapheresis largest lesion diameter", "Age", "HCT-CI", "Preleukapheresis ALC","tisacel versus axicel","JCAR014 versus axicel")
# newhead <- c("Odds Ratio", "Lower CI", "Upper CI" )
# p <- forestplot(summary(fit_mi),digits=2,title="OR for complete response (best response by PET or CT)",row.names.y = newnames,header = newhead)
# grid.arrange(p,top = grid::textGrob("Figure 3.", x = 0, hjust = 0,gp=grid::gpar(fontface="bold")))
# dev.off()

Sensitivity analysis: univariate and multivariable logistic regression for CR (best response) excluding JCAR014

Show code

options(datadist='dd2')
#Univariate
fit_resp_sub <- lrm(bestCR_b ~ product,x=T,y=T,data=df2)
fit_resp_subCT <- lrm(bestCR_CT_b ~ product,x=T,y=T,data=df2)
fit_resp_subPET <- lrm(bestCR_PET_b ~ product,x=T,y=T,data=df2)
#Multivariable
fit_mi_sub <-
  fit.mult.impute(
    bestCR_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df2,
    xtrans = mi2,
    x = T,
    y = T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       1.09            1.17            1.02            1.89 
        age       hct_score  preleuka_ALC_l 
       1.12            1.10            1.09

Rate of Missing Information:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       0.08            0.15            0.02            0.47 
        age       hct_score  preleuka_ALC_l 
       0.10            0.09            0.08

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
    1326.89          417.29        22210.14           40.53 
        age       hct_score  preleuka_ALC_l 
     836.48         1070.79         1332.40

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

fit_mi_subCT <-
  fit.mult.impute(
    bestCR_CT_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df2,
    xtrans = mi2,
    x = T,
    y = T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       1.01            1.08            1.06            1.64 
        age       hct_score  preleuka_ALC_l 
       1.02            1.06            1.03

Rate of Missing Information:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       0.01            0.07            0.06            0.39 
        age       hct_score  preleuka_ALC_l 
       0.02            0.06            0.03

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
   48727.43         1634.76         2964.14           59.38 
        age       hct_score  preleuka_ALC_l 
   37019.39         2618.45        11489.67

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

fit_mi_subPET <-
  fit.mult.impute(
    bestCR_PET_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df2,
    xtrans = mi2,
    x = T,
    y = T
  )

Variance Inflation Factors Due to Imputation:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       1.02            1.09            1.04            2.11 
        age       hct_score  preleuka_ALC_l 
       1.03            1.17            1.05

Rate of Missing Information:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
       0.02            0.08            0.03            0.53 
        age       hct_score  preleuka_ALC_l 
       0.03            0.14            0.05

d.f. for t-distribution for Tests of Single Coefficients:

  Intercept product=tisacel  preleuka_LDH_l preleuka_bulk_s 
   16015.41         1312.26         7651.04           32.64 
        age       hct_score  preleuka_ALC_l 
    8398.90          430.87         4272.16

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

#Build table for resp PET or CT
#Parse ORs and 95% CIs from univariate model
d <- summary(fit_resp_sub)
d <- data.frame(d)
resp_unadj <- cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-"))
colnames(resp_unadj) <- c("unadjusted OR","95% CI")
rownames(resp_unadj) <- c("Tisacel versus axicel")
#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_resp_sub)
#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi_sub)
#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi_sub))
resp_adj <- cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-"))
colnames(resp_adj) <- c("adjusted OR","95% CI")
rownames(resp_adj) <- c("Tisacel versus axicel")
#Create table object
CR_sub_t <- cbind(resp_unadj,unadjstats_d$coefs[2,5],resp_adj,adjstats_d$coefs[2,5])

#Build table for CR by CT only
#Parse ORs and 95% CIs from univariate model
d <- summary(fit_resp_subCT)
d <- data.frame(d)
resp_unadj <- cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-"))
colnames(resp_unadj) <- c("unadjusted OR","95% CI")
rownames(resp_unadj) <- c("Tisacel versus axicel")
#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_resp_subCT)
#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi_subCT)
#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi_subCT))
resp_adj <- cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-"))
colnames(resp_adj) <- c("adjusted OR","95% CI")
rownames(resp_adj) <- c("Tisacel versus axicel")

CR_CT_sub_t <- cbind(resp_unadj,unadjstats_d$coefs[2,5],resp_adj,adjstats_d$coefs[2,5])

#Build table for CR by PET only
#Parse ORs and 95% CIs from univariate model
d <- summary(fit_resp_subPET)
d <- data.frame(d)
resp_unadj <- cbind(round(d[2,4],2),paste(round(d[2,6],2),round(d[2,7],2),sep="-"))
colnames(resp_unadj) <- c("unadjusted OR","95% CI")
rownames(resp_unadj) <- c("Tisacel versus axicel")
#Parse p values from univariate model
unadjstats_d <- get_model_stats(fit_resp_subPET)
#Parse p values from multivariable model
adjstats_d <- get_model_stats(fit_mi_subPET)
#Parse ORs and 95% CI from multivariable model
d <- data.frame(summary(fit_mi_subPET))
resp_adj <- cbind(round(d[12,4],2),paste(round(d[12,6],2),round(d[12,7],2),sep="-"))
colnames(resp_adj) <- c("adjusted OR","95% CI")
rownames(resp_adj) <- c("Tisacel versus axicel")

CR_PET_sub_t <- cbind(resp_unadj,unadjstats_d$coefs[2,5],resp_adj,adjstats_d$coefs[2,5])

Sensitivity analysis: univariate and multivariable logistic regression for CR (best response) including bridging therapy

Show code

options(datadist='dd')
#Univariate
fit_resp <- lrm(bestCR_b ~ product,x=T,y=T,data=df)
fit_resp_CT <- lrm(bestCR_CT_b ~ product,x=T,y=T,data=df)
fit_resp_PET <- lrm(bestCR_PET_b ~ product,x=T,y=T,data=df)
#Multivariate
fit_mi <-
  fit.mult.impute(
    bestCR_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.04            1.20            1.10            1.02 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           1.84            1.06            1.01            1.09 
     bridge=yes 
           1.02 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.04            0.17            0.09            0.02 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           0.46            0.06            0.01            0.08 
     bridge=yes 
           0.02 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
        5106.59          313.23         1110.74        21678.25 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
          43.39         2553.87        64623.46         1353.93 
     bridge=yes 
       31258.08 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

fit_mi_CT <-
  fit.mult.impute(
    bestCR_CT_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.02            1.09            1.03            1.06 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           2.29            1.03            1.06            1.08 
     bridge=yes 
           1.02 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.02            0.08            0.03            0.06 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           0.56            0.02            0.06            0.07 
     bridge=yes 
           0.02 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       16244.05         1433.86        11439.88         2555.45 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
          28.43        14845.87         2710.58         1693.66 
     bridge=yes 
       26016.62 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

fit_mi_PET <-
  fit.mult.impute(
    bestCR_PET_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.01            1.05            1.05            1.03 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           1.42            1.02            1.02            1.02 
     bridge=yes 
           1.01 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.01            0.05            0.04            0.03 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           0.30            0.02            0.02            0.02 
     bridge=yes 
           0.01 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
      113205.35         3763.88         4606.43        10649.11 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
         103.40        32365.33        18331.27        16591.95 
     bridge=yes 
      159634.33 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

resp_bridge_t <- resp_table_sparser(fit_resp,fit_mi)
respCT_bridge_t <- resp_table_sparser(fit_resp_CT,fit_mi_CT)
respPET_bridge_t <- resp_table_sparser(fit_resp_PET,fit_mi_PET)

Sensitivity analysis: univariate and multivariable logistic regression for CR (best response) including preLD ALC instead of preleukapheresis ALC

Show code

options(datadist='dd')
#Univariate
fit_resp <- lrm(bestCR_b ~ product,x=T,y=T,data=df)
fit_resp_CT <- lrm(bestCR_CT_b ~ product,x=T,y=T,data=df)
fit_resp_PET <- lrm(bestCR_PET_b ~ product,x=T,y=T,data=df)
#Multivariate
fit_mi <-
  fit.mult.impute(
    bestCR_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.03            1.16            1.06            1.03 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           1.65            1.03            1.02            1.06 
     bridge=yes 
           1.02 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.03            0.14            0.05            0.03 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           0.39            0.03            0.02            0.05 
     bridge=yes 
           0.02 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       11862.87          470.17         3159.99        12088.25 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
          57.77         9194.02        25967.40         3197.65 
     bridge=yes 
       22825.73 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

fit_mi_CT <-
  fit.mult.impute(
    bestCR_CT_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.03            1.12            1.04            1.07 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           2.23            1.03            1.07            1.12 
     bridge=yes 
           1.05 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.03            0.10            0.04            0.07 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           0.55            0.03            0.06            0.11 
     bridge=yes 
           0.04 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
        8651.72          843.65         5552.07         2074.23 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
          29.58        14337.68         2389.11          747.58 
     bridge=yes 
        4527.55 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

fit_mi_PET <-
  fit.mult.impute(
    bestCR_PET_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preLD_ALC_l + bridge,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.01            1.04            1.03            1.03 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           1.33            1.01            1.02            1.03 
     bridge=yes 
           1.01 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.01            0.04            0.03            0.03 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
           0.25            0.01            0.02            0.03 
     bridge=yes 
           0.01 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
       48157.34         7079.25        10832.01         7991.64 
preleuka_bulk_s             age       hct_score     preLD_ALC_l 
         145.75       130554.37        17635.44        12719.07 
     bridge=yes 
       88352.05 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

resp_preLDALC_t <- resp_table_sparser(fit_resp,fit_mi)
respCT_preLDALC_t <- resp_table_sparser(fit_resp_CT,fit_mi_CT)
respPET_preLDALC_t <- resp_table_sparser(fit_resp_PET,fit_mi_PET)

Sensitivity analysis: univariate and multivariable logistic regression model for ORR (best response - CT criteria only)

Show code

options(datadist='dd')
fit_uni <-
  fit.mult.impute(
    bestresp_CT_b ~ product,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014 
              1               1               1 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014 
              0               0               0 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014 
   1.333025e+31    1.164304e+32    1.333328e+32 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

print(fit_uni)

Frequencies of Missing Values Due to Each Variable
bestresp_CT_b       product 
            5             0 

Logistic Regression Model
 
 fit.mult.impute(formula = bestresp_CT_b ~ product, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)
 
 
                       Model Likelihood    Discrimination    Rank Discrim.    
                             Ratio Test           Indexes          Indexes    
 Obs           124    LR chi2      1.91    R2       0.021    C       0.564    
  0             53    d.f.            2    g        0.260    Dxy     0.128    
  1             71    Pr(> chi2) 0.3851    gr       1.296    gamma   0.206    
 max |deriv| 3e-11                         gp       0.063    tau-a   0.063    
                                           Brier    0.241                     
 
                 Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept        0.5355 0.2570  2.08  0.0372  
 product=tisacel -0.4710 0.4418 -1.07  0.2864  
 product=jcar014 -0.5355 0.4571 -1.17  0.2413

Show code

fit_mi <-
  fit.mult.impute(
    bestresp_CT_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.03            1.06            1.02            1.03 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           1.28            1.06            1.01            1.01 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.03            0.05            0.02            0.03 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           0.22            0.06            0.01            0.01 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
        8765.85         3022.36        26363.77         9885.07 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
         184.87         2597.69        88191.68        44275.44 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for CR (best response - CT criteria only)"
)

Multivariable ordinal proportional odds logistic regression model for CR (best response - CT criteria only)

 fit.mult.impute(formula = bestresp_CT_b ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

   bestresp_CT_b         product  preleuka_LDH_l preleuka_bulk_s 
               5               0               0               0 
             age       hct_score  preleuka_ALC_l 
               0               0               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 124	LR χ² 21.01	R² 0.209	C 0.715
0 53	d.f. 7	g 1.076	D_xy 0.430
1 71	Pr(>χ²) 0.0052	g_r 2.940	γ 0.430
max \|∂log L/∂β\| 2×10^-6		g_p 0.225	τ_a 0.212
		Brier 0.209

	β	S.E.	Wald Z	Pr(>\|Z\|)
Intercept	5.10	2.52	2.02	0.0431
product=tisacel	-0.53	0.51	-1.04	0.2970
product=jcar014	-0.93	0.51	-1.80	0.0712
preleuka_LDH_l	-1.30	0.91	-1.43	0.1526
preleuka_bulk_s	-0.16	0.09	-1.82	0.0686
age	0.00	0.02	-0.29	0.7725
hct_score	0.06	0.12	0.49	0.6240
preleuka_ALC_l	1.90	0.73	2.59	0.0097

Show code

print(summary(fit_mi,age=c(50,60),preleuka_LDH_l=c(2,2.3),preleuka_bulk_s=c(1,2),preleuka_ALC_l=c(2,2.3)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `bestresp_CT_b`
preleuka_LDH_l	2	2.3	0.3	-0.390	0.270	-0.92	0.140
Odds Ratio	2	2.3	0.3	0.680		0.40	1.200
preleuka_bulk_s	1	2.0	1.0	-0.160	0.085	-0.32	0.012
Odds Ratio	1	2.0	1.0	0.860		0.72	1.000
age	50	60.0	10.0	-0.049	0.170	-0.38	0.280
Odds Ratio	50	60.0	10.0	0.950		0.69	1.300
hct_score	0	3.0	3.0	0.180	0.360	-0.54	0.890
Odds Ratio	0	3.0	3.0	1.200		0.59	2.400
preleuka_ALC_l	2	2.3	0.3	0.570	0.220	0.14	1.000
Odds Ratio	2	2.3	0.3	1.800		1.10	2.700
product --- tisacel:axicel	1	2.0		-0.530	0.510	-1.50	0.460
Odds Ratio	1	2.0		0.590		0.22	1.600
product --- jcar014:axicel	1	3.0		-0.930	0.510	-1.90	0.080
Odds Ratio	1	3.0		0.400		0.14	1.100

Show code

print(anova(fit_mi),digits = 2)

	χ²	d.f.	P
Wald Statistics for `bestresp_CT_b`
product	3.50	2	0.1741
preleuka_LDH_l	2.05	1	0.1526
preleuka_bulk_s	3.32	1	0.0686
age	0.08	1	0.7725
hct_score	0.24	1	0.6240
preleuka_ALC_l	6.70	1	0.0097
TOTAL	14.90	7	0.0372

Show code

options(prType=NULL)

ORR_CT_t <- resp_table_sparser(fit_uni,fit_mi)
ORR_CT_t

                      unadjusted OR    95% CI Pr(>|Z|) adjusted OR
Tisacel versus axicel          0.62 0.26-1.48   0.2864        0.59
JCAR014 versus axicel          0.59 0.24-1.43   0.2413         0.4
                         95% CI Pr(>|Z|)
Tisacel versus axicel 0.22-1.59   0.2970
JCAR014 versus axicel 0.14-1.08   0.0712

Sensitivity analysis: univariate and multivariable logistic regression model for ORR (best response - PET criteria only)

Show code

options(datadist='dd')
fit <-
  fit.mult.impute(
    bestresp_PET_b ~ product,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014 
              1               1               1 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014 
              0               0               0 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014 
   1.305893e+29             Inf    6.396243e+30 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

print(fit)

Frequencies of Missing Values Due to Each Variable
bestresp_PET_b        product 
            25              0 

Logistic Regression Model
 
 fit.mult.impute(formula = bestresp_PET_b ~ product, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)
 
 
                       Model Likelihood    Discrimination    Rank Discrim.    
                             Ratio Test           Indexes          Indexes    
 Obs           104    LR chi2      2.27    R2       0.032    C       0.587    
  0             27    d.f.            2    g        0.345    Dxy     0.173    
  1             77    Pr(> chi2) 0.3208    gr       1.413    gamma   0.280    
 max |deriv| 1e-07                         gp       0.067    tau-a   0.067    
                                           Brier    0.188                     
 
                 Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept        1.3652 0.3234  4.22  <0.0001 
 product=tisacel -0.5921 0.5901 -1.00  0.3157  
 product=jcar014 -0.7293 0.5240 -1.39  0.1640

Show code

fit_mi <-
  fit.mult.impute(
    bestresp_PET_b ~ product + preleuka_LDH_l + preleuka_bulk_s + age + hct_score + preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )


Variance Inflation Factors Due to Imputation:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           1.00            1.01            1.01            1.01 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           1.16            1.01            1.01            1.00 

Rate of Missing Information:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
           0.00            0.01            0.01            0.01 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
           0.14            0.01            0.01            0.00 

d.f. for t-distribution for Tests of Single Coefficients:

      Intercept product=tisacel product=jcar014  preleuka_LDH_l 
     1865128.84       122517.38       163454.85        50141.61 
preleuka_bulk_s             age       hct_score  preleuka_ALC_l 
         458.03       279304.57        79503.41      1314882.86 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T,
  title = "Multivariable ordinal proportional odds logistic regression model for CR (best response - PET criteria only)"
)

Multivariable ordinal proportional odds logistic regression model for CR (best response - PET criteria only)

 fit.mult.impute(formula = bestresp_PET_b ~ product + preleuka_LDH_l + 
     preleuka_bulk_s + age + hct_score + preleuka_ALC_l, fitter = lrm, 
     xtrans = mi, data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

  bestresp_PET_b         product  preleuka_LDH_l preleuka_bulk_s 
              25               0               0               0 
             age       hct_score  preleuka_ALC_l 
               0               0               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 104	LR χ² 22.23	R² 0.282	C 0.781
0 27	d.f. 7	g 1.420	D_xy 0.563
1 77	Pr(>χ²) 0.0025	g_r 4.139	γ 0.563
max \|∂log L/∂β\| 3×10^-5		g_p 0.223	τ_a 0.218
		Brier 0.150

	β	S.E.	Wald Z	Pr(>\|Z\|)
Intercept	8.56	3.32	2.58	0.0099
product=tisacel	-1.14	0.71	-1.62	0.1060
product=jcar014	-1.41	0.63	-2.25	0.0243
preleuka_LDH_l	-2.02	1.15	-1.77	0.0774
preleuka_bulk_s	-0.11	0.11	-0.99	0.3203
age	-0.01	0.02	-0.42	0.6751
hct_score	-0.17	0.14	-1.15	0.2504
preleuka_ALC_l	3.10	1.03	3.00	0.0027

Show code

print(summary(fit_mi,age=c(50,60),preleuka_LDH_l=c(2,2.3),preleuka_bulk_s=c(1,2),preleuka_ALC_l=c(2,2.3)),digits = 2)

	Low	High	Δ	Effect	S.E.	Lower 0.95	Upper 0.95
Effects Response: `bestresp_PET_b`
preleuka_LDH_l	2	2.3	0.3	-0.610	0.34	-1.300	0.067
Odds Ratio	2	2.3	0.3	0.540		0.280	1.100
preleuka_bulk_s	1	2.0	1.0	-0.110	0.11	-0.320	0.100
Odds Ratio	1	2.0	1.0	0.900		0.730	1.100
age	50	60.0	10.0	-0.096	0.23	-0.540	0.350
Odds Ratio	50	60.0	10.0	0.910		0.580	1.400
hct_score	0	3.0	3.0	-0.500	0.43	-1.300	0.350
Odds Ratio	0	3.0	3.0	0.610		0.260	1.400
preleuka_ALC_l	2	2.3	0.3	0.930	0.31	0.320	1.500
Odds Ratio	2	2.3	0.3	2.500		1.400	4.600
product --- tisacel:axicel	1	2.0		-1.100	0.71	-2.500	0.240
Odds Ratio	1	2.0		0.320		0.080	1.300
product --- jcar014:axicel	1	3.0		-1.400	0.63	-2.600	-0.180
Odds Ratio	1	3.0		0.240		0.071	0.830

Show code

print(anova(fit_mi),digits = 2)

	χ²	d.f.	P
Wald Statistics for `bestresp_PET_b`
product	5.65	2	0.0593
preleuka_LDH_l	3.12	1	0.0774
preleuka_bulk_s	0.99	1	0.3203
age	0.18	1	0.6751
hct_score	1.32	1	0.2504
preleuka_ALC_l	9.00	1	0.0027
TOTAL	15.88	7	0.0262

Show code

options(prType=NULL)

ORR_PET_t <- resp_table_sparser(fit_uni,fit_mi)
ORR_PET_t

                      unadjusted OR    95% CI Pr(>|Z|) adjusted OR
Tisacel versus axicel          0.62 0.26-1.48   0.2864        0.32
JCAR014 versus axicel          0.59 0.24-1.43   0.2413        0.24
                         95% CI Pr(>|Z|)
Tisacel versus axicel 0.08-1.28   0.1060
JCAR014 versus axicel 0.07-0.83   0.0243

Table 3. CAR T-cell product type and outcome prediction

Show code

rbind(CRS_t, NT_t, CR_t, CR_CT_t, CR_PET_t) %>%
  addHtmlTableStyle(pos.caption = "bottom") %>%
  htmlTable(
    rgroup = c(
      "CRS grade*",
      "ICANS grade*",
      "CR (Best Response)^",
      "CR (Best Response - CT criteria only)^",
      "CR (Best Response - PET criteria only)^"
    ),
    n.rgroup = c(2, 2, 2, 2, 2),
    caption = "*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC. P values per the Wald test."
  )

	unadjusted OR	95% CI	Pr(>\|Z\|)	adjusted OR	95% CI	Pr(>\|Z\|)
CRS grade*
Tisacel versus axicel	0.43	0.2-0.95	0.0375	0.47	0.21-1.06	0.0698
JCAR014 versus axicel	0.19	0.08-0.46	0.0002	0.19	0.08-0.46	2e-04
ICANS grade*
Tisacel versus axicel1	0.22	0.08-0.56	0.0017	0.17	0.06-0.48	8e-04
JCAR014 versus axicel1	0.18	0.07-0.48	7e-04	0.17	0.06-0.47	6e-04
CR (Best Response)^
Tisacel versus axicel2	0.38	0.16-0.94	0.0365	0.23	0.06-0.78	0.0191
JCAR014 versus axicel2	0.54	0.22-1.29	0.1657	0.21	0.06-0.73	0.0138
CR (Best Response - CT criteria only)^
Tisacel versus axicel3	0.89	0.32-2.46	0.8272	0.78	0.24-2.52	0.6781
JCAR014 versus axicel3	0.37	0.1-1.38	0.1383	0.25	0.06-1.06	0.0601
CR (Best Response - PET criteria only)^
Tisacel versus axicel4	0.66	0.23-1.87	0.4357	0.42	0.11-1.54	0.1886
JCAR014 versus axicel4	0.63	0.25-1.59	0.3296	0.3	0.09-1	0.05
*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC. P values per the Wald test.

Table S4. Comparison of preleukapheresis and prelymphodepletion tumor burden and ALC

Show code

#Descriptive statistics
df %>%
  select(patient_id,
         `Preleukapheresis` = preleuka_LDH,
         `Prelymphodepletion` = preLD_LDH) %>%
  pivot_longer(
    cols = c(`Preleukapheresis`, `Prelymphodepletion`),
    names_to = "TP",
    values_to = "LDH"
  ) %>%
  left_join(
    .,
    df %>%
      select(
        patient_id,
        `Preleukapheresis` = preleuka_bulk_s,
        `Prelymphodepletion` = preLD_bulk_s
      ) %>%
      pivot_longer(
        cols = c(`Preleukapheresis`, `Prelymphodepletion`),
        names_to = "TP",
        values_to = "Largest lesion diameter"
      ),
    by=c("patient_id","TP")) %>%
  left_join(
    .,
    df %>%
      select(
        patient_id,
        `Preleukapheresis` = preleuka_ALC,
        `Prelymphodepletion` = preLD_ALC
      ) %>%
      pivot_longer(
        cols = c(`Preleukapheresis`, `Prelymphodepletion`),
        names_to = "TP",
        values_to = "ALC"
      ),
    by=c("patient_id","TP")) %>% 
  select(-patient_id) %>% 
  tbl_summary(
    by="TP",
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c("{median} ({p25}, {p75})", 
                                     "{min}, {max}"),
    missing_text = "Missing data",
    sort = all_categorical() ~ "frequency") %>%
  add_p(all_continuous() ~ "wilcox.test",
             simulate.p.value=TRUE) %>% 
  bold_labels()

Characteristic	Preleukapheresis, N = 129	Prelymphodepletion, N = 129	p-value¹
LDH			>0.9
Median (IQR)	227 (174, 366)	238 (164, 344)
Range	105, 2,812	103, 2,371
Missing data	0	2
Largest lesion diameter			0.2
Median (IQR)	4.3 (2.5, 6.7)	4.7 (2.8, 7.5)
Range	0.8, 24.8	0.8, 23.5
Missing data	24	16
ALC			0.010
Median (IQR)	0.71 (0.43, 1.04)	0.59 (0.31, 0.94)
Range	0.11, 3.82	0.00, 2.22
Missing data	0	1
¹ Wilcoxon rank sum test

Table S5. Sensitivity analysis including bridging therapy in addition to the minimal adjustment set of covariates

Show code

rbind(CRS_bridge_t,
      NT_bridge_t,
      resp_bridge_t,
      respCT_bridge_t,
      respPET_bridge_t) %>%
  addHtmlTableStyle(pos.caption = "bottom") %>%
  htmlTable(
    rgroup = c("CRS grade*", "ICANS grade*", "CR (Best Response)^","CR (Best Response - CT criteria only)^","CR (Best Response - PET criteria only)^"),
    n.rgroup = c(2, 2, 2, 2, 2),
    caption = "*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC, bridging therapy post leukapheresis. P values per the Wald test."
  )

	unadjusted OR	95% CI	Pr(>\|Z\|)	adjusted OR	95% CI	Pr(>\|Z\|)
CRS grade*
Tisacel versus axicel	0.43	0.2-0.95	0.0375	0.46	0.2-1.06	0.0681
JCAR014 versus axicel	0.19	0.08-0.46	0.0002	0.2	0.08-0.49	5e-04
ICANS grade*
Tisacel versus axicel1	0.22	0.08-0.56	0.0017	0.16	0.06-0.46	6e-04
JCAR014 versus axicel1	0.18	0.07-0.48	7e-04	0.19	0.07-0.55	0.0022
CR (Best Response)^
Tisacel versus axicel2	0.38	0.16-0.94	0.0365	0.22	0.06-0.77	0.018
JCAR014 versus axicel2	0.54	0.22-1.29	0.1657	0.23	0.06-0.81	0.0228
CR (Best Response - CT criteria only)^
Tisacel versus axicel3	0.89	0.32-2.46	0.8272	0.78	0.24-2.57	0.6803
JCAR014 versus axicel3	0.37	0.1-1.38	0.1383	0.25	0.06-1.1	0.0673
CR (Best Response - PET criteria only)^
Tisacel versus axicel4	0.66	0.23-1.87	0.4357	0.42	0.11-1.56	0.1963
JCAR014 versus axicel4	0.63	0.25-1.59	0.3296	0.28	0.08-1.01	0.0513
*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC, bridging therapy post leukapheresis. P values per the Wald test.

Table S6. Sensitivity analysis excluding JCAR014 patients

Show code

rbind(CRS_sub_t, NT_sub_t, CR_sub_t, CR_CT_sub_t, CR_PET_sub_t) %>%
  addHtmlTableStyle(pos.caption = "bottom") %>%
  htmlTable(
    rgroup = c(
      "CRS grade*",
      "ICANS grade*",
      "CR (Best Response)^",
      "CR (Best Response - CT criteria only)^",
      "CR (Best Response - PET criteria only)^"
    ),
    n.rgroup = c(1, 1, 1, 1, 1),
    caption = "*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC, bridging therapy post leukapheresis. P values per the Wald test."
  )

	unadjusted OR	95% CI	Pr(>\|Z\|)	adjusted OR	95% CI	Pr(>\|Z\|)
CRS grade*
Tisacel versus axicel	0.41	0.18-0.92	0.0306	0.43	0.19-0.97	0.0433
ICANS grade*
Tisacel versus axicel1	0.21	0.08-0.56	0.0016	0.16	0.06-0.45	5e-04
CR (Best Response)^
Tisacel versus axicel2	0.38	0.16-0.94	0.0365	0.16	0.04-0.62	0.0083
CR (Best Response - CT criteria only)^
Tisacel versus axicel3	0.89	0.32-2.46	0.8272	0.68	0.21-2.2	0.5216
CR (Best Response - PET criteria only)^
Tisacel versus axicel4	0.66	0.23-1.87	0.4357	0.3	0.07-1.29	0.1059
*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, preleukapheresis ALC, bridging therapy post leukapheresis. P values per the Wald test.

Table S7. Sensitivity analysis including prelymphodepletion ALC instead of preleukapheresis ALC

Show code

rbind(CRS_preLDALC_t,
      NT_preLDALC_t,
      resp_preLDALC_t,
      respCT_preLDALC_t,
      respPET_preLDALC_t) %>%
  addHtmlTableStyle(pos.caption = "bottom") %>%
  htmlTable(
    rgroup = c(
      "CRS grade*",
      "ICANS grade*",
      "CR (Best Response)^",
      "CR (Best Response - CT criteria only)^",
      "CR (Best Response - PET criteria only)^"
    ),
    n.rgroup = c(2, 2, 2, 2, 2),
    caption = "*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, prelymphodepletion ALC, bridging therapy post leukapheresis. P values per the Wald test."
  )

	unadjusted OR	95% CI	Pr(>\|Z\|)	adjusted OR	95% CI	Pr(>\|Z\|)
CRS grade*
Tisacel versus axicel	0.43	0.2-0.95	0.0375	0.47	0.2-1.08	0.0755
JCAR014 versus axicel	0.19	0.08-0.46	0.0002	0.2	0.08-0.51	7e-04
ICANS grade*
Tisacel versus axicel1	0.22	0.08-0.56	0.0017	0.18	0.06-0.49	0.001
JCAR014 versus axicel1	0.18	0.07-0.48	7e-04	0.2	0.07-0.58	0.0028
CR (Best Response)^
Tisacel versus axicel2	0.38	0.16-0.94	0.0365	0.26	0.08-0.89	0.0316
JCAR014 versus axicel2	0.54	0.22-1.29	0.1657	0.28	0.09-0.91	0.0343
CR (Best Response - CT criteria only)^
Tisacel versus axicel3	0.89	0.32-2.46	0.8272	0.8	0.24-2.71	0.7202
JCAR014 versus axicel3	0.37	0.1-1.38	0.1383	0.24	0.05-1.1	0.0664
CR (Best Response - PET criteria only)^
Tisacel versus axicel4	0.66	0.23-1.87	0.4357	0.5	0.14-1.79	0.2865
JCAR014 versus axicel4	0.63	0.25-1.59	0.3296	0.37	0.11-1.18	0.0935
*From a multivariable proportional odds logistic regression model or ^logistic regression model including the following variables: CAR T-cell product type, preleukapheresis LDH, preleukapheresis largest lesion diameter, age, HCT-CI, prelymphodepletion ALC, bridging therapy post leukapheresis. P values per the Wald test.

Multivariable ordinal logistic regression model for CRS grade including interaction effects

Show code

#Merge tisacel and JCAR014
df <- df %>% 
  mutate(product_s=ifelse(product=="axicel","axicel","tisacel/JCAR014"))
dd <- datadist(df)
options(datadist='dd')

#### Refit model for CRS with interaction terms
fit_mi <-
  fit.mult.impute(
    crs_g_r ~ product_s*preleuka_LDH_l + product_s*preleuka_bulk_s + product_s*age + product_s*hct_score + product_s*preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,
    y=T
  )

Variance Inflation Factors Due to Imputation:

                                   y>=1 
                                   1.00 
                                 y>=2-5 
                                   1.00 
              product_s=tisacel/JCAR014 
                                   1.02 
                         preleuka_LDH_l 
                                   1.01 
                        preleuka_bulk_s 
                                   1.09 
                                    age 
                                   1.00 
                              hct_score 
                                   1.00 
                         preleuka_ALC_l 
                                   1.00

product_s=tisacel/JCAR014 * preleuka_LDH_l 1.02 product_s=tisacel/JCAR014 * preleuka_bulk_s 1.11 product_s=tisacel/JCAR014 * age 1.01 product_s=tisacel/JCAR014 * hct_score 1.00 product_s=tisacel/JCAR014 * preleuka_ALC_l 1.01

Rate of Missing Information:

                                   y>=1 
                                   0.00 
                                 y>=2-5 
                                   0.00 
              product_s=tisacel/JCAR014 
                                   0.02 
                         preleuka_LDH_l 
                                   0.01 
                        preleuka_bulk_s 
                                   0.08 
                                    age 
                                   0.00 
                              hct_score 
                                   0.00 
                         preleuka_ALC_l 
                                   0.00

product_s=tisacel/JCAR014 * preleuka_LDH_l 0.02 product_s=tisacel/JCAR014 * preleuka_bulk_s 0.10 product_s=tisacel/JCAR014 * age 0.01 product_s=tisacel/JCAR014 * hct_score 0.00 product_s=tisacel/JCAR014 * preleuka_ALC_l 0.01

d.f. for t-distribution for Tests of Single Coefficients:

                                   y>=1 
                           1.304235e+09 
                                 y>=2-5 
                           1.289710e+09 
              product_s=tisacel/JCAR014 
                           3.982849e+04 
                         preleuka_LDH_l 
                           2.146590e+05 
                        preleuka_bulk_s 
                           1.411540e+03 
                                    age 
                           6.159073e+05 
                              hct_score 
                           6.849014e+05 
                         preleuka_ALC_l 
                           1.765304e+09

product_s=tisacel/JCAR014 * preleuka_LDH_l 1.637265e+04 product_s=tisacel/JCAR014 * preleuka_bulk_s 9.072100e+02 product_s=tisacel/JCAR014 * age 3.208591e+05 product_s=tisacel/JCAR014 * hct_score 2.359576e+06 product_s=tisacel/JCAR014 * preleuka_ALC_l 1.795134e+05

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = crs_g_r ~ product_s * preleuka_LDH_l + 
     product_s * preleuka_bulk_s + product_s * age + product_s * 
     hct_score + product_s * preleuka_ALC_l, fitter = lrm, xtrans = mi, 
     data = df, x = T, y = T)

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 20.15	R² 0.163	C 0.691
0 34	d.f. 11	g 0.888	D_xy 0.381
1 46	Pr(>χ²) 0.0436	g_r 2.429	γ 0.381
2-5 49		g_p 0.160	τ_a 0.253
max \|∂log L/∂β\| 3×10^-9		Brier 0.170

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-1.02	2.74	-0.37	0.7106
y≥2-5	-2.74	2.74	-1.00	0.3180
product_s=tisacel/JCAR014	-0.44	4.55	-0.10	0.9231
preleuka_LDH_l	1.27	0.96	1.32	0.1863
preleuka_bulk_s	0.02	0.09	0.18	0.8602
age	0.00	0.02	-0.21	0.8323
hct_score	-0.06	0.13	-0.47	0.6387
preleuka_ALC_l	0.44	0.71	0.63	0.5318
product_s=tisacel/JCAR014 × preleuka_LDH_l	-0.02	1.76	-0.01	0.9907
product_s=tisacel/JCAR014 × preleuka_bulk_s	-0.05	0.12	-0.46	0.6472
product_s=tisacel/JCAR014 × age	-0.01	0.03	-0.48	0.6310
product_s=tisacel/JCAR014 × hct_score	0.25	0.20	1.21	0.2248
product_s=tisacel/JCAR014 × preleuka_ALC_l	0.17	1.27	0.13	0.8953

Show code

print(anova(fit_mi))

	χ²	d.f.	P
Wald Statistics for `crs_g_r`
product_s (Factor+Higher Order Factors)	13.09	6	0.0417
All Interactions	2.21	5	0.8189
preleuka_LDH_l (Factor+Higher Order Factors)	2.45	2	0.2942
All Interactions	0.00	1	0.9907
preleuka_bulk_s (Factor+Higher Order Factors)	0.25	2	0.8826
All Interactions	0.21	1	0.6472
age (Factor+Higher Order Factors)	0.80	2	0.6688
All Interactions	0.23	1	0.6310
hct_score (Factor+Higher Order Factors)	1.66	2	0.4359
All Interactions	1.47	1	0.2248
preleuka_ALC_l (Factor+Higher Order Factors)	0.73	2	0.6953
All Interactions	0.02	1	0.8953
product_s × preleuka_LDH_l (Factor+Higher Order Factors)	0.00	1	0.9907
product_s × preleuka_bulk_s (Factor+Higher Order Factors)	0.21	1	0.6472
product_s × age (Factor+Higher Order Factors)	0.23	1	0.6310
product_s × hct_score (Factor+Higher Order Factors)	1.47	1	0.2248
product_s × preleuka_ALC_l (Factor+Higher Order Factors)	0.02	1	0.8953
TOTAL INTERACTION	2.21	5	0.8189
TOTAL	18.02	11	0.0810

Show code

plot(anova(fit_mi))

Show code

fit_check_ord2(fit_mi)

n=129 Mean absolute error=0.027 Mean squared error=0.00149 0.9 Quantile of absolute error=0.058

n=129 Mean absolute error=0.05 Mean squared error=0.00403 0.9 Quantile of absolute error=0.094

Show code

options(prType=NULL)

Multivariable ordinal logistic regression model for ICANS grade including interaction terms

Show code

fit_mi <-
  fit.mult.impute(
    nt_g_r ~ product_s*preleuka_LDH_l + product_s*preleuka_bulk_s + product_s*age + product_s*hct_score + product_s*preleuka_ALC_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=TRUE,y=TRUE
  )

Variance Inflation Factors Due to Imputation:

                                   y>=1 
                                   1.00 
                                   y>=2 
                                   1.00 
                                 y>=3-5 
                                   1.00 
              product_s=tisacel/JCAR014 
                                   1.01 
                         preleuka_LDH_l 
                                   1.01 
                        preleuka_bulk_s 
                                   1.11 
                                    age 
                                   1.00 
                              hct_score 
                                   1.00 
                         preleuka_ALC_l 
                                   1.00

product_s=tisacel/JCAR014 * preleuka_LDH_l 1.02 product_s=tisacel/JCAR014 * preleuka_bulk_s 1.11 product_s=tisacel/JCAR014 * age 1.00 product_s=tisacel/JCAR014 * hct_score 1.01 product_s=tisacel/JCAR014 * preleuka_ALC_l 1.01

Rate of Missing Information:

                                   y>=1 
                                   0.00 
                                   y>=2 
                                   0.00 
                                 y>=3-5 
                                   0.00 
              product_s=tisacel/JCAR014 
                                   0.01 
                         preleuka_LDH_l 
                                   0.01 
                        preleuka_bulk_s 
                                   0.10 
                                    age 
                                   0.00 
                              hct_score 
                                   0.00 
                         preleuka_ALC_l 
                                   0.00

product_s=tisacel/JCAR014 * preleuka_LDH_l 0.02 product_s=tisacel/JCAR014 * preleuka_bulk_s 0.10 product_s=tisacel/JCAR014 * age 0.00 product_s=tisacel/JCAR014 * hct_score 0.01 product_s=tisacel/JCAR014 * preleuka_ALC_l 0.01

d.f. for t-distribution for Tests of Single Coefficients:

                                   y>=1 
                           265800609.25 
                                   y>=2 
                           271553554.60 
                                 y>=3-5 
                           218370964.08 
              product_s=tisacel/JCAR014 
                               77078.95 
                         preleuka_LDH_l 
                              106572.94 
                        preleuka_bulk_s 
                                 849.60 
                                    age 
                              409824.29 
                              hct_score 
                              403376.89 
                         preleuka_ALC_l 
                           428300149.88

product_s=tisacel/JCAR014 * preleuka_LDH_l 32826.65 product_s=tisacel/JCAR014 * preleuka_bulk_s 919.23 product_s=tisacel/JCAR014 * age 1027273.81 product_s=tisacel/JCAR014 * hct_score 344006.92 product_s=tisacel/JCAR014 * preleuka_ALC_l 230303.63

The following fit components were averaged over the 10 model fits:

stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = nt_g_r ~ product_s * preleuka_LDH_l + 
     product_s * preleuka_bulk_s + product_s * age + product_s * 
     hct_score + product_s * preleuka_ALC_l, fitter = lrm, xtrans = mi, 
     data = df, x = TRUE, y = TRUE)

Frequencies of Responses

   0   1   2 3-5 
  74   8  20  27

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 129	LR χ² 27.91	R² 0.218	C 0.720
max \|∂log L/∂β\| 2×10^-6	d.f. 11	g 1.125	D_xy 0.441
	Pr(>χ²) 0.0034	g_r 3.081	γ 0.441
		g_p 0.234	τ_a 0.266
		Brier 0.191

	β	S.E.	Wald Z	Pr(>\|Z\|)
y≥1	-4.55	2.60	-1.75	0.0805
y≥2	-4.87	2.61	-1.87	0.0616
y≥3-5	-5.79	2.63	-2.20	0.0275
product_s=tisacel/JCAR014	1.89	5.24	0.36	0.7180
preleuka_LDH_l	1.30	0.84	1.54	0.1233
preleuka_bulk_s	0.00	0.09	0.02	0.9836
age	0.04	0.02	1.78	0.0744
hct_score	-0.04	0.13	-0.31	0.7531
preleuka_ALC_l	0.85	0.68	1.26	0.2094
product_s=tisacel/JCAR014 × preleuka_LDH_l	-0.66	2.00	-0.33	0.7425
product_s=tisacel/JCAR014 × preleuka_bulk_s	0.04	0.12	0.32	0.7456
product_s=tisacel/JCAR014 × age	-0.04	0.03	-1.24	0.2141
product_s=tisacel/JCAR014 × hct_score	0.05	0.23	0.21	0.8301
product_s=tisacel/JCAR014 × preleuka_ALC_l	0.04	1.56	0.03	0.9782

Show code

print(anova(fit_mi))

	χ²	d.f.	P
Wald Statistics for `nt_g_r`
product_s (Factor+Higher Order Factors)	20.58	6	0.0022
All Interactions	1.90	5	0.8625
preleuka_LDH_l (Factor+Higher Order Factors)	2.50	2	0.2868
All Interactions	0.11	1	0.7425
preleuka_bulk_s (Factor+Higher Order Factors)	0.23	2	0.8927
All Interactions	0.11	1	0.7456
age (Factor+Higher Order Factors)	3.21	2	0.2008
All Interactions	1.54	1	0.2141
hct_score (Factor+Higher Order Factors)	0.10	2	0.9507
All Interactions	0.05	1	0.8301
preleuka_ALC_l (Factor+Higher Order Factors)	1.97	2	0.3731
All Interactions	0.00	1	0.9782
product_s × preleuka_LDH_l (Factor+Higher Order Factors)	0.11	1	0.7425
product_s × preleuka_bulk_s (Factor+Higher Order Factors)	0.11	1	0.7456
product_s × age (Factor+Higher Order Factors)	1.54	1	0.2141
product_s × hct_score (Factor+Higher Order Factors)	0.05	1	0.8301
product_s × preleuka_ALC_l (Factor+Higher Order Factors)	0.00	1	0.9782
TOTAL INTERACTION	1.90	5	0.8625
TOTAL	24.39	11	0.0112

Show code

plot(anova(fit_mi))

Show code

fit_check_ord3(fit_mi)

n=129 Mean absolute error=0.041 Mean squared error=0.00222 0.9 Quantile of absolute error=0.078

n=129 Mean absolute error=0.047 Mean squared error=0.00241 0.9 Quantile of absolute error=0.064

n=129 Mean absolute error=0.051 Mean squared error=0.00311 0.9 Quantile of absolute error=0.083

Show code

options(prType=NULL)

Multivariable logistic regression model for CR including interaction terms

Show code

fit_mi <-
  fit.mult.impute(
    bestCR_b ~ product_s * preleuka_ALC_l + product_s * preleuka_bulk_s + product_s *
      age + product_s * hct_score + product_s * preleuka_LDH_l,
    fitter = lrm,
    data = df,
    xtrans = mi,
    x=T,y=T
  )


Variance Inflation Factors Due to Imputation:

                                  Intercept 
                                       1.03 
                  product_s=tisacel/JCAR014 
                                       1.05 
                             preleuka_ALC_l 
                                       1.06 
                            preleuka_bulk_s 
                                       1.21 
                                        age 
                                       1.02 
                                  hct_score 
                                       1.07 
                             preleuka_LDH_l 
                                       1.02 
 product_s=tisacel/JCAR014 * preleuka_ALC_l 
                                       1.04 
product_s=tisacel/JCAR014 * preleuka_bulk_s 
                                       1.31 
            product_s=tisacel/JCAR014 * age 
                                       1.04 
      product_s=tisacel/JCAR014 * hct_score 
                                       1.06 
 product_s=tisacel/JCAR014 * preleuka_LDH_l 
                                       1.04 

Rate of Missing Information:

                                  Intercept 
                                       0.03 
                  product_s=tisacel/JCAR014 
                                       0.05 
                             preleuka_ALC_l 
                                       0.06 
                            preleuka_bulk_s 
                                       0.18 
                                        age 
                                       0.02 
                                  hct_score 
                                       0.07 
                             preleuka_LDH_l 
                                       0.02 
 product_s=tisacel/JCAR014 * preleuka_ALC_l 
                                       0.04 
product_s=tisacel/JCAR014 * preleuka_bulk_s 
                                       0.24 
            product_s=tisacel/JCAR014 * age 
                                       0.04 
      product_s=tisacel/JCAR014 * hct_score 
                                       0.06 
 product_s=tisacel/JCAR014 * preleuka_LDH_l 
                                       0.03 

d.f. for t-distribution for Tests of Single Coefficients:

                                  Intercept 
                                    9355.75 
                  product_s=tisacel/JCAR014 
                                    3485.52 
                             preleuka_ALC_l 
                                    2468.54 
                            preleuka_bulk_s 
                                     289.03 
                                        age 
                                   31859.35 
                                  hct_score 
                                    1883.91 
                             preleuka_LDH_l 
                                   23986.80 
 product_s=tisacel/JCAR014 * preleuka_ALC_l 
                                    6760.74 
product_s=tisacel/JCAR014 * preleuka_bulk_s 
                                     158.34 
            product_s=tisacel/JCAR014 * age 
                                    6379.68 
      product_s=tisacel/JCAR014 * hct_score 
                                    2702.39 
 product_s=tisacel/JCAR014 * preleuka_LDH_l 
                                    7857.15 

The following fit components were averaged over the 10 model fits:

  stats linear.predictors

Show code

options(prType="html")
print(
  fit_mi,
  digits = 2,
  coefs = T
)

Logistic Regression Model

 fit.mult.impute(formula = bestCR_b ~ product_s * preleuka_ALC_l + 
     product_s * preleuka_bulk_s + product_s * age + product_s * 
     hct_score + product_s * preleuka_LDH_l, fitter = lrm, xtrans = mi, 
     data = df, x = T, y = T)

Frequencies of Missing Values Due to Each Variable

        bestCR_b       product_s  preleuka_ALC_l preleuka_bulk_s 
               3               0               0               0 
             age       hct_score  preleuka_LDH_l 
               0               0               0

	Model Likelihood Ratio Test	Discrimination Indexes	Rank Discrim. Indexes
Obs 126	LR χ² 61.95	R² 0.519	C 0.864
0 68	d.f. 11	g 2.623	D_xy 0.729
1 58	Pr(>χ²) <0.0001	g_r 13.933	γ 0.729
max \|∂log L/∂β\| 5×10^-8		g_p 0.368	τ_a 0.365
		Brier 0.150

	β	S.E.	Wald Z	Pr(>\|Z\|)
Intercept	18.44	6.20	2.97	0.0029
product_s=tisacel/JCAR014	-7.92	8.06	-0.98	0.3263
preleuka_ALC_l	4.45	1.59	2.80	0.0051
preleuka_bulk_s	-0.48	0.17	-2.76	0.0058
age	-0.01	0.04	-0.14	0.8852
hct_score	0.11	0.27	0.42	0.6762
preleuka_LDH_l	-6.08	2.10	-2.90	0.0038
product_s=tisacel/JCAR014 × preleuka_ALC_l	-4.14	2.04	-2.03	0.0426
product_s=tisacel/JCAR014 × preleuka_bulk_s	0.32	0.24	1.36	0.1745
product_s=tisacel/JCAR014 × age	0.00	0.04	0.10	0.9172
product_s=tisacel/JCAR014 × hct_score	-0.04	0.33	-0.12	0.9022
product_s=tisacel/JCAR014 × preleuka_LDH_l	1.67	2.97	0.56	0.5739

Show code

# options(prType=NULL)
print(anova(fit_mi))

	χ²	d.f.	P
Wald Statistics for `bestCR_b`
product_s (Factor+Higher Order Factors)	11.78	6	0.0670
All Interactions	5.18	5	0.3942
preleuka_ALC_l (Factor+Higher Order Factors)	7.89	2	0.0193
All Interactions	4.11	1	0.0426
preleuka_bulk_s (Factor+Higher Order Factors)	8.08	2	0.0176
All Interactions	1.84	1	0.1745
age (Factor+Higher Order Factors)	0.02	2	0.9892
All Interactions	0.01	1	0.9172
hct_score (Factor+Higher Order Factors)	0.32	2	0.8510
All Interactions	0.02	1	0.9022
preleuka_LDH_l (Factor+Higher Order Factors)	12.89	2	0.0016
All Interactions	0.32	1	0.5739
product_s × preleuka_ALC_l (Factor+Higher Order Factors)	4.11	1	0.0426
product_s × preleuka_bulk_s (Factor+Higher Order Factors)	1.84	1	0.1745
product_s × age (Factor+Higher Order Factors)	0.01	1	0.9172
product_s × hct_score (Factor+Higher Order Factors)	0.02	1	0.9022
product_s × preleuka_LDH_l (Factor+Higher Order Factors)	0.32	1	0.5739
TOTAL INTERACTION	5.18	5	0.3942
TOTAL	23.83	11	0.0135

Show code

plot(anova(fit_mi))

Show code

plot(calibrate(fit_mi))


n=126   Mean absolute error=0.066   Mean squared error=0.00527
0.9 Quantile of absolute error=0.109

Show code

options(prType=NULL)

Figure 4. The efficacy of distinct CAR T-cell products is differentially impacted by preleukapheresis ALC

Show code

alc_breaks <- round(c(10^-.5,10^-.25,10^0,10^.25),2)
data.frame(Predict(fit_mi,preleuka_ALC_l,product_s,fun = plogis)) %>%
  mutate(product_s=recode(product_s,`axicel`="axicel",`tisacel/jcar014`="tisacel/JCAR014")) %>% 
  ggplot() +
  geom_line(aes(x=preleuka_ALC_l,y=yhat,col=product_s)) +
  geom_ribbon(aes(x=preleuka_ALC_l,ymin=lower, ymax=upper,fill=product_s),alpha=0.10)+
  labs(x = "Preleukapheresis ALC (G/L)",
       y = "Probability of CR")+
  scale_x_continuous(labels = alc_breaks,breaks = c(-.5,-.25,0,.25))+
  ggtitle("Figure 4.")+
  theme_bw() +
  theme(title = element_text(size=16,face="bold"),legend.title=element_blank())

Show code

# Save to eps
# cairo_ps("Manuscript/Fig 4.eps", width = 9, height = 8)
# data.frame(Predict(fit_mi,preleuka_ALC_l,product_s,fun = plogis)) %>%
#   mutate(product_s=recode(product_s,`axicel`="axicel",`tisacel/jcar014`="tisacel/JCAR014")) %>% 
#   ggplot() +
#   geom_line(aes(x=preleuka_ALC_l,y=yhat,col=product_s)) +
#   geom_ribbon(aes(x=preleuka_ALC_l,ymin=lower, ymax=upper,fill=product_s),alpha=0.10)+
#   labs(x = "Preleukapheresis ALC (G/L)",
#        y = "Probability of CR")+
#   scale_x_continuous(labels = alc_breaks,breaks = c(-.5,-.25,0,.25))+
#   ggtitle("Figure 4.")+
#   theme_bw() +
#   theme(plot.title = element_text(size=16,face="bold"),legend.title=element_blank())
# dev.off()

Computational environment

Show code

devtools::session_info()

─ Session info ─────────────────────────────────────────────────────
 setting  value
 version  R version 4.2.0 (2022-04-22)
 os       macOS Big Sur/Monterey 10.16
 system   x86_64, darwin17.0
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Los_Angeles
 date     2022-05-18
 pandoc   2.17.1.1 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/ (via rmarkdown)

─ Packages ─────────────────────────────────────────────────────────
 package       * version date (UTC) lib source
 abind           1.4-5   2016-07-21 [1] CRAN (R 4.2.0)
 assertthat      0.2.1   2019-03-21 [1] CRAN (R 4.2.0)
 backports       1.4.1   2021-12-13 [1] CRAN (R 4.2.0)
 base64enc       0.1-3   2015-07-28 [1] CRAN (R 4.2.0)
 brio            1.1.3   2021-11-30 [1] CRAN (R 4.2.0)
 broom           0.8.0   2022-04-13 [1] CRAN (R 4.2.0)
 broom.helpers   1.7.0   2022-04-22 [1] CRAN (R 4.2.0)
 bslib           0.3.1   2021-10-06 [1] CRAN (R 4.2.0)
 cachem          1.0.6   2021-08-19 [1] CRAN (R 4.2.0)
 callr           3.7.0   2021-04-20 [1] CRAN (R 4.2.0)
 car             3.0-13  2022-05-02 [1] CRAN (R 4.2.0)
 carData         3.0-5   2022-01-06 [1] CRAN (R 4.2.0)
 cellranger      1.1.0   2016-07-27 [1] CRAN (R 4.2.0)
 checkmate       2.1.0   2022-04-21 [1] CRAN (R 4.2.0)
 cli             3.3.0   2022-04-25 [1] CRAN (R 4.2.0)
 cluster         2.1.3   2022-03-28 [1] CRAN (R 4.2.0)
 codetools       0.2-18  2020-11-04 [1] CRAN (R 4.2.0)
 colorspace      2.0-3   2022-02-21 [1] CRAN (R 4.2.0)
 commonmark      1.8.0   2022-03-09 [1] CRAN (R 4.2.0)
 crayon          1.5.1   2022-03-26 [1] CRAN (R 4.2.0)
 data.table      1.14.2  2021-09-27 [1] CRAN (R 4.2.0)
 DBI             1.1.2   2021-12-20 [1] CRAN (R 4.2.0)
 dbplyr          2.1.1   2021-04-06 [1] CRAN (R 4.2.0)
 desc            1.4.1   2022-03-06 [1] CRAN (R 4.2.0)
 devtools        2.4.3   2021-11-30 [1] CRAN (R 4.2.0)
 digest          0.6.29  2021-12-01 [1] CRAN (R 4.2.0)
 distill         1.3     2021-10-13 [1] CRAN (R 4.2.0)
 downlit         0.4.0   2021-10-29 [1] CRAN (R 4.2.0)
 dplyr         * 1.0.8   2022-02-08 [1] CRAN (R 4.2.0)
 ellipsis        0.3.2   2021-04-29 [1] CRAN (R 4.2.0)
 evaluate        0.15    2022-02-18 [1] CRAN (R 4.2.0)
 fansi           1.0.3   2022-03-24 [1] CRAN (R 4.2.0)
 farver          2.1.0   2021-02-28 [1] CRAN (R 4.2.0)
 fastmap         1.1.0   2021-01-25 [1] CRAN (R 4.2.0)
 forcats       * 0.5.1   2021-01-27 [1] CRAN (R 4.2.0)
 foreign         0.8-82  2022-01-16 [1] CRAN (R 4.2.0)
 Formula       * 1.2-4   2020-10-16 [1] CRAN (R 4.2.0)
 fs              1.5.2   2021-12-08 [1] CRAN (R 4.2.0)
 generics        0.1.2   2022-01-31 [1] CRAN (R 4.2.0)
 ggplot2       * 3.3.5   2021-06-25 [1] CRAN (R 4.2.0)
 ggpubr        * 0.4.0   2020-06-27 [1] CRAN (R 4.2.0)
 ggsignif        0.6.3   2021-09-09 [1] CRAN (R 4.2.0)
 glue            1.6.2   2022-02-24 [1] CRAN (R 4.2.0)
 gridExtra     * 2.3     2017-09-09 [1] CRAN (R 4.2.0)
 gt              0.5.0   2022-04-21 [1] CRAN (R 4.2.0)
 gtable          0.3.0   2019-03-25 [1] CRAN (R 4.2.0)
 gtsummary     * 1.6.0   2022-04-25 [1] CRAN (R 4.2.0)
 haven           2.5.0   2022-04-15 [1] CRAN (R 4.2.0)
 highr           0.9     2021-04-16 [1] CRAN (R 4.2.0)
 Hmisc         * 4.6-0   2021-10-07 [1] CRAN (R 4.2.0)
 hms             1.1.1   2021-09-26 [1] CRAN (R 4.2.0)
 htmlTable     * 2.4.0   2022-01-04 [1] CRAN (R 4.2.0)
 htmltools       0.5.2   2021-08-25 [1] CRAN (R 4.2.0)
 htmlwidgets     1.5.4   2021-09-08 [1] CRAN (R 4.2.0)
 httr            1.4.2   2020-07-20 [1] CRAN (R 4.2.0)
 jpeg            0.1-9   2021-07-24 [1] CRAN (R 4.2.0)
 jquerylib       0.1.4   2021-04-26 [1] CRAN (R 4.2.0)
 jsonlite        1.8.0   2022-02-22 [1] CRAN (R 4.2.0)
 knitr           1.39    2022-04-26 [1] CRAN (R 4.2.0)
 labeling        0.4.2   2020-10-20 [1] CRAN (R 4.2.0)
 lattice       * 0.20-45 2021-09-22 [1] CRAN (R 4.2.0)
 latticeExtra    0.6-29  2019-12-19 [1] CRAN (R 4.2.0)
 lifecycle       1.0.1   2021-09-24 [1] CRAN (R 4.2.0)
 lubridate       1.8.0   2021-10-07 [1] CRAN (R 4.2.0)
 magrittr        2.0.3   2022-03-30 [1] CRAN (R 4.2.0)
 MASS            7.3-56  2022-03-23 [1] CRAN (R 4.2.0)
 Matrix          1.4-1   2022-03-23 [1] CRAN (R 4.2.0)
 MatrixModels    0.5-0   2021-03-02 [1] CRAN (R 4.2.0)
 memoise         2.0.1   2021-11-26 [1] CRAN (R 4.2.0)
 modelr          0.1.8   2020-05-19 [1] CRAN (R 4.2.0)
 multcomp        1.4-19  2022-04-26 [1] CRAN (R 4.2.0)
 munsell         0.5.0   2018-06-12 [1] CRAN (R 4.2.0)
 mvtnorm         1.1-3   2021-10-08 [1] CRAN (R 4.2.0)
 nlme            3.1-157 2022-03-25 [1] CRAN (R 4.2.0)
 nnet            7.3-17  2022-01-16 [1] CRAN (R 4.2.0)
 openxlsx      * 4.2.5   2021-12-14 [1] CRAN (R 4.2.0)
 ormPlot       * 0.3.4   2020-12-17 [1] CRAN (R 4.2.0)
 patchwork     * 1.1.1   2020-12-17 [1] CRAN (R 4.2.0)
 pillar          1.7.0   2022-02-01 [1] CRAN (R 4.2.0)
 pkgbuild        1.3.1   2021-12-20 [1] CRAN (R 4.2.0)
 pkgconfig       2.0.3   2019-09-22 [1] CRAN (R 4.2.0)
 pkgload         1.2.4   2021-11-30 [1] CRAN (R 4.2.0)
 png             0.1-7   2013-12-03 [1] CRAN (R 4.2.0)
 polspline       1.1.19  2020-05-15 [1] CRAN (R 4.2.0)
 prettyunits     1.1.1   2020-01-24 [1] CRAN (R 4.2.0)
 processx        3.5.3   2022-03-25 [1] CRAN (R 4.2.0)
 ps              1.7.0   2022-04-23 [1] CRAN (R 4.2.0)
 purrr         * 0.3.4   2020-04-17 [1] CRAN (R 4.2.0)
 quantreg        5.88    2022-02-05 [1] CRAN (R 4.2.0)
 R6              2.5.1   2021-08-19 [1] CRAN (R 4.2.0)
 RColorBrewer    1.1-3   2022-04-03 [1] CRAN (R 4.2.0)
 Rcpp            1.0.8.3 2022-03-17 [1] CRAN (R 4.2.0)
 readr         * 2.1.2   2022-01-30 [1] CRAN (R 4.2.0)
 readxl          1.4.0   2022-03-28 [1] CRAN (R 4.2.0)
 remotes         2.4.2   2021-11-30 [1] CRAN (R 4.2.0)
 repr            1.1.4   2022-01-04 [1] CRAN (R 4.2.0)
 reprex          2.0.1   2021-08-05 [1] CRAN (R 4.2.0)
 rlang           1.0.2   2022-03-04 [1] CRAN (R 4.2.0)
 rmarkdown       2.14    2022-04-25 [1] CRAN (R 4.2.0)
 rms           * 6.2-0   2021-03-18 [1] CRAN (R 4.2.0)
 rpart           4.1.16  2022-01-24 [1] CRAN (R 4.2.0)
 rprojroot       2.0.3   2022-04-02 [1] CRAN (R 4.2.0)
 rstatix         0.7.0   2021-02-13 [1] CRAN (R 4.2.0)
 rstudioapi      0.13    2020-11-12 [1] CRAN (R 4.2.0)
 rvest           1.0.2   2021-10-16 [1] CRAN (R 4.2.0)
 sandwich        3.0-1   2021-05-18 [1] CRAN (R 4.2.0)
 sass            0.4.1   2022-03-23 [1] CRAN (R 4.2.0)
 scales          1.2.0   2022-04-13 [1] CRAN (R 4.2.0)
 sessioninfo     1.2.2   2021-12-06 [1] CRAN (R 4.2.0)
 skimr         * 2.1.4   2022-04-15 [1] CRAN (R 4.2.0)
 SparseM       * 1.81    2021-02-18 [1] CRAN (R 4.2.0)
 stringi         1.7.6   2021-11-29 [1] CRAN (R 4.2.0)
 stringr       * 1.4.0   2019-02-10 [1] CRAN (R 4.2.0)
 survival      * 3.3-1   2022-03-03 [1] CRAN (R 4.2.0)
 testthat        3.1.3   2022-03-29 [1] CRAN (R 4.2.0)
 TH.data         1.1-1   2022-04-26 [1] CRAN (R 4.2.0)
 tibble        * 3.1.6   2021-11-07 [1] CRAN (R 4.2.0)
 tidyr         * 1.2.0   2022-02-01 [1] CRAN (R 4.2.0)
 tidyselect      1.1.2   2022-02-21 [1] CRAN (R 4.2.0)
 tidyverse     * 1.3.1   2021-04-15 [1] CRAN (R 4.2.0)
 tzdb            0.3.0   2022-03-28 [1] CRAN (R 4.2.0)
 usethis         2.1.5   2021-12-09 [1] CRAN (R 4.2.0)
 utf8            1.2.2   2021-07-24 [1] CRAN (R 4.2.0)
 vctrs           0.4.1   2022-04-13 [1] CRAN (R 4.2.0)
 withr           2.5.0   2022-03-03 [1] CRAN (R 4.2.0)
 xfun            0.30    2022-03-02 [1] CRAN (R 4.2.0)
 xml2            1.3.3   2021-11-30 [1] CRAN (R 4.2.0)
 yaml            2.3.5   2022-02-21 [1] CRAN (R 4.2.0)
 zip             2.2.0   2021-05-31 [1] CRAN (R 4.2.0)
 zoo             1.8-10  2022-04-15 [1] CRAN (R 4.2.0)

 [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library

────────────────────────────────────────────────────────────────────