1 Ill-conditioned regression problems

In the companion vignette Pure Inverse Problems we used inverse_ce() to recover a distribution from an ill-posed problem in which the number of unknowns exceeds the number of observations. That is not the only situation in which traditional estimation methods run into trouble. In many applied settings — even in cases where the sample is relatively large compared to the number of parameters — traditional methods may still behave poorly.

Such problems are quite common in the social sciences, and especially in economics, where the data are often observational rather than the outcome of a carefully designed and replicable experiment. When the data are collected passively we cannot rule out possibilities of strong dependence, including in the form of multicollinearity. The consequences of the latter are well known: standard errors are inflated, coefficients lose statistical significance, point estimates become erratic and sensitive to small changes in the data, and the numerical procedure becomes less stable. These situations are often encountered in impact evaluations in which the analyst adds polynomial and interaction terms. It is a sensible strategy to capture non-linearities, but it may also backfire by creating regressors with near-exact linear dependencies.

This vignette introduces the Generalized Maximum Entropy (GME) and Generalized Cross-Entropy (GCE) estimators (Golan et al., 1996; Golan, 2008), which are designed to deal with these ill-conditioned problems, and shows how to fit them with the linreg() function of the infometrics package. At the end of this vignette, we also demonstrate its sibling, the linreg_iv() function, which deals with instrumental variables.

But before demonstrating a working example, let’s highlight the underlying concepts and core mechanics of the GME and GCE estimators, and start with the general linear model (GLM)

\[\begin{equation} \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e} \tag{1.1} \end{equation}\]

