Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update test_cop.R #58

Merged
merged 1 commit into from
Oct 26, 2023
Merged

Update test_cop.R #58

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
164 changes: 82 additions & 82 deletions R/test_cop.R
View file
Open in desktop
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")
}

}
}