SMHS PartialCorrelation - Revision history

Dinov: /* Problems and Exercises */

2025-12-09T19:27:40Z

‎Problems and Exercises

Dinov: /* Problems and Exercises */

2025-12-09T19:25:01Z

‎Problems and Exercises

Dinov: /* Problems and Exercises */

2025-12-09T19:20:28Z

‎Problems and Exercises

Dinov: /* Multicollinearity Among Controls */

2025-12-09T19:19:43Z

‎Multicollinearity Among Controls

Dinov: /* Common Issues and Solutions */

2025-12-09T19:16:23Z

‎Common Issues and Solutions

Dinov: /* Software Implementation */

2025-12-09T19:16:03Z

‎Software Implementation

Dinov: /* Software Implementation */

2025-12-09T19:15:28Z

‎Software Implementation

Dinov: /* Causal Inference and Partial Correlation */

2025-12-09T19:10:36Z

‎Causal Inference and Partial Correlation

Show changes

Dinov: /* Causal Inference and Partial Correlation */

2025-12-09T19:07:28Z

‎Causal Inference and Partial Correlation

Dinov: /* Causal Inference and Partial Correlation */

2025-12-09T19:04:19Z

‎Causal Inference and Partial Correlation

@@ Line 1,267: / Line 1,267: @@
    cat("Coverage of 95% CI (excluding failed CI):", coverage_rate, "\n")
    return(data.frame(
      estimate = estimates,
@@ Line 1,277: / Line 1,290: @@
 sim_results <- simulate_partial_correlation(n_sim = 500, n = 50,
                                            true_rho = 0.3, n_controls = 2)
 </pre>

@@ Line 1,211: / Line 1,211: @@
 . '''Applied Problems''':
-   a) Analyze the relationship between `mpg` and `qsec` in the `mtcars` dataset, controlling for `wt` and `hp`
-   b) Compute the partial correlation matrix for the `iris` dataset (sepal and petal measurements)
+: a) Analyze the relationship between `mpg` and `qsec` in the `mtcars` dataset, controlling for `wt` and `hp`
-   c) Conduct a power analysis for detecting partial correlation ρ = 0.3 with 3 control variables, α = 0.05, power = 0.80
+: b) Compute the partial correlation matrix for the `iris` dataset (sepal and petal measurements)
 . '''Simulation Study''':
@@ Line 1,223: / Line 1,224: @@
    estimates <- numeric(n_sim)
    p_values <- numeric(n_sim)
-   coverage <- numeric(n_sim)
+   coverage <- logical(n_sim)  # better as logical
    for (i in 1:n_sim) {
      # Generate control variables
      Z <- matrix(rnorm(n * n_controls), n, n_controls)
@@ Line 1,246: / Line 1,246: @@
      p_values[i] <- result$p.value
-     # Check coverage of 95% CI
+     # Check coverage of 95% CI — safely
      ci <- result$conf.int
-     coverage[i] <- (true_rho >= ci[1] & true_rho <= ci[2])
+     if (length(ci) >= 2) {
+   }
    cat("Simulation Results (n =", n_sim, "):\n")
@@ Line 1,257: / Line 1,265: @@
    cat("RMSE:", sqrt(mean((estimates - true_rho)^2)), "\n")
    cat("Type I error rate (for true_rho = 0):", mean(p_values < 0.05), "\n")
-   cat("Coverage of 95% CI:", mean(coverage), "\n")
+   cat("Coverage of 95% CI (excluding failed CI):", coverage_rate, "\n")
    return(data.frame(
@@ Line 1,269: / Line 1,277: @@
 sim_results <- simulate_partial_correlation(n_sim = 500, n = 50,
                                            true_rho = 0.3, n_controls = 2)
 </pre>

@@ Line 1,205: / Line 1,205: @@
 . '''Conceptual Problems''':
-   a) Prove that for jointly Gaussian variables, zero partial correlation implies conditional independence
-   b) Derive the relationship between partial correlation and regression coefficients
+: a) Prove that for jointly Gaussian variables, zero partial correlation implies conditional independence
-   c) Show that \(|r_{X(Y\cdot Z)}| \leq |r_{XY\cdot Z}| \leq |r_{XY}|\)
+: b) Derive the relationship between partial correlation and regression coefficients
 . '''Applied Problems''':

