Misstatements, misperceptions, and mistakes in controlling for covariates in observational research

Roger S. Zoh; Annie Yu
rszoh@iu.edu; yu78@iu.edu

School of Public Health, Epidemiology and Biostatistics

Introduction

  • Essential to reduce bias in observational research using methods like statistical models, matching, and propensity score analysis(Yu et al. 2024).
  • Researchers often inaccurately claim to control for covariates, leading to incorrect conclusions about bias elimination.
  • Discusses 12 common misperceptions and mistakes in covariate control based on extensive experience across various life sciences.
  • Emphasizes the importance of accurate covariate control to reduce bias and improve research quality.

Misperception 1 - Construct Validity

  • Believing an observed variable (e.g., highest educational degree) accurately measures a construct (e.g., socioeconomic status) can be misleading.
  • Construct validity involves accurately measuring what a test or measure is intended to capture; proxy variables often fail to meet this criterion.

Examples of Issues

  • Psychology: Modifying the Patient Health Questionnaire-9 (PHQ-9) can reduce its effectiveness in measuring depression(Bell 1994; Indu et al. 2018).
  • Nutritional Epidemiology: Debate on the validity of ‘food addiction’ as a measurable construct due to limited empirical evidence(Hanley-Cook et al. 2022).

Why This Misperception Occurs

  • Innovations Outpace Testing: New ways to measure complex behaviors emerge faster than the validation of instruments.
  • Reliance on Old Instruments: Established instruments may not be validated for new scenarios.
  • Difficulty in Testing: Ensuring construct validity of new or adapted instruments is time-consuming and tedious.

Avoiding the Misperception

  • Use Established Instruments Fully: Avoid altering instruments unless absolutely necessary.
  • Empirical Testing for Adaptations: Use psychometric analyses like confirmatory factor analysis and test/retest reliability to ensure validity and reliability of new or adapted measures.

Misperception 2 - Measurement error in a covariate only attenuates associations or effect estimates and does not create apparent effects

  • Measurement Errors(Fuller 1987; Carroll et al. 2006):
    • Errors can be random, correlated with true values, or correlated with other variables/errors.
    • The distribution and variance of measurement errors significantly impact model bias.
    • Measurement errors can reduce, increase, or have no impact on bias in model estimation.
    • Common belief that classical measurement errors only attenuate effects is incorrect, especially in complex models.
    • Often ignored due to lack of awareness or insufficient information on error variance.
    • Essential to include measurement error in analyses to accurately control for bias.

Why This Misperception Occurs

  • Model-Dependent Influence: Measurement error effects vary by regression model; it can either attenuate or inflate covariate effects.
  • Nature of Random Measurement Error: Random errors are independent and add to true variance; high error variance can turn covariates into random noise, failing to control for potentially biasing covariates (PBCs).
  • Non-Random Errors: Errors can be correlated with important factors in the model. Such correlations make it unpredictable whether including the covariate will reduce or increase bias without understanding the error structure.

Avoiding the Misperception

  • Use measurement error correction methods (Fuller 2009).
  • Strive for better measurements to minimize error(Westfall and Yarkoni 2016).
  • Consider the error structure relative to the model to understand potential biases.

Misperception 3a - Continuous covariates divided into polychotomous categories for better interpretation are still well-controlled

Why This Misperception Occurs:

  • Investigators recognize nonlinearity and break covariates into categories for simpler interpretation(Blas Achic et al. 2018).
  • Using quantile-defined categories introduces measurement error unless the true relationship is a step function at cut points.
  • Coarse categorization increases residual bias and fails to eliminate bias from PBCs.

Avoiding the Misperception:

  • Researchers should avoid dichotomizing continuous covariates to prevent suboptimal analysis and unnecessary measurement error.

Misperception 3b - Covariates categorized in coarse rather than fine categories are more reliable in the presence of measurement error

  • Coarse categorization of covariates is believed to improve reliability in the presence of measurement error.

Why This Misperception Occurs:

  • Belief that coarse categories make data more reliable when fine-grained distinctions are unsupported due to measurement error.

  • Reality: Dichotomization or coarse categorization reduces measurement reliability and attenuates correlations, leading to negative statistical consequences(MacCallum et al. 2002a).

Avoiding the Misperception:

  • Refrain from dichotomizing low-reliability covariates as it negatively impacts analysis(MacCallum et al. 2002b).

  • Improving reliability through categorization is rarely supported by classical measurement theory and is difficult to verify in real application.

Misperception 4 - Controlling for a covariate reduces the power to detect an association of the IV of interest with the DV of interest

Why This Misperception Occurs:

  • Investigators think that including covariates will reduce power, especially if the covariate’s unadjusted correlation with the IV is zero.

  • Reality: Controlling for covariates that are strongly correlated with the DV but uncorrelated with the IV reduces residual variance, increasing power by lowering the denominator of test statistics.

Avoiding the Misperception:

  • Refrain from dichotomizing low-reliability covariates as it negatively impacts analysis.

  • Improving reliability through categorization is rarely supported by classical measurement theory and difficult to verify in practice.

Misperception 5a - If when controlling for X and Z simultaneously in a statistical model as predictors of an outcome Y, X is significant with Z in the model, but Z is not significant with X in the model, then X is a ‘better’ predictor than Z

Why This Misperception Occurs:

  • Investigators incorrectly conclude X has a true causal effect and Z does not, or that X has a stronger causal effect than Z.
  • Example: Muscle size vs. muscle strength on longevity (Attia, 2022). Conclusions drawn from statistical significance without considering collinearity, measurement error, and model fit.
  • Strong multicollinearity between predictors can lead to misleading conclusions.

Avoiding the Misperception:

  • Ensure variables in the model are not too correlated to avoid multicollinearity.
  • Use Variance Inflation Factor (VIF) tests to diagnose multicollinearity.
  • Consider measurement error correction if low reliability is suspected.

Sample R code for mispercetion 5a

X and Z are independent (n=500)

set.seed(123)
n <- 1000
n_simulations <- 500

# Initialize vectors to store coefficients
beta_X_values <- numeric(n_simulations)
beta_Z_values <- numeric(n_simulations)
beta_Z_value_model2 <- numeric(n_simulations)

for (i in 1:n_simulations) {
  # covariance matrix between X and Z
  sigma1 <- matrix(data = c(1, 0, 0, 1), nrow = 2, ncol = 2)
  mu <- c(0, 0)
  beta <- matrix(data = c(0.2, 0.6, 0.7), nrow = 1, ncol = 3)
  d1 <- matrix(data = c(1), nrow = n, ncol = 1)
  XZ1 <- mvrnorm(n, mu, sigma1)
  XZ1 <- cbind(d1, XZ1)
  e1 <- rnorm(n, 0, 0.5)
  
  # Predict Y value
  Y <- XZ1 %*% t(beta) + e1
  
  df1 <- cbind(XZ1, Y) %>% as.data.frame()
  colnames(df1) <- c('x0', 'X', 'Z', 'Y')
  
  # Fit the models
  model1 <- lm(Y ~ X + Z, data = df1)
  model2 <- lm(Y ~ Z, data = df1)
  
  # coefficients for beta_X and beta_Z
  coeff1 <- model1$coefficients
  coeff2 <- model2$coefficients
  beta_X_values[i] <- coeff1['X']
  beta_Z_values[i] <- coeff1['Z']
  beta_Z_value_model2[i]<- coeff2['Z']
  
}

# Calculate the mean values for beta_X and beta_Z in each model
mean_beta_X <- mean(beta_X_values)
mean_beta_Z <- mean(beta_Z_values)

mean_beta_Z_model2 <- mean(beta_Z_value_model2)

# Print the mean values
cat("Mean value for beta_X in model 1:", mean_beta_X, "\n")
cat("Mean value for beta_Z in model 1:", mean_beta_Z, "\n")
cat("Mean value for beta_Z in model 2:", mean_beta_Z_model2, "\n")

mean_values_df1 <- data.frame(
  Variable = c("beta_X in model 1", "beta_Z in model 1", "beta_Z in model 2"),
  Mean_Value = c(mean_beta_X, mean_beta_Z, mean_beta_Z_model2)
)

X and Z are correlate (n=500)

set.seed(123)
n <- 1000
n_simulations <- 500

# Initialize vectors to store coefficients
beta_X_values <- numeric(n_simulations)
beta_Z_values <- numeric(n_simulations)
beta_Z_value_model2 <- numeric(n_simulations)

for (i in 1:n_simulations) {
  # covariance matrix between X and Z (rho=0.8)
  sigma1 <- matrix(data = c(1, 0.8, 0.8, 1), nrow = 2, ncol = 2)
  mu <- c(0, 0)
  XZ1 <- mvrnorm(n, mu, sigma1)
  e1 <- rnorm(n, 0, 2)
  
  x <- XZ1[,1]
  z <- XZ1[,2]
 
# Simulate response
  beta0 <- 1. 
  betaz <- 0 # z has non-zero coefficient
  betax <- 2.5 # x has zero coefficient
  
  # Predict Y value
  y <- beta0 + betax*x + betaz*z + e1  
  
  #df1 <- cbind(XZ1, Y) %>% 
    #as.data.frame()
  
  #colnames(df1) <- c('x0', 'X', 'Z', 'Y')
  
  # Fit the models
  model1 <- lm(y ~ x + z)
  model2 <- lm(y ~ z)
  
  # coefficients for beta_X and beta_Z
  coeff1 <- model1$coefficients
  coeff2 <- model2$coefficients
  beta_X_values[i] <- coeff1['x']
  beta_Z_values[i] <- coeff1['z']
  beta_Z_value_model2[i]<- coeff2['z']
  
}

# Calculate the mean values for beta_X and beta_Z in each model
mean_beta_X <- mean(beta_X_values)
mean_beta_Z <- mean(beta_Z_values)

mean_beta_Z_model2 <- mean(beta_Z_value_model2)

# Print the mean values
cat("Mean value for beta_X in model 1:", mean_beta_X, "\n")
cat("Mean value for beta_Z in model 1:", mean_beta_Z, "\n")
cat("Mean value for beta_Z in model 2:", mean_beta_Z_model2, "\n")

