# 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 few questions (and no answer helping me out) referring to the adjustment of variance estimators in the case of imputed regressors because I guess 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?

This seems to be pretty 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 become soon pretty technical (at least for me) for more difficult 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 quiet 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 a 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 predictor 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 using hierarchical data everyone can access.

How would you adjust standard error here to account for the two generated regressors? Would you advise to do parametric bootstraping (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 that 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_republican<dat$share_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