From 98361701ebecbeaa87f00145ba222b06104d6906 Mon Sep 17 00:00:00 2001 From: wwang93 <142034236+wwang93@users.noreply.github.com> Date: Sun, 19 Nov 2023 19:47:02 -0500 Subject: [PATCH] Update cop_pse_auxiliary.R Fix those lines of code in cop_pse_auxiliary.R that are longer than 80 characters. Except for this part: cal_pse <- function(thr, kryx) # calculations for preserving standard error ---- (the formula for calculations stays the same) --- R/cop_pse_auxiliary.R | 134 +++++++++++++++++++++++++++--------------- 1 file changed, 87 insertions(+), 47 deletions(-) diff --git a/R/cop_pse_auxiliary.R b/R/cop_pse_auxiliary.R index cc958b81..1adca840 100644 --- a/R/cop_pse_auxiliary.R +++ b/R/cop_pse_auxiliary.R @@ -21,7 +21,16 @@ cal_rxy <- function(ryxGz, 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){ +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 < R2 + inci) {FR2max = R2 + inci} check with Ken what this means if (FR2max > R2) {D = sqrt(FR2max - R2)} @@ -35,7 +44,9 @@ cal_delta_star <- function(FR2max, R2, R2_uncond, est_eff, eff_thr, var_x, 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 + ## 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) @@ -128,19 +139,29 @@ verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv, } ) #if model can be run to verify true delta, then run it can save results - if (class(flag_cov) == "lavaan") { - fit <- lavaan::sem(model, - sample.cov = cov.matrix, - sample.nobs = n_obs) - ## the R2 extracted from summary is NOT right, do the calculation below - R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[4,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se - betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est - seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se - betaCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$est - seCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$se - } + if (class(flag_cov) == "lavaan") { + fit <- lavaan::sem( + model, + sample.cov = cov.matrix, + sample.nobs = n_obs + ) + ## the R2 extracted from summary is NOT right, do the calculation below + R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[4, ]$est) / sdy^2 + betaX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta1", ]$est + seX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta1", ]$se + betaZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta2", ]$est + seZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta2", ]$se + betaCV <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta3", ]$est + seCV <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta3", ]$se +} + + #get regression based on true delta in terms of standardized coefficent cor.matrix <- matrix(c(1,rxy, rzy, rcvy, @@ -167,18 +188,25 @@ verify_reg_Gzcv = function(n_obs, sdx, sdy, sdz, sdcv, ) # if model can be run, then run it - if (class(flag_cor) == "lavaan") { - fit <- lavaan::sem(model, - sample.cov = cor.matrix, - sample.nobs = n_obs) - std_R2 <- 1 - lavaan::parameterEstimates(fit)[4,]$est - std_betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est - std_seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se - std_betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est - std_seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se - std_betaCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$est - std_seCV <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta3',]$se - } + if (class(flag_cor) == "lavaan") { + fit <- lavaan::sem(model, + sample.cov = cor.matrix, + sample.nobs = n_obs) + std_R2 <- 1 - lavaan::parameterEstimates(fit)[4, ]$est + std_betaX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta1', ]$est + std_seX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta1', ]$se + std_betaZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta2', ]$est + std_seZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta2', ]$se + std_betaCV <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta3', ]$est + std_seCV <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta3', ]$se +} + if (class(flag_cor) == "lavaan" && class(flag_cov) == "lavaan") { result = list(R2, betaX, seX, betaZ, seZ, betaCV, seCV, @@ -223,7 +251,8 @@ cal_pse <- function(thr, kryx){ } verify_manual <- function(rxy, rxz, rxcv, ryz, rycv, rzcv, sdy, sdx){ - beta <- (rxy + rycv * rxz * rzcv + ryz * rxcv * rzcv - rxy * rzcv^2 - rycv * rxcv - ryz * rxz) / + beta <- (rxy + rycv * rxz * rzcv + ryz * rxcv * rzcv - + rxy * rzcv^2 - rycv * rxcv - ryz * rxz) / (1 + 2 * rxcv * rzcv * rxz - rxcv^2 - rzcv^2 - rxz^2) eff <- beta * sdy / sdx return(beta) @@ -261,16 +290,22 @@ verify_reg_Gz = function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){ ) #if model can be run to verify true delta, then run it can save results if (class(flag_cov) == "lavaan") { - fit <- lavaan::sem(model, - sample.cov = cov.matrix, - sample.nobs = n_obs) - ## the R2 extracted from summary is NOT right, do the calculation below - R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[3,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se - betaZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$est - seZ <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta2',]$se - } + fit <- lavaan::sem(model, + sample.cov = cov.matrix, + sample.nobs = n_obs + ) + ## the R2 extracted from summary is NOT right, do the calculation below + R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[3, ]$est) / sdy^2 + betaX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta1", ]$est + seX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta1", ]$se + betaZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta2", ]$est + seZ <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == "beta2", ]$se +} + if (class(flag_cov) == "lavaan") { result = list(R2, betaX, seX, betaZ, seZ) @@ -307,15 +342,20 @@ verify_reg_uncond = function(n_obs, sdx, sdy, rxy){ } ) #if model can be run to verify true delta, then run it can save results - if (class(flag_cov) == "lavaan") { - fit <- lavaan::sem(model, - sample.cov = cov.matrix, - sample.nobs = n_obs) - ## the R2 extracted from summary is NOT right, do the calculation below - R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[2,]$est) / sdy^2 - betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est - seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se - } + if (class(flag_cov) == "lavaan") { + fit <- lavaan::sem( + model, + sample.cov = cov.matrix, + sample.nobs = n_obs + ) + ## the R2 extracted from summary is NOT right, do the calculation below + R2 <- (sdy^2 - lavaan::parameterEstimates(fit)[2, ]$est) / sdy^2 + betaX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta1', ]$est + seX <- lavaan::parameterEstimates(fit)[ + lavaan::parameterEstimates(fit)$label == 'beta1', ]$se +} + if (class(flag_cov) == "lavaan") { result = list(R2, betaX, seX)