Restricted regression with vimpute

VIM contributors

vimpute(method = "restricted") imputes numeric variables with a least-squares model whose predictions are constrained by edit rules. The rules are written with validate::validator() and solved with ECOSolveR.

This is useful when the completed data must satisfy accounting identities, bounded totals, or row-wise consistency rules.

Minimal Example

The smallest useful example is a lower bound. Here the observed data suggest a line through (x, y), but the missing value is constrained to be at least 4. The imputation indicator y_imp marks the cell filled by vimpute().

library(VIM)

small_data <- data.frame(
  y = c(1, 2, NA),
  x = c(1, 2, 3)
)

rules <- validate::validator(y >= 4)

small_imp <- vimpute(
  small_data,
  method = list(y = "restricted"),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(
    restricted = list(
      rules = rules,
      save_optimization_problem = TRUE
    )
  )
)

small_imp
#>   y x y_imp
#> 1 1 1 FALSE
#> 2 2 2 FALSE
#> 3 4 3  TRUE

The completed value satisfies the rule, and the saved optimization problem can be inspected when save_optimization_problem = TRUE.

all(validate::values(validate::confront(small_imp, rules)))
#> [1] FALSE

names(attr(small_imp, "restricted_optimization_problems"))
#> [1] "y"

Interval Rules

Restricted regression is most natural when variables have row-wise lower and upper bounds. In this example amount is imputed while every completed value must remain between its corresponding lower and upper columns.

set.seed(42)
n <- 100L

bounded_data <- data.frame(
  x1 = stats::runif(n, 0, 10),
  x2 = stats::runif(n, -2, 2)
)
bounded_data$amount <- 20 + 1.5 * bounded_data$x1 - 0.75 * bounded_data$x2
bounded_data$lower <- bounded_data$amount - 0.5
bounded_data$upper <- bounded_data$amount + 0.5
bounded_data <- bounded_data[, c("amount", "x1", "x2", "lower", "upper")]

missing_idx <- sample.int(n, 20L)
bounded_missing <- bounded_data
bounded_missing$amount[missing_idx] <- NA_real_

bounded_rules <- validate::validator(
  amount >= lower,
  amount <= upper,
  amount >= 0,
  lower <= upper
)

bounded_imp <- vimpute(
  bounded_missing,
  method = list(amount = "restricted"),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(amount = list(rules = bounded_rules))
)

head(bounded_imp[bounded_imp$amount_imp, c("amount", "lower", "upper")])
#>      amount    lower    upper
#> 5  28.29882 27.79882 28.79882
#> 10 32.06915 31.56915 32.56915
#> 14 23.97324 23.47324 24.47324
#> 24 33.93099 33.43099 34.43099
#> 25 20.57759 20.07759 21.07759
#> 27 24.59544 24.09544 25.09544
completed <- as.data.frame(bounded_imp)[, names(bounded_data), drop = FALSE]

sum(bounded_imp$amount_imp)
#> [1] 20
all(validate::values(validate::confront(completed, bounded_rules)))
#> [1] TRUE

Rules can also involve categorical predictors and conditional logic. The restricted learner uses the edit rules only to constrain the imputed target; the remaining columns can be numeric or categorical predictors.

categorical_data <- bounded_missing
categorical_data$c1 <- factor(rep(c("A", "B", "C"), length.out = n))
categorical_data$c2 <- factor(rep(c("B", "C", "A"), length.out = n))
categorical_data$c3 <- factor(rep(c("C", "A", "B"), length.out = n))

conditional_rules <- validate::validator(
  amount >= lower,
  amount <= upper,
  amount >= 0,
  lower <= upper,
  (c1 == "A") + (c2 == "A") + (c3 == "A") <= 2,
  if (c1 == "A") amount >= 2
)

categorical_imp <- vimpute(
  categorical_data,
  method = list(amount = "restricted"),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(amount = list(rules = conditional_rules))
)

all(validate::values(validate::confront(categorical_imp, conditional_rules)))
#> [1] TRUE

Formula Control

formula can be used with method = "restricted" to choose the predictors in the regression model while the edit rules still constrain the final predictions. The left-hand side must be the untransformed target variable.

n_formula <- 80L
formula_data <- data.frame(
  y = 10 + 2 * seq_len(n_formula),
  x_signal = seq_len(n_formula),
  x_noise = 10 + 2 * seq_len(n_formula),
  lower = 0,
  upper = 200
)

formula_missing_idx <- c(12L, 24L, 36L)
expected_y <- formula_data$y[formula_missing_idx]
formula_data$y[formula_missing_idx] <- NA_real_
formula_data$x_noise[formula_missing_idx] <- -1000

formula_rules <- validate::validator(y >= lower, y <= upper, lower <= upper)

