Distance term for ergm (Minkowski and spherical distances)

[CarterButts]

This update adds a term called distance that can be used for binary ERGMs (directed or undirected). The term includes support for

• Minkowski metrics with arbitrary exponents in arbitrary numbers of dimensions
• Spherical (great circle) distances on spheres of arbitrary size (defaulting to the IUGG geosphere, for lat/lon coordinates)
• Logarithmic distance scaling (necessary for most social network applications)
• Offsets and "safety" thresholding for log distances
• Coordinates passed either directly in matrix form, or via an embedded network attribute

This functionality was previously in an ergm userterms package, but given recent demand from Steve, Sam, and Martina, I dusted the package off, upgraded it, and then integrated the functionality into the ergm package.

The function is fully documented (with examples), and includes some basic tests. I have confirmed that the package w/the term passes R CMD check, the tests pass, and the examples work correctly. I think all is copacetic, but let me know if there are things that require fixin'.

[krivit]

Thanks, Carter! Definitely worth including. My thoughts so far:

1. This could probably use a c_ function rather than d_.
2. Positions are vertex attributes, so semantically it would make more sense to store them there. This would also work better with ergm.multi and tergm.
3. ergm_get_vattr() API can fetch multiple vertex attributes. For example, nodecov(c("Lat", "Long")) or nodecov(~cbind(Lat, Long)) will internally fetch a 2-column matrix.
   1. At this point, list-form vertex attributes aren't supported, I think, but that's easily fixed.

[CarterButts]

Hi, Pavel. To your points:

1. Although we have obviously discussed it in the past, I forget the specs for c_* functions versus d_* functions. Is there much of an advantage for independence terms? Also, where are they currently documented? They postdate the ergmuserterms package, and the ergm 4 paper doesn't discuss them. I can port the d_* term if you think it is needed, but I need the relevant specs. (I always use the d_* form, being Olde Skoole.)
2. I don't store the coordinates on the vertices because this is poses various annoying issues. We need to support an arbitrary number of dimensions (one cannot assume that the data is two-dimensional), so storing the coordinates in a vertex attribute would require using vector-valued attributes. network certainly supports that, but it's much more work for the user to have to decompose coordinates into individual vectors and then write them into place using the set.vertex.attribute command than to just pass a matrix (or stuff it in a network attribute). It's also more of a pain for us to retrieve, and creates checks that have to be done for consistency, etc. I actually started to implement that, but then realized that it was inferior to just using a network attribute, so stopped. One could name every dimension, as you suggest, but then the user has to (1) store every dimension under its own name, and then (2) pass the whole vector of names to the function. This seems to me unlikely to be easier than just putting the data in a network attribute. If we really wanted to support this, it could be implemented, but then there's the question of how to signal whether an attribute name (if only one is passed) refers to a vertex or a network attribute. I suppose one way to go would be to first check to see if there is e.g. a vertex attribute, and then go look for a network attribute if that fails. Seems inelegant, and could lead to unexpected behavior if the user doesn't read the docs, but better than having to add yet another argument to the function to tell it what kind of attribute we want....
3. I wasn't familiar with that, but can use it if you think it is important to have the functionality in (2). Let me know if you have a better alternative to the failover idea (i.e., that if a single character string is passed, first we look for a vertex attribute with that name - if none exists, we then look for a network attribute with that name and use it instead).

[krivit]

I can only answer briefly at the moment:

The c_ API is documented in the Terms API vignette, and most of the terms have been ported over. It reduces boilerplate since it handles 1 toggle at a time, and it receives the current edge state as a part of the call rather than having to look it up, so it's somewhat faster.

There might be a simpler way to do it, but apply(M, 1, identity, simplify = FALSE) creates a list of rows suitable for assignment to a vertex attribute.

If you really want to access the network attribute via a nodal attribute API, you can: for example, nodecov(~.%n%"M")) will grab the network attribute matrix.

[CarterButts]

OK, making those changes (and adding a few more features). Will update....

