Advanced Statistical Data Analysis

Advanced Regression Modelling¶

Best Practices and Key Takeaways¶

1. Always check model assumptions using diagnostic plots (Tukey-Anscombe, Normal Q-Q, Scale-Location, Residual vs Leverage). The enhanced diagnostic plots with simulation envelopes (plot.glmSim()) provide rigorous checks against model assumptions.

2. Transform variables when appropriate using Tukey’s First-Aid Transformations (log for positive continuous, square-root for counts, logit for proportions). These often lead to simpler, more interpretable models.

3. Use AIC for model selection, not just p-values from marginal tests. Remember that AIC is for model selection, not hypothesis testing.

4. Check for overdispersion in count data and use quasi-models if needed. The quasi-Poisson and quasi-binomial families account for extra variation.

5. Validate models using cross-validation, PRESS, or other validation techniques. The cv.glm() function provides k-fold cross-validation.

6. Interpret coefficients in the context of the link function (e.g., odds ratios for logistic, multiplicative effects for log link). For logistic regression, always predict on the link scale and transform back to ensure valid confidence intervals.

7. Make predictions on the appropriate scale (response vs. link). For GLMs with non-linear link functions, the link scale is often more stable for prediction and confidence intervals.

Foundations of Regression Analysis¶

Regression modelling is the art of explaining or modelling the relationship between a response variable $Y$ and one or more explanatory variables (predictors) $x^{(1)}, \ldots, x^{(m)}$, where the relationship is masked by random noise. The goal is to identify a function $f\langle x^{(1)}, \ldots, x^{(m)}\rangle$ such that $Y \approx f\langle x^{(1)}, \ldots, x^{(m)}\rangle + E$, where $E$ represents random variation.

Multiple Linear Regression Model:

Where:

• $Y_i$ is the response (must be continuous)

• $x_i^{(j)}$ are predictors (can be continuous, discrete, binary, or categorical)

• $E_i \sim \mathcal{N}(0, \sigma^2)$ are independent error terms

• $\beta_0, \beta_1, \ldots, \beta_m$ are unknown regression coefficients

Matrix Form: $\underline{Y} = \boldsymbol{X} \underline{\beta} + \underline{E}$

Model Assumptions:

1. Linearity: The relationship between predictors and response is linear in parameters

2. Independence: Errors are stochastically independent

3. Homoscedasticity: Constant variance $\sigma^2$ across all levels of predictors

4. Normality: Errors are normally distributed

5. No Perfect Multicollinearity: Predictors are not linearly dependent

Data Transformation: Tukey’s First-Aid Transformations¶

When the relationship between variables is non-linear or the distribution is skewed, transformations can often linearize the relationship and stabilize variance. Mosteller and Tukey (1977) recommend these first-aid transformations as a starting point:

Why Transform?

• Linearize relationships: Makes non-linear patterns linear

• Stabilize variance: Addresses heteroscedasticity

• Normalize residuals: Helps meet normality assumption

• Improve interpretability: Coefficients represent multiplicative effects

Model Fitting and Inference¶

Least Squares Estimation: The ordinary least squares (OLS) estimator minimizes the sum of squared residuals:

This is equivalent to the maximum likelihood estimator when errors are normally distributed.

Variable Selection with AIC¶

When many predictors appear non-significant, model selection helps find the simplest adequate model. Akaike’s Information Criterion (AIC) balances model fit and complexity:

Where $p^\diamond$ is the total number of estimated parameters (including $\hat{\sigma}$).

Interpretation: Lower AIC = better model (rewards fit, penalizes complexity)

Handling Categorical Predictors¶

Categorical variables (factors) require special handling in regression models.

Dummy Coding (Treatment Contrasts):

• Creates k-1 binary variables for a factor with k levels

• Each coefficient represents the difference from the reference level

Other Contrast Types:

• Sum-to-zero: contr.sum() - coefficients sum to zero

• Helmert: contr.helmert() - compares each level to mean of subsequent levels

• Polynomial: contr.poly() - for ordered factors

Detecting Multicollinearity¶

Multicollinearity occurs when predictors are highly correlated, making it difficult to estimate individual coefficients reliably.

Variance Inflation Factor (VIF):

Where $R_j^2$ is the coefficient of determination from regressing predictor j on all other predictors.

Interpretation:

• VIF ≈ 1: No collinearity

• VIF > 5: Moderate collinearity (concerning)

• VIF > 10: Severe collinearity (problematic)

Solutions for Multicollinearity:

1. Remove one of the correlated predictors (if theoretically justified)

2. Combine predictors (e.g., create an index)

3. Use ridge regression or other regularization methods

