diff --git a/R/test_cop.R b/R/test_cop.R index 5f94ff95..7bd6529c 100644 --- a/R/test_cop.R +++ b/R/test_cop.R @@ -29,13 +29,13 @@ test_cop <- function(est_eff, # unstandardized # sdx = 2, sdy = 6, R2 = .7, eff_thr = 0, FR2max = .8, to_return = "short") ## prepare input - df = n_obs - n_covariates - 3 ## df of M3 - var_x = sdx^2 - var_y = sdy^2 + df <- n_obs - n_covariates - 3 ## df of M3 + var_x <- sdx^2 + var_y <- sdy^2 ### if the user specifies R2max directly then we use the specified R2max if (FR2max == 0){FR2max = FR2max_multiplier * R2} - var_z = sdz = 1 + 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.")} @@ -53,135 +53,135 @@ test_cop <- function(est_eff, # unstandardized } ## now standardize - beta_thr = eff_thr * sdx / sdy - beta = est_eff * sdx / sdy - SE = std_err * sdx / sdy + beta_thr <- eff_thr * sdx / sdy + beta <- est_eff * sdx / sdy + SE <- std_err * sdx / sdy ## observed regression, reg y on x Given z - tyxGz = beta / SE ### this should be equal to est_eff / std_err - ryxGz = tyxGz / sqrt(df + 1 + tyxGz^2) + tyxGz <- beta / SE ### this should be equal to est_eff / std_err + ryxGz <- tyxGz / sqrt(df + 1 + tyxGz^2) ## df + 1 because omitted variable is NOT included in M2 - ryxGz_M2 = tyxGz / sqrt(n_obs + tyxGz^2) + ryxGz_M2 <- tyxGz / sqrt(n_obs + tyxGz^2) ## 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 = T} else {illcond_ryxGz = F} + if (ryxGz^2 > R2) {illcnd_ryxGz = TRUE} else {illcond_ryxGz = FALSE} ## calculate ryz, rxz, rxy - ryz = rzy = cal_ryz(ryxGz, R2) - rxz = rzx = cal_rxz(var_x, var_y, R2, df + 1, std_err) + ryz <- rzy <- cal_ryz(ryxGz, R2) + rxz <- rzx <- cal_rxz(var_x, var_y, R2, df + 1, std_err) ## 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 - rxy = ryx = cal_rxy(ryxGz, rxz, ryz) - rxy_M2 = cal_rxy(ryxGz_M2, rxz, ryz) + 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 ## baseline regression model, no z (unconditional) - eff_uncond = sqrt((var_y / var_x)) * rxy - t_uncond = rxy * sqrt(n_obs - 2)/sqrt(1 - rxy^2) + eff_uncond <- sqrt((var_y / var_x)) * rxy + t_uncond <- rxy * sqrt(n_obs - 2)/sqrt(1 - rxy^2) ## n_obs - 2 - adjust the df in the M1 - std_err_uncond = eff_uncond / t_uncond - R2_uncond = rxy^2 + std_err_uncond <- eff_uncond / t_uncond + 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 - v = 1 - rxz^2 # this is to simplify calculation later - D = sqrt(FR2max - R2) # same above + sdcv <- var_cv <- 1 + rcvz <- rzcv <- 0 + v <- 1 - rxz^2 # this is to simplify calculation later + D <- sqrt(FR2max - R2) # same above ## calculate rxcv & rycv implied by Oster from delta_star (assumes rcvz=0) - rxcv_oster = rcvx_oster = delta_star * rxz * (sdcv / sdz) + rxcv_oster <- rcvx_oster <- delta_star * rxz * (sdcv / sdz) if (abs(rcvx_oster) <1 && (rcvx_oster^2/v)<1) - {rcvy_oster = rycv_oster = + {rcvy_oster <- rycv_oster <- D * sqrt(1 - (rcvx_oster^2 / v)) + (ryx * rcvx_oster) / (v) - (ryz * rcvx_oster * rxz) / (v)} # Checking beta, R2, se generated by delta_star with a regression - verify_oster = verify_reg_Gzcv(n_obs, sdx, sdy, sdz, sdcv, + verify_oster <- verify_reg_Gzcv(n_obs, sdx, sdy, sdz, sdcv, rxy, rxz, rzy, rycv_oster, rxcv_oster, rcvz) # 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 - 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 - 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]) - eff_cv_M3_oster = as.numeric(verify_oster[6]) - se_cv_M3_oster = as.numeric(verify_oster[7]) - cov_oster = verify_oster[[11]] - cor_oster = verify_oster[[12]] + 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 + 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 + 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]) + eff_cv_M3_oster <- as.numeric(verify_oster[6]) + se_cv_M3_oster <- as.numeric(verify_oster[7]) + cov_oster <- verify_oster[[11]] + cor_oster <- verify_oster[[12]] ## 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) - ryxcvGz_exact_sq = (FR2max - ryz^2) / (1 - ryz^2) + sdxGz <- sdx * sqrt(1 - rxz^2) + sdyGz <- sdy * sqrt(1 - ryz^2) + ryxcvGz_exact_sq <- (FR2max - ryz^2) / (1 - ryz^2) ### equation 7 in the manuscript - rxcvGz_exact = (ryxGz - sdxGz / sdyGz * beta_thr) / + rxcvGz_exact <- (ryxGz - sdxGz / sdyGz * beta_thr) / sqrt((sdxGz^2) / (sdyGz^2) * (beta_thr^2) - 2 * ryxGz * sdxGz / sdyGz * beta_thr + ryxcvGz_exact_sq) ### equation 6 in the manuscript - rycvGz_exact = ryxGz * rxcvGz_exact + + rycvGz_exact <- ryxGz * rxcvGz_exact + sqrt((ryxcvGz_exact_sq - ryxGz^2) * (1 - rxcvGz_exact^2)) ### now get unconditional exact rxcv and rycv - rycv_exact = sqrt(1 - ryz^2) * rycvGz_exact - rxcv_exact = sqrt(1 - rxz^2) * rxcvGz_exact - delta_exact = rxcv_exact / rxz + rycv_exact <- sqrt(1 - ryz^2) * rycvGz_exact + rxcv_exact <- sqrt(1 - rxz^2) * rxcvGz_exact + delta_exact <- rxcv_exact / rxz ## previous approach - comment out, but could find in cop_pse_auxiliary - ## exact_result = cal_delta_exact(ryx, ryz, rxz, beta_thr, FR2max, R2, sdx, sdz) - ## rxcv_exact = rcvx_exact = as.numeric(exact_result[1]) - ## rycv_exact = rcvy_exact = as.numeric(exact_result[2]) - ## delta_exact = as.numeric(exact_result[3]) + ## exact_result <- cal_delta_exact(ryx, ryz, rxz, beta_thr, FR2max, R2, sdx, sdz) + ## rxcv_exact <- rcvx_exact <- as.numeric(exact_result[1]) + ## rycv_exact <- rcvy_exact <- as.numeric(exact_result[2]) + ## delta_exact <- as.numeric(exact_result[3]) # Checking beta, R2, se generated by True/Exact Delta with a regression - verify_exact = verify_reg_Gzcv(n_obs, sdx, sdy, sdz, sdcv, + verify_exact <- verify_reg_Gzcv(n_obs, sdx, sdy, sdz, sdcv, rxy, rxz, rzy, rycv_exact, rxcv_exact, rcvz) ### verify_exact[1] == verify_exact[4] == FR2max ### verify_exact[2] == eff_thr ### verify_exact[5] == beta_thr # calculate % bias in delta comparing oster's delta_star with true delta - delta_pctbias = 100 * (delta_star - delta_exact) / delta_exact + delta_pctbias <- 100 * (delta_star - delta_exact) / delta_exact # 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 - 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 - 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]) - eff_cv_M3 = as.numeric(verify_exact[6]) - se_cv_M3 = as.numeric(verify_exact[7]) - cov_exact = verify_exact[[11]] - cor_exact = verify_exact[[12]] - - 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_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 - se_x_M1 = as.numeric(verify_pse_reg_M1[3]) - t_x_M1 = eff_x_M1 / se_x_M1 - - delta_star_restricted = ((est_eff - eff_thr)/(eff_x_M1 - est_eff))* + R2_M3 <- as.numeric(verify_exact[1]) + 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 + 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]) + eff_cv_M3 <- as.numeric(verify_exact[6]) + se_cv_M3 <- as.numeric(verify_exact[7]) + cov_exact <- verify_exact[[11]] + cor_exact <- verify_exact[[12]] + + 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_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 + se_x_M1 <- as.numeric(verify_pse_reg_M1[3]) + t_x_M1 <- eff_x_M1 / se_x_M1 + + delta_star_restricted <- ((est_eff - eff_thr)/(eff_x_M1 - est_eff))* ((R2_M2 - R2_M1)/(R2_M3 - R2_M2)) fTable <- matrix(c(R2_M1, R2_M2, R2_M3, R2_M3_oster, # R2 for three reg models @@ -210,7 +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 = T) + "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) @@ -226,7 +226,7 @@ test_cop <- function(est_eff, # unstandardized figTable$coef_X <- as.numeric(figTable$coef_X) figTable$R2 <- as.numeric(figTable$R2) -scale = 1/(round(max(figTable$coef_X)*10)/10) +scale <- 1/(round(max(figTable$coef_X)*10)/10) fig <- ggplot2::ggplot(figTable, ggplot2::aes(x = ModelLabel)) + ggplot2::geom_point(ggplot2::aes(y = coef_X, group = cat, shape = cat), color = "blue", size = 3) + @@ -333,4 +333,4 @@ fig <- ggplot2::ggplot(figTable, ggplot2::aes(x = ModelLabel)) + cat("\n") } -} \ No newline at end of file +}