mean_values_df2 <- data.frame(
  Variable = c("beta_X in model 1", "beta_Z in model 1", "beta_Z in model 2"),
  Mean_Value = c(mean_beta_X, mean_beta_Z, mean_beta_Z_model2)
)

Misperception 5b - Controlling for the linear effect of a covariate is equivalent to controlling for the covariate

Why This Misperception Occurs:

  • Nonlinear relationships can exist between variables, not captured by controlling only the linear term.
  • Example: U-shaped relationship between BMI and mortality, where linear modeling can lead to biased estimates.
  • Residual bias depends on the degree of nonlinearity in the relationship.

Avoiding the Misperception:

  • Assess for residual relationships or allow for nonlinearity from the onset(Andersen 2009).
  • Use nonlinear modeling techniques: splines, Box-Tidwell transformations, knotted regressions, polynomials, and fractional polynomials.

Misperception 6 - One should check whether covariates are normally distributed and take corrective action if not

Why This Misperception Occurs:

  • Common misunderstanding that variables in parametric models must be normally distributed.

  • In reality, normality is required for the residuals of the model, not the predictor variables or covariates.

  • Concern should focus on residuals and potential outliers, not the normality of covariates

Avoiding the Misperception:

  • Ensure residuals are uncorrelated with IVs and focus on residual analysis.

Misperception 7 - If the relation between a plausible confounder and the IV of interest is not statistically significant, the plausible confounder can be excluded with no concern for bias

  • Emphasis on statistical significance can lead to exclusion of important confounders using stepwise regression.
  • This approach may cause bias in coefficient estimates and significance tests by not accounting for the actual degree of confounding.

Why This Misperception Occurs:

  • Confusion between statistical significance and the actual degree of confounding.

  • Non-significant variables may still introduce bias if they are true confounders.

  • Motivation to create a parsimonious model, especially with many covariates and a modest sample size(VanderWeele 2019).

Avoiding the Misperception:

  • Include plausible confounders based on substantive knowledge, regardless of statistical significance.

  • Use sensitivity analysis: conduct the analysis with and without the confounder and report both results(VanderWeele and Ding 2017; Rosenbaum 2002).

  • Adopt the approach of seeking “an approximate answer to the right question”(Tukey1962?) rather than a precise answer to the wrong question.

Misperception 8 - Analyzing the residuals of an analysis in which a DV is regressed on the PBC is equivalent to including the PBC in an overall statistical model with the IV of interest

Why This Misperception Occurs:

  • Incorrect assumption that residual analysis (ANOVA on residuals) is mathematically equivalent to ANCOVA.

  • Residuals calculated separately for different IV levels or for the overall sample can introduce bias and incorrect significance tests.

Avoiding the Misperception:

  • Best to include PBCs directly in the model to avoid bias.

  • Residualizing for covariate effects beforehand can be used in complex models, but should be done with caution to avoid additional concerns.

Misperception 9 - Excluding a covariate that is not associated with the outcome of interest does not affect the association of the IV with the outcome

Why This Misperception Occurs:

  • The suppressor effect: A variable not correlated with the outcome can still impact the IV-outcome association.

  • Example: Adenovirus 36 (Ad36) infection affects adiposity and improves glucoregulatory function and lipid profiles, despite no direct association with the outcome(Akheruzzaman, Hegde, and Dhurandhar 2018).

  • Excluding such suppressor variables can lead to biased estimates of the IV’s effect on the outcome.

Avoiding the Misperception:

Sample R code for misperception 9

  • More details in supplementary file

The following function is used to generate data \((Y, X_A, X_B)\) with \(\beta_0\).

simDat <- function(n, ba, bb, la, lb, sa, sb, sz, se = 1) {
  Z <- rnorm(n, 0, sz)
  Xa <- la * Z + rnorm(n, 0, sa)
  Xb <- lb * Z + rnorm(n, 0, sb)
  Y <- ba * Xa + bb * Xb + rnorm(n, 0, se)
  data.frame(Y = Y, Xa = Xa, Xb = Xb)
}

The arguments to be supplied in simDat() are:

  • n is an integer value to specify the sample size
  • ba is a numerical value to specify \(\beta_A\)
  • bb is a numerical value to specify \(\beta_B\)
  • la is a numerical value to specify \(\lambda_A\)
  • lb is a numerical value to specify \(\lambda_B\)
  • sa is a numerical value to specify \(\sigma_\eta\)
  • sb is a numerical value to specify \(\sigma_\gamma\)
  • sz is a numerical value to specify \(\sigma_Z\)
  • se is a numerical value to specify \(\sigma_\epsilon\)
set.seed(0); dat <- simDat(500, .4, -.3, sqrt(.75), sqrt(.75), .5, .5, 1, .75)
dim(dat)
head(dat)

The following fits the full model and the reduced model based on dat.

## Full model
lm(Y ~ Xa + Xb, data = dat)
## Reduced model
lm(Y ~ Xa, data = dat)

Simulation output

Table 1: Parameters Used to Generate Simulated Data for the Simulation Studies Under Misperception 9

Scenario \(\beta_A\) \(\beta_B\) \(\lambda_A\) \(\lambda_B\) \(\sigma_\epsilon^2\) \(\sigma_\eta^2\) \(\sigma_\gamma^2\)
I -0.4 0.3 \(\frac{\sqrt{3}}{2}\) \(\frac{\sqrt{3}}{2}\) 0.93 0.25 0.25
II 0.4 -0.3 \(\frac{\sqrt{3}}{2}\) \(\frac{\sqrt{3}}{2}\) 0.93 0.25 0.25
III -0.5 0.24 0.8 0.6 0.8076 0.36 0.64
IV 0.5 -0.24 0.8 0.6 0.8076 0.36 0.64

Table 2: Summary of Bias When Fitting the Full Model (\(M_F\)) and the Reduced Model (\(M_R\)). The bias is defined as \(\hat{\beta}_A - \beta_A\), where \(\hat{\beta}_A\) is the least-squares estimate under the corresponding model.

\[\begin{array}{|c|cc|cc|cc|} \hline \text{Scenario} & {n=500} & & {n=1000} & &{n=2000} \\ & M_F & M_R & M_F & M_R & M_F & M_R \\ \hline \text{I} & -0.0007 & 0.2248 & 0.0001 & 0.2251 & -0.0001 & 0.2249 \\ \text{II} & 0.0005 & -0.2249 & 0.0003 & -0.2249 & 0.0002 & -0.2248 \\ \text{III} & -0.0001 & 0.2400 & 0.0003 & 0.2405 & -0.0003 & 0.2399 \\ \text{IV} & 0.0004 & -0.2396 & -0.0002 & -0.2400 & -0.0005 & -0.2402 \\ \hline \end{array}\]

Misperception 10 - If a plausible confounding variable is one that has a bivariate unadjusted correlation of zero with the IV, then it does not create bias in the association of the IV with the outcome

Why This Misperception Occurs:

  • Even if a PBC has zero correlation with the IV, it can still create bias if not properly included in the analysis.

  • Example: Collider bias can occur, altering associations and potentially masking true relationships(Munafò et al. 2017).

    • In Figure 4, both \(X\) and \(Z\) have a causal effect on \(Y_1\). \(Y_1\) can then be referred to as a ‘collider’

Avoiding the Misperception:

  • Use the backdoor criterion and causal DAGs to determine which variables to adjust for.

    • In Figure 4, for example, \(Y_1\) does not meet the backdoor criterion from \(Y_2\) to \(X\), and adjusting for it or selecting on it will bias the estimate of the effect estimate.

  • Recognize that zero correlation does not imply the absence of bias, and informed assumptions about causal structures are essential.

Sample R code for Misperception 10

  • More details in supplementary file
library(rethinking)
library(jtools)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
library(ppcor) ## To compute partial correlation
library(MTS)
set.seed(20211118)
bet.a = -.25
bet.b = .75  
lb.a =  .5
lb.binv = -bet.a/(bet.b*lb.a) ## Inverse of lb.b
lb.b <- lb.binv
rho.yb = bet.b + lb.binv*lb.a*bet.a
rho.ya = bet.a + lb.binv*lb.a*bet.b
N = 500
## Condition for lb.a 
# lb.aval =seq(0.5,1, length = 1000)
# lb.bval <- -bet.b/(bet.a*lb.aval)
# plot(lb.aval, 1/lb.bval)

SimMisconcept11 <- function(N, bet.a, bet.b, lb.a, lb.b){
  
  X  = rnorm(n=N) #belief Dieatary fat is ok 
  X <- X/sd(X)
  FatA <- rnorm(n=N, mean=lb.a*X, sd = sqrt(1 - lb.a^2)); #FatA = FatA/sd(FatA) #
  FatB <- rnorm(n=N, mean = lb.b*X, sd = sqrt(1 - lb.b^2)); #FatB = FatB/sd(FatB) #sqrt(1 - lb.b^2)
  
  var.y = bet.a^2 + bet.b^2 + 2*lb.a*lb.b*bet.a*(bet.b) 
  mu = bet.a*FatA + bet.b*FatB
  Y <- rnorm(N, mean = mu, sd = sqrt(1 - var.y)) # 1 );  #Y <- Y/sd(Y)
  
  #cat("\n---  Correlation values ---\n")
  
  #cat("Cor(Y, FATA)[empirical] - Cor(Y, FATA)[Theoritical] \n")
  c(cor(Y, FatA), rho.ya) 
  
  #cat("Cor(Y, FATB)[empirical] - Cor(Y, FATB)[Theoritical]\n")  
  c(cor(Y, FatB), rho.yb)
  
  #c(cor(X, FatA), cor(X, FatB))
  #c(cor(Y, FatA), rho.ya, cor(Y, FatB), rho.yb)
  #All_Models= list()
  mAB <- (lm(Y ~ FatA + FatB))
  mB <- (lm(Y ~ FatB))
  mA <- (lm(Y ~ FatA))
  
  # Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           cor(FatA, Y), # crude correlations
  #           cor(FatB, Y), # crude correlations
  #           pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  #           )
  Parm = c(0.0, bet.a, bet.b)
  Res <- c( c( c(summary(mAB)$coefficients[,1], (confint.lm(mAB)[,2] > Parm)*(confint.lm(mAB)[,1] < Parm)) ), #vector of length 6 (mean - cover in the 95%(1) or mot (0) )
            c(summary(mA)$coefficients[,1], (confint.lm(mA)[,2] > Parm[-3])*(confint.lm(mA)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            c(summary(mB)$coefficients[,1], (confint.lm(mB)[,2] > Parm[-2])*(confint.lm(mB)[,1] < Parm[-2])), #(mean - cover in the 95%(1) or mot (0) )
            cor(FatA, Y), # crude correlations
            cor(FatB, Y), # crude correlations
            pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB))$estimate[1,-1] # partial correlations vector of length 2
            #pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB, X = X))$estimate[1,-1]
  )
  
  (Res)
}