4. Accept it if the goal is prediction, not inference

Model Validation Techniques¶

Cross-Validation:

• Split data into training and test sets

• Fit model on training, evaluate on test

• k-fold CV: Split into k folds, use each as test set once

PRESS (Predicted Residual Sum of Squares):

Where $\hat{y}_{(i)}$ is the predicted value when observation i is left out.

Advanced Topics in Linear Regression¶

Weighted Least Squares (WLS)¶

When variance is not constant (heteroscedasticity), WLS assigns different weights to observations:

Where $w_i = 1/\text{Var}(E_i)$ if the variance structure is known.

Robust Fitting¶

Robust regression reduces the influence of outliers and leverage points using robust loss functions:

• Huber: Less sensitive to outliers than squared loss

• Tukey bisquare: Even more robust, completely ignores extreme outliers

• Huber with M-estimation: Iteratively reweighted least squares

Fitting Smooth Functions¶

Non-parametric smoothing captures non-linear relationships without specifying a functional form:

LOESS (Locally Estimated Scatterplot Smoothing):

• Fits local polynomials using weighted regression

• Span parameter: Controls smoothness (larger = smoother)

Splines:

• Natural cubic splines: Smooth, flexible functions

• B-splines: Local support, computationally efficient

Additive Models¶

Additive models extend linear regression by allowing non-linear relationships:

Where $f_j$ are smooth functions (can be linear, splines, LOESS, etc.).

Advantages:

• Flexible non-linear relationships

• Additive structure (easy to interpret each term)

• Automatic smoothness selection

Binary Response: Logistic Regression¶

When the response is binary (success/failure, yes/no), ordinary linear regression is inappropriate because:

• Predictions can fall outside [0, 1]

• The linear model assumes constant variance, but binary data has variance that depends on the mean

• The distribution is binomial, not normal

Logistic Regression Model

Where $\eta_i = \beta_0 + \beta_1 x_i^{(1)} + \cdots + \beta_m x_i^{(m)}$ is the linear predictor.

Logit Link Function

Interpretation of Coefficients

• $\beta_j$ = change in log-odds per unit change in $x^{(j)}$

• $\exp(\beta_j)$ = odds ratio (multiplicative change in odds)

• If $\beta_j > 0$: Higher $x^{(j)}$ increases probability of success

• If $\beta_j < 0$: Higher $x^{(j)}$ decreases probability of success

Generalized Linear Models (GLM)¶

GLMs unify various regression models under a single framework (Nelder & Wedderburn, 1972). They consist of three components:

1. Random Component: Response $Y_i$ follows a distribution from the exponential family

2. Systematic Component: Linear predictor $\eta_i = \sum_{j=0}^m \beta_j x_i^{(j)}$

3. Link Function: $g(\mu_i) = \eta_i$, where $\mu_i = E[Y_i]$

Exponential Family¶

• The dispersion parameter $\phi$ is connected to the variance

• The value $w_i$ is a fixed known number (weighting)

• The function $d(y_i; \phi, w_i)$ scales the probability to 1.

• The functions $b(\mu_i)$ and $c(\mu_i)$ determine the distribution.

Expected Value:

where $b'(\mu)$ and $c'(\mu)$ are the first derivatives with respect to $\mu$.

Variance:

• Distribution: The probability distribution of the response variable $Y$;Defines the random component of the GLM

• Range of $Y$: The support (set of possible values) for the response variable; Determines which distribution is appropriate for your data

• $\mathbb{E}\langle Y\rangle = \mu$: The expected value (mean) of the response; The target of your regression model

• $\text{var}\langle Y\rangle$: The variance of the response variable; Determines precision of estimates

• $b\langle \mu\rangle$: The canonical link function from the exponential family density; $\eta = b'(\mu)$

• $V\langle \mu\rangle$: The variance function; Connects mean to variance: $var(Y) = \frac{\phi}{w} V(\mu)$

• $\phi$: The dispersion parameter; Scales the variance: larger $\phi$ means more spread

• $w$: The prior weight for each observation; Known constant for each observation (e.g., $m$ for binomial)

• Temperature coefficient: -0.1077 (p = 0.0158)

  • Interpretation: Each 1F decrease in temperature multiplies the odds of failure by $\exp(0.1077) \approx 1.114$ (11.4% increase)

  • Engineering implication: At the launch temperature of 31F (well below previous launches), the failure probability was unacceptably high

• Pressure coefficient: 0.0074 (p = 0.2523)

  • Not statistically significant at 5% level

Hypothesis Test:

Does pressure affect failure probability?

$H_0: \beta_1 = 0$ (no effect), $H_1: \beta_1 \neq 0$ (effect exists)

