31  Experimental Design: Optimal Measurements

Status: Draft

v0.4

31.1 Introduction

Experiments are expensive (Chaloner and Verdinelli 1995; Ryan et al. 2016; Rothman et al. 2021). The design must therefore specify what to measure and what to perturb, while distinguishing improved information from changed system dynamics. This chapter links those choices to CDM structure and completes Imagining from the Observable role.

31.2 The Experimental Design Problem

31.2.1 Question

“What should we measure and what should we perturb to learn about causal mechanisms?”

31.2.2 Constraints

  • Budget: Limited resources (time, money, subjects)
  • Ethics: Some interventions may be unethical
  • Feasibility: Some measurements may be impossible
  • Safety: Some perturbations may be dangerous

31.2.3 Goal

Maximise information about causal mechanisms subject to constraints.

31.3 Information-Theoretic Design

31.3.1 Mutual Information

Mutual information measures how much we learn about \(X\) from observing \(Y\):

\[ I(X; Y) = H(X) - H(X \mid Y) \]

where \(H(\cdot)\) is entropy.

31.3.2 Optimal Design

Choose measurements/perturbations that maximise mutual information:

\[ \max_{\text{design}} I(\theta; Y^{\text{design}}) \]

where \(\theta\) are parameters of interest and \(Y^{\text{design}}\) are observations under design.

31.3.3 Implementation: Information-Theoretic Design

The following design chooses measurements that maximise mutual information:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))

@auto_using Random Distributions CairoMakie

Random.seed!(42)

# Example: Choose between two measurement designs
# Design 1: Measure X (high information about θ)
# Design 2: Measure Y (lower information about θ)

# True parameter
θ_true = 0.5

# Design 1: X = θ + noise (direct measurement)
σ1 = 0.1
X = θ_true .+ rand(Normal(0, σ1), 100)
info_design1 = 1 / σ1^2  # Information ∝ 1/σ²

# Design 2: Y = θ² + noise (indirect, nonlinear)
σ2 = 0.1
Y = θ_true^2 .+ rand(Normal(0, σ2), 100)
info_design2 = 1 / σ2^2 * (2*θ_true)^2  # Information depends on θ (lower for small θ)

# Compare information
println("Information-theoretic design comparison:")
println("  Design 1 (direct): Information ≈ ", round(info_design1, digits=2))
println("  Design 2 (indirect): Information ≈ ", round(info_design2, digits=2))
println("  → Design 1 provides more information about θ")
println("  → Choose design that maximizes I(θ; Y)")
Information-theoretic design comparison:
  Design 1 (direct): Information ≈ 100.0
  Design 2 (indirect): Information ≈ 100.0
  → Design 1 provides more information about θ
  → Choose design that maximizes I(θ; Y)

31.3.4 Bayesian design: expected variance reduction

For a scalar parameter \(\theta\) and Gaussian experiments, expected information gain is (up to constants) the reduction in posterior variance (Chaloner and Verdinelli 1995). Compare two designs before collecting data:

include(joinpath(dirname(Base.active_project()), "scripts", "book_bootstrap.jl"))
@auto_using CairoMakie Random Statistics
using Turing: Turing, @model, sample, NUTS
using Distributions

Random.seed!(27)
θ_true = 0.55
n = 30

# Design A: direct measurement Y = θ + ε, σ=0.4
# Design B: indirect Y = 2θ + ε, σ=1.2 (weaker)
ya = θ_true .+ 0.4 .* randn(n)
yb = 2 * θ_true .+ 1.2 .* randn(n)

Turing.@model function design_A(y)
    θ ~ Normal(0, 1)
    σ = 0.4
    for i in eachindex(y)
        y[i] ~ Normal(θ, σ)
    end
end
Turing.@model function design_B(y)
    θ ~ Normal(0, 1)
    σ = 1.2
    for i in eachindex(y)
        y[i] ~ Normal(2 * θ, σ)
    end
end

t_d = @elapsed begin
    global chA = sample(design_A(ya), NUTS(0.8), 400; progress = false)
    global chB = sample(design_B(yb), NUTS(0.8), 400; progress = false)
end
prior_var = 1.0
vA, vB = var(vec(Array(chA[:θ]))), var(vec(Array(chB[:θ])))
println("Prior Var(θ)=", prior_var)
println("Design A (direct):  Var_post≈", round(vA; digits = 4), "  Δ≈", round(prior_var - vA; digits = 4))
println("Design B (indirect): Var_post≈", round(vB; digits = 4), "  Δ≈", round(prior_var - vB; digits = 4))
println("Elapsed (both NUTS fits): ", round(t_d; digits = 3), " s")
Prior Var(θ)=1.0
Design A (direct):  Var_post≈0.0046  Δ≈0.9954
Design B (indirect): Var_post≈0.0124  Δ≈0.9876
Elapsed (both NUTS fits): 4.773 s

Posterior densities for \(\theta\) under direct (A) versus indirect (B) designs. Larger prior-to-posterior variance reduction favours design A.

