Skip to content

Commit

Permalink
Update test_cop.R
Browse files Browse the repository at this point in the history 
Maintain "=" for function parameter assignments and "<-" for variable assignments.
  • Loading branch information
@wwang93
wwang93 authored Oct 12, 2023
1 parent 8f07241 commit c5444b8
Showing 1 changed file with 82 additions and 82 deletions.
164 changes: 82 additions & 82 deletions R/test_cop.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.")}
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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) +
Expand Down Expand Up @@ -333,4 +333,4 @@ fig <- ggplot2::ggplot(figTable, ggplot2::aes(x = ModelLabel)) +
cat("\n")
}

}
}

0 comments on commit c5444b8

Please sign in to comment.