formula_imp <- vimpute(
  formula_data,
  method = list(y = "restricted"),
  formula = list(y = y ~ x_signal),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(y = list(rules = formula_rules))
)

formula_imp[formula_missing_idx, c("y", "x_signal", "x_noise", "y_imp")]
#>     y x_signal x_noise y_imp
#> 12 34       12   -1000  TRUE
#> 24 58       24   -1000  TRUE
#> 36 82       36   -1000  TRUE
max(abs(formula_imp$y[formula_missing_idx] - expected_y))
#> [1] 0

Robust Restricted Regression

Set robust = TRUE to use a Huber loss instead of ordinary least squares. This is helpful when the observed data contain influential outliers but the final imputed values must still obey the edit rules.

robust_rules <- validate::validator(y >= 0)

outlier_data <- data.frame(
  x = 0:21,
  y = 1 + 2 * (0:21)
)
outlier_data$y[21L] <- 500
outlier_data$y[22L] <- NA_real_

ordinary_imp <- vimpute(
  outlier_data,
  method = list(y = "restricted"),
  formula = list(y = y ~ x),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(y = list(rules = robust_rules))
)

robust_imp <- vimpute(
  outlier_data,
  method = list(y = "restricted"),
  formula = list(y = y ~ x),
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(y = list(
    rules = robust_rules,
    robust = TRUE,
    huber_k = 1.345
  ))
)

c(
  ordinary = ordinary_imp$y[22L],
  robust = robust_imp$y[22L],
  expected_without_outlier = 43
)
#>                 ordinary                   robust expected_without_outlier 
#>                      130                       57                       43

Synthetic LSE Data

The package includes lse_synthetic and lse_synthetic_rules, a compact synthetic business-statistics example with many linear edit rules. The rules ship as plain data frames of rule text (so the data set does not depend on validate); validate::validator(.data = ) turns a rule set into a validator. The test suite uses this data at larger scale; the vignette uses a smaller slice so it can be knitted quickly.

data("lse_synthetic", package = "VIM")
data("lse_synthetic_rules", package = "VIM")

lse_rules <- validate::validator(.data = lse_synthetic_rules$edit)

numeric_cols <- c(
  "persons_employed",
  "employees_paid",
  "self_employed",
  "employees_male",
  "employees_female",
  "employees_blue_collar",
  "employees_white_collar",
  "apprentices",
  "marginal_employees",
  "turnover_total",
  "turnover_domestic",
  "turnover_exports",
  "e_commerce_turnover",
  "material_costs",
  "purchased_services",
  "rents_leasing",
  "other_operating_expense",
  "intermediate_consumption",
  "gross_value_added",
  "personnel_costs",
  "wages_salaries",
  "social_security_costs",
  "other_personnel_costs",
  "gross_operating_surplus",
  "investments_tangible",
  "investment_machinery",
  "investment_buildings",
  "investment_software"
)

edit_cols <- c(
  "reporting_year",
  "onace_section",
  "onace_group",
  "nuts2",
  "data_source",
  "survey_mode",
  "employment_size_class",
  "turnover_size_class",
  numeric_cols
)

lse_complete <- lse_synthetic[seq_len(80L), edit_cols]

missing_map <- list(
  turnover_total = c(2L, 15L),
  intermediate_consumption = c(13L, 18L),
  gross_value_added = c(24L, 31L),
  personnel_costs = c(39L, 45L),
  investments_tangible = c(52L, 60L)
)

lse_missing <- lse_complete
for (var in names(missing_map)) {
  lse_missing[missing_map[[var]], var] <- NA_real_
}

restricted_formulas <- setNames(
  lapply(names(missing_map), function(var) {
    stats::reformulate(setdiff(numeric_cols, var), response = var)
  }),
  names(missing_map)
)

lse_imp <- vimpute(
  lse_missing,
  method = "restricted",
  formula = restricted_formulas,
  pmm = FALSE,
  sequential = FALSE,
  learner_params = list(
    restricted = list(
      rules = lse_rules,
      save_optimization_problem = TRUE
    )
  )
)

lse_completed <- as.data.frame(lse_imp)[, edit_cols, drop = FALSE]

sum(as.data.frame(lse_imp)[paste0(names(missing_map), "_imp")])
#> [1] 10
all(validate::values(validate::confront(lse_completed, lse_rules)))
#> [1] TRUE
sort(names(attr(lse_imp, "restricted_optimization_problems")))
#> [1] "gross_value_added"        "intermediate_consumption"
#> [3] "investments_tangible"     "personnel_costs"         
#> [5] "turnover_total"

The same pattern scales to many variables: define a validate::validator object for the edits, provide one formula per imputed target when you need predictor control, and pass the rules through learner_params for either the method name (restricted) or the individual variables.