Working weights are not updated properly in bayesglm.R

[NikTuzov]

Hello:

I believe there is a bug in bayesglm.R. I attached bayesglm_A_1.R where I inserted a few cat statements. If you run it using the code from the other file, you'll see from the output that the last two weights that correspond to pseudo-data are initialized but never updated.

The problem is that in .bayesglm.fit.loop.updateState() the working weights are computed as:

w.star <- c(w, sqrt(state$dispersion)/priors$scale)

but to update them one needs to use state$prior.sd in the denominator, which should be equal to priors$scale only in the beginning. It looks like you only need to change one line to fix it.

Regards,
Nik Tuzov

bayesglm_issue.zip

[gfinak]

Thanks for this report. @amcdavid can you confirm?

[gfinak]

Running this I also found it differs from the output of arm::bayesglm.fit in version 1.10.1, and the change of
w.star <- c(w, sqrt(state$dispersion)/priors$scale)
to
w.star <- c(w, sqrt(state$dispersion)/state$prior.sd)
leads to a non-convergence warning.
@amcdavid this raises the question whether we should import the function from arm, perhaps have an option to run the old code for back-compatible reproducibility.

[amcdavid]

Thanks Nik for the report. I think I agree that this may be a bug, but I definitely don't understand this code well (because we didn't write it). As such, I am hesitant to make changes without a failing unit test. Do we have an example, and an expected result to show what exactly is wrong about this, and verify that changing this line fixes it?

This is code from Gelman et al's arm package that we directly included into MAST. (We did this rather than Importing arm because there were changes being made to arm by one of the authors without regression testing and this was breaking things in MAST.) As such, I can only guess at what it's doing. On line 434 we seem to be augmenting the weights for a weighted least squares with pseudo observations that represent the prior. Then on lines 473-475 we seem to be doing something to represent the t-prior as a scale mixture. As far as I can tell, state$prior.sd is never used again in the main loop, and I suppose it should be.

[gfinak]

The main loop calls .bayesglm.fit.loop.main.ideal which calls another function, passing in priors$scale to update the state$prior.sd variable.
The latest arm code seems to be entirely in Fortran, Neve mind, I seem to have had a stroke.I don't find these utility functions there at all.

[amcdavid]

Where do you see fortran in the repo? I don't see any. FWIW, here's the diff between our version of arm and the latest: cran/arm@82ebb3c

I would rather move away from arm altogether and implement the logistic model in stan or rstanarm. Arm is legacy code, and doesn't have any unit or regression testing, whereas stan and rstanarm are tested and more modern. They may be faster anyways, since they are in C. If we do that we can add an option to use the "old version".

[gfinak]

Sorry, you're right, not sure what the heck I saw.

[amcdavid]

And actually this is the diff between our version of bayesglm and the latest version. It has been refactored to the point of being unrecognizable. Presumably for the better, but I am hesitant to mess with things given that lack of unit tests.

[gfinak]

Okay, I think we have to update things. I like the idea of using rstanarm.

[amcdavid]

👍 It's been on my todo for a while. I probably won't be able to get to it until middle of next week, so feel free to take a stab. In any case, we wouldn't want to push anything to bioconductor until after the bioconductor release.

[NikTuzov]

Hello Greg and Andrew:

Thanks for responding.
Can you say whether that bug or any other major bug was present in MAST around spring 2017 when Soneson and Robinson worked on their paper in which they evaluated MAST performance.

Nik

[gfinak]

It's not clear to me yet what the impact of this particular bug was on that paper, but it was present at the time. That said, it was a bug in arm as well, when MAST was released. It's obviously a priority to address it.

[NikTuzov]

Please make sure a similar issue is not present in rstanarm or whatever package you are going to use.

[gfinak]

@NikTuzov Let me be clear, that this was a bug in arm.

[suyusung]

@NikTuzov We have decided to revert bayesglm back to an older version (when Andrew just wrote it from scratch) because the new one was buggy. We made several mistakes when we tried to clean the codes up using many functions embedded in the bayesglm(). In any case, bayesglm() and bayesglm.fit() in arm_1.10.1 should be a good (if not buggy) version to use. Having said this, I use the zip file in your post to test it, the current bayesglm() is working. I'm not sure if I overlooked anything mentioned in your post. If you can be more clear about the bug and whether or not the bug still affects the current version, you can help us make arm better. Thanks!

[NikTuzov]

@gfinak
Hello Greg:

What is the present state of things? Are the weights updated properly?

Regards,
Nik