where \(\mathbf{y} = (y_1, \dots, y_T)'\) is the \(T\)-dimensional response, \(\mathbf{X}\) is a known \(T \times K\) matrix of explanatory variables, \(\boldsymbol{\beta} = (\beta_1, \dots, \beta_K)'\) is the \(K\)-vector of unknown parameters, and \(\mathbf{e} = (e_1, \dots, e_T)'\) is the noise, which absorbs sampling error, non-sampling error, and model misspecification.

2 The GME formulation

The idea behind GME is to treat the unknown coefficients and errors as random variables and to reparameterize them as expected values over bounded support spaces. Each estimated coefficient is written as the mean of a discrete distribution \(\mathbf{p}_k\) on a support \(\mathbf{z}_k = (z_{k1}, \dots, z_{kM})\), and each error as the mean of a distribution \(\mathbf{w}_t\) on a support \(\mathbf{v} = (v_1, \dots, v_J)\):

\[\begin{equation} \beta_k = \sum_{m=1}^M z_{km}\, p_{km}, \qquad e_t = \sum_{j=1}^J v_j\, w_{tj}. \tag{2.1} \end{equation}\]

Because \(\mathbf{P}\) and \(\mathbf{W}\) are probability distributions, they must satisfy the additivity (normalization) constraints

\[\begin{equation} \sum_{m=1}^M p_{km} = 1 \;\; (k = 1, \dots, K), \qquad \sum_{j=1}^J w_{tj} = 1 \;\; (t = 1, \dots, T), \tag{2.2} \end{equation}\]

with \(p_{km} \ge 0\) and \(w_{tj} \ge 0\). The support \(\mathbf{z}_k\) must be wide enough to contain the true value of \(\beta_k\), and the error support \(\mathbf{v}\) — symmetric around zero so that the noise is mean-zero — must be wide enough to contain every realized error.

The traditional choice is the three-sigma rule (Pukelsheim, 1994): a half-width of three times an estimate of \(\mathbf{y}\)’s standard deviation. It is worth being careful with it here. Pukelsheim’s rule bounds a single draw (at least 95% of the mass of any unimodal distribution lies within three standard deviations), but the error support has to bound all \(T\) errors at once. The largest of \(T\) draws grows like \(\sigma\sqrt{2\log T}\), which passes \(3\sigma\) at around \(T \approx 90\) and keeps growing, so a fixed three-sigma support becomes too narrow in even moderately large samples. linreg() therefore defaults the error half-width to \(\max(3, \sqrt{2\log T})\,\mathrm{sd}(y)\), and we use the same rule in the examples below. Section 5.4 shows what goes wrong when the support does not fit.

Among the infinitely many distributions consistent with the data we select the one that is maximally uninformative, i.e. that maximizes the Shannon entropy (Shannon, 1948)

\[\begin{equation} H(\mathbf{P}, \mathbf{W}) = - \sum_{k}\sum_{m} p_{km} \log p_{km} - \sum_{t}\sum_{j} w_{tj} \log w_{tj}. \tag{2.3} \end{equation}\]

Maximizing (2.3) subject to the data constraint (1.1), the reparameterization (2.1) and the normalizations (2.2) yields the following solutions (see Golan et al., 1996, Ch. 6; Golan, 2008, Section 6.1)

\[\begin{equation} \hat{p}_{km} = \frac{\exp\!\left(-z_{km} \sum_t \hat{\lambda}_t x_{tk}\right)} {\Omega_k(\hat{\boldsymbol{\lambda}})}, \qquad \hat{w}_{tj} = \frac{\exp\!\left(-\hat{\lambda}_t v_{j}\right)} {\Psi_t(\hat{\boldsymbol{\lambda}})}, \tag{2.4} \end{equation}\]

where \(\Omega_k\) and \(\Psi_t\) are the normalizing (partition) functions and \(\hat{\boldsymbol{\lambda}}\) is a \(T\)-vector of Lagrange multipliers on the data constraint. As with the pure inverse problem, \(\lambda_t\) measures the information content of observation \(t\): a multiplier near zero means that row carries little information.

3 From the primal to the concentrated dual

Substituting the solution (2.4) back into the Lagrangian collapses the primal problem into an unconstrained concentrated (dual) problem in \(\boldsymbol{\lambda}\) alone:

\[\begin{equation} \min_{\boldsymbol{\lambda}} \left\{ \sum_t y_t \lambda_t + \sum_k \log \Omega_k(\boldsymbol{\lambda}) + \sum_t \log \Psi_t(\boldsymbol{\lambda}) \right\}. \tag{3.1} \end{equation}\]

Please note that the \(\boldsymbol{\lambda}\)’s are the Lagrange multipliers associated with equation (1.1) only; we omit the Lagrange multipliers of the normalization constraints (2.2) because the closed-form solutions in (2.4) already satisfy the normalization by construction. It is also worth noting that the gradient of (3.1) is the moment condition \(\mathbf{y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{e} = \mathbf{0}\), so at the optimum the data are reproduced exactly, \(\mathbf{y} = \mathbf{X}\hat{\boldsymbol{\beta}} + \hat{\mathbf{e}}\); the estimator simply chooses the least-informative split of \(\mathbf{y}\) into signal and noise.

Solving in \(\boldsymbol{\lambda}\) rather than in \((\mathbf{P}, \mathbf{W})\) is a large computational saving. Suppose \(T = 100\), \(K = 5\), \(M = 3\), and \(J = 3\). The primal problem requires estimation of \(K \times M + T \times J = 15 + 300 = 315\) unknowns, whereas the concentrated dual is solved in just \(T = 100\) Lagrange multipliers. linreg() always uses the dual form (BFGS method with an analytic gradient), recovering \(\hat{\boldsymbol{\lambda}}\) and then mapping it back through (2.4) and (2.1) to recover \(\hat{\boldsymbol{\beta}}\) and \(\hat{\mathbf{e}}\).

4 The GCE extension: incorporating prior information

If we have prior beliefs about the data-generating process and the noise — summarized as prior probabilities \(\left\{ \mathbf{p}^0_k, \mathbf{w}^0_t \right\}\) over their respective support spaces \(\left\{ \mathbf{z}_k, \mathbf{v} \right\}\) — we can fold them in by replacing the entropy with the cross-entropy (or Kullback–Leibler divergence) between the estimated distributions and their prior information:

\[\begin{equation} D(\mathbf{P} \,\Vert\, \mathbf{P}^0) + D(\mathbf{W} \,\Vert\, \mathbf{W}^0) = \sum_k \sum_m p_{km} \log\!\left(p_{km} / p^0_{km}\right) + \sum_t \sum_j w_{tj} \log\!\left(w_{tj} / w^0_{tj}\right). \tag{4.1} \end{equation}\]

Minimizing the divergence (4.1) subject to the same constraints is the GCE estimator. A uniform prior reproduces GME; a non-uniform prior pulls the estimates toward the prior mean. In linreg() prior information on \(\boldsymbol{\beta}\) can be summarized in the single argument p0 (a \(K \times M\) matrix whose rows sum to one), and our subjective beliefs on \(\mathbf{e}\) can be incorporated in w0 (a \(T \times J\) matrix whose rows sum to one): you may leave either one or both blank to opt for the default (uniform distribution), or supply them for GCE.

5 Worked examples

We now put GME and GCE estimation to work on a deliberately ill-conditioned data set, in the spirit of the collinearity illustration in Section 6.6 of Golan (2008). linreg() also lets you place weights on the signal and noise, which the worked example presents as well.

set.seed(10)
n  <- 200
x1 <- runif(n, 0, 20)                 # clean regressor
x2 <- runif(n, 0, 20)
x3 <- x2 + rnorm(n, 0, 0.05)          # near-collinear with x2 (extreme)

b  <- c(1, -2, 3, 0)                  # (Intercept), x1, x2, x3
y  <- b[1] + b[2] * x1 + b[3] * x2 + b[4] * x3 + rnorm(n, 0, sqrt(2))
dat <- data.frame(y, x1, x2, x3)

truth <- c(`(Intercept)` = 1, x1 = -2, x2 = 3, x3 = 0)
cor(x2, x3)                           # near-perfect collinearity
#> [1] 0.9999637

ols <- lm(y ~ x1 + x2 + x3, data = dat)

Zvec <- seq(-6, 6, length.out = 5)   # coefficient support (shared)
sig  <- sd(dat$y)                    # standard deviation of y
k    <- max(3, sqrt(2 * log(n)))     # sample-size-aware multiple
vvec <- seq(-k * sig, k * sig, length.out = 5)

gme <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = vvec)             # nu = 0.5, uniform priors => GME

se <- function(fit) sqrt(diag(vcov(fit)))
tab1 <- data.frame(
  Truth    = truth,
  OLS      = coef(ols),
  `OLS SE` = se(ols),
  GME      = coef(gme),
  `GME SE` = se(gme),
  check.names = FALSE
)
knitr::kable(tab1, digits = 3,
             caption = "OLS versus GME under near-perfect collinearity of x2 and x3.")

# we can recover probabilities associated with the  
# prior information, using the inverse_ce() function  
temp <- data.frame(y=2.8, s1 = Zvec[1], s2 = Zvec[2], 
                   s3= Zvec[3], s4= Zvec[4], s5= Zvec[5])
p0_x2 <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 - 1, data=temp)$p_hat
p0_other_x <- rep(1/5, 5)
p0 <- rbind(p0_other_x, p0_other_x, p0_x2, p0_other_x)

