Nothing
#' L'Abbe plot for subgroup effect size
#'
#' this function produces a L'Abbe plot showing the treatment effect size of subgroups defined by the categories of covariates. The
#' x-axis and y-axis represent the treatment effect estimate from the control and treatment group, respectively. The dashed dignal
#' indicates no effect, and the solid diagnal line corresponds to the full population effect estimate. The squares represent subgroups,
#' and the red or blue dashed lines have a length which show the magnitude of the lower or upper bound of the 95% C.I for subgroup
#' treatment effect estimates. If the solid diagonal line cross all the blue or red lines, it may indicate homogeneity across subgroup
#' with repective to the full population effect estimate. Note that the overall size of squares which represent subgroups can be
#' adjusted by setting different values on the associated input argument. In addition, the function uses log odd ratio for displaying
#' subgroup effect sizes in binary and survival data, respectively.
#'
#'
#' @param dat a data set
#' @param covari.sel a vector of indices of the two covariates
#' @param trt.sel a covariate index specifying the treatment code
#' @param resp.sel a covariate index specifying the response variable
#' @param outcome.type a string specifying the type of the response variable, it can be "continuous", or "binary".
#' @param adj.ann.subgrp a parameter controlling the distance between a square and its corresponding subgroup label; also, the line
#' gap between the legend text for subgroups.
#' @param size.shape a parameter controlling the height and width of the squares and the size is proportional to the ratio of subgroup sample sizes
#' over the full population size.
#' @param font.size a vector specifying the size of labels and text; the first element is for the main title, the second is for
#' for x-axis and y-axis labels; the thrid is for the legend text of subgroups; the fourth is for the subgroup
#' labels near the corresponding squares.
#' @param title a string specifying the main title.
#' @param lab.xy a list of two strings specifying the labels of the x and y axes.
#' @param time time for calculating the survival in each subgroup
#' @param show.ci A logical indicating whether to show an additional line for confidence intervals
#' @param effect Either "HR" of "RMST". Only when outcome.type = "survival"
#' @param legend.position where to place the legend? either "inside" or "outside"
#'
#' @examples
#' library(dplyr)
#'
#' # Load the data to be used
#' data(prca)
#' dat <- prca
#' dat %>%
#' mutate(bm = factor(ifelse(bm == 0 , "No", "Yes")),
#' hx = factor(ifelse(hx == 0 , "No", "Yes")))-> dat
#'
#' ## 8. Labbe Plot -----------------------------------------------------------
#' lab.xy = list("Control Group Estimate", "Treatment Group Estimate")
#' plot_labbe(dat = dat,
#' covari.sel = c(4,5,6,7),
#' trt.sel = 3,
#' resp.sel = c(1,2),
#' outcome.type = "survival",
#' effect = "RMST",
#' lab.xy = lab.xy,
#' size.shape = 0.2,
#' adj.ann.subgrp = 1/30,
#' time=50, show.ci = FALSE)
#'
#' @export
#' @import grid
#' @import graphics
plot_labbe <- function(dat, covari.sel, trt.sel, resp.sel, outcome.type,
effect = c("HR","RMST"), size.shape = 1/18,
adj.ann.subgrp = 1/30, font.size = c(1, 1, 0.85, 0.85),
title = NULL, lab.xy = NULL,
time = mean(dat[,resp.sel[1]]), show.ci = TRUE,
legend.position = c("inside", "outside"))
{
################################################ 0. argument validity check #################################################################
legend.position = match.arg(legend.position)
if (missing(dat)) stop("Data have not been inputed!")
if (!(is.data.frame(dat))) stop("The data set is not with a data frame!")
if (missing(covari.sel)) stop("The variables for defining subgroups have not been specified!")
if (!(is.numeric(covari.sel))) stop("The variables for defining subgroups are not numeric!")
for (i in 1 : length(covari.sel)) if (!(is.factor(dat[,covari.sel[i]]))) stop("The variables for defining subgroups are not categorical!")
if (missing(trt.sel)) stop("The variable specifying the treatment code (for treatment / control groups) has not been specified!")
if (!(length(trt.sel) == 1)) stop("The variable specifying the treatment code can not have more than one component!")
if (!(is.factor(dat[, trt.sel]))) stop("The variable specifying the treatment code is not categorical!")
if (length(names(table(dat[, trt.sel]))) > 2) stop("The variable specifying the treatment code is not binary!")
if (sum(is.element(names(table(dat[, trt.sel])), c("0","1"))) != 2) stop("The treatment code is not 1 and 0 (for treatment / control groups)!")
type.all = c("continuous", "binary", "survival")
if (is.null(outcome.type)) stop("The type of the response variable has not been specified!")
if (!(is.element(outcome.type, type.all)) == TRUE) stop("A unrecognized type has been inputed!")
if (outcome.type == "continuous"){
if (missing(resp.sel)) stop("The response variable has not been specified!")
if (!(length(resp.sel) == 1)) stop("The response variable has more than one component!")
if (!(is.numeric(dat[, resp.sel]))) stop("The response variable is not numeric!")
}else if (outcome.type == "binary"){
if (missing(resp.sel)) stop("The response variable has not been specified!")
if (!(length(resp.sel) == 1)) stop("The response variable has more than one component!")
if (!(is.factor(dat[, resp.sel]) || is.numeric(dat[, resp.sel]) )) stop("The response variable is not categorical or numerical!")
if (length(names(table(dat[, resp.sel]))) > 2) stop("The response variable is not binary!")
if (sum(is.element(names(table(dat[, resp.sel])), c("0","1"))) != 2) stop(" The response variable is not coded as 0 and 1!")
}
if (!(is.numeric(size.shape))) stop("The argument about the shape setting of diamonds is not numeric!")
if (!(length(size.shape) == 1)) stop("The shape set-up of diamonds is not specified correctly!!")
if (!(is.numeric(adj.ann.subgrp))){
stop("The argument adjusting the distance between a square and its corresponding subgroup label, or adjusting line gap between the
legend text is not numeric!")
}
if (adj.ann.subgrp < 0){
stop("The argument adjusting the distance between a square and its corresponding subgroup label, or adjusting line gap between the
legend text is not positive!")
}
if (!(is.numeric(font.size))) stop("The argument about the font sizes of labels and text is not numeric!")
if (!(length(font.size) == 4)) stop("The font size setups for labels or text should have four components only!")
################################################ 1. create subgroup data #################################################################
n.covari = length(covari.sel)
lab.vars = names(dat)[covari.sel] # set the names of the covariates which relates to the defined subgroup; if a covariate
# are considered for multiple times, we make their name identical. (otherwise, the resulsting
# names are like var, var.1, var.2 and so on.)
names(dat)[trt.sel] = "trt" # rename the variable for treatment code
if (outcome.type == "continuous"){
names(dat)[resp.sel] = "resp" # rename the response variable
}else if (outcome.type == "binary"){
names(dat)[resp.sel] = "resp" # rename the response variable
}else if (outcome.type == "survival"){
names(dat)[resp.sel[1]] = "time" # rename the response variable for survival time
names(dat)[resp.sel[2]] = "status" # rename the response variable for survival right censoring status
}
for (i in 1: length(covari.sel)){
cond = covari.sel == covari.sel[[i]]
lab.vars[cond] = rep(lab.vars[i], length(which(cond == TRUE)))
}
cats.var = list()
n.subgrp.tol = 0
for (i in 1 : length(covari.sel)){
cats.var[[i]] = names(table(dat[,covari.sel[i]]))
n.subgrp.tol = n.subgrp.tol + length(cats.var[[i]]) # the number of the subgroups (excluding the complement)
}
cond = list()
data.subgrp = list()
ss.subgrp = matrix(rep(0, n.subgrp.tol * n.subgrp.tol), nrow = n.subgrp.tol)
ss.subgrp.T = matrix(0, nrow = n.subgrp.tol + 1)
ss.subgrp.C = matrix(0, nrow = n.subgrp.tol + 1)
k = 0
for (i in 1 : length(covari.sel)) {
for (j in 1 : length(cats.var[[i]])){
k = k + 1
cond[[k]] = which((dat[, covari.sel[i]] == cats.var[[i]][j]) == T )
ss.subgrp[k, k] = length(cond[[k]])
data.subgrp[[k]] = dat[cond[[k]], ]
ss.subgrp.T[k] = length(which(data.subgrp[[k]]$trt == 1))
ss.subgrp.C[k] = length(which(data.subgrp[[k]]$trt == 0))
}
}
ss.subgrp.T[n.subgrp.tol + 1] = length(which(dat$trt == 1))
ss.subgrp.C[n.subgrp.tol + 1] = length(which(dat$trt == 0))
k = n.subgrp.tol
r.prop = diag(n.subgrp.tol)
for (i in 1 : (n.subgrp.tol - 1) ){
for (j in (i + 1) : (n.subgrp.tol) ){
k = k + 1
cond[[k]] = intersect(cond[[i]], cond[[j]])
ss.subgrp[i, j] = length(cond[[k]])
ss.subgrp[j, i] = length(cond[[k]])
}
}
# create matrices for treatment size and standard error of MLE
treatment.C.mean = matrix(0, nrow = n.subgrp.tol + 1, ncol = 1)
treatment.T.mean = matrix(0, nrow = n.subgrp.tol + 1, ncol = 1)
treatment.mean = matrix(0, nrow = n.subgrp.tol + 1, ncol = 1)
treatment.upper = matrix(0, nrow = n.subgrp.tol + 1, ncol = 1)
treatment.lower = matrix(0, nrow = n.subgrp.tol + 1, ncol = 1)
for (i in 1:n.subgrp.tol){
if (sum((data.subgrp[[i]]$trt == "1")) == 0 | sum((data.subgrp[[i]]$trt == "0")) == 0){
treatment.mean[i] = NA
treatment.upper[i] = NA
treatment.lower[i] = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = data.subgrp[[i]])
model.sum = summary(model.int)
treatment.mean[i] = model.sum$coefficients[2, 1]
treatment.upper[i] = model.sum$coefficients[2, 1] + 1.96 * model.sum$coefficients[2, 2]
treatment.lower[i] = model.sum$coefficients[2, 1] - 1.96 * model.sum$coefficients[2, 2]
if (length(which(data.subgrp[[i]]$trt == 0)) == 0){
treatment.C.mean[i] = NA
}else{
model.int = lm(resp ~ 1, data = data.subgrp[[i]][which(data.subgrp[[i]]$trt == 0), ])
model.sum = summary(model.int)
treatment.C.mean[i] = model.sum$coefficients[1, 1]
}
if (length(which(data.subgrp[[i]]$trt == 1)) == 0){
treatment.T.mean[i] = NA
}else{
model.int = lm(resp ~ 1, data = data.subgrp[[i]][which(data.subgrp[[i]]$trt == 1), ])
model.sum = summary(model.int)
treatment.T.mean[i] = model.sum$coefficients[1, 1]
}
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = data.subgrp[[i]])
model.sum = summary(model.int)
treatment.mean[i] = model.sum$coefficients[2, 1]
treatment.upper[i] = model.sum$coefficients[2, 1] + 1.96 * model.sum$coefficients[2, 2]
treatment.lower[i] = model.sum$coefficients[2, 1] - 1.96 * model.sum$coefficients[2, 2]
if (length(which(data.subgrp[[i]]$trt == 0)) == 0){
treatment.C.mean[i] = NA
}else{
model.int = glm(resp ~ 1, family = "binomial", data = data.subgrp[[i]][which(data.subgrp[[i]]$trt == 0), ])
model.sum = summary(model.int)
treatment.C.mean[i] = model.sum$coefficients[1, 1]
}
if (length(which(data.subgrp[[i]]$trt == 1)) == 0){
treatment.T.mean[i] = NA
}else{
model.int = glm(resp ~ 1, family = "binomial", data = data.subgrp[[i]][which(data.subgrp[[i]]$trt == 1), ])
model.sum = summary(model.int)
treatment.T.mean[i] = model.sum$coefficients[1, 1]
}
}else if (outcome.type == "survival"){
if (effect == "RMST"){
dat.subgr.i = data.subgrp[[i]]
rmst = survRM2::rmst2(time = dat.subgr.i$time, status = dat.subgr.i$status,
arm = dat.subgr.i$trt, tau = time)
treatment.mean[i] = rmst$unadjusted.result[1,1]
treatment.lower[i] = rmst$unadjusted.result[1,2]
treatment.upper[i] = rmst$unadjusted.result[1,3]
if (length(which(dat$trt == 0)) == 0){
treatment.C.mean[i] = NA
}else{
treatment.C.mean[i] = rmst$RMST.arm0$rmst[1]
}
if (length(which(dat$trt == 1)) == 0){
treatment.T.mean[i] = NA
}else{
treatment.T.mean[i] = rmst$RMST.arm1$rmst[1]
}
}
if (effect == "HR"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = data.subgrp[[i]])
model.sum = summary(model.int)
treatment.mean[i] = model.sum$coef[1, 1]
treatment.upper[i] = model.sum$coef[1, 1] + 1.96 * model.sum$coef[1, 3]
treatment.lower[i] = model.sum$coef[1, 1] - 1.96 * model.sum$coef[1, 3]
surv.fit = survival::survfit(survival::Surv(time, status) ~ trt, data = data.subgrp[[i]])
difference <- summary(surv.fit, time=time)
if (length(which(data.subgrp[[i]]$trt == 0)) == 0){
treatment.C.mean[i] = NA
}else{
treatment.C.mean[i] = difference$surv[1]
}
if (length(which(data.subgrp[[i]]$trt == 1)) == 0){
treatment.T.mean[i] = NA
}else{
treatment.T.mean[i] = difference$surv[2]
}
}
}
}
}
if (sum((dat$trt == "1")) == 0 | sum((dat$trt == "0")) == 0){
treatment.mean[n.subgrp.tol + 1] = NA
treatment.upper[n.subgrp.tol + 1] = NA
treatment.lower[n.subgrp.tol + 1] = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean[n.subgrp.tol + 1] = model.sum$coefficients[2, 1]
treatment.upper[n.subgrp.tol + 1] = model.sum$coefficients[2, 1] + 1.96 * model.sum$coefficients[2, 2]
treatment.lower[n.subgrp.tol + 1] = model.sum$coefficients[2, 1] - 1.96 * model.sum$coefficients[2, 2]
if (length(which(dat$trt == 0)) == 0){
treatment.C.mean[n.subgrp.tol + 1] = NA
}else{
model.int = lm(resp ~ 1, data = dat[which(dat$trt == 0), ])
model.sum = summary(model.int)
treatment.C.mean[n.subgrp.tol + 1] = model.sum$coefficients[1, 1]
}
if (length(which(dat$trt == 1)) == 0){
treatment.T.mean[n.subgrp.tol + 1] = NA
}else{
model.int = lm(resp ~ 1, data = dat[which(dat$trt == 1), ])
model.sum = summary(model.int)
treatment.T.mean[n.subgrp.tol + 1] = model.sum$coefficients[1, 1]
}
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = dat)
model.sum = summary(model.int)
treatment.mean[n.subgrp.tol + 1] = model.sum$coefficients[2, 1]
treatment.upper[n.subgrp.tol + 1] = model.sum$coefficients[2, 1] + 1.96 * model.sum$coefficients[2, 2]
treatment.lower[n.subgrp.tol + 1] = model.sum$coefficients[2, 1] - 1.96 * model.sum$coefficients[2, 2]
if (length(which(dat$trt == 0)) == 0){
treatment.C.mean[n.subgrp.tol + 1] = NA
}else{
model.int = glm(resp ~ 1, family = "binomial", data = dat[which(dat$trt == 0), ])
model.sum = summary(model.int)
treatment.C.mean[n.subgrp.tol + 1] = model.sum$coefficients[1, 1]
}
if (length(which(dat$trt == 1)) == 0){
treatment.T.mean[n.subgrp.tol + 1] = NA
}else{
model.int = glm(resp ~ 1, family = "binomial", data = dat[which(dat$trt == 1), ])
model.sum = summary(model.int)
treatment.T.mean[n.subgrp.tol + 1] = model.sum$coefficients[1, 1]
}
}else if (outcome.type == "survival"){
if (effect == "RMST"){
rmst = survRM2::rmst2(time = dat$time, status = dat$status,
arm = dat$trt, tau = time)
treatment.mean[n.subgrp.tol + 1] = rmst$unadjusted.result[1,1]
treatment.lower[n.subgrp.tol + 1] = rmst$unadjusted.result[1,2]
treatment.upper[n.subgrp.tol + 1] = rmst$unadjusted.result[1,3]
if (length(which(dat$trt == 0)) == 0){
treatment.C.mean[n.subgrp.tol + 1] = NA
}else{
treatment.C.mean[n.subgrp.tol + 1] = rmst$RMST.arm0$rmst[1]
}
if (length(which(dat$trt == 1)) == 0){
treatment.T.mean[n.subgrp.tol + 1] = NA
}else{
treatment.T.mean[n.subgrp.tol + 1] = rmst$RMST.arm1$rmst[1]
}
}
if (effect == "HR"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean[n.subgrp.tol + 1] = model.sum$coef[1, 1]
treatment.upper[n.subgrp.tol + 1] = model.sum$coef[1, 1] + 1.96 * model.sum$coef[1, 3]
treatment.lower[n.subgrp.tol + 1] = model.sum$coef[1, 1] - 1.96 * model.sum$coef[1, 3]
surv.fit = survival::survfit(survival::Surv(time, status) ~ trt, data = dat)
difference <- summary(surv.fit, time=time)
if (length(which(dat$trt == 0)) == 0){
treatment.C.mean[n.subgrp.tol + 1] = NA
}else{
treatment.C.mean[n.subgrp.tol + 1] = difference$surv[1]
}
if (length(which(dat$trt == 1)) == 0){
treatment.T.mean[n.subgrp.tol + 1] = NA
}else{
treatment.T.mean[n.subgrp.tol + 1] = difference$surv[2]
}
}
}
}
lab.subgrp = vector()
lab.subgrp2 = vector()
k = 0
for (i in 1: length(covari.sel)){
for (j in 1 : length(cats.var[[i]])){
k = k + 1
lab.subgrp[k] = paste("(", LETTERS[i], j, ") ", lab.vars[i], " = ", cats.var[[i]][j], sep = "")
lab.subgrp2[k] = paste(LETTERS[i], j, sep = "")
}
}
lab.subgrp[n.subgrp.tol + 1] = "Full"
lab.subgrp2[n.subgrp.tol + 1] = "Full"
dimnames(treatment.C.mean) = list(c(lab.subgrp), c("C.mean"))
dimnames(treatment.T.mean) = list(c(lab.subgrp), c("T.mean"))
dimnames(treatment.mean) = list(c(lab.subgrp), c("mean diff"))
dimnames(treatment.upper) = list(c(lab.subgrp), c("upper"))
dimnames(treatment.lower) = list(c(lab.subgrp), c("lower"))
################################################ 2. produce a graph #################################################################
if (is.null(title)){
old.par = par(mar = c(3,3,0,0)+0.1)
} else{
old.par = par(mar = c(3,3,1,0)+0.1)
}
if(legend.position == "outside"){
layout(matrix(c(1, 2), nrow=1, ncol=2), widths=c(.7,.3))
}
x.lim.min = min(min(treatment.T.mean, na.rm = TRUE), min(treatment.C.mean, na.rm = TRUE)) - 0.2
x.lim.max = max(max(treatment.T.mean, na.rm = TRUE), max(treatment.C.mean, na.rm = TRUE)) + 0.2
eff.col = c("#2c75c9","#fb4b4b")
if (treatment.mean[n.subgrp.tol + 1] > 0 ){
CI.bar.y.pos = treatment.C.mean + abs(treatment.upper)
}else{
CI.bar.y.pos = treatment.C.mean - abs(treatment.lower)
}
if ((outcome.type == "survival" & effect == "HR")){
lab.x=paste0(lab.xy[[1]], " (Survival probability at time=", time,")")
lab.y=paste0(lab.xy[[2]], " (Survival probability at time=", time,")")
} else if ((outcome.type == "survival" & effect == "RMST")){
lab.x=paste0(lab.xy[[1]], " (RMST for tau=", time,")")
lab.y=paste0(lab.xy[[2]], " (RMST for tau=", time,")")
} else {
lab.x=paste0(lab.xy[[1]])
lab.y=paste0(lab.xy[[2]])
}
plot(treatment.C.mean, treatment.T.mean, type='n',
xlab = "",
ylab = "",
main = title,
ylim = c(x.lim.min, x.lim.max),
xlim = c(x.lim.min, x.lim.max),
cex.main = font.size[1],
cex.lab = font.size[2],
bty = "o")
abline(a=0, b = 1, lty = 2)
mtext(text = lab.x, side = 1, line = 2)
mtext(text = lab.y, side = 2, line = 2)
text(max(max(treatment.T.mean, na.rm = TRUE), max(treatment.C.mean, na.rm = TRUE)),
max(max(treatment.T.mean, na.rm = TRUE), max(treatment.C.mean, na.rm = TRUE)),
labels = "x = y")
if (!(outcome.type == "survival" & effect == "HR")){
eff.col = rev(eff.col)
abline(a = treatment.mean[n.subgrp.tol + 1] , b = 1, lty = 1)
}
rectangle = list()
w = vector()
h = vector()
x.range = (max(treatment.C.mean, na.rm = TRUE)+ 0.2) - (min(treatment.C.mean, na.rm = TRUE)-0.2)
y.range = (max(treatment.T.mean, na.rm = TRUE)+ 0.2) - (min(treatment.T.mean, na.rm = TRUE)-0.2)
data.C.size = length(which(dat$trt == 0))
data.T.size = length(which(dat$trt == 1))
col.seg = rep(eff.col[2], n.subgrp.tol + 1)
col.seg[which(treatment.mean < 0 )] = rep(eff.col[1], length(which(treatment.mean < 0 )))
prop.w = (x.lim.max - x.lim.min) * size.shape #0.8 # the shrinkage in width
prop.h = (x.lim.max - x.lim.min) * size.shape #0.8 # the shrinkage in height
for (i in 1: (n.subgrp.tol + 1)){
w[i] = ss.subgrp.C[i]/data.C.size * prop.w
h[i] = ss.subgrp.T[i]/data.T.size * prop.h
rectangle[[i]] = rectangleSP_test(treatment.C.mean[i], treatment.T.mean[i], w[i], h[i])
sp::plot(rectangle[[i]], add=TRUE, col= adjustcolor("white", alpha.f = 0))
segments(treatment.C.mean[i], treatment.C.mean[i], treatment.C.mean[i], treatment.T.mean[i], col = col.seg[i], lwd = 2)
if (show.ci){
segments(treatment.C.mean[i], treatment.C.mean[i], treatment.C.mean[i], CI.bar.y.pos[i], col = "violet", lwd = 3, lty = 2)
}
}
text(treatment.C.mean, treatment.T.mean - (x.lim.max - x.lim.min) * adj.ann.subgrp, labels = lab.subgrp2, pch = 4, cex = font.size[4])
x.pos = vector(); y.pos = vector()
# Draw legend ----------------------------------------------------------------
if(legend.position == "outside"){
par(mar = c(0,0,0,0) + 0.1, plt = c(0,1,0,1), mgp = c(0,0,0))
plot.new()
}
xy.current.pos = par("usr")
lab.pos.adj = (xy.current.pos[2] - xy.current.pos[1]) * adj.ann.subgrp
legend(xy.current.pos[1], xy.current.pos[4],
c("",lab.subgrp[-length(lab.subgrp)]),
bty = "n",
cex = font.size[3])
par(old.par)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.