[CarterButts]

@krivit OK, I have made the suggested changes, and added some new functionality. Among other things, I switched from d_* to c_* form for the changescore, and one can now either use vertex or network attributes (in addition to matrices) to supply coordinates. It is now possible to mix spherical and Minkowski effects, and I added a scaling option. I updated the tests and docs accordingly. Everything seems to be passing checks, and my own examples/tests are working fine. So it seems to be copacetic, as far as I can tell.

[sgoodreau]

@CarterButts This is superb, thank you! So much good functionality, plus the base case that is most useful for our students (~ distance(c('X', 'Y'), log=FALSE)) is syntactically straightforward and fits in well with their level of knowledge at the point we want to use this. I'm very glad you added the vertex attributes; that fits with the logic I would use here and how we train our students to think about this. And I love the ease to switch between Euclidean and Manhattan, the ease of logging or not, and of course the great circle functionality (something I'm unlikely to use, but great it's there).

Forgive me if you already mentioned this. Clearly log distance is commonly used in practice given its superior properties over distance in many practical applications. And I see that you have a scale factor feature in there also. Do folks ever do other transformations in practice, particularly powers? E.g. taking the sqrt of distance? (I realize 2D Euclidean distance already involves taking a sqrt, but that doesn't preclude one from doing it again). If so, is this functionality in there in a way I'm not seeing? Or is it just not done in practice?

Thanks!

[sgoodreau]

I should add that I pulled your branch and played with it a bit, including trying out some made-up examples just to ensure I understood the syntax fully. Not a deep set of tests, but at least one more pair of eyes saying the basic parts work as expected.

[CarterButts]

@sgoodreau Glad that this can be useful! To your question, we have usually (so far) seen spatial interaction functions (i.e., marginal tie probabilities, as a function of distance) with forms like pb/(1+a*d^c) or pb/(1+a*d)^c, where pb is the base tie probability at distance zero, a is a scaling factor, d is distance, and c is an exponent. In a sparse network, we have in the second case

ln Pr(y+_ij)/Pr(y^-_ij) apx log Pr(y^+_ij)  
                                 = log pb -  c log(1+a*d_ij),

so the conditional log odds of the i,j edge are approximately linear in the log of the rescaled distance added to the offset.

Now, consider an ERGM with an edge term and an edgecov of the form z_ij = log(1+a*d_ij). This gives the conditional log odds

ln Pr(y+_ij)/Pr(y^-_ij) = log θ_e + θ+z log(1+a*d_ij).  

We can immediately recognize this as equivalent to the above expression. So what this tells is is that adding a term for the log of the scaled base distance plus an offset yields a power law distance dependence, with the ERGM coefficient θ_z being the power law exponent. The power law here arises from the fact that we are logging both sides (approximately as in an OLS, where logging x and logging y gives you a power law), so we don't need to specify a power law coefficient in our distance term; that's what the ergm coefficient controls.

(Note that this is approximating what we have called an "attenuated power law" in our work. It has power law tails, but is locally concave down near the origin. The other obvious family, pb/(1+a*d^c) , can't be implemented this way. This is unfortunate, but on the bright side the two will agree when a*d>>1. So, when your goal is mainly to control for long-range behavior, this should work OK.)

BTW, we can contrast this with what happens when we have a term like d^c. Then we get

ln Pr(y+_ij)/Pr(y^-_ij) = log θ_e + θ_z d_ij^c  

and since, in the sparse regime, the log odds is approximately the log probability,

Pr(y^+_ij) apx exp(θ_e) exp(θ_z d__ij^c)

Very different: if θ_z<0, we are saying that tie probability falls off exponentially fast in d^c. This pretty well damps out long-range ties: we can put a circle around each ego, and be pretty sure that they won't have any ties outside that circle. This might well be a good model for some systems, though most interpersonal networks I have seen have too many long-range ties for that model.

On reflection, though, for this exponential model, I was more concerned in the code with changing the metric (the power that goes into how we define the distance), and not the c parameter above. I don't know if it will be useful, but it's a quick addition - I'll add it now.

[CarterButts]

@sgoodreau OK, added the pow argument. It allows one to raise the rescaled distance to an arbitrary power before further calculations. (This is distinct from the Minkowski exponent.) Let me know if there is other critical functionality that you need....

[sgoodreau]

@CarterButts, thank you much. That is all the functionality I was thinking of (and more), so we're good to go on my end. @martinamorris or @smjenness did you have anything else you were thinking of?

[krivit]

Thanks, everyone, but particularly @CarterButts!

I've rebased this over some changes to master that enable nodal attributes API to handle vector-valued vertex attributes.

I am going to make a few more changes, including putting it into the template that makes it work for valued ERGMs as well.

For future reference, when proposing pull requests, please do not make changes to the master branch of your fork, but rather create a unique branch name, since it's confusing and potentially dangerous to have multiple master/main branches floating around my clone.

A few more items for discussion in the meantime,

1. Are scale and radius redundant?
2. Would it be more user-friendly to have two separate terms, one for sphere distance and the other for Minkowski? If not, what about aliases? What should they be?
3. For Minkowski distances, would it be useful to specify a vector with a scale for each dimension?
4. For Minkowski distances, do the cases for metric of 0 (Hamming) and +Inf (Chebyshev) work? These could be very useful.

[CarterButts]

Hi, Pavel -

A few more items for discussion in the meantime,

1. Are `scale` and `radius` redundant?

No; they overlap in one special case, but are not in general redundant. scale applies to the final distance calculated by any metric, while radius only affects the sphere used for calculating geospherical distances. If you are only using great circle distances, then you could use the radius argument to change the scale of your distances, but that wouldn't carry over if you had either additional Minkowski dimensions, or if you weren't using great circle distances at all. Also, from a user-level standpoint, it is helpful to allow radius to keep a very specific default value, so that the typical user (who is going to be using lat/lon coordinates) does not have to look up the Earth radius that they want to use; by turns, it is pretty handy to be able have a one-stop shopping point for altering the scale.

2. Would it be more user-friendly to have two separate terms, one for sphere distance and the other for Minkowski? If not, what about aliases? What should they be?

It seems to me more natural to just have a single term, for several reasons. (1) in geospatial work it is not uncommon for users to switch between lat/lon and projected geometries in different settings, and having a single term makes it clear where one should go for the functionality in either case. (2) a lot of the basic functionality involved in getting the attributes and processing the distances is shared, so this minimizes code replication. And perhaps more importantly, (3) you may in some cases want to merge spherical and non-spherical distance elements (as in the case of elevation, but perhaps also Blau spaces), and the current distance term does this natively with minimal user effort and computational overhead. Splitting would thus seem to add overhead and complexity, and I'm not sure what the benefit would be.

Relatedly, it's hard to suggest an alias, because I'm not sure what the motivating use case is. I have no particular objection if it serves a purpose (other than not multiplying terms beyond necessity, and potentially setting up name conflicts).

3. For Minkowski distances, would it be useful to specify a vector with a `scale` for each dimension?

The scale parameter applies to the distance, not to the underlying dimensions. One could also introduce an argument that changes the metric by altering the dimensions themselves, but that would be something different.

4. For Minkowski distances, do the cases for `metric` of `0` (Hamming) and `+Inf` (Chebyshev) work? These could be very useful.

Not currently, though they could be added. I can indeed see these being useful for non-physical spaces....

[krivit]

I am in the process of converting the code to work for both binary and valued ERGMs. Please don't make any further edits until I push.

In the meantime, more discussion items:

5. ergm_get_vattr() can now gracefully handle vector-valued nodal attributes to return a matrix, and there is a way to access network attributes if you really want to. Do we still need the fallback?

6. Usually, vertex attribute specifiers that implement the API are named attr. Any objections to renaming coord to that?

[CarterButts]

I've rebased this over some changes to master that enable nodal attributes API to handle vector-valued vertex attributes.

BTW, can you clarify?

[sgoodreau]

5. ergm_get_vattr() can now gracefully handle vector-valued nodal attributes to return a matrix, and there is a way to access network attributes if you really want to. Do we still need the fallback?

I don't understand what this means, sorry.

Let me clarify my perspective of EpiModel's user base's needs. Typical use case woud be to have two nodal attributes representing coordinates in 2D space. Their names are arbitrary, but let's call them X and Y. Ideal code would be something like ~distance(c('X', 'Y')), which is as it currently is (with other optional arguments). I would prefer to not have to see an NxN matrix of distances constructed every time one wants to determine a change statistic that only requires knowing the distance of a single dyad, as that seems inefficient.

Is what you are saying here consistent with that?

[CarterButts]

@sgoodreau indeed, this is what the current term does. Actually, the current term lets you do any of the following:

• Pass a vector of vertex attribute names, each of which contains a single coordinate dimension (e.g., c("x","y"), c("lat","lon","elevation"), etc.);
• Pass an attribute name for a network attribute containing an n x d coordinate matrix; or
• Directly pass a coordinate matrix or data.frame.

So the user can employ whichever is most convenient. (I usually want to work with coordinate matrices, for instance, but you've said that your folks would rather do it via vertex attributes. Everyone wins!)

That's why I'm a bit confused....

[CarterButts]

Overlooked this before, but:

6. Usually, vertex attribute specifiers that implement the API are named attr. Any objections to renaming coord to that?

Please don't do that, because coord can take either a coordinate matrix, an attribute name, or a vector of attribute names. And the latter must refer to coordinates. The current argument name works for all cases, because it focuses on the content of what is being provided, not the way in which it is being provided.

[krivit]

BTW, can you clarify?

The most idiomatic way to pass each vertex's coordinates to is via vertex attributes. ergm_get_vattr() API previously supported passing multiple scalar-type vertex attributes that are then combined into a matrix. I have added support for vector-type vertex attributes that are then bound into rows of a matrix. I have pushed these changes to statnet/ergm@master.

There are two ways I could have integrated this change into this branch, merging or rebasing. Since it's a development branch, I have done the latter to keep the log a bit cleaner.

[krivit]

1. and 2.: I am not necessarily against these explanations, but I want to explore possibilities before we finalise this, particularly since this term will hopefully be used outside of the geospatial context. Do we want to assume a nontrivial value for radius? For example, what if someone is using miles. What if we instead allowed a number or a character placeholder, e.g., "earth.km", "earth.mi", "unit", etc.?

3. Unless I am misreading the current implementation, scale is currently applied after taking metricth root, so applying it to each individual difference (and to the sphere distance) separately would have exactly the same effect. Proliferating arguments similar to each other in name and functionality but subtly different strikes me as a bad idea; what's the harm of moving scale down to individual dimensions, except perhaps for a small performance hit?

4. This would also help address some of #481 as well. Now that I think about it, do we want to have a user interface more like R's dist() function's, in which the second argument (method) could be "spherical", "euclidean", etc., with an option to specify, e.g., "spherical+euclidean" or c("spherical", "euclidean")?

5. and 6. To recap (and this is also for @sgoodreau), nodal attributes UI (i.e., how the end-user can use it) lets you:
* grab network attributes by passing ~.nw%n%"coord";
* values to be used verbatim by passing I(coord) (where coord is a matrix variable in the formula's environment);
* c("X", "Y", "Z") where X, Y, and Z are scalar-valued vertex attributes; and
* (as of a few days ago), the latter can also be vector-valued vertex attributes (i.e., each vertex has a vector associated with it), in which case these vectors are concatenated, then stacked into a matrix.
This already works for other terms that can take n*k matrix input, such as nodecov, so there would be a consistent UI across all terms that do. attr (shortened from attrname, once we could pass things other than names) is the conventional name for this type of argument.