gce <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = vvec, p0 = p0)    # GCE: informative signal prior

tab2 <- data.frame(
  Truth          = truth,
  `OLS`      = coef(ols),
  `GME` = coef(gme),
  `GCE`   = coef(gce),
  check.names     = FALSE
)
knitr::kable(tab2, digits = 3,
             caption = "OLS, GME (uniform prior), and GCE (informative prior).")

nus <- c(0.1, 0.5, 0.9)
fit_nu <- lapply(nus, function(nu)
  linreg(y ~ x1 + x2 + x3, data = dat, Z = Zvec, v = vvec, p0=p0, nu = nu))

tab3 <- rbind(sapply(fit_nu, coef),
              `R-squared` = sapply(fit_nu, `[[`, "r.squared"),
              `Signal S`  = sapply(fit_nu, `[[`, "S"))
tab3 <- cbind(Truth = c(truth, `R-squared` = NA, `Signal S` = NA), round(tab3, 4))
colnames(tab3) <- c("Truth", paste0("nu = ", nus))
knitr::kable(tab3, digits = 3,
             caption = "Coefficients, ordinary R-squared, and signal normalized entropy across nu.")

c(foc_residual = gme$foc_residual, converged = gme$converged)
#> foc_residual    converged 
#> 7.520817e-05 1.000000e+00

bad <- linreg(y ~ x1 + x2 + x3, data = dat,
              Z = Zvec, v = c(-0.01, 0, 0.01))   # far too narrow
#> Warning in linreg(y ~ x1 + x2 + x3, data = dat, Z = Zvec, v = c(-0.01, 0, :
#> optim did not converge (code 1).
#> Warning in linreg(y ~ x1 + x2 + x3, data = dat, Z = Zvec, v = c(-0.01, 0, : The
#> GME/GCE dual appears unbounded: the supports are too narrow to represent y
#> (max|y - X*beta - e| = 131). Estimates are unreliable and beta may sit on the
#> boundary of the signal support. Widen the error support 'v' (and/or 'Z').
round(coef(bad), 3)
#> (Intercept)          x1          x2          x3 
#>           6          -6          -6           6
c(foc_residual = bad$foc_residual, converged = bad$converged)
#> foc_residual    converged 
#>     130.5888       0.0000

