3.1 Linear specification: confidence intervals
hdcce_inference returns a desparsified estimate for each
coefficient in COEF_INDEX_VEC. The debiasing runs a
nodewise lasso of the target regressor on the others and corrects the
lasso estimate along the resulting residual direction, which restores
asymptotic normality.
dat0 <- list(x = data_inference$x, y = data_inference$y[, 1]) # beta_1 = 0
inf0 <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T,
COEF_INDEX_VEC = 1, NFOLDS = 5)
r <- inf0$results[["1"]]
c(estimate = r$coef_despar, se = r$se, p_value = r$p_value)
#> estimate se p_value
#> -0.01843391 0.04668838 0.69296943
r$confidence_band
#> conf_band_min conf_band_max
#> alpha=0.01 -0.13869521 0.10182740
#> alpha=0.05 -0.10994146 0.07307364
#> alpha=0.1 -0.09522946 0.05836165
Rows of confidence_band follow the order of
alpha, which defaults to c(0.01, 0.05, 0.10).
The interval covers zero, as it should.
With a non-zero coefficient the interval moves away from zero:
dat2 <- list(x = data_inference$x, y = data_inference$y[, 3]) # beta_1 = 0.2
inf2 <- hdcce_inference(dat2, obs_N = obs_N, obs_T = obs_T,
COEF_INDEX_VEC = 1, NFOLDS = 5)
inf2$results[["1"]]$confidence_band
#> conf_band_min conf_band_max
#> alpha=0.01 0.04034862 0.2810507
#> alpha=0.05 0.06912383 0.2522755
#> alpha=0.1 0.08384680 0.2375525
The reported \(p\)-value tests \(H_0 \colon \beta_j = 0\) and is exactly
dual to the band: \(p \le \alpha\) if
and only if zero lies outside the \((1-\alpha)\) interval.
cb <- r$confidence_band
cbind(alpha = inf0$alpha,
p_le_alpha = r$p_value <= inf0$alpha,
zero_excluded = !(cb[, 1] <= 0 & 0 <= cb[, 2]))
#> alpha p_le_alpha zero_excluded
#> alpha=0.01 0.01 0 0
#> alpha=0.05 0.05 0 0
#> alpha=0.1 0.10 0 0
Several coefficients can be treated in one call; the result carries
one entry per index, named by it. Coefficients 1 and 5 are truly zero
here, coefficient 2 is truly one.
inf_multi <- hdcce_inference(dat0, obs_N = obs_N, obs_T = obs_T,
COEF_INDEX_VEC = c(1, 2, 5), NFOLDS = 5)
t(sapply(inf_multi$results,
function(z) c(estimate = z$coef_despar, se = z$se, p = z$p_value)))
#> estimate se p
#> 1 -0.01843391 0.04668838 6.929694e-01
#> 2 0.93925548 0.06315704 5.028456e-50
#> 5 -0.02748294 0.05404289 6.110751e-01
The HAC argument selects the variance estimator:
1 assumes homoscedastic errors without serial correlation,
2 (the default) allows heteroscedasticity, and
3 allows serial correlation as well. The choice matters
when \(T\) is small, since the
unit-wise variance estimates underlying HAC = 2 rest on few
degrees of freedom.
sapply(1:3, function(k)
hdcce_inference(dat2, obs_N, obs_T, COEF_INDEX_VEC = 1,
NFOLDS = 5, HAC = k)$results[["1"]]$se)
#> [1] 0.04718331 0.04672322 0.03764749
3.2 Dictionary specification: significance test
With a dictionary the null of interest is \(H_0 \colon m_j = 0\), a functional null.
There is no scalar to invert, so the output is a test rather than a
confidence interval: se and confidence_band
are NULL, and a statistic,
critical_values, p_value and
profile take their place.
The statistic is a maximum, over a grid of locations \(w\), of self-normalised correlations
between the residuals under the null and kernel weights centred at \(w\). Its critical values come from a
Gaussian coupling simulated B times.
A covariate with no effect. In the first response
column \(\beta_1 = 0\), so \(m_1 \equiv 0\) and the null is true.
tst1 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]),
obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 1,
dictionaries = dict_inf, NFOLDS = 5, B = 500)
q1 <- tst1$results[["1"]]
c(statistic = q1$statistic, p_value = q1$p_value)
#> statistic p_value
#> 1.7278021 0.3193613
q1$critical_values
#> alpha=0.01 alpha=0.05 alpha=0.1
#> 3.097412 2.528749 2.281566
A covariate with an effect. The second covariate has
coefficient one in the same response, so the null is false.
tst2 <- hdcce_inference(list(x = data_inference$x, y = data_inference$y[, 1]),
obs_N = obs_N, obs_T = obs_T, COEF_INDEX_VEC = 2,
dictionaries = dict_inf, NFOLDS = 5, B = 500)
q2 <- tst2$results[["2"]]
c(statistic = q2$statistic, p_value = q2$p_value)
#> statistic p_value
#> 7.440765494 0.001996008
q2$critical_values
#> alpha=0.01 alpha=0.05 alpha=0.1
#> 3.022649 2.506008 2.185410
Because the statistic is a maximum it is informative to see where the
evidence sits. The profile records \(\Psi_{w,h}\) across the grid.
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2.4, 1))
for (z in list(list(q = q1, main = "covariate 1 (no effect)"),
list(q = q2, main = "covariate 2 (effect)"))) {
pr <- z$q$profile
plot(pr$w, abs(pr$Psi_w), type = "b", pch = 16, main = z$main,
xlab = "location w", ylab = expression(group("|", Psi[list(w, h)], "|")),
ylim = range(0, abs(pr$Psi_w), z$q$critical_values))
abline(h = z$q$critical_values[2], col = "red", lty = 2)
}
For the first covariate the curve stays well below the 5% critical
value at every location; for the second it crosses it, and the crossing
identifies the region of the nodewise residuals in which the deviation
is detected.
The profile also reports n_eff, the number
of nodewise residuals falling in each bump \([w-h, w+h]\). This is the diagnostic to
check before trusting the test: when the smallest count is low the
Gaussian approximation is unreliable, and the bandwidth should be
increased or the region narrowed.
q2$profile
#> w Psi_w n_eff
#> 1 -0.8445217 -6.596739 72
#> 2 -0.2411041 -4.722727 113
#> 3 0.3623135 4.690826 123
#> 4 0.9657312 7.440765 56
Note that the projection is rebuilt from \(\bar{\boldsymbol{X}}_{(-j)}\) for each
tested index, so the grid, \(\widehat{K}\) and the residuals differ
slightly between the two tests above.
