diff --git a/DESCRIPTION b/DESCRIPTION index 1ee78b72..1f5132ae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,9 @@ Imports: ggrepel, pbkrtest, lifecycle, - glue + cli, + glue, + withr Suggests: covr, devtools, diff --git a/R/concord1.R b/R/concord1.R index e03fd5d3..a406e2f1 100644 --- a/R/concord1.R +++ b/R/concord1.R @@ -4,8 +4,6 @@ #' #' @docType data #' @name concord1 -#' @references Hamilton, Lawrence C. 1983. Saving water: -#' A causal model of household conservation. -#' Sociological Perspectives 26(4):355-374. +#' @references Hamilton, Lawrence C. 1983. Saving water: A causal model of household conservation. Sociological Perspectives 26(4):355-374. #' @format A data.frame with 496 rows and 10 variables. -NULL +NULL \ No newline at end of file diff --git a/R/cop_pse_auxiliary.R b/R/cop_pse_auxiliary.R index 40b88b6b..1205cfaf 100644 --- a/R/cop_pse_auxiliary.R +++ b/R/cop_pse_auxiliary.R @@ -5,9 +5,9 @@ #' @return R2yz value #' @importFrom lavaan parameterEstimates cal_ryz <- function(ryxGz, R2){ - R2yz <- (ryxGz^2 - R2)/(ryxGz^2 - 1) + R2yz = (ryxGz^2 - R2)/(ryxGz^2 - 1) if (R2yz >= 0) { - ryz <- sqrt(R2yz) + ryz = sqrt(R2yz) } else { stop("Error! R2yz < 0!") } @@ -24,10 +24,10 @@ cal_ryz <- function(ryxGz, R2){ #' @return R2xz value #' @importFrom lavaan parameterEstimates cal_rxz <- function(var_x, var_y, R2, df, std_err){ - R2xz <- 1 - ((var_y * (1 - R2))/(var_x * df * std_err^2)) + R2xz = 1 - ((var_y * (1 - R2))/(var_x * df * std_err^2)) if (R2xz <= 0) {stop("Error! R2xz < 0!")} ## Note output the number of R2xz - rxz <- sqrt(R2xz) + rxz = sqrt(R2xz) return(rxz) } @@ -39,7 +39,7 @@ cal_rxz <- function(var_x, var_y, R2, df, std_err){ #' @return rxy value #' @importFrom lavaan parameterEstimates cal_rxy <- function(ryxGz, rxz, ryz){ - rxy <- ryxGz * sqrt((1 - rxz^2)*(1 - ryz^2)) + rxz * ryz + rxy = ryxGz * sqrt((1 - rxz^2)*(1 - ryz^2)) + rxz * ryz return(rxy) } @@ -66,33 +66,32 @@ cal_delta_star <- function(FR2max, est_uncond, rxz, n_obs) { - if (FR2max > .99) {FR2max <- .99} - #if (FR2max < R2 + inci) {FR2max = R2 + inci} check with Ken what this means - if (FR2max > R2) {D <- sqrt(FR2max - R2)} + if (FR2max > .99) {FR2max = .99} + # if (FR2max < R2 + inci) {FR2max = R2 + inci} check with Ken what this means + if (FR2max > R2) {D = sqrt(FR2max - R2)} #elements for computing Oster's delta_star - bt_m_b <- est_eff - eff_thr - rt_m_ro_t_syy <- (R2 - R2_uncond) * var_y - b0_m_b1 <- est_uncond - est_eff - rm_m_rt_t_syy <- (FR2max - R2) * var_y + bt_m_b = est_eff - eff_thr + rt_m_ro_t_syy = (R2 - R2_uncond) * var_y + b0_m_b1 = est_uncond - est_eff + rm_m_rt_t_syy = (FR2max - R2) * var_y - t_x <- var_x * (n_obs / (n_obs - 1)) * (1 - rxz^2) + t_x = var_x * (n_obs / (n_obs - 1)) * (1 - rxz^2) ## adjust df for var_x ## var_x is population variance, need sample variance from x - ## this adjustment is to get closer to what robomit generates as they - ## run regression using the sample data - num1 <- bt_m_b * rt_m_ro_t_syy * t_x - num2 <- bt_m_b * var_x * t_x * b0_m_b1^2 - num3 <- 2 * bt_m_b^2 * (t_x * b0_m_b1 * var_x) - num4 <- bt_m_b^3 * (t_x * var_x - t_x^2) - num <- num1 + num2 + num3 + num4 - den1 <- rm_m_rt_t_syy * b0_m_b1 * var_x - den2 <- bt_m_b * rm_m_rt_t_syy * (var_x - t_x) - den3 <- bt_m_b^2 * (t_x * b0_m_b1 * var_x) - den4 <- bt_m_b^3 * (t_x * var_x - t_x^2) - den <- den1 + den2 + den3 + den4 + ## this adjustment is to get closer to what robomit generates as they run regression using the sample data + num1 = bt_m_b * rt_m_ro_t_syy * t_x + num2 = bt_m_b * var_x * t_x * b0_m_b1^2 + num3 = 2 * bt_m_b^2 * (t_x * b0_m_b1 * var_x) + num4 = bt_m_b^3 * (t_x * var_x - t_x^2) + num = num1 + num2 + num3 + num4 + den1 = rm_m_rt_t_syy * b0_m_b1 * var_x + den2 = bt_m_b * rm_m_rt_t_syy * (var_x - t_x) + den3 = bt_m_b^2 * (t_x * b0_m_b1 * var_x) + den4 = bt_m_b^3 * (t_x * var_x - t_x^2) + den = den1 + den2 + den3 + den4 #obtain delta_star which is Oster's delta - delta_star <- num / den + delta_star = num / den return(delta_star) } @@ -111,17 +110,17 @@ cal_delta_star <- function(FR2max, #' @param rcvz correlation coefficient between V and Z #' @return list of model parameters #' @importFrom lavaan sem parameterEstimates -verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, +verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv, rxy, rxz, rzy, rcvy, rcvx, rcvz){ model <- 'Y ~ beta1 * X + beta2 * Z + beta3 * CV' - ccvy <- rcvy * sdcv * sdy # cov(cv, y) - ccvx <- rcvx * sdcv * sdx # cov(cv, x) - ccvz <- rcvz * sdcv * sdz - cxy <- rxy * sdx * sdy - czy <- rzy * sdz * sdy - cxz <- rxz * sdx * sdz + ccvy = rcvy * sdcv * sdy # cov(cv, y) + ccvx = rcvx * sdcv * sdx # cov(cv, x) + ccvz = rcvz * sdcv * sdz + cxy = rxy * sdx * sdy + czy = rzy * sdz * sdy + cxz = rxz * sdx * sdz # set up the covariance matrix cov.matrix <- matrix(c(sdy^2, cxy, czy, ccvy, @@ -138,11 +137,11 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, sample.nobs = n_obs) }, error = function(e){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) }, warning = function(w){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) } ) @@ -153,18 +152,12 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, sample.nobs = n_obs) ## the R2 extracted from summary is NOT right, do the calculation below R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[4,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$se - betaZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$est - seZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$se - betaCV <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta3',]$est - seCV <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta3',]$se + betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est + seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se + betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est + seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se + betaCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$est + seCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$se } #get regression based on true delta in terms of standardized coefficent @@ -182,11 +175,11 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, sample.nobs = n_obs) }, error = function(e){ - flag_cor <- FALSE + flag_cor = F return(flag_cor) }, warning = function(w){ - flag_cor <- FALSE + flag_cor = F return(flag_cor) } ) @@ -197,22 +190,16 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, sample.cov = cor.matrix, sample.nobs = n_obs) std_R2 <- 1 - lavaan::parameterEstimates(fit)[4,]$est - std_betaX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$est - std_seX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$se - std_betaZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$est - std_seZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$se - std_betaCV <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta3',]$est - std_seCV <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta3',]$se + std_betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est + std_seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se + std_betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est + std_seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se + std_betaCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$est + std_seCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$se } if ((inherits(flag_cor, "lavaan")) && (inherits(flag_cov, "lavaan"))) { - result <- list(R2, betaX, seX, betaZ, seZ, betaCV, seCV, + result = list(R2, betaX, seX, betaZ, seZ, betaCV, seCV, std_R2, std_betaX, std_seX, cov.matrix, cor.matrix) return(result) @@ -223,8 +210,7 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv, ## TO DO ## need to code the other solutions for pse as well -## then see how the different solutions perform in -## different scenarios (run the regression) +## then see how the different solutions perform in different scenarios (run the regression) cal_pse <- function(thr, kryx){ # calculations for preserving standard error i1 <- complex(real = 1, imaginary = -sqrt(3)) @@ -261,13 +247,13 @@ verify_manual <- function(rxy, rxz, rxcv, ryz, rycv, rzcv, sdy, sdx){ return(beta) } -verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ +verify_reg_Gz = function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ model <- 'Y ~ beta1 * X + beta2 * Z' - cxy <- rxy * sdx * sdy - czy <- rzy * sdz * sdy - cxz <- rxz * sdx * sdz + cxy = rxy * sdx * sdy + czy = rzy * sdz * sdy + cxz = rxz * sdx * sdz # set up the covariance matrix cov.matrix <- matrix(c(sdy^2, cxy, czy, @@ -283,11 +269,11 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ sample.nobs = n_obs) }, error = function(e){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) }, warning = function(w){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) } ) @@ -297,19 +283,15 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ sample.cov = cov.matrix, sample.nobs = n_obs) ## the R2 extracted from summary is NOT right, do the calculation below - R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[3,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$se - betaZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$est - seZ <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta2',]$se + R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[3,]$est) / sdy^2 + betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est + seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se + betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est + seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se } if (inherits(flag_cov, "lavaan")) { - result <- list(R2, betaX, seX, betaZ, seZ) + result = list(R2, betaX, seX, betaZ, seZ) return(result) } else { stop("Error!") @@ -324,10 +306,10 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ #' @param rxy correlation coefficient between X and Y #' @return list of model parameters #' @importFrom lavaan sem parameterEstimates -verify_reg_uncond <- function(n_obs, sdx, sdy, rxy){ +verify_reg_uncond = function(n_obs, sdx, sdy, rxy){ model <- 'Y ~ beta1 * X' - cxy <- rxy * sdx * sdy + cxy = rxy * sdx * sdy # set up the covariance matrix cov.matrix <- matrix(c(sdy^2, cxy, @@ -342,11 +324,11 @@ verify_reg_uncond <- function(n_obs, sdx, sdy, rxy){ sample.nobs = n_obs) }, error = function(e){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) }, warning = function(w){ - flag_cov <- FALSE + flag_cov = F return(flag_cov) } ) @@ -357,16 +339,14 @@ verify_reg_uncond <- function(n_obs, sdx, sdy, rxy){ sample.nobs = n_obs) ## the R2 extracted from summary is NOT right, do the calculation below R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[2,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[ - lavaan::parameterEstimates(fit)$label == 'beta1',]$se + betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est + seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se } if (inherits(flag_cov, "lavaan")) { - result <- list(R2, betaX, seX) + result = list(R2, betaX, seX) return(result) } else { stop("Error!") } -} +} \ No newline at end of file diff --git a/R/core-sensitivity-mkonfound.R b/R/core-sensitivity-mkonfound.R index 442a75d8..dae16b94 100644 --- a/R/core-sensitivity-mkonfound.R +++ b/R/core-sensitivity-mkonfound.R @@ -31,13 +31,11 @@ core_sensitivity_mkonfound <- function(t, df, alpha = .05, tails = 2) { r_con <- round(sqrt(abs(itcv)), 3) out <- tibble::tibble(t, df, action, inference, pct_bias, itcv, r_con) - names(out) <- c("t", "df", "action", "inference", - "pct_bias_to_change_inference", "itcv", "r_con") - out$pct_bias_to_change_inference <- round(out$pct_bias_to_change_inference, - 3) + names(out) <- c("t", "df", "action", "inference", "pct_bias_to_change_inference", "itcv", "r_con") + out$pct_bias_to_change_inference <- round(out$pct_bias_to_change_inference, 3) out$itcv <- round(out$itcv, 3) out$action <- as.character(out$action) out$inference <- as.character(out$inference) return(out) -} +} \ No newline at end of file diff --git a/R/helper_output_dataframe.R b/R/helper_output_dataframe.R index 957e88e9..ed71d6a1 100644 --- a/R/helper_output_dataframe.R +++ b/R/helper_output_dataframe.R @@ -15,15 +15,8 @@ #' @param non_linear flag for non-linear models #' @return data frame with model information #' @importFrom dplyr tibble -output_df <- function(est_eff, - beta_threshhold, - unstd_beta, - bias = NULL, - sustain = NULL, - recase, obs_r, - critical_r, - r_con, - itcv, +output_df <- function(est_eff, beta_threshhold, unstd_beta, bias = NULL, + sustain = NULL, recase, obs_r, critical_r, r_con, itcv, non_linear) { if (abs(est_eff) > abs(beta_threshhold)) { df <- dplyr::tibble( diff --git a/R/helper_output_print.R b/R/helper_output_print.R index f98fedad..ccd11d12 100644 --- a/R/helper_output_print.R +++ b/R/helper_output_print.R @@ -25,48 +25,27 @@ #' @param alpha The level of statistical significance. #' @param index A character string indicating the index for which the output is generated ('RIR' or 'IT'). #' @importFrom crayon bold underline italic -output_print <- function(eff_diff, - beta_threshhold, - bias = NULL, - sustain = NULL, - nu, - recase, - obs_r, - critical_r, - r_con, - itcv, - alpha, +output_print <- function(eff_diff, beta_threshhold, bias = NULL, sustain = NULL, + nu, recase, obs_r, critical_r, r_con, itcv, alpha, index) { - if (index == "RIR") { - cat(crayon::bold("Robustness of Inference to Replacement (RIR):\n")) - if (abs(eff_diff) > abs(beta_threshhold)) { - cat("To invalidate an inference, ", round(bias, 3), - "% of the estimate would have to be due to bias. ") - cat("\n") - cat("This is based on a threshold of ", round(beta_threshhold, 3), - " for statistical significance (alpha = ", alpha, ").\n", - sep = "") - cat("\n") - cat("To invalidate an inference, ", round(recase, 3), - " observations would have to be replaced with cases") - cat("\n") - cat("for which the effect is ", nu, - " (RIR = ", round(recase, 3), ").\n", sep = "") - cat("\n") + if (index == "RIR"){ + cat(crayon::bold("Robustness of Inference to Replacement (RIR):\n")) + if (abs(eff_diff) > abs(beta_threshhold)) { + cat("To invalidate an inference, ", round(bias, 3), "% of the estimate would have to be due to bias. ") + cat("\n") + cat("This is based on a threshold of ", round(beta_threshhold, 3), " for statistical significance (alpha = ", alpha, ").\n", sep = "") + cat("\n") + cat("To invalidate an inference, ", round(recase, 3), " observations would have to be replaced with cases") + cat("\n") + cat("for which the effect is ", nu, " (RIR = ", round(recase, 3), ").\n", sep = "") + cat("\n") } else if (abs(eff_diff) < abs(beta_threshhold)) { - cat("To sustain an inference, ", round(sustain, 3), - "% of the estimate would have to be due to bias. ") + cat("To sustain an inference, ", round(sustain, 3), "% of the estimate would have to be due to bias. ") cat("\n") - cat("This is based on a threshold of ", round(beta_threshhold, 3), - " for statistical significance (alpha = ", alpha, ").\n", - sep = "") + cat("This is based on a threshold of ", round(beta_threshhold, 3), " for statistical significance (alpha = ", alpha, ").\n", sep = "") cat("\n") - cat("To sustain an inference, ", round(recase, 3), - " of the cases with ", nu, - " effect would have to be replaced with cases at the - threshold of inference", - " (RIR = ", round(recase, 3), ").\n", sep = "") + cat("To sustain an inference, ", round(recase, 3), " of the cases with ", nu, " effect would have to be replaced with cases at the threshold of inference"," (RIR = ", round(recase, 3), ").\n", sep = "") } else if (eff_diff == beta_threshhold) { warning("The coefficient is exactly equal to the threshold.\n") @@ -74,94 +53,65 @@ output_print <- function(eff_diff, cat("See Frank et al. (2013) for a description of the method.") cat("\n") cat("\n") - cat(crayon::underline("Citation:"), - "Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).") + cat(crayon::underline("Citation:"), "Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).") cat("\n") cat("What would it take to change an inference?") cat("\n") - cat("Using Rubin's causal model to interpret the - robustness of causal inferences.") + cat("Using Rubin's causal model to interpret the robustness of causal inferences.") cat("\n") - cat(crayon::italic("Education, Evaluation and - Policy Analysis, 35"), "437-460.") + cat(crayon::italic("Education, Evaluation and Policy Analysis, 35"), "437-460.") cat("\n") } if (index == "IT") { cat(crayon::bold("Impact Threshold for a Confounding Variable:\n")) if (abs(obs_r) > abs(critical_r) & obs_r > 0) { - cat("The minimum impact of an omitted variable to invalidate an - inference for a null hypothesis of 0 effect is based on - a correlation of ", r_con) + cat("The minimum impact of an omitted variable to invalidate an inference for a null hypothesis of 0 effect is based on a correlation of ", r_con) cat("\n") cat(" with the outcome and at ", r_con, - " with the predictor of interest - (conditioning on observed covariates) based on a threshold of ") + " with the predictor of interest (conditioning on observed covariates) based on a threshold of ") cat("\n") - cat(round(critical_r, 3), " for statistical - significance (alpha = ", alpha, ").\n",sep = "") + cat(round(critical_r, 3), " for statistical significance (alpha = ", alpha, ").\n", + sep = "") cat("\n") - cat("Correspondingly the impact of an omitted variable - (as defined in Frank 2000) must be ") + cat("Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be ") cat("\n") - cat(r_con, " X ", r_con, " = ", round(r_con^2, 3), - " to invalidate an inference for a null hypothesis of 0 effect. - \n", sep = "") + cat(r_con, " X ", r_con, " = ", round(r_con^2, 3), " to invalidate an inference for a null hypothesis of 0 effect.\n", sep = "") } else if (abs(obs_r) > abs(critical_r) & obs_r < 0) { - cat("The minimum (in absolute value) impact of an omitted - variable to invalidate an inference for a null hypothesis of 0 - effect is based on a correlation of ", -r_con) + cat("The minimum (in absolute value) impact of an omitted variable to invalidate an inference for a null hypothesis of 0 effect is based on a correlation of ", -r_con) cat("\n") cat(" with the outcome and at ", r_con, - " with the predictor of interest (conditioning on observed - covariates; signs are interchangable) based on a threshold of ") + " with the predictor of interest (conditioning on observed covariates; signs are interchangable) based on a threshold of ") cat("\n") - cat(round(critical_r, 3), " for statistical significance - (alpha = ", alpha, ").\n", + cat(round(critical_r, 3), " for statistical significance (alpha = ", alpha, ").\n", sep = "") cat("\n") - cat("Correspondingly the impact of an omitted variable - (as defined in Frank 2000) must be ") + cat("Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be ") cat("\n") - cat(-r_con, " X ", r_con, " = ", -round(r_con^2, 3), - " to invalidate an inference for a null hypothesis - of 0 effect.\n", sep = "") + cat(-r_con, " X ", r_con, " = ", -round(r_con^2, 3), " to invalidate an inference for a null hypothesis of 0 effect.\n", sep = "") } else if (abs(obs_r) < abs(critical_r) & obs_r > 0) { - cat("The maximum impact (in absolute value) of an omitted variable - to sustain an inference for a null hypothesis of 0 effect is - based on a correlation of ", -r_con) + cat("The maximum impact (in absolute value) of an omitted variable to sustain an inference for a null hypothesis of 0 effect is based on a correlation of ", -r_con) cat("\n") cat(" with the outcome and at ", r_con, - " with the predictor of interest (conditioning on observed - covariates; signs are interchangable) based on a threshold of ", - round(beta_threshhold, 3)) + " with the predictor of interest (conditioning on observed covariates; signs are interchangable) based on a threshold of ", round(beta_threshhold, 3)) cat("\n") cat(" for statistical significance (alpha = ", alpha, ").\n", sep = "") cat("\n") - cat("Correspondingly the impact of an omitted variable - (as defined in Frank 2000) must be ") + cat("Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be ") cat("\n") - cat(-r_con, " X ", r_con, " = ", -round(r_con^2, 3), - " to sustain an inference for a null hypothesis of 0 effect.\n", - sep = "") + cat(-r_con, " X ", r_con, " = ", -round(r_con^2, 3), " to sustain an inference for a null hypothesis of 0 effect.\n", sep = "") } else if (abs(obs_r) < abs(critical_r) & obs_r < 0) { - cat("The maximum impact of an omitted variable to sustain an - inference for a null hypothesis of 0 effect is based on a - correlation of ", r_con) + cat("The maximum impact of an omitted variable to sustain an inference for a null hypothesis of 0 effect is based on a correlation of ", r_con) cat("\n") cat(" with the outcome and at ", r_con, - " with the predictor of interest (conditioning on observed - covariates) based on a threshold of ", round(beta_threshhold, 3)) + " with the predictor of interest (conditioning on observed covariates) based on a threshold of ", round(beta_threshhold, 3)) cat("\n") cat(" for statistical significance (alpha = ", alpha, ").\n", sep = "") cat("\n") - cat("Correspondingly the impact of an omitted variable - (as defined in Frank 2000) must be ") + cat("Correspondingly the impact of an omitted variable (as defined in Frank 2000) must be ") cat("\n") - cat(r_con, " X ", r_con, " = ", round(r_con^2, 3), - " to sustain an inference for a null hypothesis of 0 effect.\n", - sep = "") + cat(r_con, " X ", r_con, " = ", round(r_con^2, 3), " to sustain an inference for a null hypothesis of 0 effect.\n", sep = "") } else if (obs_r == critical_r) { warning("The correlation is exactly equal to the threshold.\n") } @@ -172,10 +122,9 @@ output_print <- function(eff_diff, cat("\n") cat("Frank, K. (2000). Impact of a confounding variable on the") cat("\n") - cat("inference of a regression coefficient.", - crayon::italic("Sociological Methods and Research, 29"), - "(2), 147-194") + cat("inference of a regression coefficient.", crayon::italic("Sociological Methods and Research, 29"), "(2), 147-194") cat("\n") } } + diff --git a/R/helper_output_table.R b/R/helper_output_table.R index 479e6601..bd3892f0 100644 --- a/R/helper_output_table.R +++ b/R/helper_output_table.R @@ -24,27 +24,19 @@ output_table <- function(model_object, tested_variable) { model_output$itcv <- NA var_row <- model_output$term == tested_variable - model_output$itcv[var_row] <- abs(konfound(model_object, - !!tested_variable, - to_return = "raw_output")$itcv) + model_output$itcv[var_row] <- abs(konfound(model_object, !!tested_variable, to_return = "raw_output")$itcv) - covariate_names <- model_output$term[ - !(model_output$term %in% c("(Intercept)", tested_variable))] + covariate_names <- model_output$term[!(model_output$term %in% c("(Intercept)", tested_variable))] for (i in seq(covariate_names)) { cov_row <- model_output$term == covariate_names[i] d <- model_object$model cor_df <- as.data.frame(stats::cor(d)) - model_output$itcv[cov_row] <- round( - abs(cor_df[cov_row, 1]) * abs(cor_df[cov_row, tested_variable]), - 3) # r_zy * r_zx + model_output$itcv[cov_row] <- round(abs(cor_df[cov_row, 1]) * abs(cor_df[cov_row, tested_variable]), 3) # r_zy * r_zx } - model_output <- purrr::modify_if(model_output, - is.numeric, - round, - digits = 3) + model_output <- purrr::modify_if(model_output, is.numeric, round, digits = 3) return(model_output) } diff --git a/R/helper_plot_correlation.R b/R/helper_plot_correlation.R index 9410d956..b953d893 100644 --- a/R/helper_plot_correlation.R +++ b/R/helper_plot_correlation.R @@ -25,111 +25,28 @@ plot_correlation <- function(r_con, obs_r, critical_r) { p <- ggplot2::ggplot(d, ggplot2::aes_string(x = "x", y = "y")) + -# main arrows -ggplot2::geom_segment( - ggplot2::aes(y = .1), - xend = 0, - yend = .9, - arrow = ggplot2::arrow(), - size = 2.5, - color = "#A6CEE3" -) + # straight up -ggplot2::geom_segment( - ggplot2::aes(x = .1), - xend = 1, - yend = .9, - arrow = ggplot2::arrow(), - size = 2.5, - color = "#A6CEE3" -) + # hypotenuse -ggplot2::geom_segment( - ggplot2::aes(x = .15, y = 1), - xend = .9, - yend = 1, - arrow = ggplot2::arrow(), - size = 2.5, - color = "#A6CEE3" -) + # straight across - - -# label annotations -ggplot2::annotate( - "text", - x = 0, - y = 0, - label = paste0("Confounding\nVariable"), - fontface = 3 -) + -ggplot2::annotate( - "text", - x = 0, - y = 1, - label = paste0("Predictor of Interest"), - fontface = 3 -) + -ggplot2::annotate( - "text", - x = 1, - y = 1, - label = paste0("Outcome"), - fontface = 3 -) + - - -# CV arrows -# ggplot2::geom_segment(ggplot2::aes(x = .05, y = .25), xend = .275, -#yend = .65, arrow = ggplot2::arrow(), size = 2.5, color = "#1F78B4") + -# straight across -# ggplot2::geom_segment(ggplot2::aes(x = .175, y = .15), xend = .3, -#yend = .625, arrow = ggplot2::arrow(), size = 2.5, color = "#1F78B4") + - # straight across -ggplot2::geom_curve( - ggplot2::aes(x = .04, y = .325, xend = .35, yend = .825), - curvature = -.35, - size = 2.5, - color = "#1F78B4", - arrow = ggplot2::arrow() -) + -ggplot2::geom_curve( - ggplot2::aes(x = .225, y = .23, xend = .4, yend = .8), - curvature = .35, - size = 2.5, - color = "#1F78B4", - arrow = ggplot2::arrow() -) + -ggplot2::geom_segment( - ggplot2::aes(x = .37, y = .81), - xend = .465, - yend = .925, - arrow = ggplot2::arrow(), - size = 2.5, - color = "#1F78B4" -) + # straight across - - -# corr. annotations -ggplot2::annotate( - "text", - x = -.125, - y = .5, - label = paste0("Rx.cv | Z =\n ", r_con), - fontface = 1 -) + -ggplot2::annotate( - "text", - x = .575, - y = .35, - label = paste0("Ry.cv | Z =\n", r_con), - fontface = 1 -) + -ggplot2::annotate( - "text", - x = .25, - y = .525, - label = paste0("Rx.cv | Z * Ry.cv | Z =\n", round(r_con^2, 3)), - fontface = 1 -) + - + # main arrows + ggplot2::geom_segment(ggplot2::aes(y = .1), xend = 0, yend = .9, arrow = ggplot2::arrow(), linewidth = 2.5, color = "#A6CEE3") + # straight up + ggplot2::geom_segment(ggplot2::aes(x = .1), xend = 1, yend = .9, arrow = ggplot2::arrow(), linewidth = 2.5, color = "#A6CEE3") + # hypotenuse + ggplot2::geom_segment(ggplot2::aes(x = .15, y = 1), xend = .9, yend = 1, arrow = ggplot2::arrow(), linewidth = 2.5, color = "#A6CEE3") + # straight across + + # lagel annotations + ggplot2::annotate("text", x = 0, y = 0, label = paste0("Confounding\nVariable"), fontface = 3) + + ggplot2::annotate("text", x = 0, y = 1, label = paste0("Predictor of Interest"), fontface = 3) + + ggplot2::annotate("text", x = 1, y = 1, label = paste0("Outcome"), fontface = 3) + + + # CV arrows + # ggplot2::geom_segment(ggplot2::aes(x = .05, y = .25), xend = .275, yend = .65, arrow = ggplot2::arrow(), size = 2.5, color = "#1F78B4") + # straight across + # ggplot2::geom_segment(ggplot2::aes(x = .175, y = .15), xend = .3, yend = .625, arrow = ggplot2::arrow(), size = 2.5, color = "#1F78B4") + # straight across + ggplot2::geom_curve(ggplot2::aes(x = .04, y = .325, xend = .35, yend = .825), curvature = -.35, linewidth = 2.5, color = "#1F78B4", arrow = ggplot2::arrow()) + + ggplot2::geom_curve(ggplot2::aes(x = .225, y = .23, xend = .4, yend = .8), curvature = .35, linewidth = 2.5, color = "#1F78B4", arrow = ggplot2::arrow()) + + + ggplot2::geom_segment(ggplot2::aes(x = .37, y = .81), xend = .465, yend = .925, arrow = ggplot2::arrow(), linewidth = 2.5, color = "#1F78B4") + # straight across + + # corr. annotations + ggplot2::annotate("text", x = -.125, y = .5, label = paste0("Rx.cv | Z =\n ", r_con), fontface = 1) + + ggplot2::annotate("text", x = .575, y = .35, label = paste0("Ry.cv | Z =\n ", r_con), fontface = 1) + + ggplot2::annotate("text", x = .25, y = .525, label = paste0("Rx.cv | Z * Ry.cv | Z =\n ", round(r_con^2, 3)), fontface = 1) + # plot modifications ggplot2::xlim(-.15, 1.1) + diff --git a/R/helper_plot_threshold.R b/R/helper_plot_threshold.R index 7dcff941..764d45d1 100644 --- a/R/helper_plot_threshold.R +++ b/R/helper_plot_threshold.R @@ -40,8 +40,7 @@ plot_threshold <- function(beta_threshold, est_eff) { } else if (est_eff < beta_threshold) { # beta is below threshold dd <- dplyr::tibble(est_eff = est_eff, beta_threshold = beta_threshold) - dd <- dplyr::mutate(dd, `Above Estimated Effect, - Below Threshold` = abs(est_eff - beta_threshold)) + dd <- dplyr::mutate(dd, `Above Estimated Effect, Below Threshold` = abs(est_eff - beta_threshold)) dd <- dplyr::mutate(dd, `Below Threshold` = est_eff) dd <- dplyr::select(dd, -beta_threshold) @@ -57,26 +56,22 @@ plot_threshold <- function(beta_threshold, est_eff) { y_thresh_text <- y_thresh + sqrt(.005 * abs(y_thresh)) effect_text <- est_eff + sqrt(.025 * abs(est_eff)) # y-value of text - cols <- c("lightgray", "#1F78B4") - # dark blue and green (green not used right now) + cols <- c("lightgray", "#1F78B4") # dark blue and green (green not used right now) } p <- ggplot2::ggplot(dd, ggplot2::aes(x = inference, y = val, fill = key)) + ggplot2::geom_col(position = "stack") + ggplot2::geom_hline(yintercept = est_eff, color = "black") + - ggplot2::annotate("text", x = 1, y = effect_text, - label = "Estimated Effect") + + ggplot2::annotate("text", x = 1, y = effect_text, label = "Estimated Effect") + + ggplot2::geom_hline(yintercept = y_thresh, color = "red") + - ggplot2::annotate("text", x = 1, y = y_thresh_text, - label = "Threshold") + - # ggplot2::geom_text(aes(label = "Effect"), vjust = -.5) + - #this is discussed here: https://github.com/jrosen48/konfound/issues/5 + ggplot2::annotate("text", x = 1, y = y_thresh_text, label = "Threshold") + + # ggplot2::geom_text(aes(label = "Effect"), vjust = -.5) + this is discussed here: https://github.com/jrosen48/konfound/issues/5 ggplot2::scale_fill_manual("", values = cols) + ggplot2::theme_bw() + - ggplot2::theme(axis.text.x = ggplot2::element_blank(), - axis.ticks = ggplot2::element_blank()) + + ggplot2::theme(axis.text.x = ggplot2::element_blank(), axis.ticks = ggplot2::element_blank()) + ggplot2::xlab(NULL) + ggplot2::ylab("Effect (abs. value)") + ggplot2::theme(legend.position = "top") diff --git a/R/konfound-glm-dichotomous.R b/R/konfound-glm-dichotomous.R index 09fef4ef..f3662a90 100644 --- a/R/konfound-glm-dichotomous.R +++ b/R/konfound-glm-dichotomous.R @@ -17,13 +17,8 @@ #' @return The results of the konfound analysis. #' @importFrom broom tidy glance #' @importFrom stats glm -konfound_glm_dichotomous <- function(model_object, - tested_variable_string, - alpha, tails, - to_return, - n_treat, - switch_trm, - replace) { +konfound_glm_dichotomous <- function(model_object, tested_variable_string, alpha, tails, + to_return, n_treat, switch_trm, replace) { tidy_output <- broom::tidy(model_object) # tidying output glance_output <- broom::glance(model_object) diff --git a/R/konfound-glm.R b/R/konfound-glm.R index 60c97d32..66c0e32f 100644 --- a/R/konfound-glm.R +++ b/R/konfound-glm.R @@ -18,19 +18,13 @@ #' @importFrom dplyr select filter bind_cols #' @importFrom stats glm #' @importFrom margins margins -konfound_glm <- function(model_object, - tested_variable_string, - alpha, tails, - index = "RIR", - to_return) { +konfound_glm <- function(model_object, tested_variable_string, alpha, tails, index = "RIR", to_return) { tidy_output <- broom::tidy(model_object) # tidying output glance_output <- broom::glance(model_object) coef_df <- tidy_output[tidy_output$term == tested_variable_string, ] est_eff <- coef_df$estimate - est_eff <- suppressWarnings( - summary(margins::margins(model_object))$AME[names(summary( - margins::margins(model_object))$AME) == tested_variable_string]) + est_eff <- suppressWarnings(summary(margins::margins(model_object))$AME[names(summary(margins::margins(model_object))$AME) == tested_variable_string]) std_err <- coef_df$std.error n_obs <- glance_output$nobs n_covariates <- glance_output$df.null - glance_output$df.residual diff --git a/R/konfound-lm.R b/R/konfound-lm.R index bc57726b..3dde1166 100644 --- a/R/konfound-lm.R +++ b/R/konfound-lm.R @@ -16,12 +16,7 @@ #' @return The results of the konfound analysis for the specified variable(s). #' @importFrom broom tidy glance #' @importFrom dplyr select filter bind_cols -konfound_lm <- function(model_object, - tested_variable_string, - alpha, - tails, - index, - to_return) { +konfound_lm <- function(model_object, tested_variable_string, alpha, tails, index, to_return) { tidy_output <- broom::tidy(model_object) # tidying output glance_output <- broom::glance(model_object) diff --git a/R/konfound-lmer.R b/R/konfound-lmer.R index ce3652c8..942957ee 100644 --- a/R/konfound-lmer.R +++ b/R/konfound-lmer.R @@ -10,8 +10,7 @@ get_kr_df <- function(model_object) { L <- diag(rep(1, length(lme4::fixef(model_object)))) L <- as.data.frame(L) - out <- suppressWarnings(purrr::map_dbl(L, pbkrtest::get_Lb_ddf, - object = model_object)) + out <- suppressWarnings(purrr::map_dbl(L, pbkrtest::get_Lb_ddf, object = model_object)) names(out) <- names(lme4::fixef(model_object)) out } @@ -32,13 +31,7 @@ get_kr_df <- function(model_object) { #' @return The results of the konfound analysis for the specified fixed effect(s). #' @importFrom broom.mixed tidy #' @importFrom dplyr filter bind_cols -konfound_lmer <- function(model_object, - tested_variable_string, - test_all, - alpha, - tails, - index, - to_return) { +konfound_lmer <- function(model_object, tested_variable_string, test_all, alpha, tails, index, to_return) { tidy_output <- broom.mixed::tidy(model_object) # tidying output coef_df <- tidy_output[tidy_output$term == tested_variable_string, ] diff --git a/R/konfound.R b/R/konfound.R index b514d493..63cfbcfc 100644 --- a/R/konfound.R +++ b/R/konfound.R @@ -69,14 +69,11 @@ konfound <- function(model_object, # Stop messages if (!(class(model_object)[1] %in% c("lm", "glm", "lmerMod"))) { - stop("konfound() is currently implemented for models estimated with - lm(), glm(), and lme4::lmer(); consider using pkonfound() instead") + stop("konfound() is currently implemented for models estimated with lm(), glm(), and lme4::lmer(); consider using pkonfound() instead") } # Dealing with non-standard evaluation - tested_variable_enquo <- rlang::enquo(tested_variable) - # dealing with non-standard evaluation - #(so unquoted names for tested_variable can be used) + tested_variable_enquo <- rlang::enquo(tested_variable) # dealing with non-standard evaluation (so unquoted names for tested_variable can be used) tested_variable_string <- rlang::quo_name(tested_variable_enquo) # Dispatching based on class @@ -94,14 +91,9 @@ konfound <- function(model_object, } if (inherits(model_object, "glm") & two_by_two == FALSE) { - message("Note that for a non-linear model, - impact threshold should not be interpreted.") - message("Note that this is only an approximation. For exact results - in terms of the number of cases that must be switched from treatment - success to treatment failure to invalidate the inference see: - https://msu.edu/~kenfrank/non%20linear%20replacement%20treatment.xlsm") - message("If a dichotomous independent variable is used, consider using - the 2X2 table approach enabled with the argument two_by_two = TRUE") + message("Note that for a non-linear model, impact threshold should not be interpreted.") + message("Note that this is only an approximation. For exact results in terms of the number of cases that must be switched from treatment success to treatment failure to invalidate the inference see: https://msu.edu/~kenfrank/non%20linear%20replacement%20treatment.xlsm") + message("If a dichotomous independent variable is used, consider using the 2X2 table approach enabled with the argument two_by_two = TRUE") output <- konfound_glm( model_object = model_object, tested_variable_string = tested_variable_string, @@ -115,8 +107,7 @@ konfound <- function(model_object, if (inherits(model_object, "glm") & two_by_two == TRUE) { - if(is.null(n_treat)) stop("Please provide a value for n_treat to use - this functionality with a dichotomous predictor") + if(is.null(n_treat)) stop("Please provide a value for n_treat to use this functionality with a dichotomous predictor") output <- konfound_glm_dichotomous( model_object = model_object, @@ -143,15 +134,7 @@ konfound <- function(model_object, to_return = to_return ) - message("Note that the Kenward-Roger approximation is used to - estimate degrees of freedom for the predictor(s) of interest. - We are presently working to add other methods for calculating - the degrees of freedom for the predictor(s) of interest. - If you wish to use other methods now, consider those detailed here: - https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html - #why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have. - You can then enter degrees of freedom obtained from another method along with the coefficient, - number of observations, and number of covariates to the pkonfound() function to quantify the robustness of the inference.") + message("Note that the Kenward-Roger approximation is used to estimate degrees of freedom for the predictor(s) of interest. We are presently working to add other methods for calculating the degrees of freedom for the predictor(s) of interest. If you wish to use other methods now, consider those detailed here: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have. You can then enter degrees of freedom obtained from another method along with the coefficient, number of observations, and number of covariates to the pkonfound() function to quantify the robustness of the inference.") return(output) } @@ -160,4 +143,4 @@ konfound <- function(model_object, message("For more detailed output, consider setting `to_return` to table") } -} +} \ No newline at end of file diff --git a/R/mkonfound-data.R b/R/mkonfound-data.R index 299c2faf..db063821 100644 --- a/R/mkonfound-data.R +++ b/R/mkonfound-data.R @@ -1,8 +1,7 @@ #' Example data for the mkonfound function #' #' A dataset containing t and df values from example studies from Educational Evaluation -#' and Policy Analysis (as detailed in Frank et al., 2013): -#' https://drive.google.com/file/d/1aGhxGjvMvEPVAgOA8rrxvA97uUO5TTMe/view +#' and Policy Analysis (as detailed in Frank et al., 2013): https://drive.google.com/file/d/1aGhxGjvMvEPVAgOA8rrxvA97uUO5TTMe/view #' #' @format A data frame with 30 rows and 2 variables: #' \describe{ @@ -11,4 +10,4 @@ #' ... #' } #' @source \url{https://drive.google.com/file/d/1aGhxGjvMvEPVAgOA8rrxvA97uUO5TTMe/view} -"mkonfound_ex" +"mkonfound_ex" \ No newline at end of file diff --git a/R/mkonfound.R b/R/mkonfound.R index 2ab2425b..bbf62beb 100644 --- a/R/mkonfound.R +++ b/R/mkonfound.R @@ -37,12 +37,10 @@ mkonfound <- function(d, t, df, alpha = .05, tails = 2, return_plot = FALSE) { df_vec <- pull(select(d, !!df_enquo)) if (length(t_vec) <= 1) { - stop("To carry out sensitivity analysis - for a single study, use pkonfound()") + stop("To carry out sensitivity analysis for a single study, use pkonfound()") } - results_df <- purrr::map2_dfr(.x = t_vec, .y = df_vec, - .f = core_sensitivity_mkonfound) + results_df <- purrr::map2_dfr(.x = t_vec, .y = df_vec, .f = core_sensitivity_mkonfound) if (return_plot == TRUE) { results_df$action <- dplyr::case_when( @@ -50,10 +48,7 @@ mkonfound <- function(d, t, df, alpha = .05, tails = 2, return_plot = FALSE) { results_df$action == "to_sustain" ~ "To Sustain" ) - p <- ggplot2::ggplot(results_df, - ggplot2::aes( - x = results_df$pct_bias_to_change_inference, - fill = results_df$action)) + + p <- ggplot2::ggplot(results_df, ggplot2::aes(x = results_df$pct_bias_to_change_inference, fill = results_df$action)) + ggplot2::geom_histogram() + ggplot2::scale_fill_manual("", values = c("#1F78B4", "#A6CEE3")) + ggplot2::theme_bw() + diff --git a/R/nonlinear_auxiliary.R b/R/nonlinear_auxiliary.R index d2290955..92a90886 100644 --- a/R/nonlinear_auxiliary.R +++ b/R/nonlinear_auxiliary.R @@ -683,18 +683,16 @@ if (allnotenough) { total_switch <- final + allnotenough*final_extra -result <- list(final_switch = final, User_enter_value = table_start, - Transfer_Table = table_final, +result <- list(final_switch = final, User_enter_value = table_start, Transfer_Table = table_final, p_final = p_final, chisq_final = chisq_final, needtworows = allnotenough, taylor_pred = taylor_pred, perc_bias_pred = perc_bias_pred, final_extra = final_extra, - dcroddsratio_ob = dcroddsratio_ob, total_switch = total_switch, - isinvalidate_ob = isinvalidate_ob) + dcroddsratio_ob = dcroddsratio_ob, total_switch = total_switch, isinvalidate_ob = isinvalidate_ob) return(result) } -getswitch_fisher <- function(a, b, c, d, thr_p = 0.05, switch_trm = TRUE){ +getswitch_fisher <- function(a, b, c, d, thr_p = 0.05, switch_trm = T){ if (a > 0 & b > 0 & c > 0 & d > 0){ odds_ratio <- fisher_oddsratio(a, b, c, d) } else { diff --git a/R/pkonfound.R b/R/pkonfound.R index d18e5045..ec6e2ac5 100644 --- a/R/pkonfound.R +++ b/R/pkonfound.R @@ -89,8 +89,7 @@ FR2max, FR2max_multiplier = 1.3, to_return = "print") { - if ("table" %in% to_return) stop("a table can only be - output when using konfound") + if ("table" %in% to_return) stop("a table can only be output when using konfound") if (index == "COP") { @@ -138,8 +137,7 @@ } else if(!is.null(a)) { # error handling if (is.null(a) | is.null(b) | is.null(c) | is.null(d)) { - stop("Please enter values for a, b, c, - and d to use the 2 x 2 table functionality") + stop("Please enter values for a, b, c, and d to use the 2 x 2 table functionality") } out <- tkonfound(a = a, @@ -190,8 +188,7 @@ if (!is.null(out)) { # dealing with a strange print issue } if (to_return == "print") { - message("For other forms of output, run - ?pkonfound and inspect the to_return argument") + message("For other forms of output, run ?pkonfound and inspect the to_return argument") } message("For models fit in R, consider use of konfound().") diff --git a/R/test_cop.R b/R/test_cop.R index d9b18b69..835d70b7 100644 --- a/R/test_cop.R +++ b/R/test_cop.R @@ -34,29 +34,19 @@ test_cop <- function(est_eff, # unstandardized var_y <- sdy^2 ### if the user specifies R2max directly then we use the specified R2max - if (FR2max == 0){FR2max <- FR2max_multiplier * R2} + if (FR2max == 0){FR2max = FR2max_multiplier * R2} var_z <- sdz <- 1 ## error message if input is inappropriate - if (!(std_err > 0)) {stop("Did not run! Standard error needs to be - greater than zero.")} - if (!(sdx > 0)) {stop("Did not run! Standard deviation of x needs to be - greater than zero.")} - if (!(sdy > 0)) {stop("Did not run! Standard deviation of y needs to be - greater than zero.")} - if (!(n_obs > n_covariates + 3)) { - stop("Did not run! There are too few observations relative to the - number of observations and covariates. Please specify a less - complex model to use KonFound-It.")} + if (!(std_err > 0)) {stop("Did not run! Standard error needs to be greater than zero.")} + if (!(sdx > 0)) {stop("Did not run! Standard deviation of x needs to be greater than zero.")} + if (!(sdy > 0)) {stop("Did not run! Standard deviation of y needs to be greater than zero.")} + if (!(n_obs > n_covariates + 3)) {stop("Did not run! There are too few observations relative to the number of observations and covariates. Please specify a less complex model to use KonFound-It.")} if (!(R2 < FR2max)) {stop("Did not run! R2 Max needs to be greater than R2.")} if (!(FR2max < 1)) {stop("Did not run! R2 Max needs to be less than 1.")} - if (!(1-((sdy^2/sdx^2)*(1-R2)/((df+1) * std_err^2))>0)) { - stop("Did not run! Entered values produced Rxz^2 <=0, - consider adding more significant digits to your entered values.")} + if (!(1-((sdy^2/sdx^2)*(1-R2)/((df+1) * std_err^2))>0)) {stop("Did not run! Entered values produced Rxz^2 <=0, consider adding more significant digits to your entered values.")} - # an indicator of whether the use specified est_eff is negative, - # 1 is yes negative - negest <- 0 + negest <- 0 # an indicator of whether the use specified est_eff is negative, 1 is yes negative if (est_eff < 0) { est_eff <- abs(est_eff) negest <- 1 @@ -75,7 +65,7 @@ test_cop <- function(est_eff, # unstandardized ## ryxGz_M2 is only for simulation to recover the exact number ## make sure R2 due to x alone is not larger than overall or observed R2 - if (ryxGz^2 > R2) {illcnd_ryxGz <- TRUE} else {illcond_ryxGz <- FALSE} + if (ryxGz^2 > R2) {illcnd_ryxGz = TRUE} else {illcond_ryxGz = FALSE} ## calculate ryz, rxz, rxy ryz <- rzy <- cal_ryz(ryxGz, R2) @@ -83,8 +73,7 @@ test_cop <- function(est_eff, # unstandardized ## df + 1 because omitted variable is NOT included in M2 #### we change the n_obs to df to recover the rxz as in the particular sample - ## note that in the updated approach rxy is not necessary to calculate - ## rxcv_exact, ryxcv_exact and delta_exact + ## note that in the updated approach rxy is not necessary to calculate rxcv_exact, ryxcv_exact and delta_exact rxy <- ryx <- cal_rxy(ryxGz, rxz, ryz) rxy_M2 <- cal_rxy(ryxGz_M2, rxz, ryz) ## rxy_M2 is only for simulation to recover the exact number @@ -97,9 +86,8 @@ test_cop <- function(est_eff, # unstandardized R2_uncond <- rxy^2 ## calculate delta_star - delta_star <- cal_delta_star(FR2max, R2, R2_uncond, est_eff, - eff_thr, var_x, var_y, eff_uncond, - rxz, n_obs) + delta_star <- cal_delta_star(FR2max, R2, R2_uncond, est_eff, eff_thr, var_x, var_y, eff_uncond, rxz, n_obs) + ## now introduce cv sdcv <- var_cv <- 1 rcvz <- rzcv <- 0 @@ -120,11 +108,9 @@ test_cop <- function(est_eff, # unstandardized # prepare some other values in the final Table (long output) R2_M3_oster <- as.numeric(verify_oster[1]) - eff_x_M3_oster <- as.numeric(verify_oster[2]) - # should be equivalent or very close to eff_thr + eff_x_M3_oster <- as.numeric(verify_oster[2]) # should be equivalent or very close to eff_thr se_x_M3_oster <- as.numeric(verify_oster[3]) - beta_x_M3_oster <- as.numeric(verify_oster[9]) - # should be equivalent or very close to beta_thr + beta_x_M3_oster <- as.numeric(verify_oster[9]) # should be equivalent or very close to beta_thr t_x_M3_oster <- eff_x_M3_oster / se_x_M3_oster eff_z_M3_oster <- as.numeric(verify_oster[4]) se_z_M3_oster <- as.numeric(verify_oster[5]) @@ -133,8 +119,7 @@ test_cop <- function(est_eff, # unstandardized cov_oster <- verify_oster[[11]] cor_oster <- verify_oster[[12]] - ## calculate the exact/true rcvx, rcvy, delta - ## (updated version that does not need rxy) + ## calculate the exact/true rcvx, rcvy, delta (updated version that does not need rxy) ### the idea is to calculate everything conditional on z sdxGz <- sdx * sqrt(1 - rxz^2) sdyGz <- sdy * sqrt(1 - ryz^2) @@ -171,11 +156,9 @@ test_cop <- function(est_eff, # unstandardized # prepare some other values in the final Table (long output) R2_M3 <- as.numeric(verify_exact[1]) - eff_x_M3 <- as.numeric(verify_exact[2]) - # should be equivalent or very close to eff_thr + eff_x_M3 <- as.numeric(verify_exact[2]) # should be equivalent or very close to eff_thr se_x_M3 <- as.numeric(verify_exact[3]) - beta_x_M3 <- as.numeric(verify_exact[9]) - # should be equivalent or very close to beta_thr + beta_x_M3 <- as.numeric(verify_exact[9]) # should be equivalent or very close to beta_thr t_x_M3 <- eff_x_M3 / se_x_M3 eff_z_M3 <- as.numeric(verify_exact[4]) se_z_M3 <- as.numeric(verify_exact[5]) @@ -186,10 +169,8 @@ test_cop <- function(est_eff, # unstandardized verify_pse_reg_M2 <- verify_reg_Gz(n_obs, sdx, sdy, sdz, rxy_M2, rxz, rzy) R2_M2 <- as.numeric(verify_pse_reg_M2[1]) - eff_x_M2 <- as.numeric(verify_pse_reg_M2[2]) - # should be equivalent or very close to est_eff - se_x_M2 <- as.numeric(verify_pse_reg_M2[3]) - # should be equivalent or very close to std_err + eff_x_M2 <- as.numeric(verify_pse_reg_M2[2]) # should be equivalent or very close to est_eff + se_x_M2 <- as.numeric(verify_pse_reg_M2[3]) # should be equivalent or very close to std_err eff_z_M2 <- as.numeric(verify_pse_reg_M2[4]) se_z_M2 <- as.numeric(verify_pse_reg_M2[5]) t_x_M2 <- eff_x_M2 / se_x_M2 @@ -214,7 +195,7 @@ test_cop <- function(est_eff, # unstandardized NA, NA, eff_cv_M3, eff_cv_M3_oster, # reg coef for CV in three reg models NA, NA, se_cv_M3, se_cv_M3_oster, # se for CV in three reg models NA, NA, eff_cv_M3 / se_cv_M3, eff_cv_M3_oster / se_cv_M3_oster), # t for CV in three reg models - nrow = 8, ncol = 4, byrow = TRUE) + nrow = 8, ncol = 4, byrow = T) rownames(fTable) <- c("R2", "coef_X", "SE_X", "std_coef_X", "t_X", # "coef_Z", "SE_Z", "t_Z", @@ -229,8 +210,7 @@ test_cop <- function(est_eff, # unstandardized "Intermediate(M2)", eff_x_M2, R2, "exact", "Final(M3)", eff_x_M3, FR2max, "exact", "Intermediate(M2)", eff_x_M2, R2, "star", - "Final(M3)", eff_x_M3_oster, FR2max, "star"), - nrow = 5, ncol = 4, byrow =TRUE) + "Final(M3)", eff_x_M3_oster, FR2max, "star"), nrow = 5, ncol = 4, byrow =TRUE) colnames(figTable) <- c("ModelLabel", "coef_X", "R2", "cat") figTable <- as.data.frame(figTable) @@ -301,8 +281,7 @@ fig <- ggplot2::ggplot(figTable, ggplot2::aes(x = figTable$ModelLabel)) + if (to_return == "raw_output") { if (negest == 1) { - cat("Using the absolute value of the estimated effect, - results can be interpreted by symmetry.") + cat("Using the absolute value of the estimated effect, results can be interpreted by symmetry.") cat("\n") } output <- list("delta*" = delta_star, @@ -332,29 +311,22 @@ fig <- ggplot2::ggplot(figTable, ggplot2::aes(x = figTable$ModelLabel)) + cat("This function calculates delta* and the exact value of delta.") cat("\n") if (negest == 1) { - cat("Using the absolute value of the estimated effect, - results can be interpreted by symmetry.") + cat("Using the absolute value of the estimated effect, results can be interpreted by symmetry.") cat("\n") } - cat(sprintf("delta* is %.3f (assuming no covariates in the baseline model M1), - the exact delta is %.3f, with a bias of %.3f%%.", - delta_star, delta_exact, delta_pctbias)) + cat(sprintf("delta* is %.3f (assuming no covariates in the baseline model M1), the exact delta is %.3f, with a bias of %.3f%%.", delta_star, delta_exact, delta_pctbias)) cat("\n") - cat(sprintf("With delta*, the coefficient in the final model will be %.3f. - With the exact delta, the coefficient will be %.3f.", + cat(sprintf("With delta*, the coefficient in the final model will be %.3f. With the exact delta, the coefficient will be %.3f.", eff_x_M3_oster,eff_x_M3)) cat("\n") - cat("Use to_return = raw_ouput to see more specific results - and graphic presentation of the result.") + cat("Use to_return = raw_ouput to see more specific results and graphic presentation of the result.") cat("\n") cat("\n") - cat("This function also calculates conditional RIR that - invalidates the statistical inference.") + cat("This function also calculates conditional RIR that invalidates the statistical inference.") cat("\n") cat("If the replacement cases have a fixed value, then RIR =", cond_RIR_fixedY, ".") cat("\n") - cat("If the replacement cases follow a null distribution, - then RIR =", cond_RIR_null, ".") + cat("If the replacement cases follow a null distribution, then RIR =", cond_RIR_null, ".") cat("\n") cat("If the replacement cases satisfy rxy|Z = 0, then RIR =", cond_RIR_rxyz, ".") cat("\n") diff --git a/R/test_pse.R b/R/test_pse.R index 76ca3e0f..a1502257 100644 --- a/R/test_pse.R +++ b/R/test_pse.R @@ -14,24 +14,17 @@ test_pse <- function(est_eff, var_x <- sdx^2 var_y <- sdy^2 var_z <- sdz <- 1 - df <- n_obs - n_covariates - 3 + df = n_obs - n_covariates - 3 ## error message if input is inappropriate - if (!(std_err > 0)) {stop("Did not run! Standard error needs - to be greater than zero.")} - if (!(sdx > 0)) {stop("Did not run! Standard deviation of - x needs to be greater than zero.")} - if (!(sdy > 0)) {stop("Did not run! Standard deviation of - y needs to be greater than zero.")} - if (!(n_obs > n_covariates + 3)) { - stop("Did not run! There are too few observations relative - to the number of observations and covariates. - Please specify a less complex model to use KonFound-It.")} + if (!(std_err > 0)) {stop("Did not run! Standard error needs to be greater than zero.")} + if (!(sdx > 0)) {stop("Did not run! Standard deviation of x needs to be greater than zero.")} + if (!(sdy > 0)) {stop("Did not run! Standard deviation of y needs to be greater than zero.")} + if (!(n_obs > n_covariates + 3)) {stop("Did not run! There are too few observations relative to the number of observations and covariates. Please specify a less complex model to use KonFound-It.")} if (!(0 < R2)) {stop("Did not run! R2 needs to be greater than zero.")} if (!(R2 < 1)) {stop("Did not run! R2 needs to be less than one.")} - if (!(1-((sdy^2/sdx^2)*(1-R2)/((df+1) * std_err^2))>0)) { - stop("Did not run! Entered values produced Rxz^2 <=0, - consider adding more significant digits to your entered values.")} + if (!(1-((sdy^2/sdx^2)*(1-R2)/((df+1) * std_err^2))>0)) {stop("Did not run! Entered values produced Rxz^2 <=0, consider adding more significant digits to your entered values.")} + ## now standardize beta_thr <- eff_thr * sdx / sdy @@ -58,27 +51,19 @@ test_pse <- function(est_eff, rxcvGz <- as.numeric(Gz_pse[1]) rycvGz <- as.numeric(Gz_pse[2]) - # convert conditional correlations to unconditional correlations - # to be used in new regression + # convert conditional correlations to unconditional correlations to be used in new regression rxcv <- rxcvGz * sqrt((1 - rcvz^2) * (1 - rxz^2)) + rxz * rcvz rycv <- rycvGz * sqrt((1 - rcvz^2) * (1 - rzy^2)) + rzy * rcvz - verify_pse_reg_M3 <- verify_reg_Gzcv(n_obs, sdx, sdy, - sdz, sdcv, rxy, - rxz, rzy, rycv, - rxcv, rcvz) - verfiy_pse_manual_thr <- verify_manual(rxy, rxz, rxcv, - ryz, rycv, rzcv, - sdy, sdx) + verify_pse_reg_M3 <- verify_reg_Gzcv(n_obs, sdx, sdy, sdz, sdcv, rxy, rxz, rzy, rycv, rxcv, rcvz) + verfiy_pse_manual_thr <- verify_manual(rxy, rxz, rxcv, ryz, rycv, rzcv, sdy, sdx) cov_pse <- verify_pse_reg_M3[[11]] # prepare some other values in the final Table (long output) R2_M3 <- as.numeric(verify_pse_reg_M3[1]) - eff_x_M3 <- as.numeric(verify_pse_reg_M3[2]) - # should be equivalent or very close to eff_thr + eff_x_M3 <- as.numeric(verify_pse_reg_M3[2]) # should be equivalent or very close to eff_thr se_x_M3 <- as.numeric(verify_pse_reg_M3[3]) - beta_x_M3 <- as.numeric(verify_pse_reg_M3[9]) - # should be equivalent or very close to thr + beta_x_M3 <- as.numeric(verify_pse_reg_M3[9]) # should be equivalent or very close to thr t_x_M3 <- eff_x_M3 / se_x_M3 eff_z_M3 <- as.numeric(verify_pse_reg_M3[4]) se_z_M3 <- as.numeric(verify_pse_reg_M3[5]) @@ -87,44 +72,30 @@ test_pse <- function(est_eff, verify_pse_reg_M2 <- verify_reg_Gz(n_obs, sdx, sdy, sdz, rxy, rxz, rzy) R2_M2 <- as.numeric(verify_pse_reg_M2[1]) - eff_x_M2 <- as.numeric(verify_pse_reg_M2[2]) - # should be equivalent or very close to est_eff - se_x_M2 <- as.numeric(verify_pse_reg_M2[3]) - # should be equivalent or very close to std_err + eff_x_M2 <- as.numeric(verify_pse_reg_M2[2]) # should be equivalent or very close to est_eff + se_x_M2 <- as.numeric(verify_pse_reg_M2[3]) # should be equivalent or very close to std_err eff_z_M2 <- as.numeric(verify_pse_reg_M2[4]) se_z_M2 <- as.numeric(verify_pse_reg_M2[5]) t_x_M2 <- eff_x_M2 / se_x_M2 verify_pse_reg_M1 <- verify_reg_uncond(n_obs, sdx, sdy, rxy) - R2_M1 <- as.numeric(verify_pse_reg_M1[1]) - # should be equivalent or very close to rxy^2 - eff_x_M1 <- as.numeric(verify_pse_reg_M1[2]) - # should be equivalent or very close to rxy*sdy/sdx + R2_M1 <- as.numeric(verify_pse_reg_M1[1]) # should be equivalent or very close to rxy^2 + eff_x_M1 <- as.numeric(verify_pse_reg_M1[2]) # should be equivalent or very close to rxy*sdy/sdx se_x_M1 <- as.numeric(verify_pse_reg_M1[3]) t_x_M1 <- eff_x_M1 / se_x_M1 fTable <- matrix(c(R2_M1, R2_M2, R2_M3, # R2 for three reg models - eff_x_M1, eff_x_M2, eff_x_M3, - # unstd reg coef for X in three reg models - se_x_M1, se_x_M2, se_x_M3, - # unstd reg se for X in three reg models - rxy, ryxGz, beta_x_M3, - # std reg coef for X in three reg models - t_x_M1, t_x_M2, t_x_M3, - # t values for X in three reg models - NA, eff_z_M2, eff_z_M3, - # reg coef for Z in three reg models - NA, se_z_M2, se_z_M3, - # se for Z in three reg models - NA, eff_z_M2 / se_z_M2, eff_z_M3 / se_z_M3, - # t for Z in three reg models, - NA, NA, eff_cv_M3, - # reg coef for CV in three reg models - NA, NA, se_cv_M3, - # se for CV in three reg models - NA, NA, eff_cv_M3 / se_cv_M3), - # t for CV in three reg models - nrow = 11, ncol = 3, byrow = TRUE) + eff_x_M1, eff_x_M2, eff_x_M3, # unstd reg coef for X in three reg models + se_x_M1, se_x_M2, se_x_M3, # unstd reg se for X in three reg models + rxy, ryxGz, beta_x_M3, # std reg coef for X in three reg models + t_x_M1, t_x_M2, t_x_M3, # t values for X in three reg models + NA, eff_z_M2, eff_z_M3, # reg coef for Z in three reg models + NA, se_z_M2, se_z_M3, # se for Z in three reg models + NA, eff_z_M2 / se_z_M2, eff_z_M3 / se_z_M3, # t for Z in three reg models, + NA, NA, eff_cv_M3, # reg coef for CV in three reg models + NA, NA, se_cv_M3, # se for CV in three reg models + NA, NA, eff_cv_M3 / se_cv_M3), # t for CV in three reg models + nrow = 11, ncol = 3, byrow = T) rownames(fTable) <- c("R2", "coef_X", "SE_X", "std_coef_X", "t_X", "coef_Z", "SE_Z", "t_Z", @@ -143,19 +114,13 @@ test_pse <- function(est_eff, } if (to_return == "print") { - cat("This function calculates the conditions that set the - estimated effect approximately equal to the threshold - while preserving the standard error.") + cat("This function calculates the conditions that set the estimated effect approximately equal to the threshold while preserving the standard error.") cat("\n") - cat(sprintf("The correlation between X and CV is %.3f, - and the correlation between Y and CV is %.3f.", rxcv, rycv)) + cat(sprintf("The correlation between X and CV is %.3f, and the correlation between Y and CV is %.3f.", rxcv, rycv)) cat("\n") - cat(sprintf("Conditional on the covariates, the correlation between X - and CV is %.3f, and the correlation between Y and CV - is %.3f.", rxcvGz, rycvGz)) + cat(sprintf("Conditional on the covariates, the correlation between X and CV is %.3f, and the correlation between Y and CV is %.3f.", rxcvGz, rycvGz)) cat("\n") - cat(sprintf("Including such CV, the coefficient changes to %.3f, - and standard error is %.3f.", eff_x_M3, se_x_M3)) + cat(sprintf("Including such CV, the coefficient changes to %.3f, and standard error is %.3f.", eff_x_M3, se_x_M3)) cat("\n") cat("Use to_return = raw_ouput to see more specific results.") } diff --git a/R/test_sensitivity.R b/R/test_sensitivity.R index 8eca8db5..184a6ddb 100644 --- a/R/test_sensitivity.R +++ b/R/test_sensitivity.R @@ -34,17 +34,11 @@ test_sensitivity <- function(est_eff, to_return, model_object, tested_variable) { - if (nu != 0) warning("You entered a non-zero null hypothesis about an effect; - this is being interpreted in terms of a partial correlation. - Sampling variability is not accounted for.") + if (nu != 0) warning("You entered a non-zero null hypothesis about an effect; this is being interpreted in terms of a partial correlation. Sampling variability is not accounted for.") ## error message if input is inappropriate - if (!(std_err > 0)) { - stop("Did not run! Standard error needs to be greater than zero.")} - if (!(n_obs > n_covariates + 3)) { - stop("Did not run! There are too few observations relative to the number - of observations and covariates. Please specify a less complex model - to use KonFound-It.")} + if (!(std_err > 0)) {stop("Did not run! Standard error needs to be greater than zero.")} + if (!(n_obs > n_covariates + 3)) {stop("Did not run! There are too few observations relative to the number of observations and covariates. Please specify a less complex model to use KonFound-It.")} # calculating statistics used in every case if (est_eff < 0) { @@ -94,22 +88,16 @@ test_sensitivity <- function(est_eff, # output dispatch if (to_return == "raw_output") { - return(output_df(est_eff, beta_threshold, est_eff, bias, sustain, - recase, obs_r, critical_r, r_con, itcv)) - } else if (to_return == "thresh_plot") { - # this still makes sense for NLMs (just not quite as accurate) - return(plot_threshold(beta_threshold = beta_threshold, - est_eff = est_eff)) + return(output_df(est_eff, beta_threshold, est_eff, bias, sustain, recase, obs_r, critical_r, r_con, itcv)) + } else if (to_return == "thresh_plot") { # this still makes sense for NLMs (just not quite as accurate) + return(plot_threshold(beta_threshold = beta_threshold, est_eff = est_eff)) } else if (to_return == "corr_plot") { - return(plot_correlation(r_con = r_con, obs_r = obs_r, - critical_r = critical_r)) + return(plot_correlation(r_con = r_con, obs_r = obs_r, critical_r = critical_r)) } else if (to_return == "print") { - return(output_print(est_eff, beta_threshold, bias, sustain, nu, recase, - obs_r, critical_r, r_con, itcv, alpha, index)) + return(output_print(est_eff, beta_threshold, bias, sustain, nu, recase, obs_r, critical_r, r_con, itcv, alpha, index)) } else if (to_return == "table") { return(output_table(model_object, tested_variable)) } else { - stop("to_return must be set to 'raw_output', 'print', 'table', - 'thresh_plot', or 'corr_plot' or some combination thereof") + stop("to_return must be set to 'raw_output', 'print', 'table', 'thresh_plot', or 'corr_plot' or some combination thereof") } } diff --git a/R/test_sensitivity_ln.R b/R/test_sensitivity_ln.R index bd3caf87..13fe7756 100644 --- a/R/test_sensitivity_ln.R +++ b/R/test_sensitivity_ln.R @@ -1,5 +1,4 @@ -# Main function to test sensitivity for non-linear models to be wrapped -# with pkonfound(), konfound(), and mkonfound() +# Main function to test sensitivity for non-linear models to be wrapped with pkonfound(), konfound(), and mkonfound() test_sensitivity_ln <- function(est_eff, std_err, @@ -16,23 +15,17 @@ test_sensitivity_ln <- function(est_eff, tested_variable) { ## error message if input is inappropriate - if (!(std_err > 0)) { - stop("Did not run! Standard error needs to be greater than zero.")} - if (!(n_obs > n_covariates + 3)) { - stop("Did not run! There are too few observations relative to - the number of observations and covariates. Please specify a - less complex model to use KonFound-It.")} + if (!(std_err > 0)) {stop("Did not run! Standard error needs to be greater than zero.")} + if (!(n_obs > n_covariates + 3)) {stop("Did not run! There are too few observations relative to the number of observations and covariates. Please specify a less complex model to use KonFound-It.")} thr_t <- cal_thr_t(est_eff, alpha, tails, n_obs, n_covariates) # stop message if (n_obs <= 0 || n_treat <= 0) { - stop("Please enter positive integers for sample size - and number of treatment group cases.") + stop("Please enter positive integers for sample size and number of treatment group cases.") } if (n_obs <= n_treat) { - stop("The total sample size should be larger than - the number of treatment group cases.") + stop("The total sample size should be larger than the number of treatment group cases.") } odds_ratio <- exp(est_eff) @@ -63,8 +56,8 @@ test_sensitivity_ln <- function(est_eff, # to record the original treatment cases in case we need to adjust it # user_ntrm <- n_treat # check if the implied table solution may contain imaginary numbers - # haveimaginary <- FALSE - # changepi <- FALSE + # haveimaginary <- F + # changepi <- F # set the default value for whether we need and can adjust pi (ratio of treatment cases) # to remove the imaginary part # keyimagin <- (4 + 4 * odds_ratio^2 + odds_ratio * @@ -72,9 +65,9 @@ test_sensitivity_ln <- function(est_eff, # minimgain <- 4 + 4 * odds_ratio^2 + odds_ratio * (-8 + n_obs * std_err^2 * (4 - 0.25 * n_obs * std_err^2)) # keyx1 <- 4 + 4 * odds_ratio^2 + odds_ratio * (-8 + 4 * n_obs * std_err^2) # if (keyimagin > 0) { - # haveimaginary <- TRUE + # haveimaginary <- T # if (minimgain <= 0 && keyx1 > 0) { - # changepi <- TRUE + # changepi <- T # n_treat <- n_obs * get_pi(odds_ratio, std_err, n_obs, n_treat) # n_cnt <- n_obs - n_treat # } else { @@ -204,11 +197,9 @@ test_sensitivity_ln <- function(est_eff, change, sprintf("you would need to replace %d", RIR), RIRway, "cases") if (replace == "control") { - conclusion1a <- sprintf("with cases for which the probability of - failure in the control group applies (RIR = %d).", RIR) + conclusion1a <- sprintf("with cases for which the probability of failure in the control group applies (RIR = %d).", RIR) } else { - conclusion1a <- sprintf("with cases for which the probability of failure - in the entire sample applies (RIR = %d).", RIR) + conclusion1a <- sprintf("with cases for which the probability of failure in the entire sample applies (RIR = %d).", RIR) } conclusion1b <- paste( @@ -222,11 +213,9 @@ test_sensitivity_ln <- function(est_eff, change, sprintf("you would need to replace %d", RIR), RIRway, "cases") if (replace == "control") { - conclusion1a <- sprintf("with cases for which the probability of - failure in the control group applies (RIR = %d).", RIR) + conclusion1a <- sprintf("with cases for which the probability of failure in the control group applies (RIR = %d).", RIR) } else { - conclusion1a <- sprintf("with cases for which the probability of failure - in the entire sample applies (RIR = %d).", RIR) + conclusion1a <- sprintf("with cases for which the probability of failure in the entire sample applies (RIR = %d).", RIR) } conclusion1b <- paste( @@ -239,27 +228,22 @@ test_sensitivity_ln <- function(est_eff, conclusion1 <- paste( change, c("only transferring cases from"), transferway, "is not enough.") - conclusion1b <- paste(sprintf("We also need to transfer %d cases from", - final_solution$final_extra), - transferway_extra, c("as shown, from the User-entered - Table to the Transfer Table. ")) + conclusion1b <- paste(sprintf("We also need to transfer %d cases from", final_solution$final_extra), + transferway_extra, c("as shown, from the User-entered Table to the Transfer Table. ")) conclusion1c <- paste(sprintf("This means we need to replace %d of", RIR), RIRway, - sprintf("with null hypothesis cases; and replace %d", - RIR_extra), RIRway_extra, + sprintf("with null hypothesis cases; and replace %d", RIR_extra), RIRway_extra, c("with null hypothesis cases to change the inference.") ) } conclusion2 <- sprintf( - "For the Implied Table, we have an estimate of %.3f, - with a SE of %.3f and a t-ratio of %.3f.", + "For the Implied Table, we have an estimate of %.3f, with a SE of %.3f and a t-ratio of %.3f.", final_solution$est_eff_start, final_solution$std_err_start, final_solution$t_start ) conclusion3 <- sprintf( - "For the Transfer Table, we have an estimate of %.3f, - with a SE of %.3f and a t-ratio of %.3f.", + "For the Transfer Table, we have an estimate of %.3f, with a SE of %.3f and a t-ratio of %.3f.", final_solution$est_eff_final, final_solution$std_err_final, final_solution$t_final ) @@ -268,8 +252,7 @@ test_sensitivity_ln <- function(est_eff, if (changeSE) { notice_SE <- sprintf( - "In order to generate a usable implied contingency table, - we increased the standard error to %.3f (the original one is %.3f).", + "In order to generate a usable implied contingency table, we increased the standard error to %.3f (the original one is %.3f).", std_err, user_std_err) } diff --git a/R/tkonfound.R b/R/tkonfound.R index 0fcd1e53..ff9e6db8 100644 --- a/R/tkonfound.R +++ b/R/tkonfound.R @@ -151,12 +151,10 @@ tkonfound <- function(a, b, c, d, if (replace == "control") { conclusion1b <- paste0( - sprintf(" cases for which the probability of failure in - the control group applies (RIR = %d). ", RIR)) + sprintf(" cases for which the probability of failure in the control group applies (RIR = %d). ", RIR)) } else { conclusion1b <- paste0( - sprintf(" cases for which the probability of failure in - the entire group applies (RIR = %d). ", RIR)) + sprintf(" cases for which the probability of failure in the entire group applies (RIR = %d). ", RIR)) } conclusion1c <- paste0( @@ -171,12 +169,10 @@ tkonfound <- function(a, b, c, d, if (replace == "control") { conclusion1b <- paste0( - sprintf(" cases for which the probability of failure in - the control group applies (RIR = %d). ", RIR)) + sprintf(" cases for which the probability of failure in the control group applies (RIR = %d). ", RIR)) } else { conclusion1b <- paste0( - sprintf(" cases for which the probability of failure in - the entire group applies (RIR = %d). ", RIR)) + sprintf(" cases for which the probability of failure in the entire group applies (RIR = %d). ", RIR)) } conclusion1c <- paste0( @@ -187,35 +183,28 @@ tkonfound <- function(a, b, c, d, if (allnotenough){ conclusion1 <- paste( change, c("only transferring cases from" ), transferway, - sprintf(" is not enough. We also need to - transfer %d cases from ", final_extra)) + sprintf(" is not enough. We also need to transfer %d cases from ", final_extra)) conclusion1b <- paste0( - transferway_extra, c("as shown, from the User-entered Table - to the Transfer Table.")) + transferway_extra, c("as shown, from the User-entered Table to the Transfer Table.")) conclusion1c <- paste0(sprintf(" This means we need to replace %d of ", RIR), RIRway, - sprintf( "with null hypothesis cases; - and replace %d ", RIR_extra), RIRway_extra, + sprintf( "with null hypothesis cases; and replace %d ", RIR_extra), RIRway_extra, c(" with null hypothesis cases to change the inference.")) } if (test == "chisq"){ conclusion2 <- sprintf( - "For the User-entered Table, we have a Pearson's chi square of %.3f, - with p-value of %.3f:", chisq_ob, p_ob) + "For the User-entered Table, we have a Pearson's chi square of %.3f, with p-value of %.3f:", chisq_ob, p_ob) conclusion3 <- sprintf( - "For the Transfer Table, we have a Pearson's chi square of %.3f, - with p-value of %.3f:", chisq_final, p_final) + "For the Transfer Table, we have a Pearson's chi square of %.3f, with p-value of %.3f:", chisq_final, p_final) } if (test == "fisher"){ conclusion2 <- sprintf( - "For the User-entered Table, we have an estimated odds ratio of %.3f, - with p-value of %.3f:", fisher_ob, p_ob) + "For the User-entered Table, we have an estimated odds ratio of %.3f, with p-value of %.3f:", fisher_ob, p_ob) conclusion3 <- sprintf( - "For the Transfer Table, we have an estimated odds ratio of %.3f, - with p-value of %.3f:", fisher_final, p_final) + "For the Transfer Table, we have an estimated odds ratio of %.3f, with p-value of %.3f:", fisher_final, p_final) } info1 <- "This function calculates the number of cases that would have to be replaced " diff --git a/R/tkonfound_fig.R b/R/tkonfound_fig.R index a424181b..48aafb51 100644 --- a/R/tkonfound_fig.R +++ b/R/tkonfound_fig.R @@ -24,14 +24,12 @@ #' #' @return Returns two plots showing the effect of hypothetical case switches #' on the effect size in a 2x2 table. -tkonfound_fig <- function(a, b, c, d, thr_p = 0.05, switch_trm = TRUE, - test = "fisher", replace = "control"){ +tkonfound_fig <- function(a, b, c, d, thr_p = 0.05, switch_trm = TRUE, test = "fisher", replace = "control"){ n_obs <- a + b + c + d ###***generate the log odds for each step of switch meta <- data.frame(matrix(ncol = 10, nrow = n_obs-3)) -colnames(meta) <- c("a", "b", "c", "d", "nobs", "switch", - "logodds","cntrl_p","tr_p","pdif") +colnames(meta) <- c("a", "b", "c", "d", "nobs", "switch", "logodds","cntrl_p","tr_p","pdif") if (switch_trm == TRUE) { for (i in 1:(n_obs-3)){ if (i <= a){ @@ -125,8 +123,7 @@ if (test == "chisq"){ if (test == "fisher"){ solution <- getswitch_fisher(a, b, c, d, thr_p, switch_trm) for (i in 1:(n_obs-3)){ - meta$logodds[i] <- log(fisher_oddsratio(meta$a[i], meta$b[i], - meta$c[i], meta$d[i])) + meta$logodds[i] <- log(fisher_oddsratio(meta$a[i], meta$b[i], meta$c[i], meta$d[i])) } } @@ -275,32 +272,24 @@ if (!switch_trm && !dcroddsratio_ob) { } zoom$label <- ifelse(zoom$sigpoint=="positive", - paste("sig pos:RIR=", - zoom[zoom$sigpoint=="positive",]$RIR), - NA) + paste("sig pos:RIR=", zoom[zoom$sigpoint=="positive",]$RIR),NA) zoom$label <- ifelse(zoom$sigpoint=="negative", - paste("sig neg:RIR=", - zoom[zoom$sigpoint=="negative",]$RIR), - zoom$label) + paste("sig neg:RIR=", zoom[zoom$sigpoint=="negative",]$RIR),zoom$label) zoom$label <- ifelse(zoom$sigpoint=="current", paste("current"),zoom$label) fig2 <- ggplot2::ggplot(zoom, ggplot2::aes_string(x="RIR",y="pdif"))+ ggplot2::geom_line(ggplot2::aes_string(y="pdif"), size = 1) + - ggplot2::geom_point(ggplot2::aes_string(y="pdif", - shape = "current", - fill = "sigpoint"), + ggplot2::geom_point(ggplot2::aes_string(y="pdif", shape = "current",fill = "sigpoint"), size = 1)+ ggrepel::geom_label_repel(ggplot2::aes_string(label="label"))+ ggplot2::scale_fill_manual(values=fillcol)+ ggplot2::scale_shape_manual(values=pointshape)+ ggplot2::scale_y_continuous(name="Difference in probability of successful outcome (treatment - control)")+ ggplot2::scale_x_continuous(name="RIR (Fragility)", - breaks= c(zoom$RIR[1], - zoom$RIR[as.integer(length(zoom$RIR)/2)], + breaks= c(zoom$RIR[1], zoom$RIR[as.integer(length(zoom$RIR)/2)], zoom$RIR[as.integer(length(zoom$RIR))]), - labels= c(zoom$xaxis[1], - zoom$xaxis[as.integer(length(zoom$RIR)/2)], + labels= c(zoom$xaxis[1], zoom$xaxis[as.integer(length(zoom$RIR)/2)], zoom$xaxis[as.integer(length(zoom$RIR))])) + ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), @@ -309,13 +298,11 @@ fig2 <- ggplot2::ggplot(zoom, ggplot2::aes_string(x="RIR",y="pdif"))+ legend.position = "none") if (pos_thr_pdif <= max(zoom$pdif) && pos_thr_pdif >= min(zoom$pdif)) { - fig2 <- fig2 + ggplot2::geom_hline(yintercept = pos_thr_pdif, - linetype = "dashed", color="green4", size = 1) + fig2 <- fig2 + ggplot2::geom_hline(yintercept = pos_thr_pdif, linetype = "dashed", color="green4", size = 1) } if (neg_thr_pdif <= max(zoom$pdif) && neg_thr_pdif >= min(zoom$pdif)) { - fig2 <- fig2 + ggplot2::geom_hline(yintercept = neg_thr_pdif, - linetype = "dashed", color="red", size = 1) + fig2 <- fig2 + ggplot2::geom_hline(yintercept = neg_thr_pdif, linetype = "dashed", color="red", size = 1) } ###plot figure 3 RIS% as sample size gets larger, using t statistic as the criterion @@ -344,11 +331,9 @@ if (neg_thr_pdif <= max(zoom$pdif) && neg_thr_pdif >= min(zoom$pdif)) { #} if (switch_trm == TRUE) { - note <- "A bend in line indicates switches from the control - row because the treatment row was exhausted." + note <- "A bend in line indicates switches from the control row because the treatment row was exhausted." } else { - note <- "A bend in line indicates switches from the treatment row - because the control row was exhausted." + note <- "A bend in line indicates switches from the treatment row because the control row was exhausted." } result <- list(fig1, note, fig2) diff --git a/R/zzz.R b/R/zzz.R index ca5a237c..57e584c0 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -7,12 +7,8 @@ #' @aliases zzz #' @rdname zzz #' @importFrom utils globalVariables -if (getRversion() >= "2.15.1") utils::globalVariables(c("inference", "key", - "replace_null_cases", - "percent_bias", "val")) +if (getRversion() >= "2.15.1") utils::globalVariables(c("inference", "key", "replace_null_cases", "percent_bias", "val")) .onAttach <- function(libname, pkgname) { - packageStartupMessage("Sensitivity analysis as described in Frank, -Maroulis, Duong, and Kelcey (2013) and in -Frank (2000).\nFor more information visit http://konfound-it.com.") + packageStartupMessage("Sensitivity analysis as described in Frank, Maroulis, Duong, and Kelcey (2013) and in Frank (2000).\nFor more information visit http://konfound-it.com.") }