6 Inference: entropy diagnostics and hypothesis tests

Section 6.4 of Golan (2008) equips the GME/GCE estimator with information-theoretic analogues of the usual regression diagnostics. All of them are by-products of the fitted probability distributions, so they cost nothing extra to report.

c(signal_S    = gme$S,
  pseudo_R2   = 1 - gme$S,
  ordinary_R2 = gme$r.squared)
#>    signal_S   pseudo_R2 ordinary_R2 
#>  0.96476874  0.03523126  0.99422254

summary(gme)
#> 
#> Call:
#> linreg(formula = y ~ x1 + x2 + x3, data = dat, Z = Zvec, v = vvec)
#> 
#> Residuals:
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -3.89600 -0.94700  0.18880  0.09612  1.17100  5.44800 
#> 
#> Coefficients (entropy-ratio test, H0: beta_k = 0):
#>             Estimate Std. Error     ER Pr(>Chi)    
#> (Intercept)  0.14701    0.31082  0.001 0.969605    
#> x1          -1.89152    0.02079 15.013 0.000107 ***
#> x2           1.48305    2.25323  0.269 0.604054    
#> x3           1.49155    2.25402  0.272 0.601941    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> nu = 0.5   R-squared: 0.9942   signal norm. entropy S: 0.9648 
#> Overall ER test (H0: all coefficients = 0): ER = 56.29, df = 4, Pr(>Chi) = 1.74e-11
#> Convergence code: 0 (converged) 
#> 
#> Note: ER is the entropy-ratio statistic for H0: beta_k = 0, computed as
#> 2*(Hstar_unrestricted - Hstar_restricted_k) where the restricted fit zeroes
#> row k of the signal support and refits; under H0, ER ~ chi-squared(1) and
#> Pr(>Chi) is the upper-tail p-value (clamped ER >= 0). Std. Error is the
#> Golan (2008, p. 96) asymptotic SE, shown for reference.

wald_test <- function(fit, R, r = NULL) {
  R <- rbind(R)                                 # allow a plain vector
  if (is.null(r)) r <- rep(0, nrow(R))
  d <- as.vector(R %*% coef(fit)) - r
  W <- as.numeric(t(d) %*% solve(R %*% vcov(fit) %*% t(R), d))
  c(W = W, df = nrow(R), p.value = pchisq(W, df = nrow(R), lower.tail = FALSE))
}

# H0: beta_x3 = 0  (single restriction; equals the squared t-ratio)
wald_test(gme, c(0, 0, 0, 1))
#>         W        df   p.value 
#> 0.4378855 1.0000000 0.5081448

# H0: beta_x2 = beta_x3 = 0  (joint test of the collinear pair)
wald_test(gme, rbind(c(0, 0, 1, 0),
                     c(0, 0, 0, 1)))
#>        W       df  p.value 
#> 24079.67     2.00     0.00

7 Instrumental variables with linreg_iv()

