`residuals_from` argument for estfun_teeMod()

[benthestatistician]

We need a mechanism for creating estfun's reflecting the shape of one teeMod but the residuals of another. We can make this possible by modifying estfun_teeMod to have a new argument residuals_from=x, expecting an lm (or glm) from which working residuals and working weights can be extracted. Those working residuals and working weights are then the ones used to create the estfun value.

[jwasserman2]

@benthestatistician I went ahead and made residuals_from a formal arg so users could use the abbreviated residuals. I also did not implement the expectation that it's an lm or glm given your statement offline about possibly wanting to accommodate residuals from random effects fits, which don't inherit from lm or glm. Please review the PR when you have a chance.

[benthestatistician]

1. Nice work!
2. In estfun.teeMod() and bread.teeMod() docstubs, let's fill out the sentences describing how the model matrix from x is combined with residuals and weights from residuals_from to also carry a statement of the bigger picture: with this change it reports the estimating functions or inverse A matrix associated with the data, with replacement of unknown parameters by plug-in of estimates according to residuals_from.
3. (Note to self: circle back and review changes relating to bias corrections after I finish reading Josh's CR2 chapter.)

[jwasserman2]

I hadn't considered the changes needed for bread to incorporate the residuals_from argument for glm teeMods. After a little investigating, it doesn't seem like this can be implemented particularly cleanly. bread.glm(x) returns a scaled version of summary(x)$cov.unscaled, which is chol2inv(x$qr). I'm not sure how we'd get use residuals from residuals_from unless we write the calculations out ourselves, which would require some renovations to our bread code and also incur a hit in computation time. Thoughts @benthestatistician ?

[benthestatistician]

Interesting problem! I'm optimistic that this won't stand in the way of our ending with a clean solution, but I agree that it complicates the path to that solution: it has me thinking we should go back to values_from for the name of our new argument, not residuals_from, as I'm seeing that we'll have to pull artifacts other than the residuals it calls for some GLM / exponential family details I'm not totally on top of.

As our $A_{22}$ we want information as evaluated at the parameters of the values_from model, not the parameters of x. But we want it to be in the shape of x, not that of the values_from model (which may have fewer regressors). To get this, I think we can:

1. generate predictions from the values_from model
2. use these predictions as starting point for a single-step re-fitting of glm x
3. extract information from resulting model.

In step (2), the one-step condition is enforced by glm argument control=list(maxit=1), while the predictions are passed through argument mustart. For step (3), we can just mimic what bread.glm() does.¹ It may be that we need to write something our own to get a legit information matrix (expected or observed) in this context. To do: investigate in generalized score testing sources.

Footnotes

1. Abstractly, as $\hat{A}_{22}$ matrix we might like an inverted Hessian of the glm likelihood / Jacobian of stage 2 estimating equations, an observed information matrix. However, bread.glm() instead provides expected information, so I think we should stick with that (they should be close to one another). Ordinarily I'm happy with this substitution, but I'm wondering whether I should continue to be happy with it in special cases of the sort we'd like to accommodate here, where the estimated parameter is not necessarily an argmax of the quasi-likelihood objective. We just need to make sure it's expected information evaluated at the parameter as estimated in the values_from model, as oppposed to x. (This might be the location of a restricted maximum, restrictions reflecting treatment effect values assumed under a null hypothesis, or reflecting a stipulation that subgroup treatment effects be the same.) As long as it's expected information as evaluated at some specific value of the parameter, the information identities underlying the usual GLM variance formula $\mathcal{J} = \mathbf{X}'\mathbf{W}\mathbf{X}$, where $W$ is diagonal with $W_{ii} = \frac{(\partial \mu_i/\partial \eta_i)^2}{\mathrm{var}_\theta(y_i)}$, hold. ↩

[benthestatistician]

(Dispelled my lingering concerns about the solution recommended $\uparrow$.)