Skip to content

Commit

Permalink
Merge pull request #8 from hta-pharma/fix_checks
Browse files Browse the repository at this point in the history
Fix checks
  • Loading branch information
hoppanda authored Jul 7, 2023
2 parents 75ce2bb + 5352a70 commit 4a97fbf
Show file tree
Hide file tree
Showing 27 changed files with 518 additions and 58 deletions.
1 change: 1 addition & 0 deletions .Rbuildignore
Original file line number Diff line number Diff line change
Expand Up @@ -20,3 +20,4 @@ images
^maicplus\.Rcheck$
^maicplus.*\.tar\.gz$
^maicplus.*\.tgz$
^inst/dev/MAIC/MAIC-main/.Rbuildignore
12 changes: 4 additions & 8 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -6,21 +6,17 @@ Date: 2023-03-03
Authors@R:
person("hta-pharma", , , "hta-pharma@example.com", role = c("aut", "cre"))
Description: R package template with GitHub Actions workflows included.
License: Apache License 2.0 | file LICENSE
License: Apache License 2.0
URL: https://github.com/hta-pharma/maicplus/
BugReports: https://github.com/hta-pharma/maicplus/issues
Depends:
R (>= 3.6)
Imports:
plumber,
shiny,
survival,
stringr
graphics,
stats,
survival
Suggests:
future,
httr,
knitr,
shinytest,
testthat (>= 2.0)
VignetteBuilder:
knitr
Expand Down
23 changes: 23 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,2 +1,25 @@
# Generated by roxygen2: do not edit by hand

export(bucher)
export(center_ipd)
export(estimate_weights)
export(log_cum_haz_plot)
export(maic_tte_unanchor)
export(medSurv_makeup)
export(plot_weights)
export(process_agd)
export(process_ipd)
export(resid_plot)
export(survfit_makeup)
import(stats)
importFrom(graphics,abline)
importFrom(graphics,axis)
importFrom(graphics,hist)
importFrom(graphics,legend)
importFrom(graphics,lines)
importFrom(graphics,par)
importFrom(graphics,points)
importFrom(survival,Surv)
importFrom(survival,cox.zph)
importFrom(survival,coxph)
importFrom(survival,survfit)
36 changes: 21 additions & 15 deletions R/bucher.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,33 +4,38 @@

