---
title: "Instrumental variables in fixed-effects Poisson models"
date: "2024-08-29"
---
```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(echo = TRUE)
library(Statamarkdown)
stataexe <- find_stata()
knitr::opts_chunk$set(engine.path = list(stata = stataexe))
```
## Why a control function is needed
Suppose a panel count-data model has the conditional mean
$$E[y_{it} \mid x_{it}, c_i] = c_i \exp(x_{it} \beta)$$
If a variable in $x_{it}$ is endogenous, linear two-stage least squares does not fit this exponential conditional mean. Nonlinear IV commands such as Stata's `ivpoisson` also cannot absorb high-dimensional individual and time fixed effects.
Lin and Wooldridge (2019) propose a two-step control-function estimator for this setting.
## The two-step estimator
Let $w_{it}$ be a continuous endogenous regressor and $z_{it}$ an excluded instrument. The first stage regresses $w_{it}$ on the instrument and exogenous controls while absorbing all fixed effects:
$$w_{it} = z_{it}\pi + x_{it1}\gamma + c_i + f_t + v_{it}$$
We save the OLS residual $\hat v_{it}$ and include it in a fixed-effects Poisson model:
$$E[y_{it} \mid w_{it}, x_{it1}, c_i, f_t, \hat{v}_{it}] = c_i f_t \exp(w_{it}\beta_1 + x_{it1}\beta_2 + \rho \hat{v}_{it})$$
In Stata, the two commands are `reghdfe ..., resid` and `ppmlhdfe`. In R, we can use `fixest::feols()` and `fixest::fepois()`. Testing $H_0:\rho=0$ gives a test of exogeneity. Because $\hat v_{it}$ is estimated, the usual second-stage standard errors are not valid. I use a cluster bootstrap that repeats both stages.
## A continuous endogenous variable in Stata
I use Stata's `website.dta`, which has 500 observations. `visits` is the count outcome, `time` is the endogenous regressor, `phone` is the instrument, and `frfam` is an exogenous control. The model absorbs fixed effects for the 13 advertising campaigns and for `female`.
```{stata}
*| label: stata1
*| echo: true
*| collectcode: true
webuse website, clear
* Step 1: Linear first stage absorbing both fixed effects
reghdfe time phone frfam, absorb(ad female) resid
predict double u2h_fe, resid
* Step 2: FE Poisson with control function residual
ppmlhdfe visits time u2h_fe frfam, absorb(ad female)
```
The coefficient on `time` is $0.0448$ ($p = 0.318$). The residual coefficient is $\hat{\rho} = 0.0765$ ($p = 0.096$), indicating weak positive selection.
### A comparison with linear models
In a linear model, the control-function estimator gives the same coefficient as 2SLS. I fit both as a check on the two-step code:
```{stata}
*| label: stata2
*| echo: true
*| collectcode: true
clear all
webuse website, clear
* 1. Linear 2SLS
ivreghdfe visits frfam (time=phone), absorb(ad female)
* 2. Linear control function
reghdfe time phone frfam, absorb(ad female) resid
predict double u2h_fe, resid
reghdfe visits time u2h_fe frfam, absorb(ad female)
* 3. Naive models (ignoring endogeneity)
reghdfe visits time frfam, absorb(ad female)
ppmlhdfe visits time frfam, absorb(ad female)
```
| Model | Estimator | `time` Coef. | Std. Err. | $p$-value |
|:---|:---|:---:|:---:|:---:|
| Naive Linear FE | `reghdfe` | 0.7402 | 0.0401 | $<0.001$ |
| Linear 2SLS | `ivreghdfe` | 0.5643 | 0.2124 | 0.008 |
| Linear Control Function | Manual two-step | 0.5643 | 0.2083 | 0.007 |
| Naive FE Poisson | `ppmlhdfe` (no residual) | 0.1188 | 0.0096 | $<0.001$ |
| Control-Function FE Poisson | `ppmlhdfe` (with $\hat{v}$) | 0.0448 | 0.0449 | 0.318 |
The manual linear control function reproduces the 2SLS coefficient exactly. In the Poisson model, adding the control-function residual reduces the `time` coefficient from $0.1188$ to $0.0448$.
### Cluster bootstrap in Stata
The bootstrap can draw the same cluster more than once. Each copy must get a new identifier through `idcluster()`; otherwise Stata pools the copies into one fixed effect.
```{stata}
*| label: stata3
*| echo: true
*| collectcode: true
clear all
capture program drop ppmlhdfe_cf
program ppmlhdfe_cf, rclass
reghdfe time phone frfam, absorb(newid female) resid
predict double u2h_fe, resid
ppmlhdfe visits time u2h_fe frfam, absorb(newid female)
return scalar b_time = _b[time]
return scalar b_u2h = _b[u2h_fe]
drop u2h_fe
xtset, clear
end
webuse website, clear
xtset, clear
bootstrap r(b_time) r(b_u2h), reps(1000) seed(123) cluster(ad) idcluster(newid) nodots: ppmlhdfe_cf
```
The bootstrap standard error on `time` is $0.0477$ (compared to the unadjusted analytic standard error of $0.0449$).
## The same estimator in R
In R, I use `fixest::feols()` and `fixest::fepois()`. `feols()` drops a singleton cluster here (`ad == 10`), so we use `obs()` to align the residuals with the rows kept in the first stage:
```{r}
#| label: r-demo
#| echo: true
#| message: false
#| warning: false
library(fixest)
library(haven)
df <- as.data.frame(read_dta("https://www.stata-press.com/data/r18/website.dta"))
df$ad <- as.factor(df$ad)
df$female <- as.factor(df$female)
# Step 1: Linear first stage
fs <- feols(time ~ phone + frfam | ad + female, data = df)
# Align rows (feols drops singletons)
d <- df[obs(fs), ]
d$u2h_fe <- residuals(fs)
# Step 2: FE Poisson with and without control function
m_naive <- fepois(visits ~ time + frfam | ad + female, data = d, notes = FALSE)
m_cf <- fepois(visits ~ time + u2h_fe + frfam | ad + female, data = d, notes = FALSE)
etable(m_naive, m_cf, headers = c("Naive FE Poisson", "Control Function"), se.below = TRUE)
```
The point estimates match Stata to seven digits: $0.1187612$ (naive) and $0.0448103$ (control function).
### Cluster bootstrap in R
```{r}
#| label: r-bootstrap
#| echo: true
#| message: false
#| warning: false
set.seed(123)
B <- 1000
clusters <- unique(d$ad)
boot_coefs <- matrix(NA_real_, nrow = B, ncol = 2, dimnames = list(NULL, c("time", "resid")))
for (b in seq_len(B)) {
drawn <- sample(clusters, length(clusters), replace = TRUE)
boot_data <- do.call(rbind, lapply(seq_along(drawn), function(k) {
sub <- d[d$ad == drawn[k], ]
sub$boot_id <- k
sub
}))
boot_data$boot_id <- as.factor(boot_data$boot_id)
fs_b <- feols(time ~ phone + frfam | boot_id + female, data = boot_data, notes = FALSE)
bb <- boot_data[obs(fs_b), ]
bb$u_b <- residuals(fs_b)
m_b <- try(fepois(visits ~ time + u_b + frfam | boot_id + female, data = bb, notes = FALSE), silent = TRUE)
if (!inherits(m_b, "try-error")) boot_coefs[b, ] <- coef(m_b)[c("time", "u_b")]
}
cat("Completed replications :", sum(complete.cases(boot_coefs)), "\n")
cat("Bootstrap SE, time :", round(sd(boot_coefs[, "time"], na.rm = TRUE), 4), "\n")
cat("Bootstrap SE, residual :", round(sd(boot_coefs[, "resid"], na.rm = TRUE), 4), "\n")
```
The R bootstrap yields $\text{SE} = 0.0450$. With only 12 clusters in `ad`, small differences across software packages ($0.0477$ in Stata vs. $0.0450$ in R) reflect finite-sample sensitivity in how rare degenerate resamples and singletons are handled.
## A binary endogenous regressor
When $D_{it}$ is binary, as with program participation, a linear first stage is misspecified because its error cannot satisfy the required conditional-independence restriction.
### A CRE logit first stage
A first-stage logit with a dummy for every unit has an incidental-parameters problem when $T$ is small. For example, its slopes can be about twice as large as they should be when $T=2$. I instead use pooled logit with a Mundlak adjustment, adding the unit means of the time-varying covariates:
$$P(D_{it} = 1 \mid z_{it}, \bar{z}_i) = \Lambda(z_{it}\gamma + \bar{z}_i \pi)$$
For logit, the generalized residual is simply the response residual,
$$\hat r_{it}=D_{it}-\hat p_{it}.$$
We then include $\hat r_{it}$ in the fixed-effects Poisson model.
### A simulation
Let's simulate 2,000 units with unbalanced panels, where $T_i$ ranges from 2 to 6. The treatment $D$ is binary and endogenous, and its true coefficient is $\beta=0.30$.
```{stata}
*| label: binary_cf
*| echo: true
*| collectcode: true
preserve
clear
set seed 20260731
set obs 2000
gen i = _n
gen double ci = rnormal()*0.5
gen T = 2 + int(runiform()*5) /* 2-6 observations per unit */
expand T
bysort i: gen t = _n
gen double z = rnormal() /* instrument */
gen double e = rnormal()
gen double u = 0.8*e + rnormal()*0.6 /* correlated with e -> endogeneity */
gen byte d = (-0.2 + 0.9*z + ci + e > 0)
gen double mu = exp(0.3*d + ci + 0.5*u)
gen y = rpoisson(mu) /* true beta = 0.3 */
* 1. Naive FE Poisson (ignores endogeneity)
quietly ppmlhdfe y d, absorb(i)
scalar b_naive = _b[d]
* 2. CF with CRE (Mundlak) logit first stage
bysort i: egen double z_i = mean(z)
quietly logit d z z_i
predict double phat, pr
gen double gr = d - phat
quietly ppmlhdfe y d gr, absorb(i)
scalar b_cre = _b[d]
* 3. CF with linear first stage
quietly reghdfe d z, absorb(i) resid
predict double v_lin, resid
quietly ppmlhdfe y d v_lin, absorb(i)
scalar b_lin = _b[d]
di "True coefficient on d = 0.300000"
di "Naive FE Poisson = " %8.6f b_naive
di "CF, CRE (Mundlak) Logit = " %8.6f b_cre
di "CF, Linear First Stage = " %8.6f b_lin
restore
```
The naive fixed-effects Poisson estimate is $0.747$, more than 2.4 times the true value. The CRE logit control function gives $0.294$, while the linear first-stage control function gives $0.278$.
In an application, I would report both control-function specifications as a sensitivity check.
## What to remember
The control-function approach adds the first-stage residual to the exponential outcome model. For a continuous endogenous variable, we can combine `reghdfe` with `ppmlhdfe` in Stata, or `feols()` with `fepois()` in R. Inference should bootstrap both stages by cluster. For a binary endogenous variable, a CRE logit first stage avoids the incidental-parameters problem from unit dummies, and its residual is $d-\hat p$.