Not an answer, but some doubts:
What drives the FWL theorem is that regressors and residuals are exactly orthogonal when estimating by OLS. Since that need not be the case for other estimation techniques (but see the discussion in the comments below related to this particular problem!), I am not sure if, and if so how, such a result could carry over.
What is maybe more, standard FWL involves three (sets of) OLS regressions:
- Dependent variable on regressors to be partialled out
- Other regressors on regressors to be partialled out
- Residuals of 1. on residuals of 2.
How would we carry over this logic to Poisson regression? We can sure run a Poisson regression like in 1. and collect the corresponding residuals.
However, in 2., it will typically not be the case that the other regressors also are count variables, so that a Poisson regression would not make sense here (and even not work if these regressors contain negative values). Similarly, the residuals of 1. will also contain negative values, such that a Poisson regression in 3. would again not work.
So, "nonlinear FWL" cannot be just replacing OLS regressions by Poisson regressions. We can instead run OLS regressions in 2. and 3., but then there is, as far as I can see, no reason to expect any type of equivalence.
That said, this code suggests a certain closeness, which may however also be due to all true coefficients being zero:
n <- 3000000
x <- rnorm(n)
z <- rnorm(n)
y <- rpois(n, lambda=1)
poisson.full <- glm(y ~ x+z, family = 'poisson')
poisson.short <- glm(y ~ z, family = 'poisson')
others <- lm(x ~ z)
summary(lm(resid(poisson.short) ~ resid(others)-1))
summary(poisson.full)