Conflicting n in bread and meat calculations when there are 0 weights

[benthestatistician]

We use sandwich::bread() to calculate bread matrices, and sandwich::meatCL() for some of our meat calculations. Each of them has an internal calculation of n, which should match the n calculated within a vcov-calculator function that calls them. But it seems that sandwich::bread.lm() uses summary.lm() for its n whereas sandwich::meatCL() gets its by way of sandwich:::estfun.lm(), and that the two sources needn't agree.

I'll use the nuclearplants data set for a quick example with 0 weights.

data(nuclearplants, package="optmatch")
table(nuclearplants$pt)
##'  0  1 
##' 26  6
lm0 <- lm(cost~pr, weights=pt, data=nuclearplants)
sum(summary(lm0)$df[1L:2L]) #n according to `sandwich::bread()`
##' [1] 6
NROW(sandwich::estfun(lm0)) #n according to `sandwich::meat()`, `sandwich::meatCL()`
##' [1] 32

Consistent with this, sandwich returns oddly small standard errors.

##> sqrt(diag(sandwich(lm0)))
##' (Intercept)          pr 
##'       4.04        5.03 
##> sqrt(diag(vcov(lm0)))
##' (Intercept)          pr 
##'       23.2        32.9 
##' sqrt(diag(car::hccm(lm0, type="hc0")))
##' (Intercept)          pr 
##'       21.6        26.8 

(... and our HC1's are similiarly too small in a non-public example I'm seeing.)

We can fix much or all of the problem as it affects us by writing and using our own local versions of sandwich::bread.lm() and sandwich::bread.glm(). These functions can follow the sandwich versions in calling summary() for the matrix, and for the scaling factor in the glm case, but then they should depart in their calculation of n (in both of those functions, sum(sx$df[1L:2L])). My proposal would be to look at the length of the residuals recorded by summary(), which I expect to match the row dimension of the model's model.matrix, which is what estfun() starts from.

[josherrickson]

This seems reasonable to me. Those functions are pretty minimal. The only tricky part will be to figure out the new dispatching. It appears we call sandwich::bread a few times, e.g. https://github.com/benbhansen-stats/propertee/blob/43173f61c9d045ac65d2d67a23e6f6a0beb227ff/R/SandwichLayerVariance.R#L67 and we'll have to instead special-case if the input is lm or glm, but that's not a big deal. Just need to make sure we catch all such cases.

I can work on this next week unless @jwasserman2 is eager to hop on it.

[jwasserman2]

Good catch here, Ben. We have a bread method for teeMod objects, actually: https://github.com/benbhansen-stats/propertee/blob/43173f61c9d045ac65d2d67a23e6f6a0beb227ff/R/teeMod.R#L153 .get_tilde_a22_inverse() calls .get_a22_inverse(), so I think the only change we need is in .get_a22_inverse(). https://github.com/benbhansen-stats/propertee/blob/43173f61c9d045ac65d2d67a23e6f6a0beb227ff/R/SandwichLayerVariance.R#L533-L547 https://github.com/benbhansen-stats/propertee/blob/43173f61c9d045ac65d2d67a23e6f6a0beb227ff/R/SandwichLayerVariance.R#L309-L317 I don't see why we can't use the code in .get_tilde_a22_inverse() that calculates nq, nc, and n in .get_a22_inverse() to re-scale the output of bread.lm(). Josh, you can go ahead handling this one.

EDIT: We'll also need to add the sum(sx$df[1L:2L]) call Ben mentioned in his original comment so we know what sandwich thinks n is.

[benthestatistician]

Nice point, Josh W. I had seen bread.teeMod() and started to ask whether something like what he's outlined could be a solution as well, but then held off in the interests of keeping my comment simple.

I'll be happy with either type of solution to the problem, new bread.lm()/bread.glm() methods or fix through bread.teeMod().

[benthestatistician]

Note re Josh W's note that