Deviance test (more reliable than Wald for small samples):

Conclusion: Cannot reject $H_0$ → No evidence that pressure affects failure probability

Confidence Interval for Temperature:

• 95% CI for temperature coefficient: [-0.201, -0.022]

• Does not include 0 → confirms significant effect

• Principle: confint() uses deviance test statistic (more reliable than Wald)

Prediction at 31F: The critical question: What was the predicted failure probability at the launch temperature of 31F?

Data Selection Lesson: Engineers considered using only data where O-rings showed thermal distress.

• Problem: This drops observations with no distress, which contain valuable information

• Principle: Never drop data based on the response variable - you never know if it influences the outcome

• Solution: Use all available data for analysis

Inference for GLMs¶

Deviances¶

Deviance measures the goodness of fit of a model:

• Null deviance: Measures total variation in the response (model with only intercept)

• Residual deviance: Measures variation not explained by the model

• Difference in deviances: Measures improvement from adding predictors

For nested models, the difference in residual deviances follows a $\chi^2$ distribution:

Where df = difference in number of parameters.

Overdispersion¶

Overdispersion occurs when the observed variance exceeds what the model assumes:

• Poisson: Variance > mean

• Binomial: Variance > $n\pi(1-\pi)$

Detection:

If $\hat{\phi} > 1$, there is evidence of overdispersion.

Example - Poisson GLM with Overdispersion

The warpbreaks dataset contains the number of breaks in yarn during weaving, classified by wool type (A, B) and tension (L, M, H). This is a classic example where overdispersion is present.

Model:

R Output:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.7967     0.0499  76.030  <2e-16 ***
woolB      -0.4566     0.0802  -5.694  1.24e-08 ***
tensionM   -0.6187     0.0844  -7.330  2.30e-13 ***
tensionH   -0.5958     0.0838  -7.112  1.15e-12 ***
woolB:tensionM 0.6382    0.1222   5.224  1.75e-07 ***
woolB:tensionH 0.1884    0.1299   1.450   0.147

Interaction Term Significance:

• Null: $\beta_4 = \beta_5 = 0$ (no interaction)

• Alternative: At least one interaction coefficient ≠ 0

• Test: Deviance test (drop1 with Chisq test)

• R Code:

# Test interaction term
drop1(WB.glm, test="Chisq")
# Result: p-value = 7.962e-07 < 0.05

• Conclusion: Interaction term is needed to explain the response

Overdispersion Detection:

• Definition: More dispersion in the data than the Poisson model allows (Poisson assumes mean = variance)

• Detection: Residual deviance / residual df = 182.31 / 48 = 3.798

• Test: 1 - pchisq(182.31, 48) = 0 (essentially 0) < 0.05

• Conclusion: Strong evidence of overdispersion

Implications:

• The Poisson model is not adequate to describe the data

• Statistical conclusions based on this model (including the interaction test) are not valid

• Solution: Use a quasi-Poisson model to account for overdispersion

Modified Model:

Confidence Intervals for GLMs¶

Wald-type Confidence Intervals¶

!! Assumes normal sampling distribution !!

• Based on asymptotic normality of MLEs

• Simple to compute but can be inaccurate for small samples or sparse data

• May produce invalid ranges (e.g., probabilities outside [0,1])

Deviance-based Confidence Intervals¶

!! For GLMs, always prefer this likelihood-based approach over Wald by hand - it’s more accurate !!

• Based on likelihood ratio tests

• More reliable for small samples

• Guarantees intervals stay within valid range

Diagnostics and Model Adequacy Checking¶

Model adequacy is crucial - a model is only as good as its assumptions. Diagnostic plots help verify these assumptions.

Essential Diagnostic Plots for GLMs¶

1. Tukey-Anscombe Plot (Residuals vs Fitted):

   • Check: Linear relationship, constant variance

   • Problem: Non-linear pattern → missing non-linear term

   • Problem: Funnel shape → heteroscedasticity

2. Normal Q-Q Plot (Quantiles of residuals vs theoretical):

   • Check: Normality of residuals

   • Problem: Points deviate from line → non-normal residuals

   • Problem: Heavy tails → outliers

3. Scale-Location Plot ($\sqrt{|\text{Residuals}|}$ vs Fitted):

   • Check: Constant variance

   • Problem: Non-constant spread → heteroscedasticity

4. Residual vs Leverage:

   • Check: Influential observations

   • Problem: Cook’s distance > 1 → too influential

   • Problem: High leverage points → extrapolation

Enhanced Diagnostics with Simulation¶

Simulated envelopes provide a more rigorous check by comparing observed diagnostic curves to what we’d expect under the model assumptions:

