GeoLift-SDID Matrix Dimension Handling

This document outlines the technical approaches for handling matrix dimension mismatches in the GeoLift-SDID implementation.

Core Problem: Dimension Constraints

The Synthetic Difference-in-Differences estimator requires matrix operations that impose specific dimensional constraints:

1. Fundamental matrices:

   • Y1_pre: Treatment unit pre-intervention outcomes (dimension: 1 × pre_periods)

   • Y1_post: Treatment unit post-intervention outcomes (dimension: 1 × post_periods)

   • Y0_pre: Control units pre-intervention outcomes (dimension: n_controls × pre_periods)

   • Y0_post: Control units post-intervention outcomes (dimension: n_controls × post_periods)

2. Critical matrix operations:

   • Omega weights calculation: omega = solve_omega(Y0_pre, Y1_pre)

   • Lambda weights calculation: lambda_ = solve_lambda(Y0_pre, Y1_pre)

   • Synthetic control generation: Y_synth = omega @ Y0_post

   • ATT calculation: att = Y1_post.mean() - Y_synth.mean()

Technical Challenge: Dimension Mismatch

When pre-intervention periods ≠ post-intervention periods (e.g., 60 pre-periods vs. 30 post-periods), direct calculation of standard errors fails due to matrix dimension incompatibility.

Key error pattern:

ValueError: Input operand has a mismatch in its core dimension: operand[0] has 60 but operand[1] has 30

Implementation Solutions

1. Automated Fallback Detection

def calculate_standard_error(model, Y1_pre, Y1_post, Y0_pre, Y0_post):
    """Calculate standard error with fallback to bootstrap if needed"""
    try:
        # Attempt standard matrix-based calculation
        se = model.se
        if np.isnan(se) or se <= 0:
            raise ValueError("Invalid standard error from matrix calculation")
        return se
    except (AttributeError, ValueError, np.linalg.LinAlgError) as e:
        logging.warning(f"Matrix-based SE calculation failed: {e}")
        logging.info("Falling back to bootstrap standard error calculation")
        return bootstrap_standard_error(Y1_pre, Y1_post, Y0_pre, Y0_post)

2. Bootstrap Standard Error Implementation

def bootstrap_standard_error(Y1_pre, Y1_post, Y0_pre, Y0_post, n_bootstrap=500):
    """
    Calculate standard error via bootstrap when matrix dimensions 
    prevent direct calculation
    
    Parameters:
    -----------
    Y1_pre : numpy.ndarray
        Treatment unit pre-intervention outcomes (1, pre_periods)
    Y1_post : numpy.ndarray
        Treatment unit post-intervention outcomes (1, post_periods)
    Y0_pre : numpy.ndarray
        Control units pre-intervention outcomes (n_controls, pre_periods)
    Y0_post : numpy.ndarray
        Control units post-intervention outcomes (n_controls, post_periods)
    n_bootstrap : int, optional
        Number of bootstrap samples to use, defaults to 500
        
    Returns:
    --------
    dict
        Dictionary containing standard error and bootstrap diagnostics
    """
    bootstrap_estimates = []
    bootstrap_metrics = {
        'n_iterations': n_bootstrap,
        'n_successful': 0,
        'failure_reasons': {}
    }
    
    for i in range(n_bootstrap):
        # Sample control units with replacement
        n_controls = Y0_pre.shape[0]
        try:
            # 1. Generate bootstrap sample indices
            indices = np.random.choice(n_controls, size=n_controls, replace=True)
            
            # 2. Create bootstrap samples
            Y0_pre_boot = Y0_pre[indices, :]
            Y0_post_boot = Y0_post[indices, :]
            
            # 3. Calculate unit weights (omega)
            try:
                omega = calculate_omega_weights(Y0_pre_boot, Y1_pre)
                
                # 4. Validate weights - check for NaN or extreme values
                if np.any(np.isnan(omega)) or np.any(np.isinf(omega)):
                    reason = "Invalid weights (NaN or Inf)"
                    bootstrap_metrics['failure_reasons'][reason] = bootstrap_metrics['failure_reasons'].get(reason, 0) + 1
                    continue
                    
                # 5. Calculate synthetic control
                Y_synth = Y0_post_boot @ omega  # Matrix multiplication
                
                # 6. Calculate ATT for this bootstrap iteration
                att_boot = np.mean(Y1_post) - np.mean(Y_synth)
                bootstrap_estimates.append(att_boot)
                bootstrap_metrics['n_successful'] += 1
                
            except np.linalg.LinAlgError as e:
                reason = f"Linear algebra error: {str(e)}"
                bootstrap_metrics['failure_reasons'][reason] = bootstrap_metrics['failure_reasons'].get(reason, 0) + 1
                continue
                
        except (ValueError, TypeError) as e:
            reason = f"Sampling error: {str(e)}"
            bootstrap_metrics['failure_reasons'][reason] = bootstrap_metrics['failure_reasons'].get(reason, 0) + 1
            continue
    
    # Log bootstrap performance metrics
    success_rate = bootstrap_metrics['n_successful'] / n_bootstrap * 100
    logging.info(f"Bootstrap completed with {success_rate:.1f}% success rate ({bootstrap_metrics['n_successful']}/{n_bootstrap} iterations)")
    
    if bootstrap_metrics['n_successful'] < 0.5 * n_bootstrap:
        logging.warning(f"Low bootstrap success rate. Failure reasons: {bootstrap_metrics['failure_reasons']}")
    
    # Calculate standard error as standard deviation of bootstrap estimates
    if bootstrap_estimates:
        bootstrap_metrics['standard_error'] = np.std(bootstrap_estimates)
        bootstrap_metrics['bootstrap_samples'] = len(bootstrap_estimates)
        return bootstrap_metrics['standard_error']
    else:
        logging.error("Bootstrap failed completely - no valid iterations")
        return np.nan

