---
title: "Example 7: Survival Analysis"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Example 7: Survival Analysis}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)
```
Below is another complete example showing how **sassy**
functions interact with **tidyverse** to create a survival
analysis. This is the most complicated example in the documentation.
## Program
There are several interesting aspects to this report. The first is
that the **reporter** package gives you the freedom to perform any analysis
you want. Any statistics that R can produce may be added to a **reporter**
report. You just need to get them into a data frame or plot,
and the **reporter** package will render the analysis as desired.
Second, observe the ability to append both tabular and plot content onto the
same report. This appending ability was designed into the **reporter**
package from the beginning, and gives you the flexibility you need to document
your entire analysis.
Finally, notice that the table on page one is a complex table, that is actually
composed of two tables stacked on top of one another. The first table
`tbl1` contains the column headers and some PSGA counts. The second
table `tbl2` is marked as `headerless` and contains the Kaplan-Meier analysis.
The two tables are lined up so they appear as one table.
This technique of stacking one table on top of another
allows you to create very complex tables in a straight-forward manner.
```{r eval=FALSE, echo=TRUE}
library(tidyverse)
library(sassy)
library(broom)
library(survival)
library(survminer)
options("logr.autolog" = TRUE,
"logr.notes" = FALSE)
# Get temp location for log and report output
tmp <- tempdir()
# Open Log
lf <- log_open(file.path(tmp, "example7.log"))
# Load and Filter Data --------------------------------------------------
sep("Load and Filter Data")
# Get path to sample data
pkg <- system.file("extdata", package = "sassy")
# Get adam data
libname(adam, pkg, "sas7bdat")
# Filter data
adsl <- adam$adsl |>
select(USUBJID, SEX, AGEGR1, AGE, ARM) |>
filter(ARM != "SCREEN FAILURE") |> put()
adpsga <- adam$adpsga |>
filter(PARAMCD =="PSGA" & TRTA != "" & !is.na(AVISITN)) |>
select(USUBJID, TRTA, AVISIT, AVISITN, AVAL, CRIT1FL) |> put()
# Get population counts
arm_pop <- adsl |> count(ARM) |> deframe() |> put()
# Prepare Data ------------------------------------------------------------
sep("Prepare data for analysis")
put("Determine minimum visit at which success was achieved")
adpsga_minvsuccess <-
adpsga |>
filter(CRIT1FL == 'Y') |>
group_by(USUBJID) |>
summarize(minvisit = min(AVISITN))
put("Get subjects which did not achieve success")
adpsga_nosuccess <-
anti_join(adpsga, adpsga_minvsuccess, by = ('USUBJID')) |>
group_by(USUBJID) |>
summarize(maxvisit = max(AVISITN))
put("Combine subjects cured with subjects not cured")
adslpsga_final <-
inner_join(adsl, adpsga, by = c('USUBJID')) |>
left_join(adpsga_minvsuccess, by = c('USUBJID')) |>
left_join(adpsga_nosuccess, by = c('USUBJID')) |>
filter((AVISITN == minvisit & !is.na(minvisit)) |
(AVISITN == maxvisit & !is.na(maxvisit))) |>
mutate(cured = case_when(CRIT1FL == "Y" ~ TRUE,
TRUE ~ as.logical(FALSE))) |>
select(-minvisit, -maxvisit)
# Counts ---------------------------------------------------------------
sep("Perform Counts and Statistical Tests")
put("Count patients with PSGA <= 1")
success_counts <- adslpsga_final |>
filter(cured == TRUE) |>
count(TRTA) |>
pivot_wider(names_from = TRTA,
values_from = n) |>
transmute(block = "counts",
label = "Number of patients with PSGA <= 1",
"ARM A" = as.character(`ARM A`),
"ARM B" = as.character(`ARM B`),
"ARM C" = as.character(`ARM C`),
"ARM D" = as.character(`ARM D`)) |>
put()
put("Count patients with PSGA > 1")
failed_counts <- adslpsga_final |>
filter(cured == FALSE) |>
count(TRTA) |>
pivot_wider(names_from = TRTA,
values_from = n) |>
transmute(block = "counts",
label = "Number of Censored Subjects (PSGA > 1)",
"ARM A" = as.character(`ARM A`),
"ARM B" = as.character(`ARM B`),
"ARM C" = as.character(`ARM C`),
"ARM D" = as.character(`ARM D`)) |>
put()
count_block <- bind_rows(success_counts, failed_counts)
# Kaplan-Meier estimates ----------------------------------------------------
sep("Perform Kaplan-Meier Tests")
put("Create survival vector")
surv_vct <- Surv(time = adslpsga_final$AVISITN, event = adslpsga_final$cured) |>
put()
put("Fit model on survival vector")
stats_survfit_trta <- survival::survfit(surv_vct ~ TRTA, data = adslpsga_final, ) |>
put()
put("Get model quantiles")
stats_survfit_quantiles <- quantile(stats_survfit_trta)
put("Get lower confidence intervals")
ci_lower <-
as.data.frame(stats_survfit_quantiles$lower) |>
rownames_to_column() |>
mutate(block = "surv",
TRTA = substring(rowname,6)) |>
pivot_longer(cols = c(`25`, `50`, `75`),
names_to = "Q",
values_to = "lower") |>
put()
put("Get upper confidence intervals")
ci_upper <-
as.data.frame(stats_survfit_quantiles$upper) |>
rownames_to_column() |>
mutate(block = "surv",
TRTA = substring(rowname,6)) |>
pivot_longer(cols = c(`25`, `50`, `75`),
names_to = "Q",
values_to = "upper") |>
put()
put("Get confidence intervals")
ci <-
inner_join(ci_lower, ci_upper) |>
mutate(ci = paste0("(", ifelse(is.na(lower), "-", lower)
, ", ", ifelse(is.na(upper), "-", upper), ")")) |>
pivot_wider(id_cols = c("block", "Q"),
names_from = TRTA,
values_from = ci) |>
mutate(order=2,
label1 = case_when(Q == 25 ~ "25th percentile (weeks)",
Q == 50 ~ "Median (weeks)",
Q == 75 ~ "75th percentile (weeks)"),
label2 = "95% Confidence Interval**") |>
select(block, Q, order, label1, label2,
`ARM A`, `ARM B`, `ARM C`, `ARM D`) |>
put()
put("Get quantiles")
quants <-
as.data.frame(stats_survfit_quantiles$quantile) |>
rownames_to_column() |>
mutate(block = "surv",
TRTA = substring(rowname,6)) |>
pivot_longer(cols = c(`25`, `50`, `75`),
names_to = "Q",
values_to = "value") |>
pivot_wider(id_cols = c("block", "Q"),
names_from = TRTA,
values_from = value) |>
mutate(order=1,
label1 = case_when(Q == 25 ~ "25th percentile (weeks)",
Q == 50 ~ "Median (weeks)",
Q == 75 ~ "75th percentile (weeks)"),
label2 = "",
`ARM A` = as.character(`ARM A`),
`ARM B` = as.character(`ARM B`),
`ARM C` = as.character(`ARM C`),
`ARM D` = as.character(`ARM D`)) |>
select(block, Q, order, label1, label2,
`ARM A`, `ARM B`, `ARM C`, `ARM D`) |>
put()
put("Final arrangement")
kaplan_block <-
bind_rows(quants, ci) |>
arrange(block, Q, order) |>
transmute(block,
label1,
label2 = ifelse(label2 == "", NA, label2),
`ARM A`, `ARM B`, `ARM C`, `ARM D`) |>
put()
# Cox Proportional Hazards -----------------------------------------------
sep("Perform Cox Proportional Hazards Test")
put("Run Cox tests")
stats_surv_cph <- survival::coxph(surv_vct ~ TRTA, data = adslpsga_final) |>
put()
put("Create summary statistics on Cox results")
cph_summary <-
summary(stats_surv_cph) |>
put()
put("Extract coefficients")
cph_coef <-
as.data.frame(cph_summary$coefficients) |>
rownames_to_column() |>
mutate(block = "surv",
TRTA = substring(rowname,5)) |>
put()
put("Extract confidence intervals")
cph_ci <-
cph_summary$conf.int |>
as.data.frame(cph_summary$conf) |>
rownames_to_column() |>
put()
put("Create cox statistics block")
cox_block <-
bind_cols(cph_coef, cph_ci) |>
rename(hazard = `exp(coef)...3`, pval = `Pr(>|z|)`,
lower = `lower .95`, upper = `upper .95`) |>
select(TRTA, hazard, pval, lower, upper) |>
mutate(block = "cox",
ci = paste0("(", ifelse(is.na(lower), "-", sprintf("%.2f", lower))
, ", ", ifelse(is.na(upper), "-", sprintf("%.2f", upper)), ")"),
hazard = sprintf("%.2f", hazard),
pval = sprintf("%.3f", pval)) |>
pivot_longer(cols = c("hazard", "pval", "ci"),
names_to = "stat",
values_to = "value") |>
pivot_wider(id_cols = c("block", "stat"),
names_from = TRTA,
values_from = value) |>
mutate(label1 = case_when(stat == "hazard"
~ "Hazard Ratio (Each Treatment Group - ARM A)***",
stat == "pval" ~ "P-value",
TRUE ~ "95% CI of Hazard Ratio"),
label2 = as.character(NA),
`ARM A` = as.character(NA)) |>
select(block, label1, label2, `ARM A`, `ARM B`, `ARM C`, `ARM D`) |>
put()
put("Combine statistics blocks")
stat_block <- bind_rows(kaplan_block, cox_block) |> put()
# Create Survival Plot -------------------------------------------------------
sep("Create survival plot")
put("Create data frame with zero values for each visit")
arms <- unique(adslpsga_final$ARM)
visits <- unique(adslpsga_final$AVISITN)
all_visits <- rep(arms, length(visits))
all_visits <- all_visits[order(all_visits)]
put("Create visit template")
df <- data.frame(ARM = all_visits,
AVISIT = paste("Week", visits),
AVISITN = visits,
cured = FALSE) |> put()
put("Calculate cummulative sum and percent")
adslpsga_plot <- adslpsga_final |>
select(ARM, AVISIT, AVISITN, cured) |>
bind_rows(df) |>
group_by(ARM, AVISIT, AVISITN) |>
summarize(sumc = sum(cured)) |>
arrange(ARM, AVISITN) |>
group_by(ARM) |>
mutate(AVISIT = ifelse(AVISIT == "Week 0", "Day 1 Baseline", AVISIT),
csumc = cumsum(sumc)) |>
distinct() |>
mutate(pct = case_when(ARM == "ARM A" ~ csumc / arm_pop["ARM A"],
ARM == "ARM B" ~ csumc / arm_pop["ARM B"],
ARM == "ARM C" ~ csumc / arm_pop["ARM C"],
ARM == "ARM D" ~ csumc / arm_pop["ARM D"])) |>
put()
# Add factor to ensure sort order is correct
adslpsga_plot$AVISIT <- factor(adslpsga_plot$AVISIT,
levels = c("Day 1 Baseline",
"Week 2",
"Week 4",
"Week 6",
"Week 8",
"Week 12",
"Week 16"))
put("Generate plot")
surv_gg <- adslpsga_plot |>
ggplot(mapping = aes(y = pct, x = AVISIT , group = ARM)) +
geom_point(aes(shape = ARM, color = ARM)) +
geom_step(aes(linetype = ARM, color = ARM)) +
scale_x_discrete(name = "Study Week") +
scale_y_continuous(name = "Proportion of Subjects with Initial Success")
# Print Report ------------------------------------------------------------
sep("Create and print report")
# Create Table 1 with header
tbl1 <- create_table(count_block, width = 9) |>
column_defaults(from = `ARM A`, to = `ARM D`, align = "center", width = 1.1) |>
define(block, visible = FALSE) |>
define(label, label = "", width = 4.25) |>
define(`ARM A`, n = arm_pop["ARM A"]) |>
define(`ARM B`, n = arm_pop["ARM B"]) |>
define(`ARM C`, n = arm_pop["ARM C"]) |>
define(`ARM D`, n = arm_pop["ARM D"]) |>
titles("Table 5.0", bold = TRUE, blank_row = "above") |>
titles("Analysis of Time to Initial PSGA Success* in Weeks",
"Safety Population")
label_lookup <- c(surv = "Kaplan-Meier estimates",
cox = "Results of Proportional Hazards Regression Analysis")
# Create table 2 for statistics with stub and without header
tbl2 <- create_table(stat_block, width = 9, headerless = TRUE) |>
column_defaults(from = `ARM A`, to = `ARM D`, align = "center", width = 1.1) |>
stub(c(block, label1, label2), width = 4.25) |>
define(block, label_row = TRUE, format = label_lookup, blank_after = TRUE) |>
define(label1, indent = .25) |>
define(label2, indent = .5) |>
define(`ARM A`) |>
define(`ARM B`) |>
define(`ARM C`) |>
define(`ARM D`)
# Create plot
plt <- create_plot(surv_gg, 3.5, 9) |>
titles("Figure 5.0", bold = TRUE, blank_row = "above") |>
titles("Kaplan-Meier Plot for Time to Initial PSGA Success (PSGA <= 1)",
"Safety Population", blank_row = "none")
put("Create Report")
# Add table 1, table 2, and plot content to the same report.
# Plot will be on a separate page with it's own title.
rpt <- create_report(file.path(tmp, "output/example7.rtf"), output_type = "RTF",
font = "Arial", missing = "-") |>
set_margins(top = 1, bottom = 1) |>
page_header("Sponsor: Client", "Study: ABC/BBC") |>
add_content(tbl1, page_break = FALSE) |>
add_content(tbl2) |>
add_content(plt) |>
footnotes("Program: Surv_Table.R",
"* Success: PSGA <= 1: PSGA > 1",
"** Based on R survival package survfit() function",
paste("*** Based on proportional hazards model with treatment",
"indicator variables as explanatory variables"),
paste("Note: The end-point is cure of the disease (PSGA <= 1)."),
" The probability of remaining diseased (PSGA > 1) defines the survival function",
paste("Note: A subject who is not cured by the end of 12 weeks or is lost",
"to follow provides a censored observation for the analysis."),
"\"-\" = Not Applicable") |>
page_footer("Date Produced: " %p% fapply(Sys.time(), "%d%b%y %H:%M"),
right = "Page [pg] of [tpg]")
put("Write out the report")
res <- write_report(rpt)
# Clean Up ----------------------------------------------------------------
sep("Clean Up")
# Close log
log_close()
# View report
# file.show(res$file_path)
# View log
# file.show(lf)
```
## Output
Here is the report produced by the sample program above:
## Log
The above program produces the following log:
```
=========================================================================
Log Path: C:/Users/dbosa/AppData/Local/Temp/RtmpgrpGQ1/log/example7.log
Program Path: C:\packages\Testing\example7.R
Working Directory: C:/packages/Testing
User Name: dbosa
R Version: 4.1.2 (2021-11-01)
Machine: SOCRATES x86-64
Operating System: Windows 10 x64 build 19041
Base Packages: stats graphics grDevices utils datasets methods base
Other Packages: tidylog_1.0.2 survminer_0.4.9 ggpubr_0.4.0 survival_3.2-13 broom_0.7.10
reporter_1.2.6 libr_1.2.1 fmtr_1.5.4 logr_1.2.7 sassy_1.0.5
forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.0.2
tidyr_1.1.4 tibble_3.1.5 ggplot2_3.3.5 tidyverse_1.3.1
Log Start Time: 2021-11-21 15:30:37
=========================================================================
=========================================================================
Load and Filter Data
=========================================================================
# library 'adam': 2 items
- attributes: sas7bdat not loaded
- path: C:/Users/dbosa/Documents/R/win-library/4.1/reporter/extdata
- items:
Name Extension Rows Cols Size
1 adpsga sas7bdat 1206 42 424.8 Kb
2 adsl sas7bdat 152 56 98.7 Kb
LastModified
1 2021-10-09 13:57:48
2 2021-10-09 13:57:48
select: dropped 51 variables (STUDYID, SUBJID, SITEID, AGEU, RACE, )
filter: removed 2 rows (1%), 150 rows remaining
# A tibble: 150 x 5
USUBJID SEX AGEGR1 AGE ARM
1 ABC-01-049 M 30-39 years 39 ARM D
2 ABC-01-050 M 40-49 years 47 ARM B
3 ABC-01-051 M 30-39 years 34 ARM A
4 ABC-01-052 F 40-49 years 45 ARM C
5 ABC-01-053 F 18-29 years 26 ARM B
6 ABC-01-054 M 40-49 years 44 ARM D
7 ABC-01-055 F 40-49 years 47 ARM C
8 ABC-01-056 M 30-39 years 31 ARM A
9 ABC-01-113 M >65 years 74 ARM D
10 ABC-01-114 F >65 years 72 ARM B
# ... with 140 more rows
filter: removed 256 rows (21%), 950 rows remaining
select: dropped 36 variables (STUDYID, SUBJID, SITEID, QSSEQ, TRTP, )
# A tibble: 950 x 6
USUBJID TRTA AVISIT AVISITN AVAL
1 ABC-01-049 ARM D Day 1 Bas~ 0 3
2 ABC-01-049 ARM D Week 2 2 2
3 ABC-01-049 ARM D Week 4 4 3
4 ABC-01-049 ARM D Week 6 6 3
5 ABC-01-049 ARM D Week 8 8 2
6 ABC-01-049 ARM D Week 8 8 2
7 ABC-01-049 ARM D Week 12 12 2
8 ABC-01-049 ARM D Week 16 16 3
9 ABC-01-050 ARM B Day 1 Bas~ 0 3
10 ABC-01-050 ARM B Week 2 2 3
# ... with 940 more rows, and 1 more
# variable: CRIT1FL
count: now 4 rows and 2 columns, ungrouped
ARM A ARM B ARM C ARM D
36 38 38 38
=========================================================================
Prepare data for analysis
=========================================================================
Determine minimum visit at which success was achieved
filter: removed 833 rows (88%), 117 rows remaining
group_by: one grouping variable (USUBJID)
summarize: now 43 rows and 2 columns, ungrouped
Get subjects which did not achieve success
anti_join: added no columns
> rows only in x 659
> rows only in y ( 0)
> matched rows (291)
> =====
> rows total 659
group_by: one grouping variable (USUBJID)
summarize: now 102 rows and 2 columns, ungrouped
Combine subjects cured with subjects not cured
inner_join: added 5 columns (TRTA, AVISIT, AVISITN, AVAL, CRIT1FL)
> rows only in x ( 5)
> rows only in y ( 0)
> matched rows 950 (includes duplicates)
> =====
> rows total 950
left_join: added one column (minvisit)
> rows only in x 659
> rows only in y ( 0)
> matched rows 291
> =====
> rows total 950
left_join: added one column (maxvisit)
> rows only in x 291
> rows only in y ( 0)
> matched rows 659
> =====
> rows total 950
filter: removed 803 rows (85%), 147 rows remaining
mutate: new variable 'cured' (logical) with 2 unique values and 0% NA
select: dropped 2 variables (minvisit, maxvisit)
=========================================================================
Perform Counts and Statistical Tests
=========================================================================
Count patients with PSGA <= 1
filter: removed 104 rows (71%), 43 rows remaining
count: now 4 rows and 2 columns, ungrouped
pivot_wider: reorganized (TRTA, n) into (ARM A, ARM B, ARM C, ARM D) [was 4x2, now 1x4]
transmute: new variable 'block' (character) with one unique value and 0% NA
new variable 'label' (character) with one unique value and 0% NA
converted 'ARM A' from integer to character (0 new NA)
converted 'ARM B' from integer to character (0 new NA)
converted 'ARM C' from integer to character (0 new NA)
converted 'ARM D' from integer to character (0 new NA)
# A tibble: 1 x 6
block label `ARM A` `ARM B` `ARM C`
1 counts Number of ~ 5 11 16
# ... with 1 more variable: ARM D
Count patients with PSGA > 1
filter: removed 43 rows (29%), 104 rows remaining
count: now 4 rows and 2 columns, ungrouped
pivot_wider: reorganized (TRTA, n) into (ARM A, ARM B, ARM C, ARM D) [was 4x2, now 1x4]
transmute: new variable 'block' (character) with one unique value and 0% NA
new variable 'label' (character) with one unique value and 0% NA
converted 'ARM A' from integer to character (0 new NA)
converted 'ARM B' from integer to character (0 new NA)
converted 'ARM C' from integer to character (0 new NA)
converted 'ARM D' from integer to character (0 new NA)
# A tibble: 1 x 6
block label `ARM A` `ARM B` `ARM C`
1 counts Number of ~ 31 26 21
# ... with 1 more variable: ARM D
=========================================================================
Perform Kaplan-Meier Tests
=========================================================================
Create survival vector
[1] 16+ 12+ 12+ 4 12+ 12+ 12+ 6 16+
[10] 4 16+ 16+ 16+ 16+ 16+ 2 16+ 16+
[19] 16+ 16+ 16+ 8 16 16+ 16+ 16+ 16+
[28] 16+ 16+ 16+ 16+ 16+ 16+ 2 16+ 16+
[37] 8 16+ 16+ 16+ 16+ 16+ 16+ 16+ 16+
[46] 16+ 12 16+ 16+ 16+ 12+ 4 16+ 16+
[55] 12+ 2 12 4 16+ 6 6 6 12
[64] 16+ 16+ 16 12+ 12+ 8+ 12+ 12+ 6+
[73] 12 12 16+ 16+ 0+ 6 16+ 16+ 16+
[82] 16 16+ 12+ 12+ 4 16 6 4 6
[91] 16+ 16+ 6 6+ 16+ 8 16+ 16+ 16+
[100] 16+ 16+ 16+ 6 4+ 16+ 12+ 6 16+
[109] 8 2 12+ 8 16+ 16+ 16 6 16+
[118] 12+ 2+ 16+ 4 16+ 16+ 12+ 16+ 16+
[127] 16+ 12 16 16+ 16+ 12+ 6 16+ 16+
[136] 2 4 16+ 16+ 16+ 16+ 12+ 12+ 16+
[145] 16+ 16+ 6
Fit model on survival vector
Call: survfit(formula = surv_vct ~ TRTA, data = adslpsga_final)
n events median 0.95LCL 0.95UCL
TRTA=ARM A 36 5 NA NA NA
TRTA=ARM B 37 11 NA NA NA
TRTA=ARM C 37 16 NA 16 NA
TRTA=ARM D 37 11 NA NA NA
Get model quantiles
Get lower confidence intervals
mutate: new variable 'block' (character) with one unique value and 0% NA
new variable 'TRTA' (character) with 4 unique values and 0% NA
pivot_longer: reorganized (25, 50, 75) into (Q, lower) [was 4x6, now 12x5]
# A tibble: 12 x 5
rowname block TRTA Q lower
1 TRTA=ARM A surv ARM A 25 NA
2 TRTA=ARM A surv ARM A 50 NA
3 TRTA=ARM A surv ARM A 75 NA
4 TRTA=ARM B surv ARM B 25 6
5 TRTA=ARM B surv ARM B 50 NA
6 TRTA=ARM B surv ARM B 75 NA
7 TRTA=ARM C surv ARM C 25 6
8 TRTA=ARM C surv ARM C 50 16
9 TRTA=ARM C surv ARM C 75 NA
10 TRTA=ARM D surv ARM D 25 6
11 TRTA=ARM D surv ARM D 50 NA
12 TRTA=ARM D surv ARM D 75 NA
Get upper confidence intervals
mutate: new variable 'block' (character) with one unique value and 0% NA
new variable 'TRTA' (character) with 4 unique values and 0% NA
pivot_longer: reorganized (25, 50, 75) into (Q, upper) [was 4x6, now 12x5]
# A tibble: 12 x 5
rowname block TRTA Q upper
1 TRTA=ARM A surv ARM A 25 NA
2 TRTA=ARM A surv ARM A 50 NA
3 TRTA=ARM A surv ARM A 75 NA
4 TRTA=ARM B surv ARM B 25 NA
5 TRTA=ARM B surv ARM B 50 NA
6 TRTA=ARM B surv ARM B 75 NA
7 TRTA=ARM C surv ARM C 25 16
8 TRTA=ARM C surv ARM C 50 NA
9 TRTA=ARM C surv ARM C 75 NA
10 TRTA=ARM D surv ARM D 25 NA
11 TRTA=ARM D surv ARM D 50 NA
12 TRTA=ARM D surv ARM D 75 NA
Get confidence intervals
inner_join: added one column (upper)
> rows only in x ( 0)
> rows only in y ( 0)
> matched rows 12
> ====
> rows total 12
mutate: new variable 'ci' (character) with 4 unique values and 0% NA
pivot_wider: reorganized (rowname, TRTA, lower, upper, ci) into (ARM A, ARM B, ARM C, ARM D) [was 12x7, now 3x6]
mutate: new variable 'order' (double) with one unique value and 0% NA
new variable 'label1' (character) with 3 unique values and 0% NA
new variable 'label2' (character) with one unique value and 0% NA
select: columns reordered (block, Q, order, label1, label2, )
# A tibble: 3 x 9
block Q order label1 label2 `ARM A`
1 surv 25 2 25th p~ 95% Co~ (-, -)
2 surv 50 2 Median ~ 95% Co~ (-, -)
3 surv 75 2 75th p~ 95% Co~ (-, -)
# ... with 3 more variables: ARM B ,
# ARM C , ARM D
Get quantiles
mutate: new variable 'block' (character) with one unique value and 0% NA
new variable 'TRTA' (character) with 4 unique values and 0% NA
pivot_longer: reorganized (25, 50, 75) into (Q, value) [was 4x6, now 12x5]
pivot_wider: reorganized (rowname, TRTA, value) into (ARM A, ARM B, ARM C, ARM D) [was 12x5, now 3x6]
mutate: converted 'ARM A' from double to character (0 new NA)
converted 'ARM B' from double to character (0 new NA)
converted 'ARM C' from double to character (0 new NA)
converted 'ARM D' from double to character (0 new NA)
new variable 'order' (double) with one unique value and 0% NA
new variable 'label1' (character) with 3 unique values and 0% NA
new variable 'label2' (character) with one unique value and 0% NA
select: columns reordered (block, Q, order, label1, label2, )
# A tibble: 3 x 9
block Q order label1 label2 `ARM A`
1 surv 25 1 25th pe~ ""
2 surv 50 1 Median (~ ""
3 surv 75 1 75th pe~ ""
# ... with 3 more variables: ARM B ,
# ARM C , ARM D
Final arrangement
transmute: dropped 2 variables (Q, order)
changed 3 values (50%) of 'label2' (3 new NA)
# A tibble: 6 x 7
block label1 label2 `ARM A` `ARM B`
1 surv 25th per~ 12
2 surv 25th per~ 95% Conf~ (-, -) (6, -)
3 surv Median (w~
4 surv Median (w~ 95% Conf~ (-, -) (-, -)
5 surv 75th per~
6 surv 75th per~ 95% Conf~ (-, -) (-, -)
# ... with 2 more variables: ARM C ,
# ARM D
=========================================================================
Perform Cox Proportional Hazards Test
=========================================================================
Run Cox tests
Call:
survival::coxph(formula = surv_vct ~ TRTA, data = adslpsga_final)
coef exp(coef) se(coef) z
TRTAARM B 0.9185 2.5055 0.5395 1.703
TRTAARM C 1.3111 3.7102 0.5127 2.557
TRTAARM D 0.8648 2.3744 0.5394 1.603
p
TRTAARM B 0.0886
TRTAARM C 0.0105
TRTAARM D 0.1089
Likelihood ratio test=7.87 on 3 df, p=0.04877
n= 147, number of events= 43
Create summary statistics on Cox results
Call:
survival::coxph(formula = surv_vct ~ TRTA, data = adslpsga_final)
n= 147, number of events= 43
coef exp(coef) se(coef) z
TRTAARM B 0.9185 2.5055 0.5395 1.703
TRTAARM C 1.3111 3.7102 0.5127 2.557
TRTAARM D 0.8648 2.3744 0.5394 1.603
Pr(>|z|)
TRTAARM B 0.0886 .
TRTAARM C 0.0105 *
TRTAARM D 0.1089
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 .
0.1 1
exp(coef) exp(-coef) lower .95
TRTAARM B 2.505 0.3991 0.8704
TRTAARM C 3.710 0.2695 1.3584
TRTAARM D 2.374 0.4212 0.8249
upper .95
TRTAARM B 7.212
TRTAARM C 10.134
TRTAARM D 6.834
Concordance= 0.61 (se = 0.042 )
Likelihood ratio test= 7.87 on 3 df, p=0.05
Wald test = 6.68 on 3 df, p=0.08
Score (logrank) test = 7.32 on 3 df, p=0.06
Extract coefficients
mutate: new variable 'block' (character) with one unique value and 0% NA
new variable 'TRTA' (character) with 3 unique values and 0% NA
rowname coef exp(coef) se(coef)
1 TRTAARM B 0.9184774 2.505473 0.5394543
2 TRTAARM C 1.3110748 3.710159 0.5126631
3 TRTAARM D 0.8647564 2.374428 0.5394037
z Pr(>|z|) block TRTA
1 1.702605 0.08864206 surv ARM B
2 2.557381 0.01054637 surv ARM C
3 1.603171 0.10889687 surv ARM D
Extract confidence intervals
rowname exp(coef) exp(-coef) lower .95
1 1 2.505473 0.3991263 0.8703729
2 2 3.710159 0.2695302 1.3583539
3 3 2.374428 0.4211541 0.8249311
upper .95
1 7.212304
2 10.133796
3 6.834396
Create cox statistics block
rename: renamed 4 variables (hazard, pval, lower, upper)
select: dropped 8 variables (rowname...1, coef, se(coef), z, block, )
mutate: converted 'hazard' from double to character (0 new NA)
converted 'pval' from double to character (0 new NA)
new variable 'block' (character) with one unique value and 0% NA
new variable 'ci' (character) with 3 unique values and 0% NA
pivot_longer: reorganized (hazard, pval, ci) into (stat, value) [was 3x7, now 9x6]
pivot_wider: reorganized (TRTA, lower, upper, value) into (ARM B, ARM C, ARM D) [was 9x6, now 3x5]
mutate: new variable 'label1' (character) with 3 unique values and 0% NA
new variable 'label2' (character) with one unique value and 100% NA
new variable 'ARM A' (character) with one unique value and 100% NA
select: dropped one variable (stat)
# A tibble: 3 x 7
block label1 label2 `ARM A` `ARM B`
1 cox Hazard Ratio~ 2.51
2 cox P-value 0.089
3 cox 95% CI of Ha~ (0.87,~
# ... with 2 more variables: ARM C ,
# ARM D
Combine statistics blocks
# A tibble: 9 x 7
block label1 label2 `ARM A` `ARM B`
1 surv 25th perc~ 12
2 surv 25th perc~ 95% Con~ (-, -) (6, -)
3 surv Median (we~
4 surv Median (we~ 95% Con~ (-, -) (-, -)
5 surv 75th perc~
6 surv 75th perc~ 95% Con~ (-, -) (-, -)
7 cox Hazard Rat~ 2.51
8 cox P-value 0.089
9 cox 95% CI of ~ (0.87,~
# ... with 2 more variables: ARM C ,
# ARM D
=========================================================================
Create survival plot
=========================================================================
Create data frame with zero values for each visit
Create visit template
ARM AVISIT AVISITN cured
1 ARM A Week 16 16 FALSE
2 ARM A Week 12 12 FALSE
3 ARM A Week 4 4 FALSE
4 ARM A Week 6 6 FALSE
5 ARM A Week 2 2 FALSE
6 ARM A Week 8 8 FALSE
7 ARM A Week 0 0 FALSE
8 ARM B Week 16 16 FALSE
9 ARM B Week 12 12 FALSE
10 ARM B Week 4 4 FALSE
11 ARM B Week 6 6 FALSE
12 ARM B Week 2 2 FALSE
13 ARM B Week 8 8 FALSE
14 ARM B Week 0 0 FALSE
15 ARM C Week 16 16 FALSE
16 ARM C Week 12 12 FALSE
17 ARM C Week 4 4 FALSE
18 ARM C Week 6 6 FALSE
19 ARM C Week 2 2 FALSE
20 ARM C Week 8 8 FALSE
21 ARM C Week 0 0 FALSE
22 ARM D Week 16 16 FALSE
23 ARM D Week 12 12 FALSE
24 ARM D Week 4 4 FALSE
25 ARM D Week 6 6 FALSE
26 ARM D Week 2 2 FALSE
27 ARM D Week 8 8 FALSE
28 ARM D Week 0 0 FALSE
Calculate cummulative sum and percent
select: dropped 7 variables (USUBJID, SEX, AGEGR1, AGE, TRTA, )
group_by: 3 grouping variables (ARM, AVISIT, AVISITN)
summarize: now 29 rows and 4 columns, 2 group variables remaining (ARM, AVISIT)
group_by: one grouping variable (ARM)
mutate (grouped): changed 4 values (14%) of 'AVISIT' (0 new NA)
new variable 'csumc' (integer) with 13 unique values and 0% NA
distinct (grouped): removed one row (3%), 28 rows remaining
mutate (grouped): new variable 'pct' (double) with 16 unique values and 0% NA
# A tibble: 28 x 6
# Groups: ARM [4]
ARM AVISIT AVISITN sumc csumc pct
1 ARM A Day 1 B~ 0 0 0 0
2 ARM A Week 2 2 0 0 0
3 ARM A Week 4 4 1 1 0.0278
4 ARM A Week 6 6 3 4 0.111
5 ARM A Week 8 8 1 5 0.139
6 ARM A Week 12 12 0 5 0.139
7 ARM A Week 16 16 0 5 0.139
8 ARM B Day 1 B~ 0 0 0 0
9 ARM B Week 2 2 1 1 0.0263
10 ARM B Week 4 4 4 5 0.132
# ... with 18 more rows
Generate plot
=========================================================================
Create and print report
=========================================================================
Create Report
Write out the report
# A report specification: 2 pages
- file_path: 'output/example5.rtf'
- output_type: RTF
- units: inches
- orientation: landscape
- margins: top 1 bottom 1 left 1 right 1
- line size/count: 9/40
- page_header: left=Sponsor: Client right=Study: ABC/BBC
- footnote 1: 'Program: Surv_Table.R'
- footnote 2: '* Success: PSGA <= 1: PSGA > 1'
- footnote 3: '** Based on R survival package survfit() function'
- footnote 4: '*** Based on proportional hazards model with treatment indicator variables as explanatory variables'
- footnote 5: 'Note: The end-point is cure of the disease (PSGA <= 1).'
- footnote 6: ' The probability of remaining diseased (PSGA > 1) defines the survival function'
- footnote 7: 'Note: A subject who is not cured by the end of 12 weeks or is lost to follow provides a censored observation for the analysis.'
- footnote 8: '"-" = Not Applicable'
- page_footer: left=Date Produced: 21Nov21 15:30 center= right=Page [pg] of [tpg]
- content:
# A table specification:
- data: tibble 'count_block' 2 rows 6 cols
- show_cols: all
- use_attributes: all
- width: 9
- title 1: 'Table 5.0'
- title 2: 'Analysis of Time to Initial PSGA Success* in Weeks'
- title 3: 'Safety Population'
- define: block visible='FALSE'
- define: label '' width=4.25
- define: ARM A
- define: ARM B
- define: ARM C
- define: ARM D
# A table specification:
- data: tibble 'stat_block' 9 rows 7 cols
- show_cols: all
- use_attributes: all
- width: 9
- headerless: TRUE
- stub: block label1 label2 width=4.25 align='left'
- define: block dedupe='TRUE'
- define: label1
- define: label2
- define: ARM A
- define: ARM B
- define: ARM C
- define: ARM D
# A plot specification:
- data: 28 rows, 7 cols
- layers: 2
- height: 3.5
- width: 9
- title 1: 'Figure 5.0'
- title 2: 'Kaplan-Meier Plot for Time to Initial PSGA Success (PSGA <= 1)'
- title 3: 'Safety Population'
=========================================================================
Clean Up
=========================================================================
=========================================================================
Log End Time: 2021-11-21 15:30:39
Log Elapsed Time: 0 00:00:02
=========================================================================
```
Next: [Example 8: Patient Profile](sassy-profile.html)