Interpretation:

• If observed curves fall within the gray envelope (simulated fluctuation), there is no evidence against the assumption

• If observed curves fall outside the envelope, there is evidence the assumption is violated

Example - Machine Failures Residual Analysis:

A machinery runs in two modes. The question is whether the number of failures depends on operating time in each mode. After fitting a Poisson regression model, we perform a residual and sensitivity analysis.

Diagnostic Plots Interpretation:

1. Tukey-Anscombe Plot: Smoother indicates some non-constant expectation, but within stochastic fluctuation (gray spaghetti) → Weak evidence against adequately specified expectation

2. Location-Scale Plot: Smoother decreasing but within stochastic fluctuation → No evidence against adequately specified variance

3. Normal Q-Q Plot: Points scatter around line, all black points within gray points → No evidence for outliers or distributional violations

4. Residual vs Leverage: All Cook’s distances < 1 → No observation with too large influence

Overall Conclusion: Model might be adequate, though there’s a weak hint that expectation might not be constant (but limited data prevents strong conclusions).

Statistical Soundness: Because there is no serious evidence against model assumptions, the statistical conclusions from the hypothesis test and overdispersion check are valid.

Example - Residual Analysis with Simulation:

The six diagnostic plots (Tukey-Anscombe, Normal Q-Q, Scale-Location) with simulation envelopes (gray lines and points) provide additional insights.

Purpose of Gray Lines and Points:

• Gray lines: Show the scatter range of the smoother under model assumptions

• Gray points: Show the scatter range of data points under model assumptions

• Interpretation: If black (observed) structures lie outside the gray (simulated) structures, this is evidence the model does not adequately describe the data

Residual and Sensitivity Analysis for Yarn Breaks:

1. Tukey-Anscombe Plot:

   • Observation: Smoother is a reasonably straight horizontal line

   • Comparison: Within stochastic fluctuation (gray spaghetti)

   • Conclusion: No evidence that expectation is incorrectly specified

2. Normal Q-Q Plot:

   • Observation: Black points are not within stochastic fluctuation (gray points)

   • Interpretation: Evidence of distortions of the distributional assumption

   • Art of violation: Suggests overdispersion

3. Scale-Location Plot:

   • Observation: Red smoother has a weak trend, shifted upwards

   • Comparison: Outside stochastic fluctuation

   • Conclusion: Evidence that variance is incorrectly specified

   • Art of violation: Suggests overdispersion

Overall Conclusion: There is evidence that the model does not fit the data adequately. The identified violations suggest overdispersion is the main issue, which we addressed by using a quasi-Poisson model.

Residual Types for GLMs:

• Response residuals: $y_i - \hat{y}_i$ (observed - predicted)

• Pearson residuals: $(y_i - \hat{y}_i) / \sqrt{\text{Var}(\hat{y}_i)}$ (standardized)

• Deviance residuals: Contribute to deviance statistic

Model Improvement: Finding Appropriate Transformations¶

Generalized Additive Models (GAMs) are excellent for exploring non-linear relationships in a data-driven manner:

• Rule of Thumb: Predictor needs a transformation if a straight line does not fit between the confidence band or the P-value of the coresponding coefficient in Anova for Nonparametric Effects is smaller than 0.05

• lMaxWt needs an additional transformation; Hockeystick with L at 8

• Temp needs a transformation as well: possible a factor variable with iTemp=0 if Temp < 17, iTemp=1 otherwise

• lAR: no additional transformation is needed

• Food is a factor variable and cannot be transformed in any case

Extensions of GLMs¶

Rate Models¶

Rate models analyze counts per unit exposure (e.g., accidents per kilometer driven, failures per machine-hour):

Where $t_i$ is the exposure (offset).

Interpretation: The offset ensures we’re modeling the rate (counts per unit exposure) rather than raw counts.

Quasi-Likelihood Models¶

Quasi-models extend GLMs when the full distribution is not specified, requiring only the mean-variance relationship:

Quasi-Poisson:

Where $\phi$ is the dispersion parameter estimated from data.

Causality¶

Best Practices and Key Takeaways¶

1. Distinguish association from causation - correlation does not imply causation. The shoe size and salary example demonstrates how a common cause (gender) can create spurious associations.

2. Use DAGs to represent causal assumptions explicitly. Nodes are variables, directed edges are potential causal effects. Use the gRbase and gRain packages to work with DAGs and Bayesian networks.

3. Identify confounders and adjust for them appropriately using the backdoor criterion. The Nutrition-Disease example shows how to find valid adjustment sets.