#' Bucher method for adjusted treatment effect
#'
#' Given two estimated treatment effects of A vs. C and B vs. C (i.e. two point estimates and corresponding standard errors),
#' derive the adjusted treatment of A vs. B using Bucher method, with two-sided confidence limits and Z-test p-value
#' Given two estimated treatment effects of A vs. C and B vs. C
#' (i.e. two point estimates and corresponding standard errors),
#' derive the adjusted treatment of A vs. B using Bucher method,
#' with two-sided confidence limits and Z-test p-value
#'
#' @param trt a list with two named scalars for the study with interested experimental arm, one named 'est' for the point estimate, and the other named 'se' for the standard error
#' @param trt a list with two named scalars for the study with interested experimental arm,
#' one named `'est'` for the point estimate, and the other named `'se'` for the standard error
#' @param com same as \code{trt}, but for the study with interested control arm
#' @param conf_lv a numerical scalar, prescribe confidence level to derive two-sided confidence interval for the adjusted treatment effect
#' @param conf_lv a numerical scalar, prescribe confidence level to derive two-sided
#' confidence interval for the adjusted treatment effect
#'
#' @return a list with 5 elements,
#' \describe{
#' \item est - a scalar, point estimate of the adjusted treatment effect
#' \item se - a scalar, standard error of the adjusted treatment effect (i.e. \code{est} in return)
#' \item ci_l - a scalor, lower confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}
#' \item ci_u - a scalor, upper confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}
#' \item pval - p-value of Z-test, with null hypothesis that \code{est} is zero
#' \item{est}{a scalar, point estimate of the adjusted treatment effect}
#' \item{se}{a scalar, standard error of the adjusted treatment effect (i.e. \code{est} in return)}
#' \item{ci_l}{a scalar, lower confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}}
#' \item{ci_u}{a scalar, upper confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}}
#' \item{pval}{p-value of Z-test, with null hypothesis that \code{est} is zero}
#' }
#' @export
#' @import stats
#'
#' @examples
bucher <- function(trt, com, conf_lv = 0.95) {
est <- trt$est - com$est
se <- sqrt(trt$se^2 + com$se^2)
ci_l <- est - qnorm(0.5 + conf_lv / 2) * se
ci_u <- est + qnorm(0.5 + conf_lv / 2) * se
ci_l <- est - stats::qnorm(0.5 + conf_lv / 2) * se
ci_u <- est + stats::qnorm(0.5 + conf_lv / 2) * se
if (est > 0) {
pval <- 2 * (1 - pnorm(est, 0, se))
pval <- 2 * (1 - stats::pnorm(est, 0, se))
} else {
pval <- 2 * pnorm(est, 0, se)
pval <- 2 * stats::pnorm(est, 0, se)
}

list(
Expand All @@ -44,13 +49,14 @@ bucher <- function(trt, com, conf_lv = 0.95) {

# functions NOT to be exported ---------------------------------------

#' Report-friendly output format for resutl from Bucher's method
#' Report-friendly output format for result from Bucher's method
#'
#' @param output output from \code{\link{bucher}} function
#' @param ci_digits an integer, number of decimal places for point estimate and derived confidence limits
#' @param pval_digits an integer, number of decimal places to display Z-test p-value
#'
#' @return a character vector of two elements, first element inf format of 'est (ci_l; ci_u)', second element is Z-test p-value, rounded according to \code{pval_digits}
#' @return a character vector of two elements, first element in format of 'est (ci_l; ci_u)'`,
#' second element is the Z-test p-value, rounded according to \code{pval_digits}

print_bucher <- function(output, ci_digits = 2, pval_digits = 3) {
res <- paste0(
Expand Down
16 changes: 9 additions & 7 deletions R/maic_unanchored_tte.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
#'
#' @param useWt a numeric vector of individual MAIC weights, length should the same as \code{nrow(dat)}
#' @param dat a data frame that meet format requirements in 'Details', individual patient data (IPD) of internal trial
#' @param dat_ext a data frame, pseudo IPD from digitalized KM curve of external trial
#' @param dat_ext a data frame, pseudo IPD from digitized KM curve of external trial
#' @param trt a character string, name of the interested treatment in internal trial (real IPD)
#' @param trt_ext character string, name of the interested comparator in external trial used to subset \code{dat_ext} (pseudo IPD)
#' @param endpoint_name a character string, name of the endpoint
Expand All @@ -13,11 +13,12 @@
#' \itemize{
#' \item treatment - character or factor column
#' \item status - logical column, TRUE for censored/death, FALSE for otherwise
#' \time time - numeric column, observation time of the \code{status}; unit in days
#' \item time - numeric column, observation time of the \code{status}; unit in days
#' }
#'
#' @return
#' @importFrom survival Surv survfit coxph cox.zph
#' @importFrom graphics par axis lines points legend abline
#' @export
#'
#' @examples
Expand Down Expand Up @@ -50,10 +51,10 @@ maic_tte_unanchor <- function(useWt, dat, dat_ext, trt, trt_ext,
# derive km w and w/o weights
kmobj <- survfit(Surv(time, status) ~ treatment, dat, conf.type = "log-log")
kmobj_adj <- survfit(Surv(time, status) ~ treatment, dat, weights = dat$weight, conf.type = "log-log")

par(cex.main=0.85)
km_makeup(kmobj, kmobj_adj, time_scale = time_scale,
trt = trt, trt_ext = trt_ext,
km_makeup(kmobj, kmobj_adj, time_scale = time_scale,
trt = trt, trt_ext = trt_ext,
endpoint_name = endpoint_name)
res[["plot_km"]] <- grDevices::recordPlot()

Expand Down Expand Up @@ -131,13 +132,14 @@ maic_tte_unanchor <- function(useWt, dat, dat_ext, trt, trt_ext,
}

#' helper function: sort out a nice report table to summarize survival analysis results
#' for maic_tte_unanchor
#' for `maic_tte_unanchor`
#'
#' @param coxobj returned object from \code{\link[survival]{coxph}}
#' @param medSurvobj returned object from \code{\link{medSurv_makeup}}
#' @param tag a string, by default NULL, if specified, an extra 1st column is created in the output
#'
#' @return a data frame with sample size, incidence rate, median survival time with 95%CI, hazard ratio estimate with 95%CI and wald test of hazard ratio
#' @return a data frame with sample size, incidence rate, median survival time with 95% CI, hazard ratio estimate with
#' 95% CI and Wald test of hazard ratio

report_table <- function(coxobj, medSurvobj, tag = NULL) {
hr_res <- format(round(summary(coxobj)$conf.int[-2], 2), nsmall = 2)
Expand Down
20 changes: 12 additions & 8 deletions R/matching.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,16 +10,18 @@
#' @param data a numeric matrix, centered effect modifiers of IPD, no missing value in any cell is allowed
#' @param centered_colnames a character or numeric vector, column indicators of centered effect modifiers, by default NULL meaning all columns in \code{data} are effect modifiers
#' @param startVal a scalar, the starting value for all coefficients of the propensity score regression
#' @param method a string, name of the optimization algorithm (see 'method' argument of \code{base::optim()}). The default is "BFGS", other options are "Nelder-Mead", "CG", "L-BFGS-B", "SANN", and "Brent"
#' @param method a string, name of the optimization algorithm (see 'method' argument of \code{base::optim()}).
#' The default is `"BFGS"`, other options are `"Nelder-Mead"`, `"CG"`, `"L-BFGS-B"`, `"SANN"`, and `"Brent"`
#' @param ... all other arguments from \code{base::optim()}
#'
#' @return a list with the following 4 elements,
#' \describe{
#' \item data - a data.frame, includes the input \code{data} with appended column 'weights' and 'scaled_weights'. Scaled weights has a summation to be the number of rows in \code{data} that has no missing value in any of the effect modifiers
#' \item centered_colnames - column names of centered effect modifiers in \code{data}
#' \item nr_missing - number of rows in \code{data} that has at least 1 missing value in specified centered effect modifiers
#' \item ess - effective sample size, square of sum divided by sum of squares
#' \item opt - R object returned by \code{base::optim()}, for assess convergence and other details
#' \item{data}{a data.frame, includes the input \code{data} with appended column 'weights' and 'scaled_weights'.
#' Scaled weights has a summation to be the number of rows in \code{data} that has no missing value in any of the effect modifiers}
#' \item{centered.colnames}{column names of centered effect modifiers in \code{data}}
#' \item{nr_missing}{number of rows in \code{data} that has at least 1 missing value in specified centered effect modifiers}
#' \item{ess}{effective sample size, square of sum divided by sum of squares}
#' \item{opt}{R object returned by \code{base::optim()}, for assess convergence and other details}
#' }
#' @export
#'
Expand Down Expand Up @@ -94,14 +96,16 @@ estimate_weights <- function(data, centered_colnames = NULL, startVal = 0, metho

#' Plot MAIC weights in a histogram with key statistics in legend
#'
#' Generates a plot given the individuals weights with key summary in top right legend that includes
#' median weight, effective sample size (ESS), reduction percentage (what percent ESS takes up in the original sample size),
#' Generates a plot given the individuals weights with key summary in top right legend that
#' includes median weight, effective sample size (ESS), reduction percentage (what percent
#' ESS takes up in the original sample size),
#'
#'
#' @param wt a numeric vector of individual MAIC weights (derived use in \code{\link{cal_weights}})
#' @param main_title a character string, main title of the plot
#'
#' @return a plot
#' @importFrom graphics hist
#' @export
plot_weights <- function(wt, main_title = "Unscaled Individual Weigths") {

Expand Down
16 changes: 8 additions & 8 deletions R/process_data.R
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ process_agd <- function(raw_agd) {
with(use_agd,{ use_agd[tolower(ARM)=="total",,drop=FALSE] })
}

#' Process Individual Patient data to dummize categorical effective modifiers
#' Process Individual Patient data categorical effective modifiers to dummy variables
#'
#' @param raw_ipd
#' @param dummize_cols
Expand Down Expand Up @@ -96,7 +96,7 @@ process_ipd <- function(raw_ipd, dummize_cols, dummize_ref_level){
#'
#' @return
#' @export
center_ipd <- function(ipd,agd){
center_ipd <- function(ipd, agd){
# regulaized column name patterns
must_exist <- c("STUDY","ARM", "N")
legal_suffix <- c("MEAN","MEDIAN","SD","PROP")
Expand Down Expand Up @@ -146,10 +146,10 @@ center_ipd <- function(ipd,agd){

#' helper function: transform TTE ADaM data to suitable input for survival R pkg
#'
#' @param dd data frame, ADTTE read via haven::read_sas
#' @param dd data frame, ADTTE read via `haven::read_sas()`
#' @param time_scale a character string, 'year', 'month', 'week' or 'day', time unit of median survival time
#'
#' @return a data frame can be used as input to survival::Surv
#' @return a data frame can be used as input to [survival::Surv()]
ext_tte_transfer <- function(dd, time_scale = "month", trt = NULL) {
timeUnit <- list("year" = 365.24, "month" = 30.4367, "week" = 7, "day" = 1)

Expand All @@ -175,7 +175,7 @@ ext_tte_transfer <- function(dd, time_scale = "month", trt = NULL) {

# functions NOT to be exported ---------------------------------------

#' Calculate pooled arm statistics in AgD based on arm-specific statistics
#' Calculate pooled arm statistics in Aggregated Data (AgD) based on arm-specific statistics
#'
#' @param use_agd
#'
Expand All @@ -192,8 +192,8 @@ complete_agd <- function(use_agd) {
use_agd$ARM[rowId] <- "total"

# complete N and count
NN <- use_agd$N[rowId] <- sum(use_agd$N, na.rm=TRUE)
nn <- use_agd$N[-rowId]
NN <- use_agd[["N"]][rowId] <- sum(use_agd[["N"]], na.rm=TRUE)
nn <- use_agd[["N"]][-rowId]
for(i in grep("_COUNT$",names(use_agd),value=TRUE)){
use_agd[[i]][rowId] <- sum(use_agd[[i]], na.rm=TRUE)
}
Expand All @@ -205,7 +205,7 @@ complete_agd <- function(use_agd) {

# complete SD
for(i in grep("_SD$",names(use_agd),value=TRUE)){
use_agd[[i]][rowId] <- sqrt( sum(use_agd[[i]]^2*(nn-1))/(N-1) )
use_agd[[i]][rowId] <- sqrt( sum(use_agd[[i]]^2*(nn-1))/(NN-1) )
}

# complete MEDIAN, approximately!!
Expand Down
12 changes: 6 additions & 6 deletions R/survival-helper.R
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
#' helper function: makeup to get median survival time from a survival::survfit object
#' helper function: makeup to get median survival time from a `survival::survfit` object
#'
#' extract and display median survival time with confidence interval
#'
#' @param km_fit returned object from \code{survival::survfit}
#' @param legend a character string, name used in 'type' column in returned data frame
#' @param time_scale a character string, 'year', 'month', 'week' or 'day', time unit of median survival time
#'
#' @return a data frame with a index column 'type', mdian survival time and confidence interval
#' @return a data frame with a index column 'type', median survival time and confidence interval
#' @export
medSurv_makeup <- function(km_fit, legend = "before matching", time_scale) {
timeUnit <- list("year" = 365.24, "month" = 30.4367, "week" = 7, "day" = 1)
Expand All @@ -30,7 +30,7 @@ medSurv_makeup <- function(km_fit, legend = "before matching", time_scale) {



#' helper function: makeup survival::survfit object for km plot
#' helper function: makeup `survival::survfit` object for km plot
#'
#' @param km_fit returned object from \code{survival::survfit}
#'
Expand Down Expand Up @@ -210,9 +210,9 @@ log_cum_haz_plot <- function(clldat, time_scale, log_time = TRUE, endpoint_name
}


#' Plot schoenfeld residual plot for a Cox model fit
#' Plot Schoenfeld residuals for a Cox model fit
#'
#' @param coxobj object returned from \code{\link[surival]{coxph}}
#' @param coxobj object returned from \code{\link[survival]{coxph}}
#' @param time_scale a character string, 'year', 'month', 'week' or 'day', time unit of median survival time
#' @param log_time logical, TRUE (default) or FALSE
#' @param endpoint_name a character string, name of the endpoint
Expand All @@ -233,7 +233,7 @@ resid_plot <- function(coxobj, time_scale = "month", log_time = TRUE, endpoint_n
cex = 0.9, col = "navyblue", yaxt = "n",
ylab = "Unscaled Schoenfeld Residual", xlab = paste0(ifelse(log_time, "Log-", ""), "Time in ", time_scale),
main = paste0(
"Diagnosis Plot: Unscaled Schoefeld Residual\nEndpoint: ", endpoint_name,
"Diagnosis Plot: Unscaled Schoenfeld Residual\nEndpoint: ", endpoint_name,
ifelse(subtitle == "", "", "\n"), subtitle
)
)
Expand Down
14 changes: 14 additions & 0 deletions inst/WORDLIST
Original file line number Diff line number Diff line change
Expand Up @@ -10,3 +10,17 @@ initializer
pre
repo
Forkers
comparator
Bucher
Bucher's
unadjusted
MAIC
ci
ESS
IPD
Schoenfeld
ADTTE
ADaM
unanchored
TTE
AgD
12 changes: 6 additions & 6 deletions inst/dev/maicplus/maicplus-main/bucher.R
Original file line number Diff line number Diff line change
Expand Up @@ -5,24 +5,24 @@
#'
#' @param trt a list with two named scalars for the study with interested experimental arm, one named 'est' for the point estimate, and the other named 'se' for the standard error
#' @param com same as \code{trt} for the study with interested control arm
#' @param conf.lv a numerical scalar, prescribe confidence level to derive two-sided confidence interval for the adjusted treatment effect
#' @param conf_lv a numerical scalar, prescribe confidence level to derive two-sided confidence interval for the adjusted treatment effect
#'
#' @return a list with 5 elements,
#' \describe{
#' \item est - a scalar, point estimate of the adjusted treatment effect
#' \item se - a scalar, standard error of the adjusted treatment effect (i.e. \code{est} in return)
#' \item ci_l - a scalor, lower confidence limit of a two-sided CI with prescribed nominal level by \code{conf.lv}
#' \item ci_u - a scalor, upper confidence limit of a two-sided CI with prescribed nominal level by \code{conf.lv}
#' \item ci_l - a scalor, lower confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}
#' \item ci_u - a scalor, upper confidence limit of a two-sided CI with prescribed nominal level by \code{conf_lv}
#' \item pval - p-value of Z-test, with null hypothesis that \code{est} is zero
#' }
#' @export
#'
#' @examples
bucher <- function(trt, com, conf.lv = 0.95) {
bucher <- function(trt, com, conf_lv = 0.95) {
est <- trt$est - com$est
se <- sqrt(trt$se^2 + com$se^2)
ci_l <- est - qnorm(0.5 + conf.lv / 2) * se
ci_u <- est + qnorm(0.5 + conf.lv / 2) * se
ci_l <- est - qnorm(0.5 + conf_lv / 2) * se
ci_u <- est + qnorm(0.5 + conf_lv / 2) * se
if (est > 0) {
pval <- 2 * (1 - pnorm(est, 0, se))
} else {
Expand Down
Loading

0 comments on commit 4a97fbf

Please sign in to comment.