SimMisconcept11(N, bet.a, bet.b, lb.a, lb.b)

ResOUt <- t(sapply(rep(N, 5000), SimMisconcept11,  bet.a = bet.a, bet.b =bet.b, lb.a=lb.a, lb.b=lb.b))
colNam <- colnames(ResOUt)
colNam[15:18] <- c("rho.YFatA", "rho.YFatB", "prho.YFatA", "prho.YFatB")
colnames(ResOUt) <- colNam


cbind(colMeans(ResOUt), t(apply(ResOUt,2, quantile, probs =c(0.025, .975))))



#pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA, X =X))$estimate
#pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate

#cor(data.frame(Y = Y, FatB = FatB, FatA = FatA, X =X))
#--- Generate summaries


#SimMisconceptPlot <- function(N, bet.a, bet.b, lb.a, lb.b){
  
  X  = rnorm(n=N) #belief Dieatary fat is ok 
  X <- X/sd(X)
  FatA <- rnorm(n=N, mean=lb.a*X, sd = sqrt(1 - lb.a^2)); #FatA = FatA/sd(FatA) #
  FatB <- rnorm(n=N, mean = lb.b*X, sd = sqrt(1 - lb.b^2)); #FatB = FatB/sd(FatB) #sqrt(1 - lb.b^2)
  
  var.y = bet.a^2 + bet.b^2 + 2*lb.a*lb.b*bet.a*(bet.b) 
  mu = bet.a*FatA + bet.b*FatB
  Y <- rnorm(N, mean = mu, sd = sqrt(1 - var.y)) # 1 );  #Y <- Y/sd(Y)
  
  cat("\n---  Correlation values ---\n")
  
  cat("Cor(Y, FATA)[empirical] - Cor(Y, FATA)[Theoritical] \n")
  c(cor(Y, FatA), rho.ya) 
  
  cat("Cor(Y, FATB)[empirical] - Cor(Y, FATB)[Theoritical]\n")  
  c(cor(Y, FatB), rho.yb)
  
  #c(cor(X, FatA), cor(X, FatB))
  #c(cor(Y, FatA), rho.ya, cor(Y, FatB), rho.yb)
  #All_Models= list()
  mAB <- (lm(Y ~ FatA + FatB))
  mB <- (lm(Y ~ FatB))
  mA <- (lm(Y ~ FatA))
  plot_models(mAB,mA,mB, show.values = T)
  
  Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
            c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
            c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
            cor(FatA, Y), # crude correlations
            cor(FatB, Y), # crude correlations
            pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  )
  
  Res
}


#--- No ill formeed problem
#library(rethinking)
#cbind(lm(Y ~ FatA + FatB)$coefficients, confint.lm(lm(Y ~ FatA + FatB)))
#cbind(lm(Y ~ FatB)$coefficients, confint.lm(lm(Y ~  FatB))  )
#cbind(lm(Y ~ FatA)$coefficients, confint.lm(lm(Y ~  FatA))  )
#plot_models(mAB,mA,mB, show.values = T) #, value.offset=.3)

# summary(lm(Y ~ FatA + FatB))
# summary(lm(Y ~ FatB))
# summary(lm(Y~ FatA))
# coeftab(lm(Y~ FatA), lm(Y~ FatB), lm(Y~ FatA + FatB))
#plot( coeftab( lm(Y~ FatA), lm(Y~ FatA + FatB)), pars =c("FatA"))
#plot( coeftab( lm(Y~ FatB), lm(Y~ FatA + FatB)), pars =c("FatB"))

#plot_summs(lm(Y~ FatA), lm(Y~ FatB), lm(Y~ FatA + FatB), scale=T,model.names = c("Y~FATA","Y~FATB", "Y~FATA+FATB"))
#coeftab_plot( coeftab( lm(Y~ FatA), lm(Y~ FatB), lm(Y~ FatA + FatB) ), pars ="M")



#--- case of Misconception 10

library(rethinking)
library(jtools)
SimRo <- function(taucut){
  #runif(1,0,1)
  set.seed(2008)
  N = 1000000
  bet.z <- .5     #runif(1,0,1)
  bet.x <- .45      #runif(1,0,1)
  lb.z <- .4
  al.x <- 0.0     #runif(1,0,1)
  al.z <- .6     #runif(1,0,1)
  
  #--- simulate
  Z <- rnorm(n=N)
  #--- Simulate X 
  varx.z <- 1 - lb.z^2
  X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
  #X <- X/sd(X) 
  #--- Simulate Y2  
  vary2.z = 1 - (al.z^2 + al.x^2)
  Y2 <- al.z*Z + al.x*X + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
  #Y2 <- Y2/sd(Y2)
  #-- Simulate Y1
  vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
  Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd= sqrt(vary1.xz))
  #Y1 <- Y1/sd(Y1)
  
  id <- which(Y1 > 2.366733) #taucut*sd(Y1) + mean(Y1))
  R0 <- lm(Y1[id] ~  X[id])
  
  R00 <- lm(Y1[-id] ~  X[-id])
  
  R1 <- lm(Y1[id] ~  X[id] + Z[id])
  R10 <- lm(Y1[-id] ~  X[-id] + Z[-id])
  
  R20 <- lm(Y2 ~ X)
  R21 <- lm(Y2 ~ Z)
  R22 <- lm(Y2 ~ X + Z)
  
  R2id0 <- lm(Y2[id] ~ X[id] )
  R2id1 <- lm(Y2[id] ~ Z[id])
  R2id2 <- lm(Y2[id] ~ X[id] + Z[id])
  
  
  R2idn0 <- lm(Y2[-id] ~ X[-id])
  R2idn1 <- lm(Y2[-id] ~ Z[-id])
  R2idn2 <- lm(Y2[-id] ~ X[-id] + Z[-id])
  
  Res <- c(mean(Y1), sd(Y1), cor(X, Z) ,cor(X[id], Z[id]), cor(X[-id], Z[-id]), 
           as.numeric(R0$coefficients)[2], as.numeric(R00$coefficients)[2],
           as.numeric(R1$coefficients)[2:3],as.numeric(R10$coefficients)[2:3],
           as.numeric(R20$coefficients)[2], as.numeric(R21$coefficients)[2],as.numeric(R22$coefficients)[2:3],
           as.numeric(R2id0$coefficients)[2], as.numeric(R2id1$coefficients)[2], as.numeric(R2id2$coefficients)[2:3],
           as.numeric(R2idn0$coefficients)[2], as.numeric(R2idn1$coefficients)[2], as.numeric(R2idn2$coefficients)[2:3])
  
}   
SimRo <- Vectorize(SimRo,"taucut")#

#--- Run code
SimRo(1.30)

taucut <- seq(.1, 2, length.out = 500)
Res <- t(SimRo(taucut))
colnames(Res)  <- c("meanY1","sdY1", "rhoXZ(all.)" ,"rhoXZ(abov.)","rhoXZ(blw.)",
                    "betaxhat_Y1X(abov.)", "betaxhat_Y1X(blw.)",
                    "betaxhat_Y1XZ(abov)","betazhat_Y1Xz(blw)", "betaxhat_Y1XZ(blw)","betazhat_Y1Xz(blw)",
                    "alphax_Y2X(all)", "alphaz_Y2Z(all)", "alphax_Y2XZ(all)", "alphaz_Y2XZ(all)",
                    "alphax_Y2X(abov)", "alphaz_Y2Z(abov)", "alphax_Y2XZ(abov)", "alphaz_Y2XZ(abov)",
                    "alphax_Y2X(blw)","alphaz_Y2Z(blw)", "alphax_Y2XZ(blw)" , "alphaz_Y2XZ(blw)")


#--- Find the zero values
Idp <- which(Res[,4] > 0)
IdN <- which(Res[,4] < 0)

cat("\n ----- Chosen case ------ \n")
Res[c(Idp[length(Idp)],IdN[1]),]
taucut[c(Idp[length(Idp)],IdN[1])]
taucut0 = taucut[IdN[1]]
taucut0

Res[c(IdN[1]),]
#Res0Tep <- SimRoPlot(taucut[IdN[1]])
#dim(Res0Tep$Mod)

#coeftab(lm(Res0Tep$Mod[,1] ~ Res0Tep$Mod[,2]) )
plot(taucut,Res[,4], type="p", ylab = "cor(X,Z)")
abline(h=0, lwd=4, col="orange")
abline(v=taucut0, lwd=4, col="blue")


#--- Type error estimates
CorXZ <- function(a, Y1, X, Z){
  
  id <- which(Y1 > mean(Y1) + sd(Y1)*a)
  cor(X[id], Z[id]) 
}

CorXZ <- Vectorize(CorXZ,"a")



