Nothing
R/07-ggradial-plot-function.r
In SubgrPlots: Graphical Displays for Subgroup Analysis in Clinical Trials
Defines functions
ggplot_radial
ggplot_radial2
Documented in
ggplot_radial
ggplot_radial2
#' Radial (Galbraith) plot for subgroup effect size
#'
#' This function produces a modified Galbraith's radial plot using ggplot. It shows the
#' treatment effect size of subgroups defined by the categories
#' of covariates. The x-axis represents the reciprocal of the standard error
#' of subgroup treatment effect estimates.The y-axis means
#' standardized effect size difference (the difference between subgroup effect
#' the full popultion effect is divided by the standard error
#' of the estimator for the overall population effect.
#' Points here are for subgroups.
#' The grey region indicates whether subgroup effects
#' are homogeneous to the full population effect or not. The two arcs on the
#' right side show subgroup treatment effects in the original
#' scale, where the red spots are the projection of points from the origin on
#' the left side. Note that the vertical range of display can
#' be changed by setting different values on the associated input argument.
#' In addition, the function uses log odd ratio and log hazard
#' ratio for displaying subgroup effect sizes in binary and survival data,
#' respectively.
#' A ggplot object is returned, where further modifications can be added using the theme() function.
#'
#' @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" or "survival".
#' @param range.v a vector specifying the vertical range of graphical display.
#' @param font.size a number specifying the font size for the subgroup labels near points
#' @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 ticks.length a number from 0 to 1 to specify the length of the red
#' ticks indicating the treatment effect in the radial axis. A value of 1 will draw a complete line from the (0,0) point to the axis.
#'
#' @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
#'
#' ggplot_radial(dat,
#' covari.sel = c(4, 5, 6, 7),
#' trt.sel = 3,
#' resp.sel = c(1, 2),
#' outcome.type = "survival",
#' range.v = c(-7, 6),
#' font.size = 4)
#' @export
#' @import ggplot2
#' @import grid
#' @import graphics
#' @importFrom ggrepel geom_text_repel
ggplot_radial <- function(dat, covari.sel, trt.sel, resp.sel, outcome.type,
range.v = NULL,
font.size = 4,
title = NULL,
lab.xy = "default",
ticks.length = 0.05){
################################################ 0. argument validity check #################################################################
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!")
}else if (outcome.type == "survival"){
if (missing(resp.sel)) stop("The response variablehas not been specified!")
if (!(length(resp.sel) == 2)) stop("The response variable for analysing survival data should have two components!")
if (!(is.numeric(dat[, resp.sel[1]]))) stop("The response variable specifying survival time is not numeric!")
if (!(is.numeric(dat[, resp.sel[2]]) || is.logical(dat[, resp.sel[2]]) ) ) stop("The response variable specifying indicators of right censoring should be numerical or logical!")
if (length(names(table(dat[, resp.sel[2]]))) > 2) stop("The response variable specifying indicators of right censoring is not binary!")
if (sum(is.element(names(table(dat[, resp.sel[2]])), c("0","1"))) != 2) stop("The response variable specifying indicators of right censoring is not coded as 0 and 1!")
}
if (!(length(range.v) == 2)) stop("The vertical range of graphical display should have two components (specifying the minimum and maximum!)")
if (!(is.numeric(font.size))) stop("The font.size argument is not numeric!")
if (!(length(font.size) == 1)) stop("The font.size argument should be of length 1!")
if(lab.xy == "default") {
lab.xy = list(expression(1/hat(sigma)[hat(delta)[i]]),
expression((hat(delta)[i]-hat(delta)[F])/hat(sigma)[hat(delta)[i]]))
}
################################################ 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]])
}
cond = list()
data.subgrp = list()
ss.subgrp = matrix(rep(0, n.subgrp.tol * n.subgrp.tol), nrow = n.subgrp.tol)
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]], ]
}
}
k = 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 for subgroups and the full
# population
if (sum((dat$trt == "1")) == 0 | sum((dat$trt == "0")) == 0){
treatment.mean.full = NA
treatment.std.full = NA
zscore.full = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = treatment.mean.full / treatment.std.full
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = treatment.mean.full / treatment.std.full
}else if (outcome.type == "survival"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coef[1, 1]
treatment.std.full = model.sum$coef[1, 3]
zscore.full = treatment.mean.full / treatment.std.full
}
}
treatment.mean = matrix(0, nrow = n.subgrp.tol, ncol = 1)
treatment.std = matrix(0, nrow = n.subgrp.tol, ncol = 1)
zscore = matrix(0, nrow = n.subgrp.tol, 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.std[i] = NA
zscore[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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
zscore[i] = treatment.mean[i] / treatment.std[i] # calculate the z-score of the subgroup effect size relative
# to the standardized error of the full population effect size
# estimator
}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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
zscore[i] = treatment.mean[i] / treatment.std[i] # calculate the z-score of the subgroup effect size relative
# to the standardized error of the full population effect size
# estimator
}else if (outcome.type == "survival"){
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.std[i] = model.sum$coef[1, 3]
zscore[i] = treatment.mean[i] / treatment.std[i]
}
}
}
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 = "")
}
}
dimnames(treatment.mean) = list(c(lab.subgrp), c("mean") )
dimnames(treatment.std) = list(c(lab.subgrp), c("std") )
dimnames(zscore) = list(c(lab.subgrp), c("zscore") )
zscore = (treatment.mean - treatment.mean.full)/treatment.std
################################################ 2. produce a graph #################################################################
if (is.null(range.v)){
y.lim.max = ceiling(max(zscore, na.rm = TRUE))
y.lim.min = floor(min(zscore, na.rm = TRUE))
if (y.lim.max <= 2) {y.lim.max = y.lim.max + 6}
if (y.lim.min >= -2) {y.lim.min = y.lim.min -3}
}else{
y.lim.max = range.v[2]
y.lim.min = range.v[1]
}
expand.fac = 2 # the expansion factor of the original unit on the x-axis
length.unit.xaxis = 5 # the number of the units on the x-axis
upper.bd = 2; lower.bd = -2
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 0.1 #ceiling(max(1/treatment.std))
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 1 #ceiling(max(1/treatment.std))
y.lab.upp.pos = treatment.mean.full + upper.bd
y.lab.low.pos = treatment.mean.full + lower.bd
pos.esize.mean = treatment.mean.full # the y-coordinate of the full
pos.esize.lower = pos.esize.mean + lower.bd
pos.esize.upper = pos.esize.mean + upper.bd
yy1 = seq(pos.esize.mean, y.lim.min + treatment.mean.full, -2)
yy2 = seq(pos.esize.mean, y.lim.max + treatment.mean.full, 2)
lab = c(rev(yy1), yy2)
lab = c(rep(NA, len = length(seq(y.lim.min, lower.bd - 1, 1))),
c(lower.bd, NA, 0, NA, upper.bd),
c(pos.esize.lower, NA, pos.esize.mean, NA, pos.esize.upper),
rep(NA, len = length(seq(upper.bd +1, y.lim.max, 1))))
x <- c(0, length.unit.xaxis * expand.fac, length.unit.xaxis * expand.fac, 0)
y <- c(pos.esize.lower, pos.esize.lower, pos.esize.upper, pos.esize.upper)
x = c(x, x[1])
y = c(y, y[1])
dat_rect =
data.frame(x = x,
y = y)
expand.fac.adj = expand.fac * length.unit.xaxis * 1/(x.lab.max.pos)
datplot =
data.frame(x = c(1/treatment.std * expand.fac.adj),
y = c(zscore),
labels = lab.subgrp2,
lab_legend = lab.subgrp)
them_subgrplots = theme_bw() +
theme(panel.grid = element_blank())
ggplot(datplot) +
geom_polygon(aes_string(x = "x", y = "y"),
data = dat_rect, fill = "gray", color = NA) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.mean, yend = pos.esize.mean,
size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.lower, yend = pos.esize.lower,
linetype = 2, size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.upper, yend = pos.esize.upper,
linetype = 2, size = 0.25) +
geom_point(aes_string(x = "x", y = "y", color = "lab_legend")) +
geom_point(aes_string(x = "x", y = "y")) +
geom_text_repel(aes_string(x =" x", y = "y", label = "labels"), size = font.size, seed = 987) +
coord_equal(ylim = c(y.lim.min, y.lim.max),
xlim = c(0,13.5),
expand = FALSE) +
labs(title = title, x = lab.xy[[1]], y = lab.xy[[2]]) +
them_subgrplots +
scale_y_continuous(breaks = seq(y.lab.low.pos - 1, y.lab.upp.pos + 1, 1)[c(2,4,6)],
labels = c(NA, lower.bd, NA, 0, NA, upper.bd, NA)[c(2,4,6)]) +
scale_x_continuous(breaks = seq(0, length.unit.xaxis * expand.fac, 1 * expand.fac),
labels = round(seq(0, x.lab.max.pos, len = length.unit.xaxis + 1), 2)) +
scale_color_manual(values = rep("white", nrow(datplot))) -> pp
### draw arcs for displaying effect sizes-------------------------------------
angle = vector()
zscore.diff = zscore - treatment.mean.full
len.pt = dim(zscore)[1]
for (i in 1: len.pt){
angle[i] = atan( zscore.diff[i]/ ((1/treatment.std[i] * expand.fac.adj))) # calculate the angle between two vectors
}
if (sum(angle>0) != length(angle)){ # when the zscores of subgroups are partially larger and less than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 = 0; y0 = treatment.mean.full
npts = 200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.2)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 = 0; y0 = treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.9)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis2 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis2) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = 'yend'),
data = dat_redticks2,
color = "red") -> pp
} else if (sum(angle>0) == length(angle)){ # when all zscores of subgroups are larger than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.2* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis) -> pp
# draw the cicular axis
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
}else if (sum(angle<0) == length(angle)){ # when all zscores of subgroups are less than the zscore of the full population
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.9* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks2,
color = "red") -> pp
}
pp + guides(color = guide_legend(keywidth = unit(0, "cm"),
title = "")) +
theme(legend.key = element_blank(),
legend.justification = "top")
}
#' Modified Radial (Galbraith) plot for subgroup effect size
#'
#' This function produces a modified Galbraith's radial plot using ggplot. It shows the
#' treatment effect size of subgroups defined by the categories
#' of covariates. The x-axis represents the reciprocal of the standard error
#' of the difference between the subgroup treatment effect estimates and the
#' overall treatment effect estimate.
#' The y-axis is the
#' standardized effect size difference (the difference between subgroup effect
#' the full population effect is divided by its standard error).
#' Points here are for subgroups.
#' The grey region indicates whether subgroup effects
#' are homogeneous to the full population effect or not. The two arcs on the
#' right side show subgroup treatment effects in the original
#' scale, where the red spots are the projection of points from the origin on
#' the left side. Note that the vertical range of display can
#' be changed by setting different values on the associated input argument.
#' In addition, the function uses log odd ratio and log hazard
#' ratio for displaying subgroup effect sizes in binary and survival data,
#' respectively.
#' A ggplot object is returned, where further modifications can be added using the theme() function.
#'
#' @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" or "survival".
#' @param range.v a vector specifying the vertical range of graphical display.
#' @param font.size a number specifying the font size for the subgroup labels near points
#' @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 ticks.length a number from 0 to 1 to specify the length of the red
#' ticks indicating the treatment effect in the radial axis. A value of 1 will draw a complete line from the (0,0) point to the axis.
#'
#' @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
#'
#' ggplot_radial2(dat,
#' covari.sel = c(4, 5, 6, 7),
#' trt.sel = 3,
#' resp.sel = c(1, 2),
#' outcome.type = "survival",
#' range.v = c(-7, 6),
#' font.size = 4)
#' @export
#' @import grid
#' @import ggplot2
#' @import graphics
#' @importFrom ggrepel geom_text_repel
ggplot_radial2 <- function(dat, covari.sel, trt.sel, resp.sel, outcome.type,
range.v = NULL,
font.size = 4,
title = NULL,
lab.xy = "default",
ticks.length = 0.05){
################################################ 0. argument validity check #################################################################
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!")
}else if (outcome.type == "survival"){
if (missing(resp.sel)) stop("The response variablehas not been specified!")
if (!(length(resp.sel) == 2)) stop("The response variable for analysing survival data should have two components!")
if (!(is.numeric(dat[, resp.sel[1]]))) stop("The response variable specifying survival time is not numeric!")
if (!(is.numeric(dat[, resp.sel[2]]) || is.logical(dat[, resp.sel[2]]) ) ) stop("The response variable specifying indicators of right censoring should be numerical or logical!")
if (length(names(table(dat[, resp.sel[2]]))) > 2) stop("The response variable specifying indicators of right censoring is not binary!")
if (sum(is.element(names(table(dat[, resp.sel[2]])), c("0","1"))) != 2) stop("The response variable specifying indicators of right censoring is not coded as 0 and 1!")
}
if (!(length(range.v) == 2)) stop("The vertical range of graphical display should have two components (specifying the minimum and maximum!)")
if (!(is.numeric(font.size))) stop("The font.size argument is not numeric!")
if (!(length(font.size) == 1)) stop("The font.size argument should be of length 1!")
if(lab.xy == "default") {
lab.x = expression(1/sqrt(hat(Var)(hat(delta)[i]-hat(delta)[F])))
lab.y = expression((hat(delta)[i]-hat(delta)[F])/sqrt(hat(Var)(hat(delta)[i]-hat(delta)[F])))
lab.xy = label.xy = list(lab.x, lab.y)
}
################################################ 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
}
n_t = sum(dat$trt == 1)
n_c = sum(dat$trt == 0)
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]])
}
cond = list()
data.subgrp = list()
lambda_t = lambda_c = lambda = deltaS = list()
ss.subgrp = matrix(rep(0, n.subgrp.tol * n.subgrp.tol), nrow = n.subgrp.tol)
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]], ]
lambda_t[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 1)/sum(dat$trt == 1)
lambda_c[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 0)/sum(dat$trt == 0)
lambda[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j])/nrow(dat)
if (outcome.type == "survival"){
which. = dat[, covari.sel[i]] == cats.var[[i]][j]
# dat[which.,]
lambda[[k]] = sum(dat[which., "status"] == 1)/sum(dat$status == 1)
}
if (outcome.type == "continuous"){
deltaS[[k]] = mean(dat$resp[which(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 1)]) -
mean(dat$resp[which(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 0)])
}
}
}
k = 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 for subgroups and the full
# population
if (sum((dat$trt == "1")) == 0 | sum((dat$trt == "0")) == 0){
treatment.mean.full = NA
treatment.std.full = NA
zscore.full = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = 0
sigma2 = sum((model.int$residuals^2))/model.int$df.residual
varDeltaF = (sigma2/sum(dat$trt == 1) + sigma2/sum(dat$trt == 0))
deltaF = mean(dat$resp[which(dat$trt == 1)]) -
mean(dat$resp[which(dat$trt == 0)])
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = 0
varDeltaF = treatment.std.full^2
deltaF = treatment.mean.full
}else if (outcome.type == "survival"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coef[1, 1]
treatment.std.full = model.sum$coef[1, 3]
zscore.full = 0
varDeltaF = treatment.std.full^2
deltaF = treatment.mean.full
}
}
treatment.mean = matrix(0, nrow = n.subgrp.tol, ncol = 1)
treatment.std = matrix(0, nrow = n.subgrp.tol, ncol = 1)
zscore = VarDsDf = sdDsDf = matrix(0, nrow = n.subgrp.tol, ncol = 1)
i=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.std[i] = NA
zscore[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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
# calculate the z-score of the subgroup effect size relative
# to the standardized error of the difference between
# subgroup treatment effect and population effect size
# estimator
VarDsDf[i] = ((1-lambda_t[[i]])/lambda_t[[i]]/n_t +
(1-lambda_c[[i]])/lambda_c[[i]]/n_c) * sigma2
zscore[i] = (deltaS[[i]] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
VarDsDf[i] = treatment.std[i]^2 + varDeltaF -
2 * sqrt(lambda[[i]])*treatment.std[i]*sqrt(varDeltaF)
zscore[i] = (treatment.mean[i] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}else if (outcome.type == "survival"){
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.std[i] = model.sum$coef[1, 3]
VarDsDf[i] = treatment.std[i]^2 + varDeltaF -
2 * sqrt(lambda[[i]])*treatment.std[i]*sqrt(varDeltaF)
zscore[i] = (treatment.mean[i] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}
}
}
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 = "")
}
}
dimnames(treatment.mean) = list(c(lab.subgrp), c("mean") )
dimnames(treatment.std) = list(c(lab.subgrp), c("std") )
dimnames(zscore) = list(c(lab.subgrp), c("zscore") )
dimnames(sdDsDf) = list(c(lab.subgrp), c("sdDsDf") )
treatment.std = sdDsDf
################################################ 2. produce a graph #################################################################
if (is.null(range.v)){
y.lim.max = ceiling(max(zscore, na.rm = TRUE))
y.lim.min = floor(min(zscore, na.rm = TRUE))
if (y.lim.max <= 2) {y.lim.max = y.lim.max + 6}
if (y.lim.min >= -2) {y.lim.min = y.lim.min -3}
}else{
y.lim.max = range.v[2]
y.lim.min = range.v[1]
}
expand.fac = 1.7 # the expansion factor of the original unit on the x-axis
length.unit.xaxis = 5 # the number of the units on the x-axis
upper.bd = 2; lower.bd = -2
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 0.1 #ceiling(max(1/treatment.std))
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 1 #ceiling(max(1/treatment.std))
y.lab.upp.pos = treatment.mean.full + upper.bd
y.lab.low.pos = treatment.mean.full + lower.bd
pos.esize.mean = treatment.mean.full # the y-coordinate of the full
pos.esize.lower = pos.esize.mean + lower.bd
pos.esize.upper = pos.esize.mean + upper.bd
yy1 = seq(pos.esize.mean, y.lim.min + treatment.mean.full, -2)
yy2 = seq(pos.esize.mean, y.lim.max + treatment.mean.full, 2)
lab = c(rev(yy1), yy2)
lab = c(rep(NA, len = length(seq(y.lim.min, lower.bd - 1, 1))),
c(lower.bd, NA, 0, NA, upper.bd),
c(pos.esize.lower, NA, pos.esize.mean, NA, pos.esize.upper),
rep(NA, len = length(seq(upper.bd +1, y.lim.max, 1))))
x <- c(0, length.unit.xaxis * expand.fac, length.unit.xaxis * expand.fac, 0)
y <- c(pos.esize.lower, pos.esize.lower, pos.esize.upper, pos.esize.upper)
x = c(x, x[1])
y = c(y, y[1])
dat_rect =
data.frame(x = x,
y = y)
expand.fac.adj = expand.fac * length.unit.xaxis * 1/(x.lab.max.pos)
datplot =
data.frame(x = c(1/treatment.std * expand.fac.adj),
y = c(zscore),
labels = lab.subgrp2,
lab_legend = lab.subgrp)
them_subgrplots = theme_bw() +
theme(panel.grid = element_blank())
ggplot(datplot) +
geom_polygon(aes_string(x = "x", y = "y"),
data = dat_rect, fill = "gray", color = NA) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.mean, yend = pos.esize.mean,
size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.lower, yend = pos.esize.lower,
linetype = 2, size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.upper, yend = pos.esize.upper,
linetype = 2, size = 0.25) +
geom_point(aes_string(x = "x", y = "y", color = "lab_legend")) +
geom_point(aes_string(x = "x", y = "y")) +
geom_text_repel(aes_string(x =" x", y = "y", label = "labels"), size = font.size, seed = 987) +
coord_equal(ylim = c(y.lim.min, y.lim.max),
xlim = c(0,13.5),
expand = FALSE) +
labs(title = title, x = lab.xy[[1]], y = lab.xy[[2]]) +
them_subgrplots +
scale_y_continuous(breaks = seq(y.lab.low.pos - 1, y.lab.upp.pos + 1, 1)[c(2,4,6)],
labels = c(NA, lower.bd, NA, 0, NA, upper.bd, NA)[c(2,4,6)]) +
scale_x_continuous(breaks = seq(0, length.unit.xaxis * expand.fac, 1 * expand.fac),
labels = round(seq(0, x.lab.max.pos, len = length.unit.xaxis + 1), 2)) +
scale_color_manual(values = rep("white", nrow(datplot))) -> pp
### draw arcs for displaying effect sizes-------------------------------------
angle = vector()
zscore.diff = zscore - treatment.mean.full
len.pt = dim(zscore)[1]
for (i in 1: len.pt){
angle[i] = atan( zscore.diff[i]/ ((1/treatment.std[i] * expand.fac.adj))) # calculate the angle between two vectors
}
if (sum(angle>0) != length(angle)){ # when the zscores of subgroups are partially larger and less than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 = 0; y0 = treatment.mean.full
npts = 200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.2)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 = 0; y0 = treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.9)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis2 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis2) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = 'yend'),
data = dat_redticks2,
color = "red") -> pp
} else if (sum(angle>0) == length(angle)){ # when all zscores of subgroups are larger than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.2* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis) -> pp
# draw the cicular axis
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
}else if (sum(angle<0) == length(angle)){ # when all zscores of subgroups are less than the zscore of the full population
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.9* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks2,
color = "red") -> pp
}
pp + guides(color = guide_legend(keywidth = unit(0, "cm"),
title = "")) +
theme(legend.key = element_blank(),
legend.justification = "top")
}
Try the SubgrPlots package in your browser
Any scripts or data that you put into this service are public.
SubgrPlots documentation built on Jan. 29, 2020, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ggplot_radial ggplot_radial2
Documented in ggplot_radial ggplot_radial2
#' Radial (Galbraith) plot for subgroup effect size
#'
#' This function produces a modified Galbraith's radial plot using ggplot. It shows the
#' treatment effect size of subgroups defined by the categories
#' of covariates. The x-axis represents the reciprocal of the standard error
#' of subgroup treatment effect estimates.The y-axis means
#' standardized effect size difference (the difference between subgroup effect
#' the full popultion effect is divided by the standard error
#' of the estimator for the overall population effect.
#' Points here are for subgroups.
#' The grey region indicates whether subgroup effects
#' are homogeneous to the full population effect or not. The two arcs on the
#' right side show subgroup treatment effects in the original
#' scale, where the red spots are the projection of points from the origin on
#' the left side. Note that the vertical range of display can
#' be changed by setting different values on the associated input argument.
#' In addition, the function uses log odd ratio and log hazard
#' ratio for displaying subgroup effect sizes in binary and survival data,
#' respectively.
#' A ggplot object is returned, where further modifications can be added using the theme() function.
#'
#' @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" or "survival".
#' @param range.v a vector specifying the vertical range of graphical display.
#' @param font.size a number specifying the font size for the subgroup labels near points
#' @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 ticks.length a number from 0 to 1 to specify the length of the red
#' ticks indicating the treatment effect in the radial axis. A value of 1 will draw a complete line from the (0,0) point to the axis.
#'
#' @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
#'
#' ggplot_radial(dat,
#' covari.sel = c(4, 5, 6, 7),
#' trt.sel = 3,
#' resp.sel = c(1, 2),
#' outcome.type = "survival",
#' range.v = c(-7, 6),
#' font.size = 4)
#' @export
#' @import ggplot2
#' @import grid
#' @import graphics
#' @importFrom ggrepel geom_text_repel
ggplot_radial <- function(dat, covari.sel, trt.sel, resp.sel, outcome.type,
range.v = NULL,
font.size = 4,
title = NULL,
lab.xy = "default",
ticks.length = 0.05){
################################################ 0. argument validity check #################################################################
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!")
}else if (outcome.type == "survival"){
if (missing(resp.sel)) stop("The response variablehas not been specified!")
if (!(length(resp.sel) == 2)) stop("The response variable for analysing survival data should have two components!")
if (!(is.numeric(dat[, resp.sel[1]]))) stop("The response variable specifying survival time is not numeric!")
if (!(is.numeric(dat[, resp.sel[2]]) || is.logical(dat[, resp.sel[2]]) ) ) stop("The response variable specifying indicators of right censoring should be numerical or logical!")
if (length(names(table(dat[, resp.sel[2]]))) > 2) stop("The response variable specifying indicators of right censoring is not binary!")
if (sum(is.element(names(table(dat[, resp.sel[2]])), c("0","1"))) != 2) stop("The response variable specifying indicators of right censoring is not coded as 0 and 1!")
}
if (!(length(range.v) == 2)) stop("The vertical range of graphical display should have two components (specifying the minimum and maximum!)")
if (!(is.numeric(font.size))) stop("The font.size argument is not numeric!")
if (!(length(font.size) == 1)) stop("The font.size argument should be of length 1!")
if(lab.xy == "default") {
lab.xy = list(expression(1/hat(sigma)[hat(delta)[i]]),
expression((hat(delta)[i]-hat(delta)[F])/hat(sigma)[hat(delta)[i]]))
}
################################################ 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]])
}
cond = list()
data.subgrp = list()
ss.subgrp = matrix(rep(0, n.subgrp.tol * n.subgrp.tol), nrow = n.subgrp.tol)
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]], ]
}
}
k = 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 for subgroups and the full
# population
if (sum((dat$trt == "1")) == 0 | sum((dat$trt == "0")) == 0){
treatment.mean.full = NA
treatment.std.full = NA
zscore.full = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = treatment.mean.full / treatment.std.full
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = treatment.mean.full / treatment.std.full
}else if (outcome.type == "survival"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coef[1, 1]
treatment.std.full = model.sum$coef[1, 3]
zscore.full = treatment.mean.full / treatment.std.full
}
}
treatment.mean = matrix(0, nrow = n.subgrp.tol, ncol = 1)
treatment.std = matrix(0, nrow = n.subgrp.tol, ncol = 1)
zscore = matrix(0, nrow = n.subgrp.tol, 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.std[i] = NA
zscore[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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
zscore[i] = treatment.mean[i] / treatment.std[i] # calculate the z-score of the subgroup effect size relative
# to the standardized error of the full population effect size
# estimator
}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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
zscore[i] = treatment.mean[i] / treatment.std[i] # calculate the z-score of the subgroup effect size relative
# to the standardized error of the full population effect size
# estimator
}else if (outcome.type == "survival"){
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.std[i] = model.sum$coef[1, 3]
zscore[i] = treatment.mean[i] / treatment.std[i]
}
}
}
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 = "")
}
}
dimnames(treatment.mean) = list(c(lab.subgrp), c("mean") )
dimnames(treatment.std) = list(c(lab.subgrp), c("std") )
dimnames(zscore) = list(c(lab.subgrp), c("zscore") )
zscore = (treatment.mean - treatment.mean.full)/treatment.std
################################################ 2. produce a graph #################################################################
if (is.null(range.v)){
y.lim.max = ceiling(max(zscore, na.rm = TRUE))
y.lim.min = floor(min(zscore, na.rm = TRUE))
if (y.lim.max <= 2) {y.lim.max = y.lim.max + 6}
if (y.lim.min >= -2) {y.lim.min = y.lim.min -3}
}else{
y.lim.max = range.v[2]
y.lim.min = range.v[1]
}
expand.fac = 2 # the expansion factor of the original unit on the x-axis
length.unit.xaxis = 5 # the number of the units on the x-axis
upper.bd = 2; lower.bd = -2
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 0.1 #ceiling(max(1/treatment.std))
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 1 #ceiling(max(1/treatment.std))
y.lab.upp.pos = treatment.mean.full + upper.bd
y.lab.low.pos = treatment.mean.full + lower.bd
pos.esize.mean = treatment.mean.full # the y-coordinate of the full
pos.esize.lower = pos.esize.mean + lower.bd
pos.esize.upper = pos.esize.mean + upper.bd
yy1 = seq(pos.esize.mean, y.lim.min + treatment.mean.full, -2)
yy2 = seq(pos.esize.mean, y.lim.max + treatment.mean.full, 2)
lab = c(rev(yy1), yy2)
lab = c(rep(NA, len = length(seq(y.lim.min, lower.bd - 1, 1))),
c(lower.bd, NA, 0, NA, upper.bd),
c(pos.esize.lower, NA, pos.esize.mean, NA, pos.esize.upper),
rep(NA, len = length(seq(upper.bd +1, y.lim.max, 1))))
x <- c(0, length.unit.xaxis * expand.fac, length.unit.xaxis * expand.fac, 0)
y <- c(pos.esize.lower, pos.esize.lower, pos.esize.upper, pos.esize.upper)
x = c(x, x[1])
y = c(y, y[1])
dat_rect =
data.frame(x = x,
y = y)
expand.fac.adj = expand.fac * length.unit.xaxis * 1/(x.lab.max.pos)
datplot =
data.frame(x = c(1/treatment.std * expand.fac.adj),
y = c(zscore),
labels = lab.subgrp2,
lab_legend = lab.subgrp)
them_subgrplots = theme_bw() +
theme(panel.grid = element_blank())
ggplot(datplot) +
geom_polygon(aes_string(x = "x", y = "y"),
data = dat_rect, fill = "gray", color = NA) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.mean, yend = pos.esize.mean,
size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.lower, yend = pos.esize.lower,
linetype = 2, size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.upper, yend = pos.esize.upper,
linetype = 2, size = 0.25) +
geom_point(aes_string(x = "x", y = "y", color = "lab_legend")) +
geom_point(aes_string(x = "x", y = "y")) +
geom_text_repel(aes_string(x =" x", y = "y", label = "labels"), size = font.size, seed = 987) +
coord_equal(ylim = c(y.lim.min, y.lim.max),
xlim = c(0,13.5),
expand = FALSE) +
labs(title = title, x = lab.xy[[1]], y = lab.xy[[2]]) +
them_subgrplots +
scale_y_continuous(breaks = seq(y.lab.low.pos - 1, y.lab.upp.pos + 1, 1)[c(2,4,6)],
labels = c(NA, lower.bd, NA, 0, NA, upper.bd, NA)[c(2,4,6)]) +
scale_x_continuous(breaks = seq(0, length.unit.xaxis * expand.fac, 1 * expand.fac),
labels = round(seq(0, x.lab.max.pos, len = length.unit.xaxis + 1), 2)) +
scale_color_manual(values = rep("white", nrow(datplot))) -> pp
### draw arcs for displaying effect sizes-------------------------------------
angle = vector()
zscore.diff = zscore - treatment.mean.full
len.pt = dim(zscore)[1]
for (i in 1: len.pt){
angle[i] = atan( zscore.diff[i]/ ((1/treatment.std[i] * expand.fac.adj))) # calculate the angle between two vectors
}
if (sum(angle>0) != length(angle)){ # when the zscores of subgroups are partially larger and less than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 = 0; y0 = treatment.mean.full
npts = 200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.2)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 = 0; y0 = treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.9)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis2 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis2) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = 'yend'),
data = dat_redticks2,
color = "red") -> pp
} else if (sum(angle>0) == length(angle)){ # when all zscores of subgroups are larger than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.2* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis) -> pp
# draw the cicular axis
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
}else if (sum(angle<0) == length(angle)){ # when all zscores of subgroups are less than the zscore of the full population
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.9* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks2,
color = "red") -> pp
}
pp + guides(color = guide_legend(keywidth = unit(0, "cm"),
title = "")) +
theme(legend.key = element_blank(),
legend.justification = "top")
}
#' Modified Radial (Galbraith) plot for subgroup effect size
#'
#' This function produces a modified Galbraith's radial plot using ggplot. It shows the
#' treatment effect size of subgroups defined by the categories
#' of covariates. The x-axis represents the reciprocal of the standard error
#' of the difference between the subgroup treatment effect estimates and the
#' overall treatment effect estimate.
#' The y-axis is the
#' standardized effect size difference (the difference between subgroup effect
#' the full population effect is divided by its standard error).
#' Points here are for subgroups.
#' The grey region indicates whether subgroup effects
#' are homogeneous to the full population effect or not. The two arcs on the
#' right side show subgroup treatment effects in the original
#' scale, where the red spots are the projection of points from the origin on
#' the left side. Note that the vertical range of display can
#' be changed by setting different values on the associated input argument.
#' In addition, the function uses log odd ratio and log hazard
#' ratio for displaying subgroup effect sizes in binary and survival data,
#' respectively.
#' A ggplot object is returned, where further modifications can be added using the theme() function.
#'
#' @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" or "survival".
#' @param range.v a vector specifying the vertical range of graphical display.
#' @param font.size a number specifying the font size for the subgroup labels near points
#' @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 ticks.length a number from 0 to 1 to specify the length of the red
#' ticks indicating the treatment effect in the radial axis. A value of 1 will draw a complete line from the (0,0) point to the axis.
#'
#' @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
#'
#' ggplot_radial2(dat,
#' covari.sel = c(4, 5, 6, 7),
#' trt.sel = 3,
#' resp.sel = c(1, 2),
#' outcome.type = "survival",
#' range.v = c(-7, 6),
#' font.size = 4)
#' @export
#' @import grid
#' @import ggplot2
#' @import graphics
#' @importFrom ggrepel geom_text_repel
ggplot_radial2 <- function(dat, covari.sel, trt.sel, resp.sel, outcome.type,
range.v = NULL,
font.size = 4,
title = NULL,
lab.xy = "default",
ticks.length = 0.05){
################################################ 0. argument validity check #################################################################
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!")
}else if (outcome.type == "survival"){
if (missing(resp.sel)) stop("The response variablehas not been specified!")
if (!(length(resp.sel) == 2)) stop("The response variable for analysing survival data should have two components!")
if (!(is.numeric(dat[, resp.sel[1]]))) stop("The response variable specifying survival time is not numeric!")
if (!(is.numeric(dat[, resp.sel[2]]) || is.logical(dat[, resp.sel[2]]) ) ) stop("The response variable specifying indicators of right censoring should be numerical or logical!")
if (length(names(table(dat[, resp.sel[2]]))) > 2) stop("The response variable specifying indicators of right censoring is not binary!")
if (sum(is.element(names(table(dat[, resp.sel[2]])), c("0","1"))) != 2) stop("The response variable specifying indicators of right censoring is not coded as 0 and 1!")
}
if (!(length(range.v) == 2)) stop("The vertical range of graphical display should have two components (specifying the minimum and maximum!)")
if (!(is.numeric(font.size))) stop("The font.size argument is not numeric!")
if (!(length(font.size) == 1)) stop("The font.size argument should be of length 1!")
if(lab.xy == "default") {
lab.x = expression(1/sqrt(hat(Var)(hat(delta)[i]-hat(delta)[F])))
lab.y = expression((hat(delta)[i]-hat(delta)[F])/sqrt(hat(Var)(hat(delta)[i]-hat(delta)[F])))
lab.xy = label.xy = list(lab.x, lab.y)
}
################################################ 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
}
n_t = sum(dat$trt == 1)
n_c = sum(dat$trt == 0)
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]])
}
cond = list()
data.subgrp = list()
lambda_t = lambda_c = lambda = deltaS = list()
ss.subgrp = matrix(rep(0, n.subgrp.tol * n.subgrp.tol), nrow = n.subgrp.tol)
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]], ]
lambda_t[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 1)/sum(dat$trt == 1)
lambda_c[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 0)/sum(dat$trt == 0)
lambda[[k]] = sum(dat[, covari.sel[i]] == cats.var[[i]][j])/nrow(dat)
if (outcome.type == "survival"){
which. = dat[, covari.sel[i]] == cats.var[[i]][j]
# dat[which.,]
lambda[[k]] = sum(dat[which., "status"] == 1)/sum(dat$status == 1)
}
if (outcome.type == "continuous"){
deltaS[[k]] = mean(dat$resp[which(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 1)]) -
mean(dat$resp[which(dat[, covari.sel[i]] == cats.var[[i]][j] & dat$trt == 0)])
}
}
}
k = 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 for subgroups and the full
# population
if (sum((dat$trt == "1")) == 0 | sum((dat$trt == "0")) == 0){
treatment.mean.full = NA
treatment.std.full = NA
zscore.full = NA
}else{
if (outcome.type == "continuous"){
model.int = lm(resp ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = 0
sigma2 = sum((model.int$residuals^2))/model.int$df.residual
varDeltaF = (sigma2/sum(dat$trt == 1) + sigma2/sum(dat$trt == 0))
deltaF = mean(dat$resp[which(dat$trt == 1)]) -
mean(dat$resp[which(dat$trt == 0)])
}else if (outcome.type == "binary"){
model.int = glm(resp ~ trt, family = "binomial", data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coefficients[2, 1]
treatment.std.full = model.sum$coefficients[2, 2]
zscore.full = 0
varDeltaF = treatment.std.full^2
deltaF = treatment.mean.full
}else if (outcome.type == "survival"){
model.int = survival::coxph(survival::Surv(time, status) ~ trt, data = dat)
model.sum = summary(model.int)
treatment.mean.full = model.sum$coef[1, 1]
treatment.std.full = model.sum$coef[1, 3]
zscore.full = 0
varDeltaF = treatment.std.full^2
deltaF = treatment.mean.full
}
}
treatment.mean = matrix(0, nrow = n.subgrp.tol, ncol = 1)
treatment.std = matrix(0, nrow = n.subgrp.tol, ncol = 1)
zscore = VarDsDf = sdDsDf = matrix(0, nrow = n.subgrp.tol, ncol = 1)
i=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.std[i] = NA
zscore[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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
# calculate the z-score of the subgroup effect size relative
# to the standardized error of the difference between
# subgroup treatment effect and population effect size
# estimator
VarDsDf[i] = ((1-lambda_t[[i]])/lambda_t[[i]]/n_t +
(1-lambda_c[[i]])/lambda_c[[i]]/n_c) * sigma2
zscore[i] = (deltaS[[i]] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}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] # record subgroup effect size
treatment.std[i] = model.sum$coefficients[2, 2]
VarDsDf[i] = treatment.std[i]^2 + varDeltaF -
2 * sqrt(lambda[[i]])*treatment.std[i]*sqrt(varDeltaF)
zscore[i] = (treatment.mean[i] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}else if (outcome.type == "survival"){
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.std[i] = model.sum$coef[1, 3]
VarDsDf[i] = treatment.std[i]^2 + varDeltaF -
2 * sqrt(lambda[[i]])*treatment.std[i]*sqrt(varDeltaF)
zscore[i] = (treatment.mean[i] - deltaF) / sqrt(VarDsDf[i])
sdDsDf[i] = sqrt(VarDsDf[i])
}
}
}
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 = "")
}
}
dimnames(treatment.mean) = list(c(lab.subgrp), c("mean") )
dimnames(treatment.std) = list(c(lab.subgrp), c("std") )
dimnames(zscore) = list(c(lab.subgrp), c("zscore") )
dimnames(sdDsDf) = list(c(lab.subgrp), c("sdDsDf") )
treatment.std = sdDsDf
################################################ 2. produce a graph #################################################################
if (is.null(range.v)){
y.lim.max = ceiling(max(zscore, na.rm = TRUE))
y.lim.min = floor(min(zscore, na.rm = TRUE))
if (y.lim.max <= 2) {y.lim.max = y.lim.max + 6}
if (y.lim.min >= -2) {y.lim.min = y.lim.min -3}
}else{
y.lim.max = range.v[2]
y.lim.min = range.v[1]
}
expand.fac = 1.7 # the expansion factor of the original unit on the x-axis
length.unit.xaxis = 5 # the number of the units on the x-axis
upper.bd = 2; lower.bd = -2
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 0.1 #ceiling(max(1/treatment.std))
x.lab.max.pos = round(max(1/c(treatment.std, treatment.std.full))) + 1 #ceiling(max(1/treatment.std))
y.lab.upp.pos = treatment.mean.full + upper.bd
y.lab.low.pos = treatment.mean.full + lower.bd
pos.esize.mean = treatment.mean.full # the y-coordinate of the full
pos.esize.lower = pos.esize.mean + lower.bd
pos.esize.upper = pos.esize.mean + upper.bd
yy1 = seq(pos.esize.mean, y.lim.min + treatment.mean.full, -2)
yy2 = seq(pos.esize.mean, y.lim.max + treatment.mean.full, 2)
lab = c(rev(yy1), yy2)
lab = c(rep(NA, len = length(seq(y.lim.min, lower.bd - 1, 1))),
c(lower.bd, NA, 0, NA, upper.bd),
c(pos.esize.lower, NA, pos.esize.mean, NA, pos.esize.upper),
rep(NA, len = length(seq(upper.bd +1, y.lim.max, 1))))
x <- c(0, length.unit.xaxis * expand.fac, length.unit.xaxis * expand.fac, 0)
y <- c(pos.esize.lower, pos.esize.lower, pos.esize.upper, pos.esize.upper)
x = c(x, x[1])
y = c(y, y[1])
dat_rect =
data.frame(x = x,
y = y)
expand.fac.adj = expand.fac * length.unit.xaxis * 1/(x.lab.max.pos)
datplot =
data.frame(x = c(1/treatment.std * expand.fac.adj),
y = c(zscore),
labels = lab.subgrp2,
lab_legend = lab.subgrp)
them_subgrplots = theme_bw() +
theme(panel.grid = element_blank())
ggplot(datplot) +
geom_polygon(aes_string(x = "x", y = "y"),
data = dat_rect, fill = "gray", color = NA) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.mean, yend = pos.esize.mean,
size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.lower, yend = pos.esize.lower,
linetype = 2, size = 0.25) +
geom_segment(x = 0, xend = length.unit.xaxis * expand.fac,
y = pos.esize.upper, yend = pos.esize.upper,
linetype = 2, size = 0.25) +
geom_point(aes_string(x = "x", y = "y", color = "lab_legend")) +
geom_point(aes_string(x = "x", y = "y")) +
geom_text_repel(aes_string(x =" x", y = "y", label = "labels"), size = font.size, seed = 987) +
coord_equal(ylim = c(y.lim.min, y.lim.max),
xlim = c(0,13.5),
expand = FALSE) +
labs(title = title, x = lab.xy[[1]], y = lab.xy[[2]]) +
them_subgrplots +
scale_y_continuous(breaks = seq(y.lab.low.pos - 1, y.lab.upp.pos + 1, 1)[c(2,4,6)],
labels = c(NA, lower.bd, NA, 0, NA, upper.bd, NA)[c(2,4,6)]) +
scale_x_continuous(breaks = seq(0, length.unit.xaxis * expand.fac, 1 * expand.fac),
labels = round(seq(0, x.lab.max.pos, len = length.unit.xaxis + 1), 2)) +
scale_color_manual(values = rep("white", nrow(datplot))) -> pp
### draw arcs for displaying effect sizes-------------------------------------
angle = vector()
zscore.diff = zscore - treatment.mean.full
len.pt = dim(zscore)[1]
for (i in 1: len.pt){
angle[i] = atan( zscore.diff[i]/ ((1/treatment.std[i] * expand.fac.adj))) # calculate the angle between two vectors
}
if (sum(angle>0) != length(angle)){ # when the zscores of subgroups are partially larger and less than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 = 0; y0 = treatment.mean.full
npts = 200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.2)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 = 0; y0 = treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + ((length.unit.xaxis + 0.9)* expand.fac)^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis2 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis2) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = 'yend'),
data = dat_redticks2,
color = "red") -> pp
} else if (sum(angle>0) == length(angle)){ # when all zscores of subgroups are larger than the zscore of the full population
axis.circle1.max = max(angle) # the difference between the maximal effect size and the nearest effect size unit
axis.circle1.min = 0 # the effect size for the full population
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.2* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle1.min, axis.circle1.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis) -> pp
# draw the cicular axis
r1 = r + 0.1
theta1 <- seq(axis.circle1.min, axis.circle1.max, len = 7)
esize.unit <- seq(treatment.mean.full, max(treatment.mean), len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks1,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks1[-1, ]) -> pp
r1 = r - 0.15
angle.pos = angle[angle>0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.pos)
x1 = r1 * cos(angle.pos)
y0 = r * sin(angle.pos)
y1 = r1 * sin(angle.pos)
dat_redticks1 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks1,
color = "red") -> pp
}else if (sum(angle<0) == length(angle)){ # when all zscores of subgroups are less than the zscore of the full population
axis.circle2.max = 0 # the effect size for the full population
axis.circle2.min = min(angle) # the difference between the minmal effect size and the nearest effect size unit
x0 =0; y0= treatment.mean.full
npts =200
r = ceiling(sqrt(pos.esize.mean^2 + (5.9* expand.fac) ^2)) - 0.1 # the radius of the outer circle
theta <- seq(axis.circle2.min, axis.circle2.max , len = 200)
xh <- r * cos(theta)
yh <- r * sin(theta)
x <- x0 + xh
y <- y0 + yh
dat_axis1 = data.frame(x = x,
y = y)
pp +
geom_path(aes_string(x = "x", y = "y"), data = dat_axis1) -> pp
r1 = r + 0.1
theta1 <- seq(axis.circle2.min, axis.circle2.max, len = 7)
esize.unit <- seq(min(treatment.mean), treatment.mean.full, len = length(theta1))
esize.unit = round(esize.unit,2)
r1 = r * 1.02
x0 = r * cos(theta1)
x1 = r1 * cos(theta1)
y0 = r * sin(theta1)
y1 = r1 * sin(theta1)
dat_ticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full,
label = esize.unit)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = 'xend', yend = "yend"),
data = dat_ticks2,
color = "black") +
geom_text(aes_string(x = "xend", y = "yend", label = "label"),
hjust = 0, nudge_x = .1,
data = dat_ticks2) -> pp
angle.neg = angle[angle<0]
r1 = r * (1 - ticks.length)
x0 = r * cos(angle.neg)
x1 = r1 * cos(angle.neg)
y0 = r * sin(angle.neg)
y1 = r1 * sin(angle.neg)
dat_redticks2 = data.frame(x = x0,
y = y0 + treatment.mean.full,
xend = x1,
yend = y1 + treatment.mean.full)
pp +
geom_segment(aes_string(x = "x", y = "y", xend = "xend", yend = "yend"),
data = dat_redticks2,
color = "red") -> pp
}
pp + guides(color = guide_legend(keywidth = unit(0, "cm"),
title = "")) +
theme(legend.key = element_blank(),
legend.justification = "top")
}
Try the SubgrPlots package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.