diff --git a/R/cop_pse_auxiliary.R b/R/cop_pse_auxiliary.R index cc958b81..c671072d 100644 --- a/R/cop_pse_auxiliary.R +++ b/R/cop_pse_auxiliary.R @@ -1,7 +1,7 @@ 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!") } @@ -9,45 +9,44 @@ cal_ryz <- function(ryxGz, R2){ } 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) } @@ -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, @@ -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) } ) @@ -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) } ) @@ -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) @@ -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, @@ -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) } ) @@ -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, @@ -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) } ) @@ -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!")