Advanced Statistical Data Analysis

Overview¶

The course focuses on regression modeling as a primary tool for explaining relationships between:

• Response variable (output/dependent variable): The variable we want to explain or predict

• Explanatory variables (inputs/independent variables/predictors): Variables that explain the response

Mathematical Formulation:

where $f\left\langle \cdot \right\rangle$ is the systematic component and $E$ is the random error term.

Symbols¶

Advanced Regression Modelling¶

Multiple Linear Regression¶

Objectives of Regression Analysis¶

1. General description of data structure.

2. Assessment of the effect of explanatory variables on the response.

3. Prediction of future observations.

Error Assumptions¶

The standard assumptions for the error terms $\mathcal{E}_i$ are:

• Stochastically independent.

• Expectation zero and constant variance $\sigma^2$ (homoscedasticity).

• Normally (Gaussian) distributed: $\mathcal{E}_i \sim \mathcal{N}(0, \sigma^2)$.

Matrix Representation¶

To simplify notation, the regression equation Definition 2 is written in matrix form:

where:

• $\mathbf{Y}$ is an $n \times 1$ vector of responses.

• $\mathbf{X}$ is an $n \times p$ matrix of explanatory variables (including a column of 1s for the intercept).

• $\boldsymbol{\beta}$ is a $p \times 1$ vector of unknown coefficients ($p = m+1$).

• $\boldsymbol{\mathcal{E}}$ is an $n \times 1$ vector of unobserved random variables.

Tukey’s First-Aid Transformations¶

Standard recommendations used to linearize relationships and stabilize variance when there is no specific domain theory to guide variable transformation. These should be applied to both explanatory variables and responses unless a valid reason exists to do otherwise:

Model Fitting and Diagnostics¶

Least Squares Estimation¶

The coefficients $\boldsymbol{\beta}$ are estimated by minimizing the sum of squared residuals.

The OLS estimator is given by:

Model Adequacy (Residual Analysis)¶

Model adequacy is checked using diagnostic plots:

• Tukey-Anscombe Plot: Residuals vs. Fitted values to check for non-linearity or heteroscedasticity.

• Normal Q-Q Plot: To check the normality assumption of errors.

• Scale-Location Plot: To check for constant variance.

• Residuals vs. Leverage: To identify influential observations (Cook’s Distance).

Model Formulation:

Assumptions:

1. $E_i \sim \mathcal{N}(0, \sigma^2)$ (normally distributed errors)

2. $\mathbb{E}[E_i] = 0$ (zero expectation)

3. $\text{Var}(E_i) = \sigma^2$ (homoscedasticity)

4. $E_i$ are stochastically independent

5. No perfect multicollinearity

Matrix Notation:

Gauss-Markov Theorem: Under assumptions 1-4 (without normality), OLS is BLUE (Best Linear Unbiased Estimator).

Properties:

• Unbiased: $\mathbb{E}[\hat{\underline{\beta}}] = \underline{\beta}$

• Consistent: $\hat{\underline{\beta}} \xrightarrow{P} \underline{\beta}$

• Asymptotic Normality: $\hat{\underline{\beta}} \sim \mathcal{N}(\underline{\beta}, \sigma^2 (\boldsymbol{X}^T \boldsymbol{X})^{-1})$

Advanced Topics in Linear Regression Modelling¶

Learning Objectives¶

• Weighted Least Squares (for non-constant variance)

• Robust Fitting (to handle outliers and leverage points)

• Fitting Smooth Functions (using LOESS and splines)

• Additive Regression Models

• Model Building strategies

Problem: OLS is not robust¶

Measuring Robustness

How robust an estimator is, can be investigated by two simple measures:

• influence function and gross error sensitivity: The (gross error) sensitivity is based on the influence function and measures the maximum effect of a single observation on the estimated value

• breakdown point: returns the minimum proportion of data that can be altered without causing completely unreliable estimates

Hence an estimator with good robustness properties has:

• a bounded gross error sensitivity and

• a breakdown point of (about) 0.5

Robust Regression: From OLS to MM-estimator

We need S-estimators, M-estimators, and MM-estimators because real-world data often violates the assumptions of OLS regression. Therefore, we need robust estimators:

1. OLS has 0% breakdown point - even one bad data point can ruin your analysis

2. M-estimators alone can’t handle leverage points - outliers in predictors are still problematic

3. S-estimators are too computationally expensive for routine use

4. MM-estimators give us the best of both worlds - robustness AND practicality

Use MM-estimators as your default robust regression method. They provide protection against both outliers and leverage points while maintaining good statistical efficiency.

Logistic Regression (Binary Regression)¶

Theoretical Background¶

Logistic Regression Model:

Logit Link:

Odds Ratio:

Asymptotic Properties:

• $\hat{\underline{\beta}} \sim \mathcal{N}(\underline{\beta}, \mathcal{I}^{-1})$

• Wald tests valid for large samples

Generalized Linear Models (GLM) - A Unifying Model Family¶

Learning Objectives¶

• Identify members of the exponential family

• Know the two basic elements that define the GLM

• Fit GLMs in R and know the algorithm underlying it

• Interpret R output of a GLM fit

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:

Common Distributions¶

• 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)

glm(y ~ x1 + x2, family = gaussian)
# or simply: lm(y ~ x1 + x2)

# Binary data
glm(y ~ x1 + x2, family = binomial)

# Aggregated data (successes and failures)
glm(cbind(success, failure) ~ x1 + x2, family = binomial)

glm(y ~ x1 + x2, family = poisson)

glm(y ~ x1 + x2, family = Gamma)

glm(y ~ x1 + x2, family = inverse.gaussian)

Distribution Decision Tree¶

• Check variance: For Poisson, verify mean ≈ variance. If variance > mean, consider Negative Binomial (not in exponential family)

• Check range: Ensure response values are valid for the distribution (e.g., Poisson requires y ≥ 0)

• Check skewness: Right-skewed data may need Gamma or Inverse Gaussian

• Check zeros: Excess zeros may require zero-inflated models

GLM Structure¶

The Generalized Linear Model (GLM) is defined by two fundamental elements:

1. Distributional Element
   The distribution of the response $Y_i$, given the explanatory variables $\underline{x}_i$, belongs to the two-parameter exponential family with:

• Expectation: $\mathbb{E}\langle Y_i \mid \underline{x}_i \rangle = \mu_i$

• Dispersion: $\phi$ The responses $Y_i$, $i = 1, \ldots, n$, are independently distributed.

2. Structural Element
   The expectation $\mu_i$ is related to the linear predictor $\eta_i = \underline{x}_i^T \underline{\beta}$ of the explanatory variables by applying a (possibly nonlinear) function $g\left\langle \cdot \right\rangle$ on $\mu_i$:

   $$g\left\langle \mu_i \right\rangle = \eta_i = \underline{x}_i^T \underline{\beta}.$$

   (15)

   The function $g\left\langle \cdot \right\rangle$ is called the link function.

The estimator for $\underline{\beta}$ is derived based on the maximum likelihood principle:

• It is defined by a nonlinear system of equations

• It can be solved most efficiently by the IRLS algorithm

This two-element structure (distributional + structural) is what unifies various regression models (linear, logistic, Poisson, etc.) under the GLM framework.

IRLS Algorithm for GLM Fitting¶

The Iteratively Reweighted Least Squares (IRLS) algorithm fits Generalized Linear Models through the following iterative steps:

1. Initialize
   Put $\hat{\mu}_i = Y_i$

2. Ensure Valid Weights
   If necessary, modify $\hat{\mu}_i$ so that the weights $\tilde{w}_i$ below are non-zero (e.g., as done by the empirical logits)