#--- Looking at Y1
SimMisconcept12 <- function(N){
  
  bet.z <- .5     #runif(1,0,1)
  bet.x <- .45      #runif(1,0,1)
  lb.z <- .4
  al.x <- 0.0     #runif(1,0,1)
  al.z <- .6     #runif(1,0,1)
  
  #--- simulate
  Z <- rnorm(n=N)
  #--- Simulate X 
  varx.z <- 1 - lb.z^2
  X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
  #X <- X/sd(X) 
  #--- Simulate Y2  
  vary2.z = 1 - al.z^2
  Y2 <- al.z*Z + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
  #Y2 <- Y2/sd(Y2)
  #-- Simulate Y1
  vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
  Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd= sqrt(vary1.xz))
  #Y1 <- Y1/sd(Y1)
  # CutVec <- seq(.1,2,length.out =300)
  # CoVec <- CorXZ(CutVec, Y1, X, Z)
  # Idp <- which(CoVec > 0)
  # IdN <- which(CoVec < 0)
  # 
  # id00 <- c(Idp[length(Idp)],IdN[1])[which.min(abs(CoVec[c(Idp[length(Idp)],IdN[1])]))]
  # taucut0 = CutVec[id00]
  taucut0 = 2.366733
  id <- which(Y1 >  2.366733)#taucut0*sd(Y1) + mean(Y1))
  
  #cat("\n---  Correlation values ---\n")
  
  #cat("Cor(Y, FATA)[empirical] - Cor(Y, FATA)[Theoritical] \n")
  #c(cor(Y, FatA), rho.ya) 
  
  #cat("Cor(Y, FATB)[empirical] - Cor(Y, FATB)[Theoritical]\n")  
  #c(cor(Y, FatB), rho.yb)
  
  #c(cor(X, FatA), cor(X, FatB))
  #c(cor(Y, FatA), rho.ya, cor(Y, FatB), rho.yb)
  #All_Models= list()
  mAB <- (lm(Y1[id] ~ X[id] + Z[id]))
  mB <- (lm(Y1[id] ~ X[id]))
  mA <- (lm(Y1[id] ~ Z[id]))
  
  # Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           cor(FatA, Y), # crude correlations
  #           cor(FatB, Y), # crude correlations
  #           pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  #           )
  Parm = c(0.0, bet.x, bet.z)
  Res <- c( c(cor(X[id], Z[id]),taucut0, c(summary(mAB)$coefficients[,1], (confint.lm(mAB)[,2] > Parm)*(confint.lm(mAB)[,1] < Parm)) ), #vector of length 6 (mean - cover in the 95%(1) or mot (0) )
            c(summary(mA)$coefficients[,1], (confint.lm(mA)[,2] > Parm[-3])*(confint.lm(mA)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            c(summary(mB)$coefficients[,1], (confint.lm(mB)[,2] > Parm[-2])*(confint.lm(mB)[,1] < Parm[-2])), #(mean - cover in the 95%(1) or mot (0) )
            cor(X[id], Y1[id]), # crude correlations
            cor(Z[id], Y1[id]), # crude correlations
            pcor(data.frame(Y1 = Y1[id], X = X[id], Z = Z[id]))$estimate[1,-1] # partial correlations vector of length 2
            #pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB, X = X))$estimate[1,-1]
  )
  
  (Res)
}

SimMisconcept12(500000)

#ResOUt_500 <- t(sapply(rep(500, 10000), SimMisconcept12))
#ResOUt_1000 <- t(sapply(rep(1000, 10000), SimMisconcept12))
ResOUt_2000 <- t(sapply(rep(500000, 10000), SimMisconcept12))


Temp <- ResOUt_2000 
colNam <- colnames(Temp)
colNam[c(1:2,17:20)] <- c("rho.ZX[Y1>tau]", "Tau", "rho.Y1X", "rho.Y1Z", "prho.Y1X", "prho.Y1Z")
colnames(Temp) <- colNam

#head(ResOUt)
Temp <- (Temp[abs(Temp[,1]) < 0.0006, ]);
dim(Temp)


xtable::xtable(round(cbind(colMeans(Temp), t(apply(Temp,2, quantile, probs =c(0.025, .975)))), 4))

#---------------------------------------------------
#--- Looking that Y2
#---------------------------------------------------

SimMisconcept12 <- function(N){
  
  bet.z <- .5     #runif(1,0,1)
  bet.x <- .45      #runif(1,0,1)
  lb.z <- .4
  al.x <- 0.     #runif(1,0,1)
  al.z <- .6     #runif(1,0,1)
  
  #--- simulate
  Z <- rnorm(n=N)
  #--- Simulate X 
  varx.z <- 1 - lb.z^2
  X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
  #X <- X/sd(X) 
  #--- Simulate Y2  
  vary2.z = 1 - al.z^2
  Y2 <- al.z*Z + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
  #Y2 <- Y2/sd(Y2)
  #-- Simulate Y1
  vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
  Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd= sqrt(vary1.xz))
  #Y1 <- Y1/sd(Y1)
  # CutVec <- seq(.1,2,length.out =300)
  # CoVec <- CorXZ(CutVec, Y1, X, Z)
  # Idp <- which(CoVec > 0)
  # IdN <- which(CoVec < 0)
  # 
  # id00 <- c(Idp[length(Idp)],IdN[1])[which.min(abs(CoVec[c(Idp[length(Idp)],IdN[1])]))]
  # taucut0 = CutVec[id00]
  OpRes <- optim(1.6, CorrelValZX, al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
  #taucut0 = 2.366733
  taucut0 = OpRes$par
  id <- which(Y1 > taucut0 )#taucut0*sd(Y1) + mean(Y1))
  
  #cat("\n---  Correlation values ---\n")
  
  #cat("Cor(Y, FATA)[empirical] - Cor(Y, FATA)[Theoritical] \n")
  #c(cor(Y, FatA), rho.ya) 
  
  #cat("Cor(Y, FATB)[empirical] - Cor(Y, FATB)[Theoritical]\n")  
  #c(cor(Y, FatB), rho.yb)
  
  #c(cor(X, FatA), cor(X, FatB))
  #c(cor(Y, FatA), rho.ya, cor(Y, FatB), rho.yb)
  #All_Models= list()
  mAB <- (lm(Y2[id] ~ X[id] + Z[id]))
  mB <- (lm(Y2[id] ~ X[id]))
  mA <- (lm(Y2[id] ~ Z[id]))
  
  # Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           cor(FatA, Y), # crude correlations
  #           cor(FatB, Y), # crude correlations
  #           pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  #           )
  Parm = c(0.0, al.x, al.z)
  Res <- c( c(cor(X[id], Z[id]),taucut0, c(summary(mAB)$coefficients[,1], (confint.lm(mAB)[,2] > Parm)*(confint.lm(mAB)[,1] < Parm)) ), #vector of length 6 (mean - cover in the 95%(1) or mot (0) )
            c(summary(mA)$coefficients[,1], (confint.lm(mA)[,2] > Parm[-3])*(confint.lm(mA)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            c(summary(mB)$coefficients[,1], (confint.lm(mB)[,2] > Parm[-2])*(confint.lm(mB)[,1] < Parm[-2])), #(mean - cover in the 95%(1) or mot (0) )
            cor(X[id], Y2[id]), # crude correlations
            cor(Z[id], Y2[id]), # crude correlations
            pcor(data.frame(Y2 = Y2[id], X = X[id], Z = Z[id]))$estimate[1,-1] # partial correlations vector of length 2
            #pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB, X = X))$estimate[1,-1]
  )
  
  (Res)
}


CorrelValZX <- function(al.z, bet.z, bet.x, lb.z, a){
  
  Sig <- CovMat(al.z, bet.z, bet.x, lb.z)
  m1 = MTS::msqrt(Sig[1:2,1:2])
  
  
  #---------------
  Al = c(1/sqrt( 1 - (t(Sig[1:2,3])%*%solve(Sig[1:2,1:2])%*%(Sig[1:2,3]) ) ))*solve(Sig[1:2,1:2])%*%Sig[1:2,3]
  Del = (1/sqrt(c(1 + t(Al)%*%Sig[1:2,1:2]%*%Al)))*Sig[1:2,1:2]%*%Al
  
  (1/sqrt( 1 - (t(Sig[1:2,3])%*%solve(Sig[1:2,1:2])%*%(Sig[1:2,3]) ) ))*(1.6)
  1.6*sqrt(c(1 + t(Al)%*%Sig[1:2,1:2]%*%Al))
  
  #delta = c(msqrt(1 + t(Lb1)%*%Sig[1:2,1:2]%*%(Lb1))$invsqrt)*Sig[1:2,1:2]%*%Lb1
  mu = (Del*(dnorm(-a)/pnorm(-a)))
  
  # Sig[1:2,1:2] - diag(1/c(mu))%*%(tcrossprod(mu))%*%diag(1/c(mu)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  SigF =  Sig[1:2,1:2] + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a)))*(tcrossprod(Del)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  
  Res = list(Mu = mu, Sig = SigF, CorSig = cov2cor(SigF)) 
  (SigF[1,2])^2
  
}


CovMat <- function(al.z, bet.z, bet.x, lb.z){
  
  Sig <- matrix(0, nrow=3,ncol=3)
  Sig[1,2] = Sig[2,1] =  lb.z  # cor(X,Z)
  Sig[1,3] = Sig[3,1] = bet.z + lb.z*bet.x # cor(Y1, Z)
  Sig[2,3] = Sig[3,2] = bet.x + bet.z*lb.z  #cor(Y1, X)
  diag(Sig) <- 1
  
  Sig
  
}

CovMat(al.z, bet.z, bet.x, lb.z)


optim(1.6, CorrelValZX, al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)

SimMisconcept12(50000)

#ResOUt2_500 <- t(sapply(rep(500, 10000), SimMisconcept12))
#ResOUt2_1000 <- t(sapply(rep(1000, 10000), SimMisconcept12))
ResOUt2_2000 <- t(sapply(rep(500000, 10000), SimMisconcept12))


Temp <- ResOUt2_2000 
colNam <- colnames(Temp)
colNam[c(1:2,17:20)] <- c("rho.ZX[Y1>tau]", "Tau", "rho.Y2X", "rho.Y2Z", "prho.Y2X", "prho.Y2Z")
colnames(Temp) <- colNam

#head(ResOUt)
sum(abs(Temp[,1]) < 0.0006)
Temp <- (Temp[abs(Temp[,1]) < 0.0006, ]);
dim(Temp)

round(cbind(colMeans(Temp), t(apply(Temp,2, quantile, probs =c(0.025, .975)))), 4)



#----------------------------------------------------------------------
#------ Including Y2 
#----------------------------------------------------------------------
bet.z <- .5     #runif(1,0,1)
bet.x <- .45      #runif(1,0,1)
lb.z <- .4  # old
al.x <- 0.0     #runif(1,0,1)
al.z <- .6     #runif(1,0,1)

#--- New
bet.z <- .2   #runif(1,0,1)
bet.x <- .6      #runif(1,0,1)
lb.z <- .4  # old


#--- CHoosing paramaters

ParmSelc <- function(lb.z){
  
  betax <- c(0, sqrt(1/(1+lb.z^2)), runif(n=1, 0, sqrt(1/(1+lb.z^2))) )
  betaz <- c(0, sqrt(1- betax[3]^2) - lb.z*betax[3], runif(n=1, 0, sqrt(1- betax[3]^2) - lb.z*betax[3]))
  cbind(betax, betaz, (betaz[3] + lb.z*betax[3])*(betax[3] + lb.z*betaz[3]))
}

ParmSelc(.4)

ParmSelc <- function(lb.z, a,b){
  
  betax <-  a*sqrt(1/(1+lb.z^2))
  betaz <- (b/a)* betax*(sqrt(1+lb.z*lb.z - a*a) - lb.z*a)
  c(betax, betaz, lb.z, (betaz + lb.z*betax)*(betax + lb.z*betaz) - lb.z, lb.z/((betaz + lb.z*betax)*(betax + lb.z*betaz)) )
}
ParmSelc <- Vectorize(ParmSelc,"lb.z")


ParmSelcMin <- function(lb.z, a,b){
  
  betax <-  a*sqrt(1/(1+lb.z^2))
  betaz <- (b/a)* betax*(sqrt(1+lb.z*lb.z - a*a) - lb.z*a)
  res = c(betax, betaz, lb.z, (betaz + lb.z*betax)*(betax + lb.z*betaz) - lb.z, lb.z/((betaz + lb.z*betax)*(betax + lb.z*betaz)) )
  (res[4])^2
}
ParmSelcMin <- Vectorize(ParmSelcMin,"lb.z")

a <- c(.2,.9)
b <- c(.2,.9)


lb.z0 <- seq(0, 1, length.out = 1000)
par(mfrow=c(2,2), font=12, cex.axis=1.6, font.axis=12, cex.lab=1.3, font.lab= 12)

LambdaZFin <- c(0.02, .035, 0.27,.4, 0.12, 0.19, 0.26, 0.45)  ## These are the value of lb.z we will use

cp =0
names <-c("N","betx","betz","lbz","tau",
          "corSel",
          "corFull",
          "PcorSel",
          "PcorFull",
          
          "mb_1","mb_2",
          "Fmb_1","Fmb_2",
          "mab_1","mab_2", "mab_3",
          "Fmab_1","Fmab_2","Fmab_3",
          #texting
          "mab_intp","mab_alx","mab_alz","mab_pb1","mab_pb2","mab_pb3", 
          "mb_intp","mb_alx","mb_pb1","mb_pb2",
          "Fmab_intp","Fmab_alx","Fmab_alz","Fmab_pb1","Fmab_pb2","Fmab_pb3", 
          "Fmb_intp","Fmb_alx","Fmb_pb1","Fmb_pb2")
Final_Res <- list() 
OutFinal <- NULL
for(i in 1:length(a)){
  Final_Res[[i]] <- list()
  for(j in 1:length(b)){
    Final_Res[[i]][[j]] <- list()
    #-Cat generate plots
    Re0 <- optimize(ParmSelcMin, c(0, .99) , a= a[i], b= b[j])
    lb.z00 <- Re0$minimum
    V00 <- ParmSelc(lb.z0, a[i], b[j])
    plot(lb.z0, V00[4,], xlab = expression(lambda[z]), ylab= "Objective", main = paste("a=",a[i],";","b=",b[j], sep=""), type="l", lwd=3)
    abline(h = 0,lwd=2)
    head(V00[,1:10])
    x0 <- lb.z0[lb.z0 < lb.z00]
    y0<- V00[4, lb.z0 < lb.z00] 
    polygon(c(x0, rev(x0)), c(y0 ,rep(0, length(y0))),
            col = "#6BD7AF")
    abline(v =lb.z00, col="grey",lwd=3,lty=2)
    abline(v =0, col="grey",lwd=3,lty=2)
    #- Perform Simulation Choosing two cases
    #-- case
    cp <- cp+1
    Parm0 <- ParmSelc(LambdaZFin[c((2*(cp-1)+1):(2*cp))], a[i],b[j])
    for(m in 1:2){
      cat("case i, j m===>", i,j, m,"\n")
      Final_Res[[i]][[j]][[m]] <- list()
      bet.x = Parm0[1,m]
      bet.z = Parm0[2,m]
      lb.z <- Parm0[3,m] 
      OT <- t(replicate(1000, SimMisconcept12V2(N, bet.x, bet.z, lb.z)))
      colnames(OT) <- names
      
      
      Final_Res[[i]][[j]][[m]] <- c(a=a[i], b= b[j], lb.z00,colMeans(OT))
      OutFinal <- rbind(OutFinal, c(a=a[i], b= b[j], lb.z00,colMeans(OT)))
    }
  }
}


#---- Summarise the simulation 
rbind(Final_Res[[1]][[1]][[1]],Final_Res[[1]][[1]][[2]])
OutFinal



V00 <- ParmSelc(lb.z0, .2,.2)
par(mfrow=c(2,2), font=12, cex=1.5)
plot(lb.z0, V00[4,],xlab = expression(lambda[z]))
abline(h = 0,lwd=3)
head(V00[,1:10])


plot(lb.z0, V00[5,])
abline(h = 1,lwd=3)
head(V00[,1:10])
a = seq(-6,6,length.out=1000)
Y00 <- (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a)))
plot(a,-Y00)

which.min(abs(V00[4,]))
id = 400
V00[,id]
bet.x <- V00[1,id]
bet.z <- V00[2,id]
lb.z <- V00[3,id]

#@SImulate data

SimMisconcept12V2 <- function(N, bet.x, bet.z, lb.z){
  
  #bet.z <- .5     #runif(1,0,1)
  #bet.x <- .45      #runif(1,0,1)
  #lb.z <- .4
  al.x <- 0.     #runif(1,0,1)
  al.z <- .6     #runif(1,0,1)
  
  #--- simulate
  Z <- rnorm(n=N)
  #--- Simulate X 
  #varx.z <- 1 - lb.z^2
  varx.z <- 1 - var(Z)*lb.z^2
  X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
  X <- X/sd(X) 
  #--- Simulate Y2  
  #vary2.z = 1 - al.z^2
  vary2.z = 1 - (var(Z)*al.z^2 + (al.x^2)*var(X) + 2*(al.x)*al.z*cov(X, Z))
  Y2 <- al.z*Z + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
  Y2 <- Y2/sd(Y2)
  #-- Simulate Y1
  #vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
  vary1.xz <- 1 - (var(Z)*bet.z^2 +  var(X)*(bet.x^2) + 2*bet.z*bet.x*cov(X, Z))
  Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd= sqrt(vary1.xz))
  Y1 <- Y1/sd(Y1)
  # CutVec <- seq(.1,2,length.out =300)
  # CoVec <- CorXZ(CutVec, Y1, X, Z)
  # Idp <- which(CoVec > 0)
  # IdN <- which(CoVec < 0)
  # 
  # FInd cut-off pount a for selecting on Y1
  #OpRes <- optim(1.6, CorrelValZX, al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
  OpRes = optimize(CorrelValObjV2, c(-7, 7), al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
  #taucut0 = 2.366733
  taucut0 = OpRes$minimum
  id <- which(Y1 > taucut0) #taucut0*sd(Y1) + mean(Y1))
  #All_Models= list()
  mAB <- (lm(Y2[id] ~ X[id] + Z[id])) #(4)
  mB <- (lm(Y2[id] ~ X[id])) # (3)
  ## mA <- (lm(Y2[id] ~ Z[id]))
  FmAB <-  (lm(Y2 ~ X + Z)) #(4)
  FmB <-  (lm(Y2 ~ X )) # (3)
  #FmAB <-  (lm(Y2 ~ X + Z))
  # Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           cor(FatA, Y), # crude correlations
  #           cor(FatB, Y), # crude correlations
  #           pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  #           )
  Parm = c(0.0, al.x, al.z)
  #Out <- NULL
  
  Res <- c( N, bet.x, bet.z, lb.z, taucut0,
            cor(X[id], Y2[id]), # crude correlations (1)
            cor(X, Y2), # crude correlations   (1)
            pcor(data.frame(Y2 = Y2[id], X = X[id], Z = Z[id]))$estimate[1,2], #-1], # partial correlations vector of length 2 (2)
            pcor(data.frame(Y2 = Y2, X = X, Z = Z))$estimate[1, 2], #-1] # partial correlations vector of length 2  (2)
            
            c(summary(mB)$coefficients[,1] - Parm[-3]), #(3)
            c(summary(FmB)$coefficients[,1] - Parm[-3]), #(3)
            c(summary(mAB)$coefficients[,1] - Parm), # (4)
            c(summary(FmAB)$coefficients[,1] - Parm),# (4)
            #-- Testing
            c(c(summary(mAB)$coefficients[,1], (confint.lm(mAB)[,2] > Parm)*(confint.lm(mAB)[,1] < Parm)) ), #vector of length 6 (mean - cover in the 95%(1) or mot (0) )
            #c(summary(mA)$coefficients[,1], (confint.lm(mA)[,2] > Parm[-3])*(confint.lm(mA)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            c(summary(mB)$coefficients[,1], (confint.lm(mB)[,2] > Parm[-3])*(confint.lm(mB)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            
            c(summary(FmAB)$coefficients[,1], (confint.lm(FmAB)[,2] > Parm)*(confint.lm(FmAB)[,1] < Parm)),
            c(summary(FmB)$coefficients[,1], (confint.lm(FmB)[,2] > Parm[-3])*(confint.lm(FmB)[,1] < Parm[-3])) #(mean - cover in the 95%(1) or mot (0) )
            
            
            #pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB, X = X))$estimate[1,-1]
  )
  #Out <- rbind(Out, Res)
  return(Res)
}

#@ estimate partial Slope - will take a covariance matrix with the last row and column for Y 
PartialSlope <- function(Cov){
  J = ncol(Cov)
  sdt = (sqrt(diag(Cov)))
  CorMat <- cov2cor(Cov)
  Bet <- numeric(J-1)
  Bet[1] = (sdt[J]/sdt[1])*(CorMat[1,3] - CorMat[2,3]*CorMat[1,2] )/(1 - CorMat[1,2]^2) ## Z
  Bet[2] = (sdt[J]/sdt[2])*(CorMat[2,3] - CorMat[1,3]*CorMat[1,2] )/(1 - CorMat[1,2]^2) ## for X
  
  Bet  
}

#@ Compute the inverse of the Correlation Matrix
InvMat3by3 <- function(Sig){
  InvSig <- matrix(0, nrow=3, ncol = 3)
  
  InvSig[1,1] = Sig[3,3]*Sig[2,2] - Sig[2,3]*Sig[2,3]
  InvSig[1,2] = Sig[1,3]*Sig[2,3] - Sig[3,3]*Sig[1,2]
  InvSig[1,3] = Sig[1,2]*Sig[2,3] - Sig[1,3]*Sig[2,2]
  
  InvSig[2,2] = Sig[3,3]*Sig[1,1] - Sig[1,3]*Sig[1,3]
  InvSig[2,3] = Sig[1,2]*Sig[1,3] - Sig[1,1]*Sig[2,3] ## 
  
  InvSig[3,3] = Sig[1,1]*Sig[2,2] - Sig[1,2]*Sig[1,2]
  
  
  DetVal <- Sig[1,1]*InvSig[1,1] + Sig[1,2]*InvSig[1,2] + Sig[1,3]*InvSig[1,3]
  DetVal
  
  InvSigF <- (1/DetVal)*(InvSig + t(InvSig) - diag(diag(InvSig)) )
  list(InverseCov = InvSigF, Partial_Correlation = (-1)*cov2cor(InvSigF))
}



names <-c("N","betx","betz","lbz","tau",
          "corSel",
          "corFull",
          "PcorSel",
          "PcorFull",
          
          "mb_1","mb_2",
          "Fmb_1","Fmb_2",
          "mab_1","mab_2", "mab_3",
          "Fmab_1","Fmab_2","Fmab_3",
          #texting
          "mab_intp","mab_alx","mab_alz","mab_pb1","mab_pb2","mab_pb3", 
          "mb_intp","mb_alx","mb_pb1","mb_pb2",
          "Fmab_intp","Fmab_alx","Fmab_alz","Fmab_pb1","Fmab_pb2","Fmab_pb3", 
          "Fmb_intp","Fmb_alx","Fmb_pb1","Fmb_pb2")


SimMisconcept12 <- function(N){
  
  bet.z <- .5     #runif(1,0,1)
  bet.x <- .45      #runif(1,0,1)
  lb.z <- .4
  al.x <- 0.     #runif(1,0,1)
  al.z <- .6     #runif(1,0,1)
  
  #--- simulate
  Z <- rnorm(n=N)
  #--- Simulate X 
  #varx.z <- 1 - lb.z^2
  varx.z <- 1 - var(Z)*lb.z^2
  X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
  #X <- X/sd(X) 
  #--- Simulate Y2  
  #vary2.z = 1 - al.z^2
  vary2.z = 1 - (var(Z)*al.z^2 + (al.x^2)*var(X) + 2*(al.x)*al.z*cov(X, Z))
  Y2 <- al.z*Z + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
  #Y2 <- Y2/sd(Y2)
  #-- Simulate Y1
  #vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
  vary1.xz <- 1 - (var(Z)*bet.z^2 +  var(X)*(bet.x^2) + 2*bet.z*bet.x*cov(X, Z))
  Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd= sqrt(vary1.xz))
  #Y1 <- Y1/sd(Y1)
  # CutVec <- seq(.1,2,length.out =300)
  # CoVec <- CorXZ(CutVec, Y1, X, Z)
  # Idp <- which(CoVec > 0)
  # IdN <- which(CoVec < 0)
  # 
  # FInd cut-off pount a for selecting on Y1
  OpRes <- optim(1.6, CorrelValZX, al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
  #taucut0 = 2.366733
  taucut0 = OpRes$par
  id <- which(Y1 > taucut0 )#taucut0*sd(Y1) + mean(Y1))
  #All_Models= list()
  mAB <- (lm(Y2[id] ~ X[id] + Z[id]))
  mB <- (lm(Y2[id] ~ X[id]))
  mA <- (lm(Y2[id] ~ Z[id]))
  
  # Res <- c( c((cbind(summary(mAB)$coefficients[,1], confint.lm(mAB)))), #vector of length 9 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mA)$coefficients[,1], confint.lm(mA)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           c((cbind(summary(mB)$coefficients[,1], confint.lm(mB)))), #vector of length 6 (mean - 2.5% - 97.5%)
  #           cor(FatA, Y), # crude correlations
  #           cor(FatB, Y), # crude correlations
  #           pcor(data.frame(Y = Y, FatB = FatB, FatA = FatA))$estimate[1,-1] # partial correlations vector of length 2
  #           )
  Parm = c(0.0, al.x, al.z)
  Res <- c( c(cor(X[id], Z[id]),taucut0, c(summary(mAB)$coefficients[,1], (confint.lm(mAB)[,2] > Parm)*(confint.lm(mAB)[,1] < Parm)) ), #vector of length 6 (mean - cover in the 95%(1) or mot (0) )
            c(summary(mA)$coefficients[,1], (confint.lm(mA)[,2] > Parm[-3])*(confint.lm(mA)[,1] < Parm[-3])), #(mean - cover in the 95%(1) or mot (0) )
            c(summary(mB)$coefficients[,1], (confint.lm(mB)[,2] > Parm[-2])*(confint.lm(mB)[,1] < Parm[-2])), #(mean - cover in the 95%(1) or mot (0) )
            
            cor(X[id], Y2[id]), # crude correlations
            cor(Z[id], Y2[id]), # crude correlations
            pcor(data.frame(Y2 = Y2[id], X = X[id], Z = Z[id]))$estimate[1,-1], # partial correlations vector of length 2
            pcor(data.frame(Y2 = Y2, X = X, Z = Z))$estimate[1,-1] # partial correlations vector of length 2
            #pcor(data.frame(Y = Y, FatA = FatA, FatB = FatB, X = X))$estimate[1,-1]
  )
  
  (Res)
}

#@ CovMat :  Compute the theoritical correlation  matrix (yield a covariance matrix for the vector Z, X, Y2, Y1)
CovMat <- function(al.z, bet.z, bet.x, lb.z, al.x){
  
  Sig <- matrix(0, nrow=4,ncol=4)
  Sig[1,2] = Sig[2,1] =  lb.z  # cor(X,Z)
  Sig[1,3] = Sig[3,1] =  al.z + al.x*lb.z # cor(Y2, Z)
  Sig[1,4] = Sig[4,1] = bet.z + bet.x*lb.z  #cor(Y1, Z)
  
  Sig[2,3] = Sig[3,2] = al.x + al.z*lb.z  #cor(Y2, X)
  Sig[2,4] = Sig[4,2] = bet.x + bet.z*lb.z  #cor(Y1, X)
  
  Sig[4,3] = Sig[3,4] = bet.x*(al.x + al.z*lb.z) + bet.z*(al.z + lb.z*al.x)   #(al.x*lb.z + al.z)*(bet.x*lb.z + bet.z)  #cor(Y2, Y1)
  diag(Sig) <- 1
  
  Sig
  
}

#@ EstimCov :  Compute empirical mean and covariance matrix between Z, X Y2
EstimCov<- function(Z,X,Y2){
  
  meanVal = c(mean(Z), mean(X), mean(Y2))
  
  dat <- data.frame(Z, X, Y2)
  Covm <- cov(dat)
  list(Mean = meanVal, Cov = Covm)
  
}

#@ CorrelValZX4 :  Compute theoritical means and covariance matrix between Z, X Y2, Y1 (given value of a)
CorrelValZX4 <- function(al.z, bet.z, bet.x, lb.z, a){
  
  Sig <- CovMat(al.z, bet.z, bet.x, lb.z, al.x)
  m1 = MTS::msqrt(Sig[1:3,1:3])
  
  
  #---------------
  Al = c(1/sqrt( 1 - (t(Sig[1:3,4])%*%solve(Sig[1:3,1:3])%*%(Sig[1:3,4]) ) ))*solve(Sig[1:3,1:3])%*%Sig[1:3,4]
  Del = (1/sqrt(c(1 + t(Al)%*%Sig[1:3,1:3]%*%Al)))*Sig[1:3,1:3]%*%Al
  
  # (1/sqrt( 1 - (t(Sig[1:3,4])%*%solve(Sig[1:3,1:3])%*%(Sig[1:3,4]) ) ))*(1.6)
  #  1.6*sqrt(c(1 + t(Al)%*%Sig[1:2,1:2]%*%Al))
  
  #delta = c(msqrt(1 + t(Lb1)%*%Sig[1:2,1:2]%*%(Lb1))$invsqrt)*Sig[1:2,1:2]%*%Lb1
  mu = c((Del*(dnorm(-a)/pnorm(-a))), dnorm(a)/pnorm(-a))
  
  # Sig[1:2,1:2] - diag(1/c(mu))%*%(tcrossprod(mu))%*%diag(1/c(mu)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  SigF <- matrix(0, nrow=4, ncol=4)
  SigF[1:3,1:3] =  Sig[1:3,1:3] + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a)))*(tcrossprod(Del)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  SigF[,4] = SigF[4,] = ( 1 + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a))) )*Sig[,4]
  
  
  Res = list(Mu = mu, Sig = SigF, CorSig = cov2cor(SigF)) 
  Res
  
}

