Covariance Methods for MPLE IV

[schmid86]

I added code to estimate MPLE covariance matrices. If a user requests the MPLE of a dyad-dependent model, standard errors are obtained from the glm()-function, which uses the inverse Hessian matrix to calculate the covariance matrix. However, MPLE standard errors obtained from the logistic regression output lead to MPLE-based confidence intervals with coverage rates far below the nominal confidence level (see for example van Duijn, Gile and Handcock (2009)) and are therefore suspect.

I added two new methods in the code that are simulation based approximation of standard errors. One uses the Godambe matrix (sometimes also referred to as the Sandwich estimator), the other one approximates standard errors using parametric bootstrapping. Simulation studies have shown that these standard errors are more reliable than MPLE s.e from the logistic regression output. See Schmid and Desmarais (2017) and Schmid and Hunter (2023). The file ergm_mplecov.R contains the code for these two methods.

Sample code: data(florentine) m1 <- ergm(flomarriage~ edges+triangles, estimate="MPLE", control=control.ergm(MPLE.covariance.method="Godambe")

m2 <- ergm(flomarriage~ edges+triangles, estimate="MPLE", control=control.ergm(MPLE.covariance.method="Bootstrap")

References: Schmid, C. S. and Desmarais B. A., Exponential Random Graph Models with Big Networks: Maximum pseudolikelihood estimation and the parametric bootstrap, IEEE International Conference on Big Data, 2017.

Schmid, C. S. and Hunter D. R., Computing Pseudolikelihood Estimators for Exponential-Family Random Graph Models, Journal of Data Science, 2023

[krivit]

To confirm, does this supersede #235, so it can be closed?

[drh20drh20]

To confirm, does this supersede #235, so it can be closed?

Yes. @schmid86 couldn't get the code on which the earlier pull request was based to work, so this is a fresh start.

[drh20drh20]

@krivit , a request: Since the journal editors have required that the (already accepted) code runs against the public version, would it be possible to release a version of the code that does not address each of these comments, or at the very least zero in on the comments that are non-negotiable?

[krivit]

@krivit , a request: Since the journal editors have required that the (already accepted) code runs against the public version, would it be possible to release a version of the code that does not address each of these comments, or at the very least zero in on the comments that are non-negotiable?

I'll go ahead and fix the code issues, but someone needs to work how offsets fit into this.

[krivit]

@drh20drh20 , @schmid86 , also, please add some test cases.

[drh20drh20]

Not sure I understand your comment @krivit ; the check is for the SE estimates. What sort of a test were you envisioning? Also, are offsets implemented at the moment?

[schmid86]

@drh20drh20 yes, offsets are implemented on this branch

[drh20drh20]

Oh! Thanks @schmid86 . Yes, I will write something for this then

[krivit]

Not sure I understand your comment @krivit ; the check is for the SE estimates. What sort of a test were you envisioning?

Check that vcov(fit) estimates returned by the code equal to the known-good Godambe and bootstrap estimates.

[drh20drh20]

Not sure I understand your comment @krivit ; the check is for the SE estimates. What sort of a test were you envisioning?

Check that vcov(fit) estimates returned by the code equal to the known-good Godambe and bootstrap estimates.

That's what it does, but just for the diagonal entries.

[krivit]

That's what it does, but just for the diagonal entries.

Ah, I see. I saw coef() and thought it was testing point estimates. (It might be better to test sqrt(diag(vcov(fit))) instead, though.

[krivit]

@drh20drh20 , @schmid86 , can you please update the tests to test offsets and also fix the default values to be powers of 2 (since that's what we use for the rest of the analysis)? I need to put the release out as soon as possible.

[drh20drh20]

Please give this a shot.

[drh20drh20]

If I'm reading this correctly, some of the tests failed because different platforms gave different results using the same random seeds. Could that be correct?

[drh20drh20]

@krivit I fixed the rounding error that created the only test error on my machine & pushed the change. Is it possible to run the tests again to see if that was the culprit? There was output in one of the failures that suggested it was more than this (hence my previous question regarding seeds) but how can we be sure?

[krivit]

@krivit I fixed the rounding error that created the only test error on my machine & pushed the change. Is it possible to run the tests again to see if that was the culprit? There was output in one of the failures that suggested it was more than this (hence my previous question regarding seeds) but how can we be sure?

The tests run automatically every time you push.

[krivit]

@krivit I fixed the rounding error that created the only test error on my machine & pushed the change. Is it possible to run the tests again to see if that was the culprit? There was output in one of the failures that suggested it was more than this (hence my previous question regarding seeds) but how can we be sure?

@drh20drh20 , you can also run all tests locally using

testthat::test_local(load_package="library")

The argument is needed, because otherwise weird errors result (see https://github.com/r-lib/testthat/issues/1751).

[drh20drh20]

@krivit I fixed the rounding error that created the only test error on my machine & pushed the change. Is it possible to run the tests again to see if that was the culprit? There was output in one of the failures that suggested it was more than this (hence my previous question regarding seeds) but how can we be sure?

@drh20drh20 , you can also run all tests locally using

testthat::test_local(load_package="library")

The argument is needed, because otherwise weird errors result (see r-lib/testthat#1751).

Thanks @krivit. Am on vacation so response time slower than usual, but pushed what appear to be fixes. Everything passed when I ran test_local as you suggested