#' @title Dose proportionality test (power model)
#' @description log(y)~log(y) lm fit to test for dose proportionality
#' @param df data frame
#' @param x x variable (e.g. dose) passed in untransformed
#' @param y y parameter (e.g. auc) passed in untransformed
#' @param ci.level confidence level for lm fit in plot, Default: 0.95
#' @param plot.par add table with fit parameters, Default: F
#' @param max_position upper y limit passed to coord_carteesian (relative max(y)), Default: 1.1
#' @param min_tab_position lower y-value for parameter table (relative max(y)), Default: 1.05
#' @param min_text_position position of text string (relative max(y)), Default: 1
#' @param signif number of significant numbers in parameter table, Default: 3
#' @return list with plot, power lm summary and power lm fit
#' @rdname dose_proportionality_power
#' @export
#' @importFrom gridExtra tableGrob
#' @importFrom rlang enexpr expr_text sym eval_tidy quo
#' @importFrom tibble rownames_to_column
dose_proportionality_power <-
function(df, x, y, ci.level=0.95,
max_position = 1.2,
min_tab_position = 1.05,
min_text_position = 0.90,
signif=3){
x <- rlang::enexpr(x)
y <- rlang::enexpr(y)
# Check that x and y is in dataset
if(!rlang::expr_text(x) %in% names(df)){
stop(paste0(rlang::expr_text(x), " not found in dataset"))
}
if(!rlang::expr_text(y) %in% names(df)){
stop(paste0(rlang::expr_text(y), " not found in dataset"))
}
log_x_name <- paste0("log_", rlang::expr_text(x))
log_y_name <- paste0("log_", rlang::expr_text(y))
log_x <- rlang::sym(log_x_name)
log_y <- rlang::sym(log_y_name)
# 1. Add log of x and y
newdf <- rlang::eval_tidy(
rlang::quo(df %>%
mutate(!!log_x_name := log(!!x),
!!log_y_name := log(!!y))
)
)
# Fit power model lm
lm_power_fit <- rlang::eval_tidy(
rlang::quo(
lm(formula=!!log_y~!!log_x, data=newdf)
)
)
lm_power_sum <- summary(lm_power_fit)
# find max and min values of x and y, log(x) and log(y)
x_max_log <- max(newdf[[log_x_name]], na.rm=T)
x_min_log <- min(newdf[[log_x_name]], na.rm=T)
y_max <- max(newdf[[rlang::expr_text(y)]], na.rm=T)
y_min <- min(newdf[[rlang::expr_text(y)]], na.rm=T)
x_max <- max(newdf[[rlang::expr_text(x)]], na.rm=T)
x_min <- min(newdf[[rlang::expr_text(x)]], na.rm=T)
# predict values and back-transform
x_pred <- data.frame(x=seq(x_min_log, x_max_log, length=100)) # improve. by as a par?
names(x_pred) <- log_x_name
predicted <- predict.lm(lm_power_fit, x_pred,
se.fit=T,
interval=c("confidence"),
level=ci.level)
# put in data frame and back-tranform
pred_data <- data.frame(x = exp(x_pred[,1]),
y = exp(predicted$fit[,'fit']),
lwr = exp(predicted$fit[,'lwr']),
upr = exp(predicted$fit[,'upr']))
# calculate critical region
x_ratio <- x_max / x_min
critial_region <- format(
c( 1 + ( log(0.80)/log(x_ratio) ),
1 + ( log(1.25)/log(x_ratio) ) ),
digits=signif)
critial_region_string <-
paste0("Critical Region = [",
critial_region[1],",",
critial_region[2],"]")
# Extract the data from the summary to a data.frame for plotting
# 1. lm parameter data
lm_power_data <- as.data.frame(lm_power_sum$coefficients)
lm_power_data <- format(lm_power_data, digits=signif)
lm_power_data <- tibble::rownames_to_column(lm_power_data)
names(lm_power_data)[1] <- "Parameter"
# Rename parameters to intercept and slope
lm_power_data$Parameter <-
str_replace(lm_power_data$Parameter, "\\(Intercept\\)", "Intercept")
lm_power_data$Parameter <-
str_replace(lm_power_data$Parameter, log_x_name, "Power")
power_est <- as.numeric(lm_power_data$Estimate[lm_power_data$Parameter=="Power"])
power_sd <- as.numeric(lm_power_data$`Std. Error`[lm_power_data$Parameter=="Power"])
lm_power_string <- paste0("Power [95% CI] = ",
power_est," [",
signif(power_est - 1.96*power_sd, digits=signif), "-",
signif(power_est + 1.96*power_sd, digits=signif), "]")
# plot settings
y_max_pos <- y_max*max_position
y_text <- y_max*min_text_position
y_min_tab <- y_max*min_tab_position
x_diff <- (x_max - x_min) / 4
x_min_pos <- x_min + x_diff
x_max_pos <- x_max - x_diff
# plot.par=F, # improve, output both plots so that one can choose at a later stage
# plots without parameters but w critical region
p_woparwCR <- rlang::quo(
ggplot(newdf, aes(x=!!x, y=!!y)) +
geom_point() +
# add predicted line (defaults from geom_smooth defaults)
geom_line(data = pred_data, aes(x=x, y=y), inherit.aes=F,
color="#3366FF", size=1) +
# add se (defaults from geom_smooth defaults)
geom_ribbon(data=pred_data, aes(x=x, ymin=lwr,ymax=upr),
fill="grey60", alpha="0.4", inherit.aes=F) +
coord_cartesian(ylim = c(y_min, !!y_max_pos)) +
labs(title="Power model") +
# give CI of power estimate and critical region
annotate(geom="text", x=!!x_min_pos, y=!!y_text,
label=paste(lm_power_string, critial_region_string, sep="\n"))
)
# plots without parameters or critical region
p_woparwoCR <- rlang::quo(
ggplot(newdf, aes(x=!!x, y=!!y)) +
geom_point() +
# add predicted line (defaults from geom_smooth defaults)
geom_line(data = pred_data, aes(x=x, y=y), inherit.aes=F,
color="#3366FF", size=1) +
# add se (defaults from geom_smooth defaults)
geom_ribbon(data=pred_data, aes(x=x, ymin=lwr,ymax=upr),
fill="grey60", alpha="0.4", inherit.aes=F) +
coord_cartesian(ylim = c(y_min, !!y_max_pos)) +
labs(title="Power model") +
# give CI of power estimate and critical region
annotate(geom="text", x=!!x_min_pos, y=!!y_text,
label=lm_power_string)
)
# plots with parameters
p_wparwCR <- rlang::quo(
ggplot(newdf, aes(x=!!x, y=!!y)) +
# add lm parameters (first layer to make sure not overplotting outliers)
annotation_custom(gridExtra::tableGrob(lm_power_data, rows=NULL),
xmin=!!x_min_pos, xmax=!!x_max_pos,
ymin=!!y_min_tab, ymax=!!y_max_pos) +
geom_point() +
# add predicted line (defaults from geom_smooth defaults)
geom_line(data = pred_data, aes(x=x, y=y), inherit.aes=F,
color="#3366FF", size=1) +
# add se (defaults from geom_smooth defaults)
geom_ribbon(data=pred_data, aes(x=x, ymin=lwr,ymax=upr),
fill="grey60", alpha="0.4", inherit.aes=F) +
coord_cartesian(ylim = c(y_min, !!y_max_pos)) +
labs(title="Power model") +
# give critical region
annotate(geom="text", x=!!x_min_pos, y=!!y_text,
label=critial_region_string)
)
# plots with parameters
p_wparwoCR <- rlang::quo(
ggplot(newdf, aes(x=!!x, y=!!y)) +
# add lm parameters (first layer to make sure not overplotting outliers)
annotation_custom(gridExtra::tableGrob(lm_power_data, rows=NULL),
xmin=!!x_min_pos, xmax=!!x_max_pos,
ymin=!!y_min_tab, ymax=!!y_max_pos) +
geom_point() +
# add predicted line (defaults from geom_smooth defaults)
geom_line(data = pred_data, aes(x=x, y=y), inherit.aes=F,
color="#3366FF", size=1) +
# add se (defaults from geom_smooth defaults)
geom_ribbon(data=pred_data, aes(x=x, ymin=lwr,ymax=upr),
fill="grey60", alpha="0.4", inherit.aes=F) +
coord_cartesian(ylim = c(y_min, !!y_max_pos)) +
labs(title="Power model")
)
return(list(summary=summary(lm_power_fit),
fit=lm_power_fit,
plot = rlang::eval_tidy(p_woparwoCR),
plot_cr = rlang::eval_tidy(p_woparwCR),
plot_par = rlang::eval_tidy(p_wparwoCR),
plot_par_cr = rlang::eval_tidy(p_wparwCR))
)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.