Skip to content

Update cop_pse_auxiliary.R #1

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

Open
wwang93 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Update cop_pse_auxiliary.R #1

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
93 changes: 46 additions & 47 deletions R/cop_pse_auxiliary.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,53 +1,52 @@
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!")
}
return(ryz)
}

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)
}

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)
}

cal_delta_star <- function(FR2max, R2, R2_uncond, est_eff, eff_thr, var_x, var_y, est_uncond, rxz, n_obs){
if (FR2max > .99) {FR2max = .99}
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 > 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

t_x = var_x * (n_obs / (n_obs - 1)) * (1 - rxz^2)
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)
## 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
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)
}

Expand Down Expand Up @@ -92,17 +91,17 @@ cal_delta_star <- function(FR2max, R2, R2_uncond, est_eff, eff_thr, var_x, var_y

# see test_cop for updated approach to calculate delta exact

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,
Expand All @@ -119,11 +118,11 @@ verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv,
sample.nobs = n_obs)
},
error = function(e){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
},
warning = function(w){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
}
)
Expand Down Expand Up @@ -157,11 +156,11 @@ verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv,
sample.nobs = n_obs)
},
error = function(e){
flag_cor = F
flag_cor <- FALSE
return(flag_cor)
},
warning = function(w){
flag_cor = F
flag_cor <- FALSE
return(flag_cor)
}
)
Expand All @@ -181,7 +180,7 @@ verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv,
}

if (class(flag_cor) == "lavaan" && class(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)
Expand Down Expand Up @@ -233,9 +232,9 @@ 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,
Expand All @@ -251,11 +250,11 @@ verify_reg_Gz = function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
sample.nobs = n_obs)
},
error = function(e){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
},
warning = function(w){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
}
)
Expand All @@ -273,17 +272,17 @@ verify_reg_Gz = function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
}

if (class(flag_cov) == "lavaan") {
result = list(R2, betaX, seX, betaZ, seZ)
result <- list(R2, betaX, seX, betaZ, seZ)
return(result)
} else {
stop("Error!")
}
}

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,
Expand All @@ -298,11 +297,11 @@ verify_reg_uncond = function(n_obs, sdx, sdy, rxy){
sample.nobs = n_obs)
},
error = function(e){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
},
warning = function(w){
flag_cov = F
flag_cov <- FALSE
return(flag_cov)
}
)
Expand All @@ -318,7 +317,7 @@ verify_reg_uncond = function(n_obs, sdx, sdy, rxy){
}

if (class(flag_cov) == "lavaan") {
result = list(R2, betaX, seX)
result <- list(R2, betaX, seX)
return(result)
} else {
stop("Error!")
Expand Down