#@ CorrelValObj :  objective dunction to obtain a
CorrelValObj <- function(al.z, bet.z, bet.x, lb.z, a){
  
  Sig <- CovMat(al.z, bet.z, bet.x, lb.z, al.x)
  m1 = MTS::msqrt(Sig[1:3,1:3])
  
  
  #---------------
  Al = c(1/sqrt( 1 - (t(Sig[1:3,4])%*%solve(Sig[1:3,1:3])%*%(Sig[1:3,4]) ) ))*solve(Sig[1:3,1:3])%*%Sig[1:3,4]
  Del = (1/sqrt(c(1 + t(Al)%*%Sig[1:3,1:3]%*%Al)))*Sig[1:3,1:3]%*%Al
  
  # (1/sqrt( 1 - (t(Sig[1:3,4])%*%solve(Sig[1:3,1:3])%*%(Sig[1:3,4]) ) ))*(1.6)
  #  1.6*sqrt(c(1 + t(Al)%*%Sig[1:2,1:2]%*%Al))
  
  #delta = c(msqrt(1 + t(Lb1)%*%Sig[1:2,1:2]%*%(Lb1))$invsqrt)*Sig[1:2,1:2]%*%Lb1
  mu = (Del*(dnorm(-a)/pnorm(-a)))
  
  # Sig[1:2,1:2] - diag(1/c(mu))%*%(tcrossprod(mu))%*%diag(1/c(mu)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  SigF =  Sig[1:3,1:3] + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a)))*(tcrossprod(Del)) #m1$invsqrt%*%(tcrossprod(mu))%*%m1$invsqrt
  
  Res = list(Mu = mu, Sig = SigF, CorSig = cov2cor(SigF)) 
  (SigF[1,2])^2
  
}
CorrelValObj(al.z, bet.z, bet.x, lb.z, 1.6)