We'll also need to add the sum(sx$df[1L:2L]) call Ben mentioned in his original comment so we know what sandwich thinks n is.

--Actually, I don't know that we'll need that. sandwich::sandwich() has n <- NROW(estfun(x)), matching the n used in meat(), meatCL() and recommendations in this thread. So I'm not seeing that we'd need the mistaken n figured in bread.lm() and bread.glm().

[jwasserman2]

bread() (specifically bread.lm()) doesn't use NROW(estfun(x)), though, it uses sum(sx$df[1L:2L]). So we'll need sum(sx$df[1L:2L]) to re-scale the output of bread().

[benthestatistician]

Should have been clearer: my last comment pertained to the strategy where we write new bread.lm() and bread.glm() methods, setting things up so that these rather than the sandwich package's methods are invoked for our A11 and A22 calculations. On that path, if we can get it to go, I think we get to avoid having to undo the sum(sx$df[1L:2L]). (OTOH, if we leave there methods in place while making changes to bread.teeMod(), then I'd agree that we'll need the sandwich::bread version of n.)

[josherrickson]

To clarify, something like

nq <- nrow(stats::model.frame(x)) 
nc_not_q <- sum(!ca@keys$in_Q) 
sum(sx$df[1L:2L])
out <- out * n / nq 

added to get_a22_inverse should address it?

Do we have a example I can test against?

[benthestatistician]

1. To take the route of patching within .get_tilde_a22_inverse(), I think it would have to be

   nq <- nrow(stats::model.frame(x)) 
   nc_not_q <- sum(!ca@keys$in_Q) 
   n <- nq + nc_not_q 
   nq_as_figured_by_s_b <- sum(sx$df[1L:2L])
   out <- out * n / nq_as_figured_by_s_b 

2. From the maintainability point of view, introducing our own bread methods would be quite a bit better. Even now when I have the details and relevant tech docs near to mind, it's very difficult to keep track of the constants when they're spread across functions in this way. By making the change within .get_tilde_a22_inverse(), we saddle our future selves with a .get_a22_inverse() that doesn't really get an A22 inverse, and a .get_tilde_a22_inverse() that does some of the things the spec describes as converting an $A{22}^{-1}$ into and $\tilde{A}{22}^{-1}$, as well as some other things relating to correction of the bread. In other words, I'm coming around to seeing some strong advantages for new bread methods rather than a .get_tilde_a22_inverse() patch.
3. Adding to the last point, fixing bread.teeMod() without fixing bread.lm() leaves us with the same vulnerability when we pull $A_{11}^{-1}$'s.
4. To deal with the problem of having to special-case for lm's and glm's that you noted above, @josherrickson, perhaps we could define our own bread() generic. Perhaps then, when from within a propertee function we invoked bread() we'd get dispatch on our own methods first, with sandwich methods as fallback. If this works, I don't imagine it would require propertee::bread() to be exported.
5. To answer Josh E's second question: I don't have a shareable example with teeMod, but I suppose the nuclearplants example at the top of the thread could be refashioned into one of those.

[josherrickson]

bread is already a generic in the Sandwich package, which we import. Should we create our own, bread_lmitt or somesuch, with generics for lm, glm, and maybe propertee? Then as Ben alludes to in 4, we do a tryCatch to fall back to regular bread?

[benthestatistician]

I was thinking of deliberately overloading sandwich::bread(), sort of like the lmerTest package overloads lme4::lmer(). (An example I have you to thank for introducing me to, @josherrickson.) Differences being that our bread would be a generic, and it wouldn't be exported.

If one or the other of these differences make the overloading strategy seem unlikely to you to succeed, then yes, I'm happy to have us create a fresh generic. I might prefer _bread() or .bread() to bread_lmitt(), as the only lmitt or propertee connection is that we've set things up in such a way that we'll more frequently stumble on the sandwich::bread() issue that we've identified than other users might.

