0
$\begingroup$

I want to translate the "score bootstrap" that Kline and Santos (2011) propose into R code. The algorithm is described on page 6 as follows.

STEP 1: Obtain the full sample OLS estimate $\hat\beta$ and employ it to generate the fitted score contributions $\{X_i(Y_i-X^{'}_i\hat\beta)\}^n_{i=1}$.

STEP 2: Using an i.i.d. sample of random weights $\{W_i\}_{i=1}^n$ independent of $\{Y_i,X_i\}_{i=1}^n$ and satisfying $E[W_i]=0$ and $E[W^2_i]=1$, construct a new set of perturbed score contributions $\{X_i(Y_i-X^{'}_i\hat\beta)W_i\}_{i=1}^n$.

STEP 3: Multiply the constructed perturbed score by the inverse Hessian to obtain $H_n^{-1}n^{-\frac12}\sum_i(Y_i-X_i\hat\beta)X_iW_i$ and use its distribution conditional on $\{Y_i,X_i\}_ {i=1}^n$ as an estimate of the distribution of $\sqrt{n}(\hat\beta-\beta_0)$.

STEPS 1 and 2 seem to be straightforward.

set.seed(42)
n <- 100
x <- rnorm(n)
y <- x + rnorm(n)

fit <- lm(y ~ x)

X <- model.matrix(fit)
yhat <- fit$fitted.values  # Xb

S <- X*(y - yhat)  # score matrix
H <- - solve(vcov(fit))  # hessian
w <- ifelse(runif(n) < .5, 1, -1)  # draws from Rademacher distr.

However, when I attempt to calculate STEP 3, I'm getting an error, because the number of columns of $H^{-1}$ won't match with the number of rows of $S$.

# sqrt(n) * solve(H) %*% S * w  # STEP 3
## Error in solve(H) %*% S : non-conformable arguments

Does anyone see if I'm maybe misinterpreting the math somewhere?

Edit

I obviously made an error with the multiplication in STEP 3, where I need to take account of the $\sum$, i.e. colSums of the perturbed scores. Thus, after STEP 3 I get this result:

solve(H) * sqrt(n) * colSums(S * w)  # STEP 3
#             (Intercept)           x
# (Intercept)  1.15362358 -0.03490467
# x            0.04228051 -1.30034588

To verify the result I set $w=1$, i.e. I don't apply any weights, and I get essentially zero.

w <- 1
solve(H) * sqrt(n) * colSums(S * w)
#               (Intercept)            x
# (Intercept) -4.149004e-16 1.255345e-17
# x           -3.034329e-19 9.332144e-18

This feels right, because the sum of the scores is (by definition) zero, too:

colSums(S)
#  (Intercept)             x 
# 5.023759e-15 -1.214306e-16

Where I'm stuck

However I'm not sure what I'm actually getting here. According to the algorithm I have an "estimate of the distribution of $\sqrt{n}(\hat\beta-\beta_0)$". I know this has to do with maximum likelihood, but how to derive my variance estimator from this? If the result is already my variance matrix, then taking the square root of the diagonal should give me standard errors similar to those of the model summary, but they're pretty off:

sqrt(diag(solve(H) * colSums(t(S) %*% S) * solve(H)))
# (Intercept)           x 
#  0.06225221  0.07737588

summary(fit)$coe[, 2]
# (Intercept)           x 
#  0.09087774  0.08766507

Is this the wrong interpretation of the estimate of $\sqrt{n}(\hat\beta-\beta_0)$? Did I make an error in interpreting the equations? What I finally need is the test statistic of one bootstrap iteration, and I am not sure how to get it out of the equation of STEP 3.

Note: I'm asking for statistical help (distinct from programming help) to understand how to derive a variance estimator from the estimate of $\sqrt{n}(\hat\beta-\beta_0)$ of STEP 3.

For sake of completeness below the formula for the test statistic of the score bootstrap given in the paper on page 7.

$$T_n^{*score}:=(H_n^{-1}\sum_{n}^*(\hat\beta)H^{-1}_n)^{-\frac12}H^{-1}_n\frac1{\sqrt{n}}\sum_{i=1}^nX_i\epsilon_i^*,\quad where \quad\sum_n^*(\beta) := \frac1n\sum_iX_iX_i^T(Y_i^*-X_i^T\beta)^2$$

|cite|improve this question|||||
$\endgroup$
  • $\begingroup$ Looks like $X$ should be transposed so that matrix columns and rows are conformable. Doing this, the code runs smoothly without errors. $\endgroup$ – compbiostats Oct 16 '19 at 14:31
  • 1
    $\begingroup$ Replace X*(y-yhat) by apply(X, 1, `*`, y-yhat) $\endgroup$ – AdamO Oct 16 '19 at 14:44
  • 1
    $\begingroup$ @AdamO When I run this line in conjunction with the rest of the OP's code, I still get a conformability error... $\endgroup$ – compbiostats Oct 16 '19 at 15:03
  • $\begingroup$ @compbiostats you need to use dim to debug your code. It's a good R skill. Your implementation isn't wrong from a stat side, you just need to sooth out your code. So probably off topic here. $\endgroup$ – AdamO Oct 16 '19 at 15:14
  • $\begingroup$ @AdamO Thanks. Yes, I see that on comparing S[1,] with y - yhat. $\endgroup$ – compbiostats Oct 16 '19 at 15:17
0
$\begingroup$

$X$ (the model matrix) is a $n \times 2$ matrix. $y-\hat{y}$ is a $n \times 1$ vector.

The * operator in R is not generic, meaning it can't figure out $X$ is a matrix for you, so it just recycles the multiplication over a $2n$ vector of the entries of $X$ by column. Not the desired output.

Research %*% and apply and study the output carefully.

|cite|improve this answer|||||
$\endgroup$
  • $\begingroup$ This is not right, X*(y - yhat) yields exactly what it should yield. Compare the output of sandwich::estfun(fit) which is designed to yield the scores, see @AchimZeileis's answer on Stack Overflow. $\endgroup$ – jay.sf Oct 16 '19 at 17:14
  • 1
    $\begingroup$ @jay.sf I think you're right. an earlier comment suggested transposing X and naively I had that in my workspace. However, my comment on good coding practice is still well, warranted: don't trust R's recycling for matrix operations. Use %*%, sweep, and apply. $\endgroup$ – AdamO Oct 16 '19 at 18:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.