CorrelValObjV2 <- function(al.z, bet.z, bet.x, lb.z, a){
  
  (lb.z/((bet.z+lb.z*bet.x)*(bet.x+lb.z*bet.z)) + (-(dnorm(-a)/pnorm(-a))^2 + (a*dnorm(-a)/pnorm(-a))))^2
  
}

Obj = optimize(CorrelValObjV2, c(-7, 7), al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
Obj
CorrelValObjV2(al.z, bet.z, bet.x, lb.z, 1.6)
X0 = seq(-3,3,length.out=1000)

plot(X0, CorrelValObjV2(al.z, bet.z, bet.x, lb.z, X0))



optim(1.6, CorrelValObj, al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)

Obj = optimize(CorrelValObj, c(-7, 7), al.z=al.z, bet.z = bet.z, bet.x = bet.x, lb.z = lb.z)
Obj

#--- SImulate data
#--- simulate
N = 50000
Z <- rnorm(n=N)
Z <- Z/sd(Z)
#--- Simulate X 
#----------------------------------------------------------------------
varx.z <- 1 - var(Z)*lb.z^2
X = lb.z*Z + rnorm(n=N, mean = 0, sd = sqrt(varx.z))
X <- X/sd(X) 
#--- Simulate Y2  
#----------------------------------------------------------------------
#vary2.z = 1 - (al.z^2 +  (al.x^2)*(lb.z^2) + 2*al.z*al.x*lb.z + (al.x^2))
vary2.z = 1 - (var(Z)*al.z^2 + (al.x^2)*var(X) + 2*(al.x)*al.z*cov(X, Z))
Y2 <- al.z*Z + al.x*X + rnorm(n=N, mean = 0, sd = sqrt(vary2.z))
Y2 <- Y2/sd(Y2)
#-- Simulate Y1
#----------------------------------------------------------------------
#vary1.xz <- 1 - (bet.z^2 +  (bet.x^2)*(lb.z^2) + 2*bet.z*bet.x*lb.z + (bet.x^2))
vary1.xz <- 1 - (var(Z)*bet.z^2 +  var(X)*(bet.x^2) + 2*bet.z*bet.x*cov(X, Z))
Y1 <-  bet.z*Z + bet.x*X + rnorm(n=N, mean = 0, sd = sqrt(vary1.xz))
Y1 <- Y1/sd(Y1)
id <- which(Y1 > Obj$minimum) ## id to be selected
length(id)

SigTe <- CovMat(al.z, bet.z, bet.x, lb.z, al.x) ## Covariance matrix of 4 times 4
SigTe
solve(SigTe[1:3,1:3])
InvMat3by3(SigTe[1:3,1:3])
dat=data.frame(Z=Z, X=X, Y2=Y2, Y1=Y1)
cov(dat)
PartialSlope(SigTe[1:3,1:3])

al.x*bet.x + al.x*lb.z*bet.z+al.z*bet.z + al.z*lb.z*bet.x
bet.x*(al.x + al.z*lb.z) + bet.z*(al.z + lb.z*al.x)


EstimCov(Z, X, Y2)                     ## Estimate covariance based on the data
EstimCov(Z[id], X[id], Y2[id])         ## Estimated Mean and Covriance after selection
dat=data.frame(Z=Z[id], X=X[id], Y2=Y2[id], Y1=Y1[id])
cov(dat)
cov2cor(cov(dat))
SigTea <- CorrelValZX4(al.z, bet.z, bet.x, lb.z, Obj$minimum)
SigTea
PartialSlope(SigTea$Sig[1:3,1:3])

summary(lm(Y2[id] ~  X[id]))
summary(lm(Y2[id] ~  Z[id]))

summary(lm(Y2[id] ~ X[id] + Z[id]))
summary(lm(Y2 ~ X+Z))
summary(lm(Y2 ~ X))

summary(lm(Y1[id] ~  X[id]))
summary(lm(Y1[id] ~  Z[id]))
summary(lm(Y1[id] ~ X[id] + Z[id]))
summary(lm(Y1 ~ X+Z))
summary(lm(Y1 ~ X))



a = Obj$minimum
1 + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a)))
(1 + (-(dnorm(-a)/pnorm(-a))^2 - (-a*dnorm(-a)/pnorm(-a))))*(cov(dat)[,4])

