# NOT RUN {
#======================================================================================
# Create a p-value function for an estimate using the normal distribution
#======================================================================================
res <- conf_dist(
estimate = c(-0.13)
, stderr = c(0.224494)
, type = "general_z"
, plot_type = "p_val"
, n_values = 1e4L
, est_names = c("Parameter value")
, log_yaxis = FALSE
, cut_logyaxis = 0.05
, conf_level = c(0.95)
, null_values = c(0)
, trans = "identity"
, alternative = "two_sided"
, xlab = "Var"
, xlim = c(-1, 1)
, together = TRUE
, plot_p_limit = 1 - 0.9999
, plot_counternull = TRUE
, title = NULL
, ylab = NULL
, ylab_sec = NULL
, inverted = FALSE
, x_scale = "default"
, plot = TRUE
)
#======================================================================================
# P-value function for a single regression coefficient (Agriculture in the model below)
#======================================================================================
mod <- lm(Infant.Mortality~Agriculture + Fertility + Examination, data = swiss)
summary(mod)
res <- conf_dist(
estimate = c(-0.02143)
, df = c(43)
, stderr = (0.02394)
, type = "linreg"
, plot_type = "p_val"
, n_values = 1e4L
, conf_level = c(0.95, 0.90, 0.80)
, null_values = c(0)
, trans = "identity"
, alternative = "two_sided"
, log_yaxis = TRUE
, cut_logyaxis = 0.05
, xlab = "Coefficient Agriculture"
, together = FALSE
, plot_p_limit = 1 - 0.999
, plot_counternull = FALSE
, title = NULL
, ylab = NULL
, ylab_sec = NULL
, inverted = FALSE
, x_scale = "default"
, plot = TRUE
)
#=======================================================================================
# P-value function for an odds ratio (logistic regression), plotted with inverted y-axis
#=======================================================================================
res <- conf_dist(
estimate = c(0.804037549)
, stderr = c(0.331819298)
, type = "logreg"
, plot_type = "p_val"
, n_values = 1e4L
, est_names = c("GPA")
, conf_level = c(0.95, 0.90, 0.80)
, null_values = c(log(1)) # null value on the log-odds scale
, trans = "exp"
, alternative = "two_sided"
, log_yaxis = FALSE
, cut_logyaxis = 0.05
, xlab = "Odds Ratio (GPA)"
, xlim = log(c(0.7, 5.2)) # axis limits on the log-odds scale
, together = FALSE
, plot_p_limit = 1 - 0.999
, plot_counternull = TRUE
, title = NULL
, ylab = NULL
, ylab_sec = NULL
, inverted = TRUE
, x_scale = "default"
, plot = TRUE
)
#======================================================================================
# Difference between two independent proportions: Newcombe with continuity correction
#======================================================================================
res <- conf_dist(
estimate = c(68/100, 98/150)
, n = c(100, 150)
, type = "propdiff"
, plot_type = "p_val"
, n_values = 1e4L
, conf_level = c(0.95, 0.90, 0.80)
, null_values = c(0)
, trans = "identity"
, alternative = "two_sided"
, log_yaxis = FALSE
, cut_logyaxis = 0.05
, xlab = "Difference between proportions"
, together = FALSE
, col = "#A52A2A" # Color curve in auburn
, plot_p_limit = 1 - 0.9999
, plot_counternull = FALSE
, title = NULL
, ylab = NULL
, ylab_sec = NULL
, inverted = FALSE
, x_scale = "default"
, plot = TRUE
)
#======================================================================================
# Difference between two independent proportions: Agresti & Caffo
#======================================================================================
# First proportion
x1 <- 8
n1 <- 40
# Second proportion
x2 <- 11
n2 <- 30
# Apply the correction
p1hat <- (x1 + 1)/(n1 + 2)
p2hat <- (x2 + 1)/(n2 + 2)
# The original estimator
est0 <- (x1/n1) - (x2/n2)
# The unmodified estimator and its standard error using the correction
est <- p1hat - p2hat
se <- sqrt(((p1hat*(1 - p1hat))/(n1 + 2)) + ((p2hat*(1 - p2hat))/(n2 + 2)))
res <- conf_dist(
estimate = c(est)
, stderr = c(se)
, type = "general_z"
, plot_type = "p_val"
, n_values = 1e4L
, log_yaxis = FALSE
, cut_logyaxis = 0.05
, conf_level = c(0.95, 0.99)
, null_values = c(0, 0.3)
, trans = "identity"
, alternative = "two_sided"
, xlab = "Difference of proportions"
, together = FALSE
, plot_p_limit = 1 - 0.9999
, plot_counternull = FALSE
, title = "P-value function for the difference of two independent proportions"
, ylab = NULL
, ylab_sec = NULL
, inverted = FALSE
, x_scale = "default"
, plot = TRUE
)
#========================================================================================
# P-value function and confidence distribution for the relative survival effect (1 - HR%)
# Replicating Figure 1 in Bender et al. (2005)
#========================================================================================
# Define the transformation function and its inverse for the relative survival effect
rse_fun <- function(x){ # x is the log-hazard ratio
100*(1 - exp(x))
}
rse_fun_inv <- function(x){
log(1 - (x/100))
}
res <- conf_dist(
estimate = log(0.72)
, stderr = 0.187618
, type = "coxreg"
, plot_type = "p_val"
, n_values = 1e4L
, est_names = c("RSE")
, conf_level = c(0.95, 0.8, 0.5)
, null_values = rse_fun_inv(0)
, trans = "rse_fun"
, alternative = "two_sided"
, log_yaxis = FALSE
, cut_logyaxis = 0.05
, xlab = "Relative survival effect (1 - HR%)"
, xlim = rse_fun_inv(c(-30, 60))
, together = FALSE
, plot_p_limit = 1 - 0.999
, plot_counternull = TRUE
, inverted = TRUE
, title = "Figure 1 in Bender et al. (2005)"
, x_scale = "default"
, plot = TRUE
)
# }
Run the code above in your browser using DataLab