3. Calculate Weights

   $$\tilde{w}_i \coloneqq \frac{w_i}{V\left\langle \hat{\mu}_i \right\rangle \left( g'\left\langle \hat{\mu}_i \right\rangle \right)^2}$$

   (16)

4. Form Adjusted Response

   $$z_i \coloneqq \underline{x}_i^T \hat{\underline{\beta}} + \left( Y_i - \hat{\mu}_i \right) g'\left\langle \hat{\mu}_i \right\rangle$$

   (17)

5. Weighted Least Squares
   Regress $\sqrt{\tilde{w}_i} z_i$ on $\sqrt{\tilde{w}_i} \underline{x}_i$ by least squares to obtain $\hat{\underline{\beta}}$. That is, solve the system of linear equations:

   $$\sum_{i=1}^n \tilde{w}_i \left( z_i - \underline{x}_i^T \hat{\underline{\beta}} \right) \underline{x}_i = \underline{0} \quad \text{for } \hat{\underline{\beta}}.$$

   (18)

6. Update Fitted Values
   Calculate the fitted values:

   $$\hat{\mu}_i = g^{-1}\left\langle \underline{x}_i^T \hat{\underline{\beta}} \right\rangle$$

   (19)

7. Iterate
   Return to step 2 and iterate until convergence

The dispersion parameter $\phi$ must be estimated as well:

• It can also be done by maximizing the likelihood.

• As we know from multiple linear regression, we rather use an unbiased estimator for estimating the dispersion parameter $\phi = \sigma^2$

• In the case of binomial and Poisson distribution: $\phi = 1$

• For other distributions: see table above

Key R Code Snippets¶

Week 6: Inference, Model Simplification, Variable Selection¶

Theoretical Background¶

Deviance:

Deviance Test:

Overdispersion:

Confidence Intervals:

• Wald: $\hat{\beta}_j \pm z_{\alpha/2} \cdot \text{SE}(\hat{\beta}_j)$

• Deviance-based: More accurate for small samples

• Bootstrap: Nonparametric

Learning Objectives¶

• Know what deviances are in GLM

• Know what overdispersion is and can identify it

• Know when to apply Wald-type vs deviance-based confidence intervals

• Apply methods in statistical data analysis using R

Key R Code Snippets¶

model <- glm(y ~ x1 + x2, family = poisson, data = mydata)
null_deviance <- model$null.deviance
resid_deviance <- model$deviance
anova(model_simple, model_complex, test = "LRT")

model_poisson <- glm(y ~ x1 + x2, family = poisson, data = mydata)
dispersion <- model_poisson$deviance / model_poisson$df.residual
if (dispersion > 1.5) {
  model_quasi <- glm(y ~ x1 + x2, family = quasipoisson, data = mydata)
  library(MASS)
  model_nb <- glm.nb(y ~ x1 + x2, data = mydata)
}

AIC(model1, model2, model3)
BIC(model1, model2, model3)
step(model1, direction = "both")

Week 7: Diagnostics / Model Adequacy Checking¶

Theoretical Background¶

Diagnostic Plots:

1. Residuals vs Fitted: Check linearity and homoscedasticity

2. Normal Q-Q: Check normality

3. Scale-Location: Check constant variance

4. Leverage: Identify influential points

Multicollinearity:

• High VIF (> 5 or 10) indicates problematic collinearity

• Effects: Unstable estimates, inflated standard errors

Learning Objectives¶

• Understand how AIC is generalized to GLMs

• Check model adequacy and determine which assumptions are violated

• Find appropriate transformations of predictors in a data-driven manner

• Apply methods in statistical data analysis using R

Key R Code Snippets¶

AIC(model1, model2)

par(mfrow = c(2, 2))
plot(model_poisson, which = 1:4)
library(halfnorm)
halfnorm(residuals(model_poisson, type = "deviance"))

library(car)
cooks.distance(model1)
hat.values(model1)
dfbetas(model1)

Week 8: Some Extensions of GLM¶

Theoretical Background¶

Rate Models:

For Poisson with log link:

Quasi Models:

• Consistent $\hat{\underline{\beta}}$ even if variance function misspecified

• Standard errors adjusted by $\hat{\phi}$

Instrumental Variables:

1. First Stage: $X = \pi_0 + \pi_1 Z + U$

2. Second Stage: $Y = \beta_0 + \beta_1 \hat{X} + E$

Learning Objectives¶

• Know what a rate model is and how to analyse it with GLM

• Know how a quasi model extends a GLM and when to apply it

• Fit methods to data, interpret results, make inference and predictions

Key R Code Snippets¶

model_rate <- glm(y ~ x1 + x2, family = poisson, offset = log(t), data = mydata)

model_quasipoisson <- glm(y ~ x1 + x2, family = quasipoisson, data = mydata)
model_quasibinomial <- glm(cbind(Y, m-Y) ~ x1 + x2, family = quasibinomial, data = mydata)

library(AER)
model_iv <- ivreg(Y ~ X + W | Z + W, data = mydata)
summary(model_iv)

Causality¶

Introduction¶

Theoretical Background¶

The Ladder of Causation (Pearl):

1. Association (Seeing): Observing patterns

2. Intervention (Doing): Estimating effects of actions

3. Counterfactuals (Imagining): Reasoning about hypotheticals

Notation:

• $P(Y | X = x)$: Conditional probability

• $P(Y | \text{do}(X = x))$: Causal effect

Causal Graphical Models:

• Nodes: Variables

• Edges: Causal relationships

• DAG: Directed Acyclic Graph (no cycles)

Graphical Structures:

1. Chain: $X → Z → Y$

• $X$ and $Y$ are conditionally independent given Z; $(X \perp\!\!\!\perp Y \mid Z)$

2. Fork: $X ← Z → Y$

• $X$ and $Y$ are conditionally independent given Z; $(X \perp\!\!\!\perp Y \mid Z)$

3. Collider: $X → Z ← Y$

• $X$ and $Y$ are independent $(X \perp\!\!\!\perp Y)$, but are conditionally dependent on Z and any descendants of Z; $(X \not\perp\!\!\!\perp Y \mid Z)$

Learning Objectives¶

• Understand when and why causal reasoning is important

• Be familiar with causal graphical models

• Understand the difference between association, interventions, and counterfactuals

• Know the difference between experimental and observational data

Intervention changes distribution¶

Key R Code Snippets¶

Simpson’s Paradox & Graphical Models¶

Theoretical Background¶

Simpson’s Paradox: A trend appears in different groups but disappears or reverses when combined.

Factorization:

• Chain: $P(X,Z,Y) = P(X)P(Z|X)P(Y|Z)$

• Fork: $P(X,Z,Y) = P(Z)P(X|Z)P(Y|Z)$

• Collider: $P(X,Z,Y) = P(X)P(Y)P(Z|X,Y)$

Bayesian Network Factorization:

Learning Objectives¶

• Understand Simpson’s paradox and how to analyze it

• Express joint distribution by factorization

• Understand (conditional) independence

• Identify (conditional) independences for chain, fork, and collider

Key R Code Snippets¶

D-Separation & Causal Effect Estimation¶

Theoretical Background¶

D-Separation:

X and Y are d-separated given Z if all paths between X and Y are blocked by Z if and only if:

1. a middle node of a chain (X → B → Y) or a fork (X ← B → Y) is conditioned on (i.e. it is in Z)

2. a middle node of a collider (X → B ← Y) is not conditioned on (i.e. is not in Z), and no descendant of B is in Z.

If Z blocks every path between two nodes X and Y, then X and Y are d-separated conditioning on Z and thus are conditionally independent on Z.

Adjustment Formula:

Backdoor Criterion: Z satisfies backdoor criterion if:

1. Z blocks all backdoor paths from X to Y

2. Z does not contain any descendants of X

Learning Objectives¶

• Know D-Separation

• Estimate causal effect using the adjustment formula

• Define adjustment sets with the backdoor criterion

Structural Causal Models¶

Theoretical Background¶

Linear SCM:

Direct vs. Total Effects:

• Total: Effect through all paths

• Direct: Effect through direct path only

Learning Objectives¶

• Be familiar with (linear) structural causal models

• Understand direct and total causal effects

• Estimate direct and total causal effect from linear SCM

Key R Code Snippets¶

Advanced Topics¶

Theoretical Background¶

Instrumental Variables:

1. Relevance: $\text{Cov}(Z, X) \neq 0$

2. Exclusion: Z has no direct effect on Y except through X

3. Exogeneity: $\text{Cov}(Z, E) = 0$

Two-Stage Least Squares:

1. First stage: $X = \pi_0 + \pi_1 Z + U$

2. Second stage: $Y = \beta_0 + \beta_1 \hat{X} + E$

Counterfactuals: Pearl’s three-step method:

1. Abduction: Use the observed information to determine the values of the noise variables.

2. Action: Modify the causal model M by removing the structural equation for the variable X and replacing with X = x.

3. Prediction: Use the modified model $M_x$ to compute the counterfactual.

Markov Equivalence:

• Chains and Forks are Markov equivalent, i.e. their graphs imply the same conditional independencies.

• Colliders encode a different independence.

→ Colliders are identifiable from data, whereas fork and chains are not
→ All graphs where we can flip the arrows without introducing new colliders and keeping all existing colliders are Markow Equivalent.

Markov Equivalence can be represented by a CPDAG.

Completed Partially DAG (CPDAG)

represents the Markov equivalence class of DAGs

• Directed edge: causal direction is clear (all graphs agree)

• Undirected edge: causal direction is ambiguous (graphs disagree)

• Encodes what can be identified from data alone

Learning Objectives¶

• Understand instrumental variables

• Be familiar with counterfactual reasoning

• Know the three steps in computing counterfactuals

• Understand Markov Equivalence

Key R Code Snippets¶

library(AER)
model_iv <- ivreg(Y ~ X + W | Z + W, data = mydata)
summary(model_iv)
first_stage <- lm(X ~ Z + W, data = mydata)
summary(first_stage)

library(pcalg)
dag1 <- dag(c("X", "Z", "Y"), list(c("X", "Z"), c("Z", "Y")))
dag2 <- dag(c("X", "Z", "Y"), list(c("Z", "X"), c("Z", "Y")))
meq(dag1, dag2)

Causal Structure Learning¶

Theoretical Background¶

Assumptions of Independence-Based Methods:

• Causal sufficiency, i.e., no unobserved confounding variables.

• Acyclicity: No directed cycles in the graph

• Markov: Graph → Data, i.e., if X and Y are d-separated by Z in the graph, then X and Y are conditionally independent given Z.

• Faithfulness: Data → Graph, i.e., if X and Y are conditionally independent given Z, then X and Y are d-separated by Z in the graph.

$\to$ What about non-Gaussian structural equations? LinGAM to the rescue!

4 Steps of Causal Inference¶

1. DAG: Create a causal model using expert knowledge.

2. Identification: Determine whether and how the causal effect can be identified from the observational data.

3. Estimation: Estimate the causal effect from the data.

4. Validity: Test the estimated causal effect:

Key R Code Snippets¶