# Adjusting Standard Error for Imputed/Generated Regressors

This is my first question, so I hope this is a valid question. I am surprised that I have seen only a few questions (and no answer helping me) referring to the adjustment of variance estimators in the case of imputed regressors because this is a common issue for applied statisticians in different fields.

Regarding my problem: I am running a simple linear regression model using hierarchical data, including several imputed (or generated) predictors, which I estimated in a first stage, and I want to take into account the additional uncertainty by including those imputed regressors. The second stage is estimated on a subsample of the population used in the first stage.

What is the easiest way to adjust the standard errors in the second regression?

It is clear how to do this, for example, for two-stage least squares or other methods where only one estimated parameter is included. But method descriptions soon become pretty technical (at least for me) for more complex cases

Should I do parametric bootstrapping? Other methods I have seen are, for example, the Murphy-Topel variance estimator (see this related question here at CrossValidated or two articles in the Stata Journal here and here) or influence functions (e.g. Hahn & Rider, 2013). Sometimes authors also refer to the delta method.

Perhaps some more details to my question: My estimation in the first step takes a while because my sample size is quite huge, and I use there at least non-linear models. I am calculating a two-step estimation using a complicated hierarchical longitudinal dataset with individual and firm information. In the first step, I estimate different parameters using different methods. Subsequently, I include those estimated parameters as predictors in the second regression (for example, random effects from a mixed logit model and higher-level inequality measures). For these predictors, I extract coefficients and standard errors available.

Unfortunately, due to data protection regulations, I cannot show my data here. But I demonstrate this issue using the data for the analysis of incumbency advantage in the book "Bayesian Data Analysis" by Gelman et al. Using this data, I generate two predictors, which I include later in my final regression, where I look whether the incumbent won an election. The first predictors are random effects of how republican leaning different districts are using a mixed model, and the second regression is a linear regression on the republican vote share in a given election. Regressions here do not make much sense, but this should just be an example of using hierarchical data everyone can access.

How would you adjust the standard error to account for the two generated regressors? Would you advise doing parametric bootstrapping (assuming that the generated regressors are independent)?

How to account additionally for possible intra-state correlation such as clustered-standard error?

Furthermore, should I worry that both stages are used on the same data (or, in my case, I use a subsample of the same data in the second stage)?

> rm(list=ls(all=TRUE))
> library(lme4)
> library(arm)
> #
> dat$$year<-1964 > for(i in seq(1966,1992,2)){ + eval(parse(text=paste("tmp<-read.table('http://www.stat.columbia.edu/~gelman/book/data/incumbency/",i,".asc',header=F)",sep=""))) + eval(parse(text=paste("tmp$$year<-",i,"",sep="")))
+ dat<-rbind(dat,tmp)
+ }
> #
1     1        1          1           141310            60654 1964
2     1        2          0           119530            69403 1964
3     1        3          1           126353            71393 1964
4     1        4         -1           117220           109027 1964
5     1        5          1           133072            64651 1964
6     1        6          1           115498            81105 1964

> # Some Data Preparation
> dat$$share_republican<-dat$$votes_republican/(dat$$votes_republican+dat$$votes_democratic)
> dat$$share_democratic<-dat$$votes_democratic/(dat$$votes_republican+dat$$votes_democratic)
> dat$$incumbency2<-factor(NA,levels=c("democratic","republican","open")) > dat$$incumbency2[which(dat$$incumbency==(-1))]<-"republican" > dat$$incumbency2[which(dat$$incumbency==1)]<-"democratic" > dat$$incumbency2[which(dat$$incumbency==0)]<-"open" > # Who won the election? > dat$$winner<-factor(NA,levels=c("democratic","republican"))
> dat$$winner[which(dat$$share_republican>dat$$share_democratic)]<-"republican" > dat$$winner[which(dat$$share_republicanshare_democratic)]<-"democratic"
> # Did the incumbent won the election?
> dat$$incumbent_winner<-NA > dat$$incumbent_winner[which(as.character(dat$$incumbency2)==as.character(dat$$winner))]<-1
> dat$$incumbent_winner[which(as.character(dat$$incumbency2)!=as.character(dat$$winner) & dat$$incumbency2!="open")]<-0

> # First-Stage (How republican leaning are US districts?)
> first_stage1<-lmer(share_republican~as.factor(year)+(1|district),data=subset(dat,share_democratic>0 & share_republican>0))
> first_stage1_results<-data.frame(rownames(ranef(first_stage1)$$district),ranef(first_stage1)$$district,se.ranef(first_stage1)\$district)
> colnames(first_stage1_results)<-c("district","estimate_republican_leaning","estimate_republican_leaning_set")

> # Second-Stage (How well performed republicans overall in a given election)?
> first_stage2<-lm(I(share_republican)~as.factor(year)-1,data=dat)
> first_stage2_results<-data.frame(
+ as.numeric(gsub("as.factor(year)","",rownames(summary(first_stage2)$$coefficients)[which(grepl("year",rownames(summary(first_stage2)$$coefficients)))],fixed=TRUE)),
+ summary(first_stage2)$$coefficients[which(grepl("year",names(coef(first_stage2)))),1], + summary(first_stage2)$$coefficients[which(grepl("year",names(coef(first_stage2)))),2]
+ )
> colnames(first_stage2_results)<-c("year","estimate_year","estimate_year_standard_error")

> # Merge Data
> dat<-merge(dat,first_stage1_results,"district",all.x=T)
> dat<-merge(dat,first_stage2_results,"year",all.x=T)

> # Calculated final regression
> summary(glm(incumbent_winner~as.factor(state)+winner+estimate_year+estimate_republican_leaning,data=dat))

Call:
glm(formula = incumbent_winner ~ as.factor(state) + winner +
estimate_year + estimate_republican_leaning, data = dat)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-0.99181   0.02736   0.04733   0.06850   0.21052

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                  0.656765   0.062446  10.517  < 2e-16 ***
as.factor(state)2            0.052053   0.053561   0.972  0.33117
...
as.factor(state)82           0.082348   0.057742   1.426  0.15389
winnerrepublican            -0.019976   0.006562  -3.044  0.00235 **
estimate_year                0.600357   0.131925   4.551 5.46e-06 ***
estimate_republican_leaning -0.093269   0.131099  -0.711  0.47684
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.05087851)

Null deviance: 289.34  on 5620  degrees of freedom
Residual deviance: 283.29  on 5568  degrees of freedom
(691 observations deleted due to missingness)
AIC: -734.65

Number of Fisher Scoring iterations: 2

• This remains an excellent question! What did you end up doing? Please update!
– ABC
Commented Dec 7, 2022 at 16:45
• Thank you @ABC, in the end, I used only one imputed predictor and did parametric bootstrap! But would be still interested in an answer! Commented Dec 8, 2022 at 7:26
• @Are, do you know if the Murphy-Topel variance estimator is applicable only to settings where the first stage is linear? It seems to be like it does, because it looks like it requires a 1st stage variance/covariance matrix. You wouldn't get that if your 1st stage was (say) a machine learning model.
– ABC
Commented Dec 12, 2022 at 16:29