4. Use d-separation to determine conditional independencies in the graph. The sprinkler network example illustrates the three basic structures (chain, fork, collider) and their independence properties.

5. Apply the backdoor criterion to find valid adjustment sets for causal effect estimation. This allows estimating causal effects from observational data when appropriate variables are measured.

6. Consider counterfactuals for individual-level causal questions using Pearl’s three-step method (abduction, action, prediction). The online sales example demonstrates practical counterfactual reasoning.

7. Validate causal models through sensitivity analysis, robustness checks, and external validation. Always question whether the causal assumptions are reasonable and test the stability of results.

Why Causal Inference Matters¶

Traditional statistics and machine learning excel at finding associations and making predictions, but they cannot answer causal questions without additional assumptions. Causal inference is essential when we want to understand what would happen if we changed something.

Key Distinction:

• Association: $P(Y | X)$ - What is the probability of Y given X? Do people who drink more coffee tend to be more productive?

• Intervention: $P(Y | do(X = x))$ - What would Y be if we set X to x? If I drink 20 cups of coffee, will I be more productive?

• Counterfactual: $P(Y^{X=x} | X=x', Y=y')$ - What would Y have been if X had been x, given that X was actually x’ and Y was y’? Would I have been more productive yesterday if I had drunk 20 cups of coffee instead of 2?

The first question can be answered with observational data, but the second and third require causal assumptions.

Causal Graphical Models¶

Causal graphical models (DAGs - Directed Acyclic Graphs) provide a language for representing causal assumptions explicitly.

Components:

• Nodes: Random variables

• Directed edges: Potential causal effects

• Absence of edge: No direct causal effect

Terminology:

• Parent: Node with directed edge to another node

• Child: Node with directed edge from another node

• Ancestor: Node with directed path to another node

• Descendant: Node with directed path from another node

• Adjacent: Nodes connected by an edge

Example - Bayesian Networks: Building and querying Bayesian networks using the gRain package:

Graphical Building Blocks:

1. Chain (Mediation): $X \rightarrow Z \rightarrow Y$

   • X affects Y through Z

   • Example: Education → Income → Happiness

2. Fork (Confounding): $X \leftarrow Z \rightarrow Y$

   • Z is a common cause of X and Y

   • Example: Gender → (Shoe Size, Salary)

3. Collider (Selection): $X \rightarrow Z \leftarrow Y$

   • Z is a common effect of X and Y

   • Example: Rain → Wet ← Sprinkler

Conditional Independence and D-Separation¶

Conditional Independence¶

Two variables X and Y are conditionally independent given Z if:

D-Separation: A path is d-separated (blocked) by a set of nodes Z if it contains a node that is:

• Not in Z, and

• Is a collider on that path

Two nodes X and Y are d-separated given Z if all paths between them are blocked by Z.

D-Separation Rules¶

DAG Analysis¶

Descendants of Place of Residence¶

• Answer: Medical Care, Nutrition, Disease

• Explanation: Descendants are all nodes reachable via directed paths from Place of Residence. From Place of Residence, we can reach Nutrition directly, Medical Care directly, and Disease through Medical Care.

Direct Causal Effect of Place of Residence on Disease¶

• Answer: 0

• Justification: There is no direct edge from Place of Residence to Disease in the DAG. The effect is only indirect through Medical Care.

Independence Statements¶

i. Nutrition $\perp\!\!\!\perp$ Medical Care | Place of Residence ?

• Correct: There is a fork at Place of Residence (common cause of both Nutrition and Medical Care)

• Reason: Conditioning on the common cause (Place of Residence) blocks the path between Nutrition and Medical Care, making them conditionally independent

ii. Nutrition $\perp\!\!\!\perp$ Physical Activity | Disease ?

• Incorrect: There is a collider at Disease (common effect of Nutrition and Physical Activity)

• Reason: Conditioning on the collider (Disease) opens the path between Nutrition and Physical Activity, making them conditionally dependent

iii. Nutrition $\perp\!\!\!\perp$ Genetics | Medical Care ?

• Correct: The collider is Disease (not in conditioning set)

• Reason: Medical Care is not a descendant of Disease. The collider path (Nutrition $\rightarrow$ Disease $\leftarrow$ Genetics) remains blocked because Disease is not in the conditioning set and Medical Care is not its descendant.

Regression for Total Causal Effect of Nutrition on Disease¶

To estimate the total causal effect of Nutrition on Disease, we must block the backdoor path: Nutrition $\leftarrow$ Place of Residence $\rightarrow$ Medical Care $\rightarrow$ Disease

lm(Disease ~ Nutrition + PlaceOfResidence, data = data)

Asymptotically Optimal Adjustment¶