Re what to do when our bread sees an object that's neither an lm nor a glm, yes in these cases we can go with the sandwich::bread() answer. From a glance at sandwich::bread()'s other methods, I'm not seeing this problem there. In fact, the default method uses a function sandwich:::nobs0() for this purpose, which in turn falls back to NROW(residuals(x, ...)), which is pretty close to the solution I've suggested.

[josherrickson]

As far as overloading sandwich:::bread.lm, it appears that defining our own bread.lm works for dispatches from bread, but dispatches internal to sandwich (e.g. in sandwich) don't:

> bread(lm0)
[1] "in propertee's bread"
            (Intercept)        pr
(Intercept)    10.66667 -10.66667
pr            -10.66667  21.33333
> sandwich(lm0)
            (Intercept)        pr
(Intercept)    16.34794 -16.34794
pr            -16.34794  25.33712

We can get the result we want if we are willing to muck with sandwich internals:

> utils::assignInNamespace("bread.lm", bread.lm, "sandwich")
> bread(lm0)
[1] "in our bread"
            (Intercept)        pr
(Intercept)    10.66667 -10.66667
pr            -10.66667  21.33333
> sandwich(lm0)
[1] "in our bread"
            (Intercept)        pr
(Intercept)    465.0081 -465.0081
pr            -465.0081  720.7004

I'm of mixed feelings about whether we can do this - on the one hand, reaching into another package and replacing their functions feels much more troubling than overloading. On the other hand, practically speaking, both would "break" user space in the same way, they'd both require most users to restart R with no packages loaded if they wanted to get back to vanilla Sandwich. (Technically, with overloading, I believe if you use unloadNamespace, that would remove the need to unload sandwich as well.)

[benthestatistician]

Thanks Josh! I should have been more clear, though: the intention of my proposal was to avoid interfering with users' independent uses of the sandwich package. I.e., my hope was to have bread() dispatch to our bread methods rather than sandwich's when and only when invoked from within a propertee function. It sounds like you're getting that just by overloading bread.lm() -- if so, that's great, and I wouldn't have us muck with the sandwich package's namespace or anything like that.

[josherrickson]

> library(sandwich)
> data(nuclearplants, package="optmatch")
> nuclearplants$id <- seq_len(nrow(nuclearplants))
> lm0 <- lm(cost~pr, weights=pt, data=nuclearplants)
> sqrt(diag(sandwich(lm0)))
(Intercept)          pr 
   4.043259    5.033599 
> sqrt(diag(vcov(lm0)))
(Intercept)          pr 
   23.24920    32.87933 
> sqrt(diag(car::hccm(lm0, type="hc0")))
(Intercept)          pr 
   21.56405    26.84586 
> des <- rct_design(pr ~ uoa(id), data = nuclearplants)
> mod <- lmitt(cost ~ 1, data = nuclearplants, design = des,
+ weights = nuclearplants$pt)
> sqrt(diag(sandwich(mod)))
[1] "in propertee's bread"
(Intercept)         pr. 
   21.56405    26.84586 
> sqrt(diag(vcov(mod)))
[1] "in propertee's bread"
(Intercept)         pr. 
   21.90909    27.27542 
> sqrt(diag(car::hccm(mod, type="hc0")))
[1] "in propertee's bread"
(Intercept)         pr. 
   21.56405    26.84586 

This is using length(sx$residuals) as the scaling factor. Is this what you're looking for, both in terms of size of standard errors, and which dispatches to sandwich::bread.lm vs propertee::lm?

[benthestatistician]

Yes! Precisely what I was looking for. Thanks!

[josherrickson]

Great, I've pushed it to main. For the sake of documentation, I used length(sx$deviance.resid) in the glm version.

[benthestatistician]

For sake of posterity, here's Josh's bread.{lm,glm}() fix, in the form of a patch to https://github.com/cran/sandwich/blob/master/R/bread.R: bread.patch