3. Direct Synthetic Control Generation

For cases where weight calculation fails due to dimension constraints, direct outcome-based synthetic control generation is implemented:

def generate_synthetic_outcomes_direct(Y1_pre, Y0_pre, Y0_post):
    """
    Generate synthetic outcomes directly when standard weight
    calculation fails
    
    This approach solves for weights that minimize the pre-period 
    difference between treatment and synthetic control
    """
    try:
        # Fit regression model to predict treatment outcomes from control outcomes
        model = LinearRegression(fit_intercept=False)
        model.fit(Y0_pre.T, Y1_pre.T)
        
        # Use coefficients as weights
        omega = model.coef_
        
        # Normalize weights to sum to 1
        omega = omega / np.sum(omega)
        
        # Generate synthetic control for post-period
        Y_synth = omega @ Y0_post
        
        return Y_synth
    except Exception as e:
        print(f"Direct synthetic generation failed: {e}")
        return None

Multiple Dimensions Handling Strategy

The implementation handles dimension mismatches through a cascading approach:

1. Attempt standard calculation: Try using the core SDID implementation.

2. Detect failure: Catch exceptions related to dimension mismatch.

3. First fallback: Use bootstrap resampling for standard error calculation.

4. Second fallback: Generate synthetic control directly if bootstrap fails.

5. Final fallback: Return NaN for standard error with warning if all methods fail.

Implementation in Single-Cell vs. Multi-Cell Analysis

Single-Cell Implementation

The single-cell implementation has specific handling for dimension challenges:

def _generate_results(self, model, Y1_pre, Y1_post, Y0_pre, Y0_post, omega_weights, lambda_weights):
    """Generate results with robust standard error calculation"""
    calculation_method = "model"
    
    try:
        # Try to get ATT directly from model
        att = model.att
    except AttributeError as e:
        # Calculate manually if model API fails
        calculation_method = "manual"
        logging.info(f"Model ATT retrieval failed, calculating manually: {e}")
        Y_sc = omega_weights @ Y0_post
        att = np.mean(Y1_post) - np.mean(Y_sc)
    
    # Handle standard error calculation
    se_calculation = "model"
    try:
        se = model.se
    except AttributeError as e:
        # Fallback to bootstrap calculation
        logging.info(f"Model SE retrieval failed, using bootstrap: {e}")
        se_calculation = "bootstrap"
        se = calculate_standard_error(Y1_pre, Y1_post, Y0_pre, Y0_post)
    
    # Validate SE
    if np.isnan(se) or se <= 0:
        logging.warning(f"Invalid SE value ({se}), results may be unreliable")
    
    # Calculate p-value
    z = att / se if se > 0 else np.nan
    p_value = 2 * (1 - stats.norm.cdf(abs(z))) if not np.isnan(z) else np.nan
    
    # Comprehensive results with metadata
    results = {
        'att': att,
        'se': se,
        'pvalue': p_value,
        'ci_lower': att - 1.96 * se,
        'ci_upper': att + 1.96 * se,
        'meta': {
            'att_calculation': calculation_method,
            'se_calculation': se_calculation,
            'omega_weight_sum': np.sum(omega_weights),  # Diagnostic for weight quality
            'n_control': Y0_pre.shape[0],
            'pre_periods': Y1_pre.shape[1],
            'post_periods': Y1_post.shape[1]
        }
    }
    
    return results