The asymptotically optimal adjustment set minimizes the variance of the causal effect estimator:

• $\text{cn}(X, Y)$ = all causal nodes in the directed path from X to Y

• $\text{pa}(S)$ = parents of set S

• $\text{forb}(X, Y)$ = forbidden nodes (all descendants of causal nodes including X)

Calculation:

• $\text{cn}(Nutrition, Disease) = \{Disease\}$ (Disease is the only causal neighbor)

• $\text{pa}(\text{cn}(Nutrition, Disease)) = \{Nutrition, Genetics, Medical Care, Physical Activity\}$ (parents of Disease)

• $\text{forb}(Nutrition, Disease) = \{Nutrition, Disease\}$ (the variables themselves)

Optimal adjustment: $\{ Genetics, Medical Care, Physical Activity \}$

D-Separation Testing¶

Testing conditional independence in DAGs and verifying with data:

Simpson’s Paradox¶

Simpson’s Paradox occurs when a trend appears in different groups of data but disappears or reverses when these groups are combined.

Example: Layoffs

• Overall: More layoffs in Department A than Department B

• Within each gender: More layoffs in Department B than Department A

• Resolution: Gender distribution differs between departments (more men in A, more women in B, and men are laid off more regardless of department)

Causal Explanation: The paradox arises because gender is a confounder - it affects both the department (exposure) and layoffs (outcome). When we aggregate across genders, the confounding effect creates the paradox.

Solution: Adjust for the confounder (gender) to see the true relationship within each gender group.

Factorization of Joint Distributions¶

For a DAG, the joint distribution factorizes according to the graph structure:

Where $\text{pa}(X_i)$ are the parents of node $X_i$ in the DAG.

Example: Medical Diagnosis

Where:

• D = Disease

• S = Symptom

• C = Test Result

Interpretation: Disease affects both Symptom and Test Result, but Symptom and Test Result are independent given Disease (d-separated by Disease).

Causal Effect Estimation¶

Causal effect measures how an intervention on X affects Y:

Adjustment Formula: If we can find a set Z that satisfies the backdoor criterion, we can estimate causal effects from observational data:

Backdoor Criterion: A set of variables Z satisfies the backdoor criterion relative to X and Y if:

1. Z blocks all backdoor paths from X to Y

2. Z does not contain any descendants of X

Finding Adjustment Sets:

1. Identify all backdoor paths from X to Y

2. For each path, find a set that blocks it

3. Combine these sets, ensuring no new paths are opened

Example - Adjustment Formula:

Estimating causal effects with adjustment: Use backdoor criterion and adjust for all parents of Nutrition that are not descendants.

lm(Disease ~ Nutrition + PlaceOfResidence, data = data)

Structural Causal Models (SCMs)¶

Structural Causal Models combine causal graphs with structural equations that describe the data-generating process:

Where:

• $X_j$ = endogenous variable

• $\text{pa}(X_j)$ = parents of $X_j$ in the causal graph

• $E_j$ = exogenous (noise) variable

• $f_j$ = structural function

Linear SCMs:

Direct vs. Total Causal Effects:

• Direct effect: Effect of X on Y not mediated through other variables

• Total effect: Sum of direct effect + all indirect effects

Example: Linear SCM

• Direct effect of X on Y: $\beta_{YX}$

• Direct effect of X on Z: $\beta_{ZX}$

• Direct effect of Y on Z: $\beta_{ZY}$

• Total effect of X on Z: $\beta_{ZX} + \beta_{ZY} \cdot \beta_{YX}$

Estimating Causal Effects: For linear SCMs, causal effects can be estimated using regression with appropriate adjustment:

# Estimate total causal effect of Nutrition on Disease
# Adjust for confounding variables (Place of Residence)
lm(Disease ~ Nutrition + PlaceOfResidence, data = data)

Assumptions of Linear SCMs:

1. Known causal structure

2. Linear relationships

3. Uncorrelated errors

4. No unmeasured confounders

Limitations:

• Requires prior knowledge of causal structure

• Linear assumption may not hold

• Errors may be correlated in practice

Example - Linear SCMs: Defining and working with linear structural causal models:

Interventional vs. Observational Data¶

Experimental Data:

• Researchers randomly assign treatments

• Groups are comparable (except for treatment)

• Gold standard for causal inference

• Example: Randomized controlled trial (RCT)

Observational Data:

• Researchers observe naturally occurring variation

• Groups may differ systematically (confounding)

• Requires adjustment for valid causal inference

• Example: Cohort study, case-control study

Key Difference:

• Experimental: $P(Y | X = x) = P(Y | do(X = x))$ (randomization blocks confounders)

