Correcting (or bootstrapping) the standard errrors for a two stage glm

Cross posted on StackOverflow

I want to somehow correct the standard errors of my two stage residual inclusion, where in contrast to the 2SLS, the residuals are included in addition to the instrumented variable (see Terza et al. 2008). I have put some example code and example data below.

On how correct the standard errors (From Terza et al. 2008, note 8) :

"The 2SRI estimator can be cast as a special case of the conventional generic two-stage optimization estimator. Therefore, its asymptotic properties (in particular correct asymptotic standard errors) follow directly from the discussions found in Newey and McFadden (1994), White (1994, Chapter 6), or Wooldridge (2002, Chapter 12). Alternatively, the standard errors of the 2SRI estimator can be obtained via bootstrapping."

What I tried is the following (based on this link):

library(sure)
# for residual function and sample data sets
library(MASS)
# for polr function
rm(df1)
df1 <- df1
df1$$x2 <- df2$$y
df1$$x3 <- df2$$x
df1$$y <- df3$$x/10
fit.polr <- polr(x2 ~ x + x3, data = df1, method = "probit")

# BOOTSTRAPPING
library(boot)
glm_2siv_coef <- function(data,indices){
d <- data[indices,]
stage_1 <- polr(x2 ~ x + x3, Hess=TRUE, data=d)
stage_2 <- glm(y ~ x + x2 + x3 + resids(stage_1),
family = "quasibinomial", data=d)
return(summary(stage_2)$estimate["x2",1]) } boot.results <- boot(data=df1,statistic=glm_2siv_coef,R=1) boot.results  I am not sure whether this is the appropriate way to do it, because I am not really familiar with bootstrapping techniques. In addition, I get the error Error in boot.out$t[, index] : subscript out of bounds. Cross posting this on stack for that reason.

Could anyone perhaps show me how to apply this to the example data/code below?

NOTE: Residuals are calculated with surrogate residuals from the sure package (link).

Terza JV, Basu A, Rathouz PJ. (2008) Two‐stage residual inclusion estimation: addressing endogeneity in health econometric modeling. J Health Econ. 2008;27(3):531‐543.

Bootstrapping: Yes, if you are bootstrapping, you want to bootstrap the whole process, i.e. resample your dataset (with replacement), and then estimate first stage and second stage with that data. Intuitively, the idea is that you're estimating the first stage and using the residuals in the second stage, so you want to account for possible noise in both stages.

Code: Take a closer look at your function glm_2siv_coef. You return summary(stage_2)$estimate["x2",1] but this is not defined. I think you mean to look at  summary(stage_2)$coefficients which is a matrix, and then you want to extract the relevant estimate. Looking at the matrix, it's not clear what value you exactly want to extract, but it should be easy to get it by properly indexing that matrix. So fix the return() part of your function. That's why you get the error in boot: you're feeding it nulls through that function!

This answer is based on this link, answer by @jay.sf. Although it runs, I am not completely sure whether this would be correct. I would be very happy if someone could comment (or check).

set.seed(2)
sandbox$$Group <- as.factor(sandbox$$Group)
reduced.form <- polr(Group ~ z + random_variable + year, data=sandbox)
consistent.glm <- glm(y ~ Group + resids(reduced.form) + random_variable + year, family="quasibinomial", data=sandbox)

FUN <- function(x) {
reduced.form <-  polr(Group ~ z + random_variable + year, data=x)
fit <- glm(y ~ Group + resids(reduced.form) + random_variable + year,family="quasibinomial", data=x)
fit$coefficients } set.seed(42) R <- 200 bs <- t(replicate(R, FUN(sandbox[sample(nrow(sandbox), nrow(sandbox), replace=T), ]))) # To scrape out a summary, the matrixStats package is most convenient. library(matrixStats) b <- consistent.glm$coefficients
SE <- colSds(bs)
z <- b/SE
p <- 2 * pt(-abs(z), df = Inf)
ci <- colQuantiles(bs, probs=c(.025, .975))
res <- signif(cbind(b, SE, z, p, ci), 4)
res


DATA

a    <- 2    # structural parameter of interest
b    <- 1    # strength of instrument
rho  <- 0.5  # degree of endogeneity

N    <- 1000
z    <- rnorm(N)
res1 <- rnorm(N)
res2 <- res1*rho + sqrt(1-rho*rho)*rnorm(N)
x    <- z*b + res1
ys   <- x*a + res2
d    <- (ys>0) #dummy variable
y    <- (10-(d*ys))/10
random_variable <- rnorm(100, mean = 0, sd = 1)

library(data.table)
DT_1 <- data.frame(y,x,z, random_variable)
DT_2 <- structure(list(ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50), year = c(1995, 1995, 1995, 1995, 1995,
1995, 1995, 1995, 1995, 1995, 2000, 2000, 2000, 2000, 2000, 2000,
2000, 2000, 2000, 2000, 2005, 2005, 2005, 2005, 2005, 2005, 2005,
2005, 2005, 2005, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,
2010, 2010, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015,
2015), Group = c("A", "A", "A", "A", "B", "B", "B", "B", "C",
"C", "A", "A", "A", "A", "B", "B", "B", "B", "C", "C", "A", "A",
"A", "A", "B", "B", "B", "B", "C", "C", "A", "A", "A", "A", "B",
"B", "B", "B", "C", "C", "A", "A", "A", "A", "B", "B", "B", "B",
"C", "C"), event = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), win_or_lose = c(-1,
-1, -1, -1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, 1, 1, 1, 1, 0, 0,
-1, -1, -1, -1, 1, 1, 1, 1, 0, 0)), row.names = c(NA, -50L), class = c("tbl_df",
"tbl", "data.frame"))
DT_1 <- setDT(DT_1)
DT_2 <- setDT(DT_2)
DT_2 <- rbind(DT_2 , DT_2 [rep(1:50, 19), ])
sandbox <- cbind(DT_1, DT_2)
sandbox <- setDT(sandbox)[y<0, y:=0]