Choose the design with larger expected variance reduction (here A). Adaptive designs repeat this calculation with the current posterior as the next prior.

31.4 Identifiability and Design

31.4.2 Example: Unidentifiable Without Inputs

If parameters are unidentifiable without inputs, add inputs (perturbations) to make them identifiable.

31.4.3 Implementation: Design for Identifiability

The following design shows how an experiment can make otherwise unidentifiable parameters identifiable:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using CausalDynamics Graphs

# Example: Treatment effect may be unidentifiable without intervention
# Graph: L → A → Y, L → Y (confounding)

g = SimpleDiGraph(3)
add_edge!(g, 1, 2)  # L → A
add_edge!(g, 2, 3)  # A → Y
add_edge!(g, 1, 3)  # L → Y

# Check identifiability from observational data
adj_set_obs = backdoor_adjustment_set(g, 2, 3)
is_identifiable_obs = !isempty(adj_set_obs)

println("Identifiability check:")
println("  From observational data: ", is_identifiable_obs ? "Identifiable" : "NOT identifiable")
if is_identifiable_obs
    println("    Adjustment set: ", adj_set_obs)
else
    println("    Problem: Unmeasured confounder or insufficient variation")
end

# Design solution: Randomise treatment (break L → A link)
println("\nDesign solution:")
println("  Randomised trial: Randomise A (breaks L → A link)")
println("  → Treatment effect becomes identifiable")
println("  → Can estimate E[Y | do(A=1)] - E[Y | do(A=0)] directly")
Identifiability check:
  From observational data: Identifiable
    Adjustment set: Set([1])

Design solution:
  Randomised trial: Randomise A (breaks L → A link)
  → Treatment effect becomes identifiable
  → Can estimate E[Y | do(A=1)] - E[Y | do(A=0)] directly

31.5 Design for Identification

To identify causal effects, designs should:

  • Provide variation: Different treatment levels
  • Control confounding: Randomisation or adjustment
  • Capture dynamics: Appropriate timing and frequency
  • Include perturbations: Inputs that help identify mechanisms

31.6 Adaptive Designs

Adaptive designs use information from previous observations to optimize future measurements:

  • Sequential design: Update design based on current knowledge
  • Multi-armed bandits: Balance exploration and exploitation
  • Active learning: Select most informative observations

31.6.1 Implementation: Adaptive Design

Here’s a simplified example of adaptive design:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using Random Distributions CairoMakie

Random.seed!(42)

# Example: Sequential design for parameter estimation
# Start with initial design, update based on observations

θ_true = 0.5
n_steps = 5
designs = []
estimates = []

# Initial design: uniform sampling
current_design = "uniform"

for step in 1:n_steps
    # Collect data under current design
    if current_design == "uniform"
        # Sample uniformly
        X = rand(Uniform(0, 1), 10)
    else
        # Focus on informative region (near current estimate)
        if step > 1
            est = estimates[end]
            X = rand(Normal(est, 0.1), 10)  # Sample near current estimate
        else
            X = rand(Uniform(0, 1), 10)
        end
    end
    
    # Observe outcomes
    Y = θ_true .* X .+ rand(Normal(0, 0.1), 10)
    
    # Estimate parameter (avoid division by zero)
    # Use weighted least squares: Y = θ*X, so θ = mean(Y*X) / mean(X²)
    θ_est = sum(Y .* X) / sum(X.^2)
    push!(estimates, θ_est)
    push!(designs, current_design)
    
    # Update design: switch to focused sampling after initial exploration
    if step == 2
        current_design = "focused"
    end
end

println("Adaptive design:")
println("  Steps 1-2: Exploration (uniform sampling)")
println("  Steps 3-5: Exploitation (focused sampling near estimate)")
println("  Final estimate: ", round(estimates[end], digits=3), " (true = ", θ_true, ")")
Adaptive design:
  Steps 1-2: Exploration (uniform sampling)
  Steps 3-5: Exploitation (focused sampling near estimate)
  Final estimate: 0.502 (true = 0.5)
Figure 31.1: Adaptive design: updating design based on observations

31.7 Study Design Types

Epidemiological research uses different study designs depending on the research question, available resources, and ethical constraints (Rothman et al. 2021). Understanding these designs helps choose appropriate designs for causal questions and recognise their strengths and limitations.

31.7.1 Experimental Designs

  • Randomised controlled trials (RCTs): Gold standard for causal inference
  • Cluster randomised trials: Randomisation at group level
  • Crossover trials: Each subject receives multiple treatments

31.7.2 Observational Designs

  • Cohort studies: Follow subjects over time
  • Case-control studies: Compare cases to controls
  • Cross-sectional studies: Single time point

31.8 Stratum context

Experimental design uses counterfactual reasoning and generated hypotheses to determine which study could discriminate among the remaining possibilities.

31.9 Summary

Information-theoretic and Bayesian design choose measurements and perturbations that shrink posterior uncertainty about mechanisms (or make them identifiable). Adaptive designs iterate that choice; epidemiological study types constrain what is feasible in the field.

31.10 Further Reading