Skip to content

Update cop_pse_auxiliary.R #2

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 #2

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
134 changes: 87 additions & 47 deletions R/cop_pse_auxiliary.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)}
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down