@@ Line 1,105: / Line 1,105: @@
 # Check for multicollinearity in control variables
 check_multicollinearity <- function(Z) {
-   if (is.null(dim(Z))) Z <- as.matrix(Z)
+  # Handle vector input
    # Variance Inflation Factor (VIF)
-   if (require(car)) {
+   if (requireNamespace("car", quietly = TRUE)) {
-     # Create a dummy regression to compute VIF
+     # Create a temporary data frame
-     dummy_y <- rnorm(nrow(Z))
+     df <- data.frame(dummy_y = rnorm(nrow(Z)), Z)
-     dummy_lm <- lm(dummy_y ~ Z)
-     vif_values <- vif(dummy_lm)
+    # Fit model using all columns in Z as predictors
      cat("Variance Inflation Factors (VIF):\n")
      print(vif_values)
-    # Rule of thumb: VIF > 10 indicates problematic multicollinearity
+     if (any(vif_values > 10, na.rm = TRUE)) {
        cat("\nWarning: High multicollinearity detected in control variables.\n")
        cat("Consider removing or combining variables with VIF > 10.\n")
@@ Line 1,127: / Line 1,136: @@
    # Condition number
    if (ncol(Z) > 1) {
-     kappa_value <- kappa(cor(Z), exact = TRUE)
+     # Use correlation matrix
      cat("\nCondition number of correlation matrix:", round(kappa_value, 2), "\n")
      if (kappa_value > 30) {
@@ Line 1,139: / Line 1,151: @@
 </pre>
-####2) Missing Data Handling
+====Missing Data Handling====
 <pre>
 # Multiple approaches for missing data in partial correlation

@@ Line 1,100: / Line 1,100: @@
 ===Common Issues and Solutions===
-####1) Multicollinearity Among Controls
+====Multicollinearity Among Controls====
 <pre>
 # Check for multicollinearity in control variables

@@ Line 906: / Line 906: @@
 <pre>
 # Comprehensive partial correlation analysis function
 run_partial_correlation_analysis <- function(data, var1, var2, controls,
                                              method = "standard",

@@ Line 906: / Line 906: @@
 <pre>
 # Comprehensive partial correlation analysis function
 run_partial_correlation_analysis <- function(data, var1, var2, controls,
                                              method = "standard",
@@ Line 931: / Line 932: @@
        pcor_result <- pcor.test(X, Y, Z)
        cat("\n2. Partial correlation:", round(pcor_result$estimate, 4), "\n")
        if (inference) {
          cat("   p-value:", format.pval(pcor_result$p.value, digits = 4), "\n")
          cat("   95% CI: [",
-             round(pcor_result$conf.int[1], 4), ", ",
+             round(r_lower, 4), ", ",
-             round(pcor_result$conf.int[2], 4), "]\n", sep = "")
+             round(r_upper, 4), "]\n", sep = "")
+       }
      } else {
@@ Line 964: / Line 984: @@
    # 4. Visualization
    if (visualization) {
      par(mfrow = c(2, 2))
@@ Line 973: / Line 997: @@
      if (!is.null(Z)) {
        # Residuals after controlling
-       res_X <- resid(lm(as.formula(paste(var1, "~", paste(controls, collapse = "+"))),
+       formula_x <- as.formula(paste(var1, "~", paste(controls, collapse = "+")))
-                       data = data))
+      res_X <- resid(lm(formula_x, data = data))
        res_Y <- resid(lm(formula_y, data = data))
@@ Line 993: / Line 1,017: @@
+       )
-       barplot(cor_types$Value, names.arg = cor_types$Type,
+      # Remove NA values for plotting
-               ylim = c(-1, 1), col = c("blue", "red", "green"),
+      cor_types_complete <- cor_types[!is.na(cor_types$Value), ]
                main = "Comparison of Correlation Measures",
                ylab = "Correlation Coefficient")
@@ Line 1,000: / Line 1,029: @@
+     }
-     # Network visualization if multiple variables
+     # Network visualization if multiple variables (with fallback)
      if (ncol(data) >= 3) {
-       pcor_matrix <- pcor(data)$estimate
+       tryCatch({
-      qgraph(pcor_matrix, layout = "spring",
+        pcor_matrix <- pcor(data)$estimate
-             labels = colnames(data),
-             title = "Partial Correlation Network",
+        # Try different visualization methods
-             maximum = 0.8)
+        if (require(corrplot, quietly = TRUE)) {
+     }
+   }
@@ Line 1,019: / Line 1,074: @@
      p_value = if(exists("pcor_result")) pcor_result$p.value else NA
+   )
    return(results)
+}
 # Example usage with mtcars dataset

@@ Line 707: / Line 707: @@
 <pre>
 # NEED TO INSTALL: R package pcalg and BiocManager::install(c("graph", "RBGL", "Rgraphviz"))
 # Partial correlation in causal inference