Skip to content

Commit

Permalink
Revert "Update cop_pse_auxiliary.R（F - FALSE） (#78)"
Browse files Browse the repository at this point in the history 
This reverts commit 0b39585.
  • Loading branch information
@wwang93
wwang93 authored Feb 19, 2024
1 parent 0b39585 commit 478e031
Show file tree
Hide file tree
Showing 25 changed files with 290 additions and 639 deletions.
4 changes: 3 additions & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,9 @@ Imports:
ggrepel,
pbkrtest,
lifecycle,
glue
cli,
glue,
withr
Suggests:
covr,
devtools,
Expand Down
6 changes: 2 additions & 4 deletions R/concord1.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,6 @@
#'
#' @docType data
#' @name concord1
#' @references Hamilton, Lawrence C. 1983. Saving water:
#' A causal model of household conservation.
#' Sociological Perspectives 26(4):355-374.
#' @references Hamilton, Lawrence C. 1983. Saving water: A causal model of household conservation. Sociological Perspectives 26(4):355-374.
#' @format A data.frame with 496 rows and 10 variables.
NULL
NULL
160 changes: 70 additions & 90 deletions R/cop_pse_auxiliary.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,9 @@
#' @return R2yz value
#' @importFrom lavaan parameterEstimates
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!")
}
Expand All @@ -24,10 +24,10 @@ cal_ryz <- function(ryxGz, R2){
#' @return R2xz value
#' @importFrom lavaan parameterEstimates
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)
}

Expand All @@ -39,7 +39,7 @@ cal_rxz <- function(var_x, var_y, R2, df, std_err){
#' @return rxy value
#' @importFrom lavaan parameterEstimates
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)
}

Expand All @@ -66,33 +66,32 @@ cal_delta_star <- function(FR2max,
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)}
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)}

#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
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)
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
## 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
#obtain delta_star which is Oster's delta
delta_star <- num / den
delta_star = num / den
return(delta_star)
}

Expand All @@ -111,17 +110,17 @@ cal_delta_star <- function(FR2max,
#' @param rcvz correlation coefficient between V and Z
#' @return list of model parameters
#' @importFrom lavaan sem parameterEstimates
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 @@ -138,11 +137,11 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv,
sample.nobs = n_obs)
},
error = function(e){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
},
warning = function(w){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
}
)
Expand All @@ -153,18 +152,12 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv,
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
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
Expand All @@ -182,11 +175,11 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv,
sample.nobs = n_obs)
},
error = function(e){
flag_cor <- FALSE
flag_cor = F
return(flag_cor)
},
warning = function(w){
flag_cor <- FALSE
flag_cor = F
return(flag_cor)
}
)
Expand All @@ -197,22 +190,16 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv,
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
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 ((inherits(flag_cor, "lavaan")) && (inherits(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 All @@ -223,8 +210,7 @@ verify_reg_Gzcv <- function(n_obs, sdx, sdy, sdz, sdcv,

## TO DO
## need to code the other solutions for pse as well
## then see how the different solutions perform in
## different scenarios (run the regression)
## then see how the different solutions perform in different scenarios (run the regression)
cal_pse <- function(thr, kryx){
# calculations for preserving standard error
i1 <- complex(real = 1, imaginary = -sqrt(3))
Expand Down Expand Up @@ -261,13 +247,13 @@ verify_manual <- function(rxy, rxz, rxcv, ryz, rycv, rzcv, sdy, sdx){
return(beta)
}

verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
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 @@ -283,11 +269,11 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
sample.nobs = n_obs)
},
error = function(e){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
},
warning = function(w){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
}
)
Expand All @@ -297,19 +283,15 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
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
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 (inherits(flag_cov, "lavaan")) {
result <- list(R2, betaX, seX, betaZ, seZ)
result = list(R2, betaX, seX, betaZ, seZ)
return(result)
} else {
stop("Error!")
Expand All @@ -324,10 +306,10 @@ verify_reg_Gz <- function(n_obs, sdx, sdy, sdz, rxy, rxz, rzy){
#' @param rxy correlation coefficient between X and Y
#' @return list of model parameters
#' @importFrom lavaan sem parameterEstimates
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 @@ -342,11 +324,11 @@ verify_reg_uncond <- function(n_obs, sdx, sdy, rxy){
sample.nobs = n_obs)
},
error = function(e){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
},
warning = function(w){
flag_cov <- FALSE
flag_cov = F
return(flag_cov)
}
)
Expand All @@ -357,16 +339,14 @@ verify_reg_uncond <- function(n_obs, sdx, sdy, rxy){
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
betaX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$est
seX <- lavaan::parameterEstimates(fit)[lavaan::parameterEstimates(fit)$label == 'beta1',]$se
}

if (inherits(flag_cov, "lavaan")) {
result <- list(R2, betaX, seX)
result = list(R2, betaX, seX)
return(result)
} else {
stop("Error!")
}
}
}
8 changes: 3 additions & 5 deletions R/core-sensitivity-mkonfound.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,11 @@ core_sensitivity_mkonfound <- function(t, df, alpha = .05, tails = 2) {
r_con <- round(sqrt(abs(itcv)), 3)

out <- tibble::tibble(t, df, action, inference, pct_bias, itcv, r_con)
names(out) <- c("t", "df", "action", "inference",
"pct_bias_to_change_inference", "itcv", "r_con")
out$pct_bias_to_change_inference <- round(out$pct_bias_to_change_inference,
3)
names(out) <- c("t", "df", "action", "inference", "pct_bias_to_change_inference", "itcv", "r_con")
out$pct_bias_to_change_inference <- round(out$pct_bias_to_change_inference, 3)
out$itcv <- round(out$itcv, 3)
out$action <- as.character(out$action)
out$inference <- as.character(out$inference)

return(out)
}
}
11 changes: 2 additions & 9 deletions R/helper_output_dataframe.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,8 @@
#' @param non_linear flag for non-linear models
#' @return data frame with model information
#' @importFrom dplyr tibble
output_df <- function(est_eff,
beta_threshhold,
unstd_beta,
bias = NULL,
sustain = NULL,
recase, obs_r,
critical_r,
r_con,
itcv,
output_df <- function(est_eff, beta_threshhold, unstd_beta, bias = NULL,
sustain = NULL, recase, obs_r, critical_r, r_con, itcv,
non_linear) {
if (abs(est_eff) > abs(beta_threshhold)) {
df <- dplyr::tibble(
Expand Down
Loading

0 comments on commit 478e031

Please sign in to comment.