Second Differences (Interaction Effects)

Quantifying how marginal effects vary across moderator levels

Conceptual Foundation

What Are Second Differences?

Second differences measure whether the marginal effect of one variable differs across levels of another variable. In generalized linear models, where effects are inherently nonlinear, second differences provide the natural way to quantify interaction effects on the predicted outcome scale.

Definition: For two variables X (focal variable) and Z (moderator), the second difference is defined as:

\[\Delta^2 = [E(Y|X=1,Z=1) - E(Y|X=0,Z=1)] - [E(Y|X=1,Z=0) - E(Y|X=0,Z=0)]\]

This represents the difference in marginal effects, a "difference of differences", capturing how the effect of X changes when Z changes.

Relationship to Interaction Terms

In linear models, the interaction coefficient directly represents the second difference. In nonlinear models (logit, probit, Poisson), the interaction coefficient does not directly represent the interaction on the predicted scale, necessitating explicit computation of second differences from marginal effects.

Methodological Context

Second differences extend the conceptual framework of Margins.jl by addressing a fundamental question: "Does the effect of X depend on Z?" This complements the two core questions of marginal effects analysis:

1. Effects: "How much does Y change when I change X?"
2. Predictions: "What value of Y should I expect?"

Second differences answer: "Does the answer to question 1 change depending on Z?"

Computational Framework

Two Approaches to Second Differences

Margins.jl implements second differences through two complementary approaches:

1. Discrete Contrast Approach (Population-Based)

The discrete contrast approach uses pre-computed population-level marginal effects (AMEs) via functions like second_differences(), second_differences_pairwise(), and second_differences_all_contrasts(). This approach provides:

• Population representativeness: Effects averaged over the sample distribution
• Statistical validity: Proper delta-method standard errors accounting for full covariance structure
• Flexibility: Works with binary, categorical, and continuous moderators
• Use case: Comparing AMEs across discrete moderator levels

2. Local Derivative Approach (Profile-Based)

The local derivative approach uses second_differences_at() to compute derivatives at specific evaluation points via finite differences. This approach provides:

• Profile-specific analysis: Evaluate how effects vary at representative covariate combinations
• Continuous moderation: Direct derivative ∂AME/∂modifier at specific points
• Scenario control: Hold other variables constant while varying only the modifier
• Use case: Understanding effect heterogeneity at particular policy-relevant profiles
• Requirement: Modifier must be continuous/numeric (focal variable can be continuous or categorical)

Statistical Inference

Standard errors for second differences employ rigorous delta-method computation. For a second difference comparing two AMEs (AME₁ and AME₂) with parameter gradients g₁ and g₂:

\[\text{SE}(\Delta^2) = \sqrt{(g_2 - g_1)' \Sigma (g_2 - g_1)}\]

where Σ is the model's parameter covariance matrix. This approach properly accounts for:

• Covariance between the two AME estimates
• Full uncertainty propagation from model parameters
• Appropriate test statistics and confidence intervals

The delta-method computation ensures that hypothesis tests maintain appropriate Type I error rates and confidence intervals achieve nominal coverage probabilities, meeting the statistical validity standards required for econometric research.

Function Reference

Main Interface

second_differences()

Unified interface for computing second differences across all modifier types. This is the recommended entry point that automatically handles binary, categorical, and continuous moderators.

Signature:

second_differences(
    ame_result::EffectsResult,
    variables::Union{Symbol, Vector{Symbol}},
    modifier::Symbol,
    vcov::Matrix{Float64};
    contrast::String = "derivative",
    modifier_type::Symbol = :auto,
    all_contrasts::Bool = true
)

Arguments:

• ame_result: Result from population_margins() with scenarios over the modifier
• variables: Focal variable(s) to analyze (single Symbol or Vector{Symbol})
• modifier: The moderating variable
• vcov: Parameter covariance matrix from the model (via vcov(model))

Keyword Arguments:

• contrast: Focal variable contrast for categorical variables (default: "derivative")
• modifier_type: Modifier classification (:auto, :binary, :categorical, :continuous)
• all_contrasts: Compute all focal variable contrasts when applicable (default: true)

Returns: DataFrame with columns:

• variable: Focal variable name
• modifier: Modifier variable name
• contrast: Contrast description (for categorical focal variables)
• modifier_level1, modifier_level2: Levels being compared
• second_diff: Second difference estimate (or slope for continuous)
• se: Standard error (delta method)
• z_stat: Z-statistic
• p_value: P-value for H₀: second difference = 0
• ame_at_level1, ame_at_level2: AME values at each level
• modifier_type: Detected modifier type
• significant: Boolean indicator (p < 0.05)

Specialized Functions

second_differences_pairwise()

Computes all pairwise modifier comparisons for a single focal variable contrast.

Use case: When analyzing a categorical or continuous moderator with multiple levels, and you want all pairwise comparisons for one specific contrast.

Signature:

second_differences_pairwise(
    ame_result::EffectsResult,
    variable::Symbol,
    modifier::Symbol,
    vcov::Matrix{Float64};
    contrast::String = "derivative",
    modifier_type::Symbol = :auto
)

Returns: DataFrame with one row per pairwise modifier comparison.

second_differences_all_contrasts()

Computes second differences for all contrasts of a categorical focal variable.

Use case: When your focal variable is categorical with multiple contrasts (e.g., pairwise comparisons across education levels) and you want to see how each contrast varies across modifier levels.

Signature:

second_differences_all_contrasts(
    ame_result::EffectsResult,
    variable::Symbol,
    modifier::Symbol,
    vcov::Matrix{Float64};
    modifier_type::Symbol = :auto
)

Returns: DataFrame with one row per (focal contrast × modifier pair) combination.

second_differences_at()

Computes local derivatives of marginal effects with respect to a continuous modifier at specified evaluation points using finite differences.

Use case: Profile-based analysis where you want to understand how effects vary at specific covariate combinations, holding other variables constant.

Important: The modifier must be numeric/continuous because the function uses finite differences (perturbing the modifier by ±δ). For categorical moderators, use second_differences() or second_differences_pairwise() instead.

Signature:

second_differences_at(
    model::RegressionModel,
    data,
    variables::Union{Symbol, Vector{Symbol}},
    modifier::Symbol,
    vcov;
    at::Union{Symbol, Real, Vector, NamedTuple} = :mean,
    profile::NamedTuple = NamedTuple(),
    delta::Union{Symbol, Real} = :auto,
    scale::Symbol = :response
)

Arguments:

• model: Fitted regression model
• data: Data frame used to fit the model
• variables: Focal variable(s) to compute derivatives for (can be continuous or categorical)
• modifier: Continuous modifier variable (differentiate with respect to this - must be numeric)
• vcov: Parameter covariance matrix or function

Keyword Arguments:

• at: Where to evaluate the derivative
  • :mean: At mean(modifier) [default]
  • :median: At median(modifier)
  • Numeric value: At specified modifier value
  • Vector: Multiple evaluation points
  • NamedTuple: Full explicit profile including modifier
• profile: Additional variables to hold fixed (scalar or vector values)
• delta: Step size for finite difference (:auto or numeric value)
• scale: Prediction scale (:link or :response)

Returns: DataFrame with columns:

• variable: Focal variable name
• contrast: Contrast description
• modifier: Modifier variable name
• eval_point: Evaluation point for modifier
• derivative: ∂AME/∂modifier (per unit change in modifier)
• se: Standard error of derivative
• z_stat, p_value: Statistical inference
• delta_used: Actual step size used
• Additional columns for profile variables

Examples:

# Continuous focal variable, continuous modifier
sd = second_differences_at(model, data, :income, :age, vcov(model))
# → 1 row: derivative of income effect w.r.t. age

# Categorical focal variable, continuous modifier
# (education is categorical: "hs", "college", "grad")
sd = second_differences_at(model, data, :education, :income, vcov(model))
# → Multiple rows: one per education contrast
# → e.g., "college - hs", "grad - hs", "grad - college"
# → Shows how each education gap varies with income

# At specific profile with other variables held constant
sd = second_differences_at(model, data, :education, :income, vcov(model);
                          at=50000,           # Income = 50k
                          profile=(age=40,    # Hold age constant
                                  region="north"))  # Hold region constant

# Multiple variables across evaluation points
sd = second_differences_at(model, data, [:x1, :x2], :age, vcov(model);
                          at=[30, 45, 60])
# → 2 variables × 3 eval points = 6 rows

# ERROR: Categorical modifier not allowed
sd = second_differences_at(model, data, :income, :region, vcov(model))
# → Error: "Modifier :region must be numeric for second_differences_at()"
# → Use second_differences() instead for categorical moderators

Variable Type Requirements:

• Focal variable (first argument): Can be continuous or categorical
  • Continuous → 1 row per evaluation point
  • Categorical → Multiple rows per evaluation point (one per contrast)
• Modifier (second argument): Must be continuous/numeric
  • Function uses finite differences: modifier ± δ
  • For categorical moderators, use second_differences() instead

Statistical Notes:

• Uses two-point symmetric finite difference: (AME(at + δ) - AME(at - δ)) / (2δ)
• Delta-method SE computed from gradient information
• Choice of δ (via delta parameter) trades off approximation bias vs variance

Usage Patterns

Basic Workflow

The standard workflow for computing second differences involves three steps:

1. Fit your regression model
2. Compute AMEs across modifier levels using population_margins() with scenarios
3. Calculate second differences using second_differences()

using Margins, GLM, DataFrames

# Step 1: Fit model with interaction
model = lm(@formula(y ~ x * treated + z), data)

# Step 2: Compute AMEs at both treatment levels
ames = population_margins(model, data;
                         scenarios=(treated=[0, 1],),
                         type=:effects)

# Step 3: Calculate second difference
sd = second_differences(ames, :x, :treated, vcov(model))
DataFrame(sd)

Output structure:

1×12 DataFrame
 Row │ variable  modifier  contrast    modifier_level1  modifier_level2  second_diff  se        z_stat   p_value   ame_at_level1  ame_at_level2  significant
     │ Symbol    Symbol    String      Int64           Int64            Float64      Float64   Float64  Float64   Float64        Float64        Bool
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ x         treated   derivative               0                1       0.152    0.0428    3.55    0.0004          0.298          0.450      true

Interpretation: The marginal effect of x is 0.152 units larger when treated=1 compared to treated=0 (p < 0.001), indicating a significant positive interaction.

Binary Moderators

Binary moderators produce a single second difference comparing the two levels.

Example: Treatment Effect Heterogeneity

# Does the effect of age vary by treatment status?
model = glm(@formula(outcome ~ age * treated + education),
            data, Binomial(), LogitLink())

ames = population_margins(model, data;
                         scenarios=(treated=[0, 1],),
                         type=:effects)

sd = second_differences(ames, :age, :treated, vcov(model))

Interpretation contexts:

• Positive second difference: Effect is stronger when treated
• Negative second difference: Effect is weaker when treated
• Near-zero second difference: Effect is homogeneous across treatment status

Categorical Moderators

Categorical moderators with K levels produce K(K-1)/2 pairwise comparisons.

Example: Education-Specific Effects

# Does the income effect vary across education levels?
model = lm(@formula(consumption ~ income + education), data)

ames = population_margins(model, data;
                         scenarios=(education=["hs", "college", "grad"],),
                         type=:effects)

sd = second_differences(ames, :income, :education, vcov(model))

Output structure (3 education levels → 3 pairwise comparisons):

3×12 DataFrame
 Row │ variable  modifier   modifier_level1  modifier_level2  second_diff  se      z_stat  p_value  significant
     │ Symbol    Symbol     String          String           Float64      Float64  Float64  Float64  Bool
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ income    education  hs              college               0.152    0.042    3.62    0.0003   true
   2 │ income    education  hs              grad                  0.243    0.051    4.76   <1e-5    true
   3 │ income    education  college         grad                  0.091    0.048    1.90    0.058    false

Interpretation: The marginal effect of income on consumption is 0.152 units higher for college graduates compared to high school graduates (p < 0.001). The effect continues to increase for graduate degree holders, though the college-to-grad increase (0.091) is not statistically significant.

Continuous Moderators

Continuous moderators provide slopes representing the rate of change in marginal effects per unit change in the moderator.

Example: Age-Varying Treatment Effects

# How does the treatment effect change with age?
model = lm(@formula(outcome ~ treatment + age), data)

# Evaluate at three age values
ames = population_margins(model, data;
                         scenarios=(age=[30, 45, 60],),
                         type=:effects)

sd = second_differences(ames, :treatment, :age, vcov(model);
                       modifier_type=:continuous)

Continuous modifier interpretation:

• second_diff: Change in marginal effect per unit increase in moderator
• For age comparison 30→45: slope = (AME₄₅ - AME₃₀) / 15
• Slopes indicate whether effects strengthen (positive) or weaken (negative) with age

Multiple Focal Variables

Analyze multiple variables simultaneously to compare which effects show strongest moderation.

Example: Comparative Moderation Analysis

# Which demographic effects vary most by treatment?
model = lm(@formula(outcome ~ age + education_yrs + experience + treated), data)

ames = population_margins(model, data;
                         scenarios=(treated=[0, 1],),
                         type=:effects,
                         vars=[:age, :education_yrs, :experience])

sd = second_differences(ames, [:age, :education_yrs, :experience],
                       :treated, vcov(model))

Output structure (3 variables × 1 binary contrast = 3 rows):

3×12 DataFrame
 Row │ variable       modifier  second_diff  se      z_stat  p_value  significant
     │ Symbol         Symbol    Float64      Float64  Float64  Float64  Bool
─────┼────────────────────────────────────────────────────────────────────────
   1 │ age            treated        0.045    0.018    2.50    0.012    true
   2 │ education_yrs  treated        0.123    0.031    3.97   <1e-4    true
   3 │ experience     treated        0.012    0.024    0.50    0.617    false

Interpretation: Education effects show the strongest treatment heterogeneity (second diff = 0.123, p < 0.001), followed by age effects (0.045, p = 0.012). Experience effects do not significantly vary by treatment status.

Categorical Focal Variables

When the focal variable is categorical, second differences can be computed for each contrast.

Example: Education Contrasts Across Regions

using CategoricalArrays

# Prepare categorical variable with explicit ordering
data.education = categorical(data.education,
                            levels=["hs", "college", "grad"],
                            ordered=true)

model = lm(@formula(income ~ education + region), data)

# Compute pairwise education contrasts across regions
ames = population_margins(model, data;
                         scenarios=(region=["north", "south", "west"],),
                         type=:effects,
                         vars=[:education])

# Get all education contrasts × all region pairs
sd = second_differences(ames, :education, :region, vcov(model))

Output dimensions: If education has 3 levels (→ 3 pairwise contrasts) and region has 3 levels (→ 3 pairwise comparisons), the result contains 3 × 3 = 9 rows representing the full matrix of focal contrasts by modifier comparisons.

Practical use: Identify which education gaps (e.g., college - hs) vary most across geographic regions.

Local Derivatives at Profiles (second_differences_at())

The profile-based approach computes how marginal effects change continuously with a modifier at specific evaluation points.

Example: Income Effect Variation Across Age

# Does the income effect on consumption strengthen or weaken with age?
model = lm(@formula(consumption ~ income * age + education), data)

# Evaluate at mean age
sd = second_differences_at(model, data, :income, :age, vcov(model))
# Interpretation: derivative shows whether income effect increases/decreases per year of age

# Evaluate at multiple age points
sd = second_differences_at(model, data, :income, :age, vcov(model);
                          at=[25, 35, 45, 55, 65])
# → 5 rows showing how income effect changes across life course

# With fixed profile
sd = second_differences_at(model, data, :income, :age, vcov(model);
                          at=40,
                          profile=(education="college", region="urban"))
# → "For urban college graduates, how does income effect vary with age?"

Key differences from discrete approach:

• Continuous moderation: Direct derivative rather than discrete contrasts
• Profile control: Can hold other variables constant at specific values
• Interpretation: Per-unit change in moderator rather than level-to-level comparison

When to use each approach:

Variable Type Compatibility:

Common error: Trying to use second_differences_at() with a categorical modifier will produce:

Error: Modifier :region must be numeric for second_differences_at()

Solution: Use second_differences() or second_differences_pairwise() for categorical moderators.

Advanced Analysis Patterns

Robust Standard Errors

Second differences support robust standard errors through CovarianceMatrices.jl integration.

using CovarianceMatrices

# Heteroskedasticity-robust second differences
vcov_robust = CovarianceMatrices.HC1(model)
sd = second_differences(ames, :x, :treated, vcov_robust)

# Clustered standard errors
vcov_clustered = CovarianceMatrices.Clustered(model, :firm_id)
sd = second_differences(ames, :x, :treated, vcov_clustered)

The delta-method computation automatically incorporates the robust covariance structure, ensuring appropriate uncertainty quantification under heteroskedasticity or clustering.

Multiple Testing Considerations

When computing many pairwise comparisons (categorical moderators with many levels), consider adjusting for multiple testing:

using MultipleTesting

# Compute all pairwise second differences
sd = second_differences(ames, :x, :education, vcov(model))

# Apply Bonferroni correction
α = 0.05
n_tests = nrow(sd)
sd.significant_bonferroni = sd.p_value .< (α / n_tests)

# Or use false discovery rate control
sd.significant_fdr = adjust(sd.p_value, BenjaminiHochberg()) .< α

Filtering to Specific Contrasts

For categorical focal variables with many contrasts, focus on specific contrasts of interest:

# Compute only the "college - hs" contrast across regions
sd = second_differences(ames, :education, :region, vcov(model);
                       contrast="college - hs",
                       all_contrasts=false)

Visualization Workflows

Second differences integrate naturally with visualization workflows:

using StatsPlots

# Plot second differences with confidence intervals
sd_df = DataFrame(sd)
@df sd_df scatter(:modifier_level2, :second_diff,
                  yerr=1.96 .* :se,
                  xlabel="Modifier Level",
                  ylabel="Second Difference",
                  title="Effect Heterogeneity Across Moderator",
                  legend=false)
hline!([0], linestyle=:dash, color=:gray)  # Reference line at zero

Methodological Notes

Comparison to Coefficient-Based Interactions

In linear models, second differences equal interaction coefficients:

# Linear model: coefficient = second difference
model_linear = lm(@formula(y ~ x * z), data)
interaction_coef = coef(model_linear)[end]  # Coefficient on x:z

ames = population_margins(model_linear, data; scenarios=(z=[0,1],), type=:effects)
sd = second_differences(ames, :x, :z, vcov(model_linear))
# sd.second_diff ≈ interaction_coef (within numerical precision)

In nonlinear models, they diverge:

# Logistic model: coefficient ≠ second difference
model_logit = glm(@formula(y ~ x * z), data, Binomial(), LogitLink())
interaction_coef = coef(model_logit)[end]

ames = population_margins(model_logit, data; scenarios=(z=[0,1],), type=:effects)
sd = second_differences(ames, :x, :z, vcov(model_logit))
# sd.second_diff ≠ interaction_coef (second diff is on probability scale)

This divergence motivates the second differences framework: to express interactions on the interpretable predicted outcome scale rather than the abstract coefficient scale.

Relationship to Ai & Norton (2003)

Ai & Norton (2003) demonstrated that in nonlinear models (logit, probit), the interaction effect:

1. Is not equal to the interaction coefficient
2. Varies across observations
3. Can have different signs than the coefficient

Second differences in Margins.jl operationalize the Ai & Norton framework by:

• Computing population-averaged interaction effects (second differences from AMEs)
• Providing proper standard errors via delta method
• Enabling hypothesis tests for interaction significance

Reference: Ai, C., & Norton, E. C. (2003). Interaction terms in logit and probit models. Economics Letters, 80(1), 123-129.

Population vs Profile Interpretation

Current second differences use population-averaged marginal effects:

• AMEs computed by averaging over the sample distribution at each modifier level
• Second differences reflect population-level interaction effects
• Appropriate for policy analysis requiring external validity

For profile-based local derivatives, use second_differences_at():

• Evaluate interaction effects at specific covariate combinations
• Useful for scenario-specific analysis or representative case interpretation
• See the Local Derivatives at Profiles section above for full documentation

Significance Testing

Hypothesis tests evaluate H₀: second difference = 0 (no interaction on predicted scale).

Important distinction:

• Significant second difference → interaction exists on predicted outcome scale
• Significant interaction coefficient → interaction exists on coefficient scale

In nonlinear models, these are distinct hypotheses. Second differences test the hypothesis most relevant for applied interpretation.

Integration Examples

With GLM.jl Ecosystem

using GLM, CategoricalArrays

# Logistic regression with interaction
model = glm(@formula(employed ~ education * experience + age),
            data, Binomial(), LogitLink())

# Second differences on probability scale
ames = population_margins(model, data;
                         scenarios=(experience=[0, 10, 20],),
                         type=:effects,
                         scale=:response)  # Probability scale

sd = second_differences(ames, :education, :experience, vcov(model);
                       modifier_type=:continuous)
# Interpretation: change in employment effect of education per year of experience

With MixedModels.jl

Second differences support mixed-effects models for clustered/hierarchical data:

using MixedModels

# Linear mixed model with interaction
model = fit(MixedModel,
           @formula(outcome ~ treatment * time + (1 + time | subject)),
           data)

ames = population_margins(model, data;
                         scenarios=(time=[0, 6, 12],),
                         type=:effects)

sd = second_differences(ames, :treatment, :time, vcov(model);
                       modifier_type=:continuous)
# Interpretation: change in treatment effect over time (time-varying treatment effect)

With DataFrames Ecosystem

using DataFrames, Chain, CSV

# Complete analysis pipeline
results = @chain begin
    population_margins(model, data; scenarios=(region=["north","south","west"],), type=:effects)
    second_differences(_, :income, :region, vcov(model))
    DataFrame(_)
    subset(_, :significant => x -> x .== true)  # Filter to significant interactions
    sort(_, :second_diff, rev=true)  # Sort by effect size
end

# Export for reporting
CSV.write("significant_interactions.csv", results)

Error Handling

Common Error Patterns

Insufficient modifier levels:

# Error: only one modifier level
ames = population_margins(model, data; scenarios=(treated=[1],), type=:effects)
second_differences(ames, :x, :treated, vcov(model))
# → Error: "Need at least 2 modifier levels"

Solution: Ensure scenarios include at least 2 modifier levels.

Missing variable in AME result:

# Error: variable not in AME computation
ames = population_margins(model, data; scenarios=(treated=[0,1],),
                         type=:effects, vars=[:x])
second_differences(ames, :z, :treated, vcov(model))
# → Error: No rows found for variable z

Solution: Include the focal variable in the original population_margins() call or omit vars parameter to include all continuous variables.

Dimension mismatch:

# Error: vcov dimensions don't match gradient dimensions
wrong_vcov = vcov(different_model)
second_differences(ames, :x, :treated, wrong_vcov)
# → Error: "Dimension mismatch: vcov has N parameters but gradients has M"

Solution: Ensure vcov comes from the same model used to compute AMEs.

Validation Workflow

function validate_second_differences(sd_result::DataFrame)
    # Check for numerical issues
    if any(isnan.(sd_result.second_diff)) || any(isinf.(sd_result.second_diff))
        @warn "NaN or Inf values detected in second differences"
    end

    # Check for zero standard errors (indicates no interaction)
    zero_se = sum(sd_result.se .≈ 0.0)
    if zero_se > 0
        @info "$zero_se second difference(s) have zero SE (likely no interaction)"
    end

    # Summary statistics
    n_significant = sum(sd_result.significant)
    println("Significant interactions: $n_significant / $(nrow(sd_result))")
    println("Mean absolute second difference: $(mean(abs.(sd_result.second_diff)))")

    return sd_result
end

# Usage
sd = second_differences(ames, :x, :treated, vcov(model))
validated_sd = validate_second_differences(DataFrame(sd))

Performance Considerations

Computational Complexity

Second differences computation is extremely efficient:

• Primary cost: Computing the underlying AMEs via population_margins()
• Second differences calculation: Negligible additional cost (vector operations only)
• Scalability: O(1) with respect to number of pairwise comparisons

# Performance example
@time ames = population_margins(model, large_data;
                                scenarios=(education=levels,),
                                type=:effects)
# → 95% of total computation time

@time sd = second_differences(ames, :income, :education, vcov(model))
# → <5% of total time, even with many pairwise comparisons

Memory Efficiency

Second differences leverage the gradient information already stored in EffectsResult:

• No additional model evaluations required
• No dataset traversal
• Minimal additional allocations

For large-scale applications with many moderator levels, the computational bottleneck remains the AME computation, not the second differences calculation.

Literature and References

Key Methodological Papers

1. Ai, C., & Norton, E. C. (2003). Interaction terms in logit and probit models. Economics Letters, 80(1), 123-129.

   • Seminal paper demonstrating problems with interpreting interaction coefficients in nonlinear models
   • Established need for marginal effects-based interaction analysis

2. Norton, E. C., Wang, H., & Ai, C. (2004). Computing interaction effects and standard errors in logit and probit models. The Stata Journal, 4(2), 154-167.

   • Practical implementation guidance
   • Standard error computation for interaction effects

3. Greene, W. H. (2010). Testing hypotheses about interaction terms in nonlinear models. Economics Letters, 107(2), 291-296.

   • Hypothesis testing framework for interactions
   • Multiple comparison considerations

4. Karaca-Mandic, P., Norton, E. C., & Dowd, B. (2012). Interaction terms in nonlinear models. Health Services Research, 47(1pt1), 255-274.

   • Applied examples in health economics
   • Interpretation guidance for practitioners

Software Implementation References

• Stata: margins command with dydx() operator and at() option
• R: margins package (Leeper et al.) and interactionTest package
• Python: statsmodels.discrete.discrete_model.Logit.get_margeff() with at specification

Margins.jl's second differences implementation follows this established methodological tradition while leveraging Julia's high-performance computational capabilities for efficient large-scale analysis.

This documentation provides comprehensive coverage of second differences functionality in Margins.jl. For the underlying marginal effects framework, see Mathematical Foundation. For population-based analysis that produces the AME inputs, see Population Scenarios. For integration with robust standard errors, see Advanced Features.