Multi-Cell Implementation

The multi-cell implementation handles dimension constraints through aggregation:

def _process_results(self, results_per_unit):
    """Aggregate results across multiple treatment units"""
    # Extract ATT values and standard errors
    atts = [r['att'] for r in results_per_unit.values()]
    ses = [r['se'] for r in results_per_unit.values()]
    
    # Calculate aggregate metrics
    combined_att = np.mean(atts)
    
    # Properly combine standard errors
    # Use correct formula for variance of mean (divide by n²)
    if any(np.isnan(ses)):
        # Fall back to bootstrap if any standard errors are missing
        logging.warning("Missing standard errors detected, using bootstrap for aggregate SE")
        combined_se = bootstrap_aggregate_se(results_per_unit)
    else:
        # Calculate combined standard error - CORRECTED FORMULA
        # This is the standard error of the mean of the ATTs
        # var(mean) = sum(var_i) / n² = sum(se_i²) / n²
        combined_se = np.sqrt(np.sum(np.array(ses)**2) / len(ses)**2)
        logging.info(f"Combined SE calculated from {len(ses)} unit-level SEs")
    
    # Calculate p-value and confidence intervals
    z = combined_att / combined_se if combined_se > 0 else np.nan
    p_value = 2 * (1 - stats.norm.cdf(abs(z)))
    
    # Create comprehensive result dictionary
    aggregate_results = {
        'aggregate_att': combined_att,
        'aggregate_se': combined_se,
        'aggregate_pvalue': p_value,
        'aggregate_ci_lower': combined_att - 1.96 * combined_se,
        'aggregate_ci_upper': combined_att + 1.96 * combined_se,
        'unit_results': results_per_unit,
        'meta': {
            'n_units': len(atts),
            'unit_weights': 'equal',  # Equal weighting across units
            'calculation_method': 'bootstrap' if any(np.isnan(ses)) else 'analytical',
            'heterogeneity': np.std(atts) / abs(combined_att) if combined_att != 0 else np.nan  # CV of effects
        }
    }
    
    return aggregate_results

Performance Considerations

Bootstrap standard error calculation is computationally intensive:

1. Time complexity: O(n_bootstrap × n_controls × max(pre_periods, post_periods))

2. Memory usage: Proportional to n_bootstrap × n_controls × (pre_periods + post_periods)

Optimizations Implemented

1. Vectorized Operations

   • Direct matrix multiplication instead of element-wise operations

   • Numpy array operations for performance-critical computations

2. Early Validation

   • Weight validation before expensive calculations

   • Specialized exception handling for different failure modes

3. Progressive Fallbacks

   • Tiered approach from direct calculation to bootstrap

   • Success rate tracking to identify systematic issues

4. Parallel Processing Option

   # Using joblib for parallelization
   from joblib import Parallel, delayed
   
   def parallel_bootstrap(Y1_pre, Y1_post, Y0_pre, Y0_post, n_bootstrap=500, n_jobs=-1):
       """Parallel implementation of bootstrap standard error"""
       # Run bootstrap iterations in parallel
       results = Parallel(n_jobs=n_jobs)(delayed(single_bootstrap_iteration)(
           Y1_pre, Y1_post, Y0_pre, Y0_post, i) for i in range(n_bootstrap))
       
       # Filter out failed iterations (None results)
       valid_results = [r for r in results if r is not None]
       return np.std(valid_results) if valid_results else np.nan
   

Calibration Guidelines

Dataset Size	Recommended Bootstrap Iterations	Parallelization
Small (<50 units)	500 iterations	Single process
Medium (50-200 units)	300 iterations	4 processes
Large (>200 units)	200 iterations	8+ processes

These guidelines balance statistical precision with computational efficiency. For extremely large datasets, consider subsampling control units before bootstrapping.