Input data

This method is applicable when: - Data is individual-level with both the outcome variable and auxiliary covariates. - The outcome variable follows a known distribution (e.g., normal).

Required Structure

Outcome Variable: A continuous variable subject to nonresponse (e.g., y).
Covariates for outcome model: Variables trying to model Outcome Variable (e.g., x).
Instrumental Variable: (optional) A variable affecting response probability but not directly related to the outcome.

set.seed(1108)
head(riddles_case1)
#>            x          y     y_true delta
#> 1 0.45763402 -1.7736289 -1.7736289     1
#> 2 0.51906283         NA -0.4653327     0
#> 3 0.46298539         NA -0.5328213     0
#> 4 1.59425312         NA -0.7447774     0
#> 5 1.33728505  0.6900794  0.6900794     1
#> 6 0.07122291  0.7444041  0.7444041     1

y_original <- riddles_case1$y_true # to compare with true values

Engine Configuration

The estimation is configured using exptilt_engine. Critical parameters include the assumed outcome density for respondents and the variance estimation method

The model is specified via a formula: Outcome ~ Response_Covariates | Instrument.

LHS: The outcome variable (Y) containing missing values.
RHS (Left of |): Covariates included in the response propensity model .
RHS (Right of |): The instrumental variable. The method assumes this variable predicts the outcome but is conditionally independent of the response status .

# Exptilt engine configuration
exptilt_config <- exptilt_engine(
  y_dens = 'normal', # Assume f(y|x, delta=1) is Normal
  family = 'logit', # Logit link for response probability
  variance_method = 'bootstrap',
  bootstrap_reps = 5,
  control = list(maxit = 10), # Solver control parameters
  stopping_threshold = 0.01, # Convergence threshold
  standardize = FALSE
)

formula = y ~x
res <- nmar(formula = formula, data = riddles_case1, engine = exptilt_config, trace_level = 0)

Results

The function returns an object containing the estimated coefficients for the response model and the adjusted mean estimate for the outcome.

The adjusted mean is calculated using propensity score weighting, where the weights are derived from the fitted nonignorable response model

print(res)
#> NMAR Result (Exponential tilting)
#> -------------------------------
#> y mean: -1.003197 (0.039987)
#> Converged: TRUE 
#> Variance method: bootstrap 
#> Estimator: exponential_tilting 
#> 
#> Exptilt diagnostics:
#>   Loss value: 0.000000 0.000000 
#>   Iterations: 5 
#>   Variance method: bootstrap 
#>   Bootstrap reps: 5 
#>   Stopping threshold: 0.010000

coef(res)
#> (Intercept)           y 
#>   0.8629256  -0.1722635

Comparing the results to the true population mean (known from the synthetic generation) demonstrates the bias correction.

y_true_mean <- mean(riddles_case1$y_true)
est_mean <- as.numeric(res$y_hat)
se <- res$se

cat('Estimated (NMAR) Mean:', sprintf('%.4f', est_mean),
    ' (SE:', sprintf('%.4f', se), ')\n')
#> Estimated (NMAR) Mean: -1.0032  (SE: 0.0400 )
cat('True Y Mean:          ', sprintf('%.4f', y_true_mean), '\n')
#> True Y Mean:           -1.0105
cat('Naive (MCAR) Mean:     ', sprintf('%.4f', mean(riddles_case1$y, na.rm = TRUE)), '\n')
#> Naive (MCAR) Mean:      -1.0716

Survey designs

The engine supports survey.design objects. When provided, the algorithm incorporates sampling weights into the likelihood estimation and the final propensity score weighting.

The weights affect both the estimation of the respondent density f(y∣x) and the solution to the mean score equation, ensuring population-level inference