When a regressor is correlated with the error — through omitted variables, measurement error, or simultaneity — OLS is biased and inconsistent. The classical remedy is an instrument: a variable that shifts the endogenous regressor but is unrelated to the error. linreg_iv() is the instrumental-variables sibling of linreg(). It identifies \(\boldsymbol{\beta}\) through the instrument moments \(\mathbf{W}'(\mathbf{y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{e}) = \mathbf{0}\), where \(\mathbf{W}\) is the instrument matrix (Golan, 2008, pp. 89–91), and keeps the same support-and-entropy machinery. It uses a matrix interface, linreg_iv(y, X, IV, Z), and supports both just-identified (\(\text{ncol}(IV) = \text{ncol}(X)\)) and over-identified systems.

We simulate an endogenous regressor \(x_e\) that shares a shock \(u\) with the error, and an instrument \(z\) that drives \(x_e\) but is independent of \(u\):

set.seed(125)
n  <- 200
z  <- rnorm(n)                       # instrument
u  <- rnorm(n)                       # shared shock => endogeneity
xe <- 0.7 * z + u + rnorm(n)         # endogenous regressor
y  <- 1 + 1.5 * xe + u               # true slope = 1.5

X   <- cbind(`(Intercept)` = 1, xe = xe)   # design: intercept + endogenous xe
IV  <- cbind(`(Intercept)` = 1, z  = z)    # instruments: intercept (self) + z
Ziv <- matrix(c(-10, 0, 10), nrow = 2, ncol = 3, byrow = TRUE)  # signal support

Because \(x_e\) is positively correlated with \(u\) and \(u\) also enters \(y\), OLS over-states the slope. Two-stage least squares (2SLS) and GME-IV both use \(z\) to purge the endogeneity.

ols_iv <- lm(y ~ xe)
b_2sls <- solve(crossprod(IV, X), crossprod(IV, y))   # just-identified 2SLS
gme_iv <- linreg_iv(y, X, IV, Ziv)

tab_iv <- data.frame(
  Truth    = c(`(Intercept)` = 1, xe = 1.5),
  OLS      = coef(ols_iv),
  `2SLS`   = as.vector(b_2sls),
  `GME-IV` = coef(gme_iv),
  check.names = FALSE
)
knitr::kable(tab_iv, digits = 3,
             caption = "OLS is biased by endogeneity; 2SLS and GME-IV recover the true slope.")

OLS pulls the slope above its true value of 1.5, while both instrumented estimators bring it back close to the truth. As the signal support widens the GME-IV estimate approaches the exact 2SLS solution; a narrower support regularizes it, which can help when the instruments are weak.

summary() reports a coefficient table with standard errors. Choose the method with se_method: "sandwich" (the default, robust), "delta" (classical), or "bootstrap". The robust sandwich matches a Monte-Carlo sampling SD and reduces to the 2SLS robust standard error as the support widens.

summary(gme_iv)
#> 
#> Stochastic-moments GME-IV (Golan 2008, pp. 89-91)
#> Call: linreg_iv(y = y, X = X, IV = IV, Z = Ziv)
#> 
#>   N=200  K=2  instruments=2 (just-identified)  M=3  nu=0.5
#> 
#> Residuals (y - X beta):
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> -2.219000 -0.709000 -0.005686  0.004531  0.698500  2.524000 
#> 
#> Coefficients (beta = Z p;  SE method: sandwich):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.89749    0.06923   12.96   <2e-16 ***
#> xe           1.59323    0.08526   18.69   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> normalized signal entropy S = 0.9885   max|moment resid| = 0.0001844
#> convergence: yes (0)
#> 
#> Note: asymptotic standard errors. The sampling distribution of beta is
#> support-bounded and can be skewed, so Wald intervals are approximate with
#> weak instruments or small n.

These are asymptotic standard errors: because the coefficient supports bound \(\hat{\boldsymbol{\beta}}\), its sampling distribution can be skewed, so Wald intervals are approximate with weak instruments or small \(n\).

8 Practical guidance

• Supports are the key input. Center the coefficient support on plausible values and make it wide enough to contain the truth; a wider support behaves more like OLS, a narrower one regularizes more strongly. For the symmetric noise support v, use a half-width of \(\max(3, \sqrt{2\log T})\,\mathrm{sd}(y)\) rather than a fixed three sigma — it must bound all \(T\) errors, not one.
• Check foc_residual. It is near zero for a healthy fit; a large value (or a warning about an unbounded dual) means the supports cannot represent the data, and the estimates are not trustworthy.
• Reach for GCE when you have real prior information. Encode it through p0; a uniform p0 is exactly GME.
• Leave nu = 0.5 unless you have a reason not to. It is the standard equal-weight GME.
• GME shines under collinearity and other ill-conditioned designs, where it trades a little bias for a large reduction in variance and instability.
• Report the entropy diagnostics. summary() gives the entropy-ratio test and the normalized entropy \(S\); the pseudo-\(R^2 = 1 - S\) and the Wald test (Section 6) complete the inferential toolkit.
• Endogeneity? Reach for linreg_iv(). Supply an instrument matrix (each exogenous regressor instruments itself); it accepts just- and over-identified systems and shares linreg()’s supports and priors.

See ?linreg and ?linreg_iv for the full argument lists, summary() for the entropy-ratio coefficient test, and the Pure Inverse Problems vignette for the underlying maximum-entropy machinery.