(bet.x + lb.z*bet.z)*(bet.z + bet.x*lb.z)

a <- seq(-10,10, length.out = 5000)
a= 2.36494
V0 <- ((dnorm(-a)/pnorm(-a))^2 + (-a*dnorm(-a)/pnorm(-a)))
plot(a, V0)

lb.z0 <- seq(0,1,length= 1000)
V1 <- (bet.x + lb.z0*bet.z)*(bet.z + bet.x*lb.z0)
plot(lb.z0, V1)
abline(a=0,b=1,lwd=3,col="red")

(bet.x + lb.z*bet.z)*(bet.z + bet.x*lb.z)

cor(X,Z)
cor(Z, Y2)
cor(Z, Y1)

cor(X,Y1)
lb.z*(bet.z + bet.x*lb.z)

cor(X, Y2)
cor(X[id], Y2[id])

cor(Y2, Y1)


Sig <-  CovMat(al.z, bet.z, bet.x, lb.z, al.x)[1:3,1:3] ## Check the partial correlation not conditionning on Y1
Sig <- CorrelValZX4(al.z, bet.z, bet.x, lb.z, Obj$minimum)$CorSig[1:3,1:3] ## Check the partial correlation not conditionning on Y1


#--- Compute the inverse of the Correlation Matrix

InvMat3by3 <- function(Sig){
  InvSig <- matrix(0, nrow=3, ncol = 3)
  
  InvSig[1,1] = Sig[3,3]*Sig[2,2] - Sig[2,3]*Sig[2,3]
  InvSig[1,2] = Sig[1,3]*Sig[2,3] - Sig[3,3]*Sig[1,2]
  InvSig[1,3] = Sig[1,2]*Sig[2,3] - Sig[1,3]*Sig[2,2]
  
  InvSig[2,2] = Sig[3,3]*Sig[1,1] - Sig[1,3]*Sig[1,3]
  InvSig[2,3] = Sig[1,2]*Sig[1,3] - Sig[1,1]*Sig[2,3] ## 
  
  InvSig[3,3] = Sig[1,1]*Sig[2,2] - Sig[1,2]*Sig[1,2]
  
  
  DetVal <- Sig[1,1]*InvSig[1,1] + Sig[1,2]*InvSig[1,2] + Sig[1,3]*InvSig[1,3]
  DetVal
  
  InvSigF <- (1/DetVal)*(InvSig + t(InvSig) - diag(diag(InvSig)) );
  InvSigF
}

Simulation output

Figure 5: Possible values of \(λ_z\) based on each choice of the pairs of a, b. The area shaded in green denotes the area for which a \(λ_z\) value has a value \(\tau\) that makes Equation 13 equal zero.

Table 3: Estimated Average Bias of \(\alpha_x\) Under Various Scenarios where \(\tau\), \(\beta_x\), \(\beta_z\) are selected to induce a zero correlation between \(X\) and \(Z\) after selecting on \(Y_1\). Results are based on a sample size of \(n = 50000\) and 1000 samples obtained from the data-generating model described above.

\[\begin{array}{|c|c|c|c|c|c|c|cc|cc|} \hline a & b & \lambda_z \text{(max)} & \beta_x & \beta_z & \lambda_z & \tau & {\text{All Data}} && {\text{Select on } Y_1 > \tau} \\ & & & & & & & \alpha_{x,z} & \alpha_x & \alpha_{x,z} & \alpha_x \\ \hline 0.2 & 0.2 & 0.0422 & 0.2000 & 0.1952 & 0.0200 & -0.5774 & -0.0000 & 0.0118 & -0.0002 & -0.0002 \\ 0.2 & 0.2 & 0.0422 & 0.1999 & 0.1946 & 0.0350 & 1.3809 & -0.0001 & 0.0208 & -0.0004 & -0.0008 \\ 0.2 & 0.9 & 0.5519 & 0.1931 & 0.8361 & 0.2700 & 0.4647 & -0.0001 & 0.1620 & -0.0001 & -0.0001 \\ 0.2 & 0.9 & 0.5519 & 0.1857 & 0.8175 & 0.4000 & 1.8276 & -0.0001 & 0.2398 & -0.0003 & -0.0003 \\ 0.9 & 0.2 & 0.2277 & 0.8936 & 0.0683 & 0.1200 & 0.6717 & -0.0001 & 0.0721 & -0.0000 & 0.0009 \\ 0.9 & 0.2 & 0.2277 & 0.8842 & 0.0598 & 0.1900 & 3.0584 & -0.0001 & 0.1140 & -0.0000 & -0.0118 \\ 0.9 & 0.9 & 0.5156 & 0.8710 & 0.2383 & 0.2600 & -0.1627 & -0.0002 & 0.1556 & -0.0002 & -0.0005 \\ 0.9 & 0.9 & 0.5156 & 0.8207 & 0.1818 & 0.4500 & 2.3590 & -0.0002 & 0.2699 & 0.0011 & -0.0005 \\ \hline \end{array}\]

