The objective of this simulation experiment is to provide a toy example on how to use simChef
and showcase the automated R Markdown-generated documentation. For the sake of illustration, this toy simulation experiment studies the performance of linear regression at the surface-level with the sole purpose of facilitating an easy-to-understand walkthrough.
[Typically, the objective of the simulation experiment (and this blurb) will be more scientific than instructive and will warrant additional context/background and domain knowledge.]
In the Linear Gaussian DGP, we simulate the feature/design matrix \(\mathbf{X} \in \mathbb{R}^{n \times p}\) from a normal distribution and the response vector \(\mathbf{y} \in \mathbb{R}^n\) from a linear model. Specifically,
\[\begin{gather*}\mathbf{X} \sim N\left(\mathbf{0}, \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\right), \\\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon},\\\boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I}_n)\end{gather*}\]
Default Parameters in DGP
[In practice, documentation of DGPs should answer the questions “what” and “why”. That is, “what” is the DGP, and “why” are we using/studying it? As this simulation experiment is a contrived example, we omit the “why” here.]
Function
#> function(n, beta, rho, sigma) {
#> cov_mat <- matrix(c(1, rho, rho, 1), byrow = T, nrow = 2, ncol = 2)
#> X <- MASS::mvrnorm(n = n, mu = rep(0, 2), Sigma = cov_mat)
#> y <- X %*% beta + rnorm(n, sd = sigma)
#> return(list(X = X, y = y))
#> }
#> <bytecode: 0x564b4ddbb488>
Input Parameters
#> $n
#> [1] 200
#>
#> $beta
#> [1] 1 0
#>
#> $rho
#> [1] 0
#>
#> $sigma
#> [1] 1
Given some data \(\mathbf{X}\) and \(\mathbf{y}\), we fit ordinary least squares (OLS) and examine the p-values for each coefficient in the model. The p-values are computed using a two-sided t-test (see summary.lm()
).
Elaborating further on the testing, we are interested in testing:
\[\begin{align*}H_0: \beta_i = 0 \quad vs. \quad H_1: \beta_i \neq 0.\end{align*}\]
To test this, we compute the observed T-statistic, defined as
\[\begin{align*}T = \frac{\hat{\beta}_i}{\hat{SE}(\hat{\beta_i})},\end{align*}\]
and then compute the two-sided p-value under the t distribution with \(n - p - 1\) degrees of freedom. If the p-value is lower than some significance value \(\alpha\), then there is sufficient evidence to reject the null hypothesis \(H_0\). Otherwise, there is not sufficient evidence, and we fail to reject the null hypothesis.
[In practice, documentation of methods should answer the questions “what” and “why”. That is, “what” is the method, and “why” are we using/studying it? As this simulation experiment is a contrived example, we omit the “why” here.]
Function
#> function(X, y, cols = c("X1", "X2")) {
#> lm_fit <- lm(y ~ X)
#> pvals <- summary(lm_fit)$coefficients[cols, "Pr(>|t|)"] %>%
#> setNames(paste(names(.), "p-value"))
#> return(pvals)
#> }
#> <bytecode: 0x564b4d193390>
Input Parameters
#> list()
We define the rejection probability as the proportion of repetitions in the simulation experiment that result in a p-value \(\leq \alpha\). Here, we choose to set the significance level \(\alpha = 0.1\).
[In practice, documentation of evaluation metrics should answer the questions “what” and “why”. That is, “what” is the metric, and “why” are we using/studying it? As this simulation experiment is a contrived example, we omit the “why” here.]
Function
#> function(fit_results, alpha = 0.05) {
#> group_vars <- c(".dgp_name", ".method_name")
#> eval_out <- fit_results %>%
#> dplyr::group_by(across({{group_vars}})) %>%
#> dplyr::summarise(
#> `X1 Reject Prob.` = mean(`X1 p-value` < alpha),
#> `X2 Reject Prob.` = mean(`X2 p-value` < alpha)
#> )
#> return(eval_out)
#> }
Input Parameters
#> $alpha
#> [1] 0.1
Function
#> function(fit_results, vary_params = NULL, alpha = 0.05) {
#> group_vars <- c(".dgp_name", ".method_name", vary_params)
#> eval_out <- fit_results %>%
#> dplyr::group_by(across({{group_vars}})) %>%
#> dplyr::summarise(
#> `X1 Reject Prob.` = mean(`X1 p-value` < alpha),
#> `X2 Reject Prob.` = mean(`X2 p-value` < alpha)
#> )
#> return(eval_out)
#> }
#> <bytecode: 0x564b4e919d18>
Input Parameters
#> $alpha
#> [1] 0.1
Function
#> function(fit_results, vary_params = NULL, alpha = 0.05) {
#> group_vars <- c(".dgp_name", ".method_name", vary_params)
#> eval_out <- fit_results %>%
#> dplyr::group_by(across({{group_vars}})) %>%
#> dplyr::summarise(
#> `X1 Reject Prob.` = mean(`X1 p-value` < alpha),
#> `X2 Reject Prob.` = mean(`X2 p-value` < alpha)
#> )
#> return(eval_out)
#> }
#> <bytecode: 0x564b4ec36538>
Input Parameters
#> $alpha
#> [1] 0.1
Function
#> function(fit_results, vary_params = NULL, alpha = 0.05) {
#> group_vars <- c(".dgp_name", ".method_name", vary_params)
#> eval_out <- fit_results %>%
#> dplyr::group_by(across({{group_vars}})) %>%
#> dplyr::summarise(
#> `X1 Reject Prob.` = mean(`X1 p-value` < alpha),
#> `X2 Reject Prob.` = mean(`X2 p-value` < alpha)
#> )
#> return(eval_out)
#> }
Input Parameters
#> $alpha
#> [1] 0.1
To examine the power of the test, we plot the rejection probability as a function of \(\alpha\), that is, \(\mathbb{P}(\text{p-value} \leq \alpha)\) vs. \(\alpha\). If the coefficient is non-zero in the underlying DGP, then a larger AUC would indicate better performance in terms of the power. We will primarily focus on plotting the power of the first coefficient \(\beta_1\).
[In practice, documentation of the plotters should answer the questions “what” and “why”. That is, “what” is the plot, and “why” are we using/studying it? As this simulation experiment is a contrived example, we omit the “why” here.]
Function
#> function(fit_results, col = "X1") {
#> plt <- ggplot2::ggplot(fit_results) +
#> ggplot2::aes(x = .data[[paste(col, "p-value")]],
#> color = as.factor(.method_name)) +
#> ggplot2::geom_abline(slope = 1, intercept = 0,
#> color = "darkgray", linetype = "solid", size = 1) +
#> ggplot2::stat_ecdf(size = 1) +
#> ggplot2::scale_x_continuous(limits = c(0, 1)) +
#> ggplot2::labs(x = "t", y = "P( p-value \u2264 t )",
#> linetype = "", color = "Method")
#> return(plt)
#> }
Input Parameters
#> list()
Function
#> function(fit_results, vary_params = NULL, col = "X1") {
#>
#> if (!is.null(vary_params)) {
#> # deal with the case when we vary across a parameter that is vector-valued
#> if (is.list(fit_results[[vary_params]])) {
#> fit_results[[vary_params]] <- list_col_to_chr(fit_results[[vary_params]],
#> name = vary_params,
#> verbatim = TRUE)
#> }
#> }
#>
#> plt <- ggplot2::ggplot(fit_results) +
#> ggplot2::aes(x = .data[[paste(col, "p-value")]],
#> color = as.factor(.method_name)) +
#> ggplot2::geom_abline(slope = 1, intercept = 0,
#> color = "darkgray", linetype = "solid", linewidth = 1) +
#> ggplot2::stat_ecdf(size = 1) +
#> ggplot2::scale_x_continuous(limits = c(0, 1)) +
#> ggplot2::labs(x = "t", y = "P( p-value \u2264 t )",
#> linetype = "", color = "Method")
#> if (!is.null(vary_params)) {
#> plt <- plt + ggplot2::facet_wrap(~ .data[[vary_params]])
#> }
#> return(plt)
#> }
#> <bytecode: 0x564b49065620>
Input Parameters
#> list()
Function
#> function(fit_results, vary_params = NULL, col = "X1") {
#>
#> if (!is.null(vary_params)) {
#> # deal with the case when we vary across a parameter that is vector-valued
#> if (is.list(fit_results[[vary_params]])) {
#> fit_results[[vary_params]] <- list_col_to_chr(fit_results[[vary_params]],
#> name = vary_params,
#> verbatim = TRUE)
#> }
#> }
#>
#> plt <- ggplot2::ggplot(fit_results) +
#> ggplot2::aes(x = .data[[paste(col, "p-value")]],
#> color = as.factor(.method_name)) +
#> ggplot2::geom_abline(slope = 1, intercept = 0,
#> color = "darkgray", linetype = "solid", linewidth = 1) +
#> ggplot2::stat_ecdf(size = 1) +
#> ggplot2::scale_x_continuous(limits = c(0, 1)) +
#> ggplot2::labs(x = "t", y = "P( p-value \u2264 t )",
#> linetype = "", color = "Method")
#> if (!is.null(vary_params)) {
#> plt <- plt + ggplot2::facet_wrap(~ .data[[vary_params]])
#> }
#> return(plt)
#> }
#> <bytecode: 0x564b48bc18f8>
Input Parameters
#> list()
Function
#> function(fit_results, vary_params = NULL, col = "X1") {
#>
#> if (!is.null(vary_params)) {
#> # deal with the case when we vary across a parameter that is vector-valued
#> if (is.list(fit_results[[vary_params]])) {
#> fit_results[[vary_params]] <- list_col_to_chr(fit_results[[vary_params]],
#> name = vary_params,
#> verbatim = TRUE)
#> }
#> }
#>
#> plt <- ggplot2::ggplot(fit_results) +
#> ggplot2::aes(x = .data[[paste(col, "p-value")]],
#> color = as.factor(.method_name)) +
#> ggplot2::geom_abline(slope = 1, intercept = 0,
#> color = "darkgray", linetype = "solid", linewidth = 1) +
#> ggplot2::stat_ecdf(size = 1) +
#> ggplot2::scale_x_continuous(limits = c(0, 1)) +
#> ggplot2::labs(x = "t", y = "P( p-value \u2264 t )",
#> linetype = "", color = "Method")
#> if (!is.null(vary_params)) {
#> plt <- plt + ggplot2::facet_wrap(~ .data[[vary_params]])
#> }
#> return(plt)
#> }
#> <bytecode: 0x564b47c45170>
Input Parameters
#> list()
We plot the rejection probability for \(\beta_1\) across varying parameters of the DGP to understand how characteristics of the DGP affect the test.We define the rejection probability as the proportion of repetitions in the simulation experiment that result in a p-value \(\leq \alpha\). Here, we choose to set the significance level \(\alpha = 0.1\).
[In practice, documentation of the plotters should answer the questions “what” and “why”. That is, “what” is the plot, and “why” are we using/studying it? As this simulation experiment is a contrived example, we omit the “why” here.]
Function
#> function(eval_results, vary_params = NULL,
#> alpha = 0.05) {
#> eval_results <- eval_results$`Rejection Prob. (alpha = 0.1)`
#> if (is.list(eval_results[[vary_params]])) {
#> # deal with the case when we vary across a parameter that is vector-valued
#> eval_results[[vary_params]] <- list_col_to_chr(eval_results[[vary_params]],
#> name = vary_params,
#> verbatim = TRUE)
#> }
#> plt <- ggplot2::ggplot(eval_results) +
#> ggplot2::aes(x = .data[[vary_params]], y = `X1 Reject Prob.`,
#> color = as.factor(.method_name),
#> fill = as.factor(.method_name)) +
#> ggplot2::labs(x = vary_params,
#> y = sprintf("Rejection Probability (alpha = %s)", alpha),
#> color = "Method", fill = "Method") +
#> ggplot2::scale_y_continuous(limits = c(0, 1))
#> if (is.numeric(eval_results[[vary_params]])) {
#> plt <- plt +
#> ggplot2::geom_line() +
#> ggplot2::geom_point(size = 2)
#> } else {
#> plt <- plt +
#> ggplot2::geom_bar(stat = "identity")
#> }
#> return(plotly::ggplotly(plt))
#> }
#> <bytecode: 0x564b4900dc40>
Input Parameters
#> $alpha
#> [1] 0.1
Function
#> function(eval_results, vary_params = NULL,
#> alpha = 0.05) {
#> eval_results <- eval_results$`Rejection Prob. (alpha = 0.1)`
#> if (is.list(eval_results[[vary_params]])) {
#> # deal with the case when we vary across a parameter that is vector-valued
#> eval_results[[vary_params]] <- list_col_to_chr(eval_results[[vary_params]],
#> name = vary_params,
#> verbatim = TRUE)
#> }
#> plt <- ggplot2::ggplot(eval_results) +
#> ggplot2::aes(x = .data[[vary_params]], y = `X1 Reject Prob.`,
#> color = as.factor(.method_name),
#> fill = as.factor(.method_name)) +
#> ggplot2::labs(x = vary_params,
#> y = sprintf("Rejection Probability (alpha = %s)", alpha),
#> color = "Method", fill = "Method") +
#> ggplot2::scale_y_continuous(limits = c(0, 1))
#> if (is.numeric(eval_results[[vary_params]])) {
#> plt <- plt +
#> ggplot2::geom_line() +
#> ggplot2::geom_point(size = 2)
#> } else {
#> plt <- plt +
#> ggplot2::geom_bar(stat = "identity")
#> }
#> return(plotly::ggplotly(plt))
#> }
#> <bytecode: 0x564b48b95098>
Input Parameters
#> $alpha
#> [1] 0.1
Parameter Values
#> $dgp
#> $dgp$`Linear Gaussian DGP`
#> $dgp$`Linear Gaussian DGP`$beta
#> $dgp$`Linear Gaussian DGP`$beta[[1]]
#> [1] 1 0
#>
#> $dgp$`Linear Gaussian DGP`$beta[[2]]
#> [1] 1.0 0.5
#>
#> $dgp$`Linear Gaussian DGP`$beta[[3]]
#> [1] 1 1
#>
#> $dgp$`Linear Gaussian DGP`$beta[[4]]
#> [1] 1.0 1.5
#>
#> $dgp$`Linear Gaussian DGP`$beta[[5]]
#> [1] 1 2
#>
#>
#>
#>
#> $method
#> list()
Parameter Values
#> $dgp
#> $dgp$`Linear Gaussian DGP`
#> $dgp$`Linear Gaussian DGP`$rho
#> [1] 0.0 0.2 0.5 0.9
#>
#>
#>
#> $method
#> list()
Parameter Values
#> $dgp
#> $dgp$`Linear Gaussian DGP`
#> $dgp$`Linear Gaussian DGP`$sigma
#> [1] 1 2 4 8
#>
#>
#>
#> $method
#> list()