# if survey is installed
if (requireNamespace("survey", quietly = TRUE)) {

  library('survey')
  surv_test_weights <- runif(length(riddles_case1$y), 0.9, 1.1)
  surv_test_data <- svydesign(ids = ~1, data = riddles_case1, weights = ~surv_test_weights)

  exptilt_config <- exptilt_engine(
     standardize = FALSE,
     on_failure = "error",
     bootstrap_reps = 5,
     supress_warnings = FALSE,
     control = list(maxit = 15, method = "Newton"),
     family = "logit",
     y_dens = "normal",
     variance_method = "bootstrap",
     stopping_threshold = 1
)


  formula = y ~ x
  res <- nmar(formula = formula, data = surv_test_data, engine = exptilt_config, trace_level = 3)
  print(res)

}
#> [STEP-1] ============================================================ 
#> [STEP-1]   EXPTILT ESTIMATION STARTED 
#> [STEP-1] ============================================================ 
#> [INFO] Running with trace_level = 3 
#> [INFO] Formula: y ~ x 
#> [INFO]  
#> [INFO] -- DATA SUMMARY -- 
#> [INFO]   Total observations:      500 
#> [INFO]   Respondents:             368 (26.4%) 
#> [INFO]   Non-respondents:         132 (73.6%) 
#> [INFO]  
#> [INFO] -- MODEL SPECIFICATION -- 
#> [INFO]   Outcome variable:        y 
#> [INFO]   Auxiliary variables:     x 
#> [INFO]   Missingness predictors:  (intercept only) 
#> [INFO]   Response model:          logit 
#> [INFO]   Outcome density:         normal 
#> [INFO]   Standardization:         disabled 
#> [INFO]  
#> [INFO] -- PARAMETER INITIALIZATION -- 
#> [INFO]   Number of parameters:    2 
#> [DETAIL-3]  
#> [INFO]   Initial values: 
#> [INFO]     (Intercept)               =  0.0481  [response model intercept] 
#> [INFO]     y                         =  0.0928  [outcome effect on response] 
#> [INFO]  
#> [INFO] -- CONDITIONAL DENSITY ESTIMATION -- 
#> [INFO]   Selected distribution:   normal 
#> [INFO]   Density parameters:      3 
#> [INFO]   Fitted density model summary: 
#> [INFO]  
#> 
#> Call:
#> glm(formula = formula, family = gaussian(), data = data)
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -1.06281    0.04866  -21.84   <2e-16 ***
#> x            1.11333    0.06712   16.59   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 0.8711681)
#> 
#>     Null deviance: 558.57  on 367  degrees of freedom
#> Residual deviance: 318.85  on 366  degrees of freedom
#> AIC: 997.58
#> 
#> Number of Fisher Scoring iterations: 2
#> 
#> [INFO]   Computing density matrices... 
#> [INFO]  
#> [INFO] -- NONLINEAR SOLVER (nleqslv) -- 
#> [INFO]   Early stopping threshold: 1.0000 
#> [INFO]   Method:                   Newton 
#> [INFO]   Control parameters: 
#> [INFO]     maxit = 15 
#> [INFO]   Solving... 
#> [INFO]  
#> [INFO]   OK Converged 
#> [INFO]   Iterations:               3 
#> [INFO]   Termination code:         1 
#> [INFO]   Max |score|:              0.000000 
#> [INFO]   Final parameter estimates (scaled): 
#> [INFO]     (Intercept)          =  0.865942  [response intercept] 
#> [INFO]     y                    = -0.159226  [outcome -> response] 
#> [INFO]  
#> [INFO] -- VARIANCE ESTIMATION (Bootstrap) -- 
#> [INFO]   Bootstrap replications:   5 
#> [INFO]   OK Bootstrap complete 
#> [INFO]   Standard error:           0.054672 
#> [INFO]  
#> [RESULT] ============================================================ 
#> [RESULT]   ESTIMATION COMPLETE 
#> [RESULT] ============================================================ 
#> [RESULT]   Mean estimate:            -1.004793 
#> [RESULT]   Standard error:           0.054672 
#> [RESULT]   95% CI:                   [-1.111949, -0.897637] 
#> [INFO]  
#> [INFO]   Response model coefficients: 
#> [INFO]     (Intercept)         :  0.865942  (Intercept) 
#> [INFO]     y                   : -0.159226  (Effect of y on response prob.) 
#> [RESULT] ============================================================ 
#> NMAR Result (Exponential tilting)
#> -------------------------------
#> y mean: -1.004793 (0.054672)
#> Converged: TRUE 
#> Variance method: bootstrap 
#> Estimator: exponential_tilting 
#> 
#> Exptilt diagnostics:
#>   Loss value: 0.000000 0.000000 
#>   Iterations: 3 
#>   Variance method: bootstrap 
#>   Bootstrap reps: 5 
#>   Stopping threshold: 1.000000

Exponential Tilting

Overview

Input data

Required Structure

Engine Configuration

Results

Survey designs