Misperception 11 - The method used to control for a covariate can be assumed to have been chosen appropriately and other methods would not, on average, produce substantially different results

Why This Misperception Occurs:

  • Surveys do not describe how aware authors are of the consequences of choosing models to achieve desired results.

  • Some evidence shows authors intentionally select covariates to achieve statistical significance.

  • Example: (Banks et al. 2015) found that management researchers sometimes chose covariates for significance.

Avoiding the Misperception:

  • Preregister analyses to prevent p-hacking and ensure transparency(Dal-Ré et al. 2014).

  • Disclose all model-building steps and decisions, including exploratory analyses, to allow readers to make informed judgments(Lenz and Sahn 2020).

  • Employ multiverse-style methods or specification curve analysis to explore the impact of covariate selection flexibility(Steegen et al. 2016).

Misperception 12 - p values derived from implementing statistical methods incorporating covariates mean exactly what they appear to mean and can be interpreted at face value

Why This Misperception Occurs:

  • Six Sigma Findings: Rare findings (six standard deviations from the null hypothesis) are frequently invalidated, highlighting issues with underlying assumptions(Linderman et al. 2002).

  • Assumption Violations: Measurement properties, data distributions, and test statistic performance assumptions are often not fully met, leading to inaccuracies(Greenland et al. 2016; Santosh Bangalore, Wang, and Allison 2009).

  • Sensitivity to Assumption Violations: Tests robust at conventional significance levels (e.g., 0.05) may become sensitive to even minor assumption violations at extremely low significance levels, such as those in Six Sigma findings.

Avoiding the Misperception:

  • Robustness at one significance level does not necessarily imply robustness at a different significance level.

  • Conduct independent replication studies to account for stochastic errors and detect biases from unknown factors.

Discussion

Considerations and Limitations:

  • List is not exhaustive; other useful references exist (Cinelli, Forney, and Pearl 2020; Gelman, Hill, and Vehtari 2020; Ding and Miratrix 2015)

  • Not all solutions discussed; Bayesian approaches offer alternative methods for addressing covariate selection and bias.

  • Use of DAGs and criteria from do-calculus framework for determining controlling variables.

  • Tools like dagitty.net (https://www.dagitty.net/) help specify DAGs and identify controlling variables.

  • Focus on estimating associations accurately and discussing their potential causal implications.

  • Practice intellectual humility by providing balanced considerations of whether associations represent causal effects.

Conclusion

  • Addressed 12 misperceptions involving covariate use, applicable to both linear and nonlinear models.

  • Aim to enhance understanding of effective bias control without introducing further biases.

  • Utilize tools like DAGs and Bayesian methods for better covariate selection and bias control.

  • Maintain intellectual humility and transparency in reporting and interpreting results.

Reference

Akheruzzaman, Md, Vijay Hegde, and Nikhil V. Dhurandhar. 2018. “Twenty-Five Years of Research about Adipogenic Adenoviruses: A Systematic Review.” Obesity Reviews 20 (4): 499–509. https://doi.org/10.1111/obr.12808.
Andersen, Robert. 2009. “Nonparametric Methods for Modeling Nonlinearity in Regression Analysis.” Annual Review of Sociology 35 (1): 67–85. https://doi.org/10.1146/annurev.soc.34.040507.134631.
Banks, George C., Ernest H. O’Boyle, Jeffrey M. Pollack, Charles D. White, John H. Batchelor, Christopher E. Whelpley, Kristie A. Abston, Andrew A. Bennett, and Cheryl L. Adkins. 2015. “Questions About Questionable Research Practices in the Field of Management.” Journal of Management 42 (1): 5–20. https://doi.org/10.1177/0149206315619011.
Bell, Carl C. 1994. “DSM-IV: Diagnostic and Statistical Manual of Mental Disorders.” JAMA: The Journal of the American Medical Association 272 (10): 828. https://doi.org/10.1001/jama.1994.03520100096046.
Blas Achic, Betsabé G., Tianying Wang, Ya Su, Victor Kipnis, Kevin Dodd, and Raymond J. Carroll. 2018. “Categorizing a Continuous Predictor Subject to Measurement Error.” Electronic Journal of Statistics 12 (2). https://doi.org/10.1214/18-ejs1489.
Carroll, Raymond J., David Ruppert, Leonard A. Stefanski, and Ciprian M. Crainiceanu. 2006. Measurement Error in Nonlinear Models. Chapman; Hall/CRC. https://doi.org/10.1201/9781420010138.
Cinelli, Carlos, Andrew Forney, and Judea Pearl. 2020. “A Crash Course in Good and Bad Controls.” SSRN Electronic Journal. https://doi.org/10.2139/ssrn.3689437.
Dal-Ré, Rafael, John P. Ioannidis, Michael B. Bracken, Patricia A. Buffler, An-Wen Chan, Eduardo L. Franco, Carlo La Vecchia, and Elisabete Weiderpass. 2014. “Making Prospective Registration of Observational Research a Reality.” Science Translational Medicine 6 (224). https://doi.org/10.1126/scitranslmed.3007513.
Ding, Peng, and Luke W. Miratrix. 2015. “To Adjust or Not to Adjust? Sensitivity Analysis of m-Bias and Butterfly-Bias.” Journal of Causal Inference 3 (1): 41–57. https://doi.org/10.1515/jci-2013-0021.
Fuller, Wayne A. 1987. “Measurement Error Models.” Wiley Series in Probability and Statistics, June. https://doi.org/10.1002/9780470316665.
Gelman, Andrew, Jennifer Hill, and Aki Vehtari. 2020. “Regression and Other Stories,” July. https://doi.org/10.1017/9781139161879.
Greenland, Sander, Stephen J. Senn, Kenneth J. Rothman, John B. Carlin, Charles Poole, Steven N. Goodman, and Douglas G. Altman. 2016. “Statistical Tests, P Values, Confidence Intervals, and Power: A Guide to Misinterpretations.” European Journal of Epidemiology 31 (4): 337–50. https://doi.org/10.1007/s10654-016-0149-3.
Hanley-Cook, Giles T., Aisling J. Daly, Roseline Remans, Andrew D. Jones, Kris A. Murray, Inge Huybrechts, Bernard De Baets, and Carl Lachat. 2022. “Food Biodiversity: Quantifying the Unquantifiable in Human Diets.” Critical Reviews in Food Science and Nutrition 63 (25): 7837–51. https://doi.org/10.1080/10408398.2022.2051163.
Indu, Pillaveetil Sathyadas, Thekkethayyil Viswanathan Anilkumar, Krishnapillai Vijayakumar, K. A. Kumar, P. Sankara Sarma, Saradamma Remadevi, and Chittaranjan Andrade. 2018. “Reliability and Validity of PHQ-9 When Administered by Health Workers for Depression Screening Among Women in Primary Care.” Asian Journal of Psychiatry 37 (October): 10–14. https://doi.org/10.1016/j.ajp.2018.07.021.
Lenz, Gabriel S., and Alexander Sahn. 2020. “Achieving Statistical Significance with Control Variables and Without Transparency.” Political Analysis 29 (3): 356–69. https://doi.org/10.1017/pan.2020.31.
Linderman, Kevin, Roger G. Schroeder, Srilata Zaheer, and Adrian S. Choo. 2002. “Six Sigma: A Goal-Theoretic Perspective.” Journal of Operations Management 21 (2): 193–203. https://doi.org/10.1016/s0272-6963(02)00087-6.
MacCallum, Robert C., Shaobo Zhang, Kristopher J. Preacher, and Derek D. Rucker. 2002a. “On the Practice of Dichotomization of Quantitative Variables.” Psychological Methods 7 (1): 19–40. https://doi.org/10.1037/1082-989x.7.1.19.
———. 2002b. “On the Practice of Dichotomization of Quantitative Variables.” Psychological Methods 7 (1): 19–40. https://doi.org/10.1037/1082-989x.7.1.19.
Munafò, Marcus R, Kate Tilling, Amy E Taylor, David M Evans, and George Davey Smith. 2017. “Collider Scope: When Selection Bias Can Substantially Influence Observed Associations.” International Journal of Epidemiology 47 (1): 226–35. https://doi.org/10.1093/ije/dyx206.
Rosenbaum, Paul R. 2002. Observational Studies. Springer New York. https://doi.org/10.1007/978-1-4757-3692-2.
Santosh Bangalore, Sai, Jelai Wang, and David B. Allison. 2009. “How Accurate Are the Extremely Small -Values Used in Genomic Research: An Evaluation of Numerical Libraries.” Computational Statistics & Data Analysis 53 (7): 2446–52. https://doi.org/10.1016/j.csda.2008.11.028.
Steegen, Sara, Francis Tuerlinckx, Andrew Gelman, and Wolf Vanpaemel. 2016. “Increasing Transparency Through a Multiverse Analysis.” Perspectives on Psychological Science 11 (5): 702–12. https://doi.org/10.1177/1745691616658637.
“The Book of Why: The New Science of Cause and Effect The Book of Why: The New Science of Cause and Effect Judea Pearl and Dana Mackenzie Basic Books, 2018. 429 Pp.” 2018. Science 361 (6405): 855–55. https://doi.org/10.1126/science.aau9731.
VanderWeele, Tyler J. 2019. “Principles of Confounder Selection.” European Journal of Epidemiology 34 (3): 211–19. https://doi.org/10.1007/s10654-019-00494-6.
VanderWeele, Tyler J., and Peng Ding. 2017. “Sensitivity Analysis in Observational Research: Introducing the E-Value.” Annals of Internal Medicine 167 (4): 268. https://doi.org/10.7326/m16-2607.
Westfall, Jacob, and Tal Yarkoni. 2016. “Statistically Controlling for Confounding Constructs Is Harder Than You Think.” Edited by Ulrich S Tran. PLOS ONE 11 (3): e0152719. https://doi.org/10.1371/journal.pone.0152719.
Yu, Xiaoxin, Roger S Zoh, David A Fluharty, Luis M Mestre, Danny Valdez, Carmen D Tekwe, Colby J Vorland, et al. 2024. “Misstatements, Misperceptions, and Mistakes in Controlling for Covariates in Observational Research.” eLife 13 (May). https://doi.org/10.7554/elife.82268.