• Observational: $P(Y | X = x) \neq P(Y | do(X = x))$ (confounding may be present)

Example: Shoe Size and Salary

• Observational: People with larger shoes tend to have higher salaries

• Reason: Gender affects both shoe size and salary (confounding)

• Experimental: If we randomly assign shoe sizes, no effect on salary

• Conclusion: No causal effect of shoe size on salary

Example - Ice Cream, Temperature, Emergencies:

Context: Observational data from 146 cities recording temperature (moderate/high), emergency calls (low/high), and ice cream sales (low/high).

Causal Graph:

Why Classical Model Finds Relationship:

• Observation: Ice cream sales and emergency calls are positively associated

• Reason: Temperature is a confounding variable - it affects both ice cream sales and emergency calls

• Spurious correlation: The association is due to the common cause, not a direct causal effect

Conditional Probabilities:

# P(nEmergencies = high | IceCream = high)
(14 + 2) / 103  # 0.1553

# P(nEmergencies = high | IceCream = low)
(4 + 1) / 43   # 0.1163

# Risk ratio
0.1553 / 0.1163  # 1.34 times higher

Interventional Probabilities:

do(IceCream = high) meaning: An intervention that sets ice cream sales to high for everyone, regardless of temperature or other factors.

Difference between P(Y|X) and P(Y|do(X)):

• P(Y|X): Conditional probability - likelihood of Y among those who happened to have X

• P(Y|do(X)): Interventional probability - likelihood of Y if we set X to a specific value

• In most cases: P(Y|X) ≠ P(Y|do(X)) due to confounding

How the Graph Changes with do(IceCream = high):

• The arrow from Temperature → IceCream is removed (intervention makes IceCream independent of Temperature)

• Temperature still affects Emergency Calls directly

Calculating P(nEmergency = high | do(IceCream = high)): Using the adjustment formula:

# P(Temp = high) = 88/146 = 0.6027
# P(Temp = moderate) = 58/146 = 0.3973

# P(nEmergency = high | IceCream = high, Temp = high)
14 / (56 + 14)  # 0.2000

# P(nEmergency = high | IceCream = high, Temp = moderate)
2 / (31 + 2)    # 0.0606

# Interventional probability
0.2000 * 0.6027 + 0.0606 * 0.3973  # 0.1446

Calculating P(nEmergency = high | do(IceCream = low)):

# P(nEmergency = high | IceCream = low, Temp = high)
4 / (14 + 4)  # 0.2000

# P(nEmergency = high | IceCream = low, Temp = moderate)
1 / (24 + 1)  # 0.0400

# Interventional probability
0.2000 * 0.6027 + 0.0400 * 0.3973  # 0.1498

Causal Effect of IceCream on Emergency:

0.1446 - 0.1498  # -0.0052

Conclusion: The estimated causal effect is approximately 0, suggesting no causal effect of ice cream sales on emergency calls. The observed association (1.34 times higher risk) is entirely due to confounding by temperature. This demonstrates why it’s crucial to distinguish between association and causation.

Instrumental Variables (IV)¶

When confounders are unmeasured, standard regression cannot estimate causal effects. Instrumental variables provide a solution.

Valid Instrument Z must satisfy:

1. Relevance: Z is associated with X (not too weak)

2. Exclusion restriction: Z has no direct effect on Y (only through X)

3. Independence: Z is independent of all confounders of X and Y

Two-Stage Least Squares (2SLS):

1. Regress X on Z to get predicted values $\hat{X}$

2. Regress Y on $\hat{X}$ to estimate causal effect

library(AER)

# 2SLS using ivreg
iv.model <- ivreg(Y ~ X + W | Z + W, data = data)
summary(iv.model)

Example: Education and Earnings

• X: Education (years)

• Y: Earnings

• Z: Quarter of birth (affects education due to school entry age, but not earnings directly)

• Assumption: Quarter of birth is independent of ability, family background, etc.

Counterfactual Reasoning¶

Counterfactuals answer questions like: “What would have happened if...?” for specific individuals.

Pearl’s Three-Step Method:

1. Abduction: Use observed data to infer the distribution of exogenous variables

   • Estimate the structural error terms from observed data

2. Action: Modify the structural equations to reflect the intervention

   • Set the variable of interest to its counterfactual value

3. Prediction: Use the modified model to predict the counterfactual outcome

   • Compute the new value based on the modified equations

Example - Online Sales Counterfactual:

Given observed: Visits = 179, Discount = 19%, Sales = 101 Question: What would sales be if Discount = 10%?

Estimated SCM from LiNGAM: Using the LiNGAM algorithm on historical sales data:

library(pcalg)
fit.sales <- lingam(dat_sales)
B <- round(fit.sales$Bpruned, 2)
colnames(B) <- rownames(B) <- colnames(dat_sales)
B
#        sales adv visits discount
# sales     0 0.00   0.24    0.18
# adv       0 0.00   0.00    0.00
# visits    0 1.79   0.00    0.00
# discount  0 0.00   0.00    0.00

Estimated DAG:

• Adv → Visits (coefficient: 1.79)

• Visits → Sales (coefficient: 0.24)

• Discount → Sales (coefficient: 0.18)

• No direct effect of Adv on Sales (only indirect through Visits)

Step 1 - Abduction:

# Estimate structural error for Sales
ES <- 101 - (55.65 + 0.24 * 179 + 0.18 * 19)
# ES = -1.03

Step 2 - Action: Modify SCM for intervention do(Discount = 10):

Step 3 - Prediction:

visits <- 179  # Held constant (no intervention on visits)
discount <- 10  # Intervention
Sales <- 55.65 + 0.24 * visits + 0.18 * discount + ES
# Sales = 99.38

Conclusion: If the discount had been 10% instead of 19%, they would have sold approximately 99 items (2 fewer than the observed 101).

Example - Pearl’s Three-Step Method: General template for counterfactual reasoning:

# Pearl's three-step method
# Step 1: Abduction
ES <- observed.Y - (beta0 + beta1 * observed.X)

# Step 2: Action
# Modify structural equation for intervention

# Step 3: Prediction
counterfactual.Y <- beta0 + beta1 * intervention.X + ES

Markov Equivalence¶

Markov Equivalence refers to the fact that different DAGs can represent the same conditional independence relationships.

Markov Equivalent Classes:

• Chains and forks are Markov equivalent

• CPDAGs (Completed Partially Directed Acyclic Graphs) represent equivalence classes

Example: The following DAGs are Markov equivalent (they encode the same conditional independencies):

• $X \rightarrow Y \rightarrow Z$

• $X \leftarrow Y \rightarrow Z$

• $X \leftarrow Y \leftarrow Z$

All three imply: $X \perp\!\!\!\perp Z | Y$

Implications:

• Without additional assumptions, we cannot distinguish between Markov equivalent DAGs from observational data alone

• We need experimental data or strong assumptions to identify the true causal structure

Causal Structure Learning¶

Causal structure learning is the process of discovering causal relationships from data.

PC Algorithm (Peter & Clark, 1991)¶

The PC algorithm is a constraint-based method that:

1. Starts with a fully connected undirected graph

2. Tests for conditional independencies in the data

3. Removes edges where conditional independence is found

4. Orients edges based on v-structures (colliders)

library(pcalg)

# Create skeleton (undirected graph)
skeleton <- pc(data, alpha = 0.05)

# Orient edges
pc.fit <- pc(skeleton, data, alpha = 0.05)

# Visualize
plot(pc.fit)

Assumptions:

• Causal sufficiency: All common causes are measured

• Causal Markov condition: D-separation implies conditional independence

• Faithfulness: No independence relationships are “accidental”

LiNGAM Algorithm (Linear Non-Gaussian Acyclic Model)¶

LiNGAM identifies the full causal structure (not just Markov equivalence class) by assuming:

• Linear relationships

• Non-Gaussian errors

• Acyclic structure

How it works: Uses Independent Component Analysis (ICA) to identify the causal order.

library(pcalg)

# Fit LiNGAM model
fit.lingam <- lingam(data)

# B matrix contains edge weights
B <- fit.lingam$Bpruned
colnames(B) <- rownames(B) <- colnames(data)

# Confidence intervals
ci <- fit.lingam$ci

# Visualize
plot(fit.lingam)

Example: Online Sales Data

# Fit LiNGAM to sales data
fit.sales <- lingam(dat_sales)
B <- round(fit.sales$Bpruned, 2)

# Resulting DAG:
# Adv -> Visits -> Sales
# Discount -> Sales
# (No direct effect of Adv on Sales, only through Visits)

Four Steps of Causal Inference¶

Step 1: Create a Causal Model

• Use expert knowledge to specify a DAG

• Represent all relevant variables and their causal relationships

Step 2: Identify the Causal Effect

• Determine if the causal effect can be identified from observational data

• Use d-separation and backdoor criterion to find valid adjustment sets

Step 3: Estimate the Causal Effect

• Apply adjustment formula or other identification strategies

• Use appropriate statistical methods (regression, matching, etc.)

Step 4: Test/Validate the Results

• Check sensitivity to modeling assumptions

• Test robustness to different specifications

• Validate with external data if available