Distance term for ergm (Minkowski and spherical distances)

R/InitErgmTerm.R

@@ -25,6 +25,7 @@
 #   C:   <concurrent>       <cycle>           <ctriple>=<ctriad> 
 #   D:   <degree>           <density>         <dsp>
 #        <dyadcov>          <degcrossprod>    <degcor>
+#        <distance>
 #   E:   <edgecov>          <edges>           <esp>
 #   G:   <gwb1degree>       <gwb2degree>      <gwdegree>
 #        <gwdsp>            <gwesp>           <gwidegree>
@@ -2412,6 +2413,296 @@
        )
 }
 
+################################################################################
+
+#' @templateVar name distance
+#' @title Inter-point distances
+#' @description This term adds a single statistic to the model whose value is
+#'   the sum over all edge variables of the edge variable value times either
+#'   the distance between the respective vertices, or its log (if
+#'   \code{log==TRUE}).  The \code{coord} agument must contain a vector,
+#'   matrix, or \code{data.frame} of coordinates (vertex x dimension), the
+#'   name of a network attribute with such coordinates, or a vector of vertex
+#'   attributes to be used as coordinates.
+#'
+#' @usage
+#' # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714,
+#' #             log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1)
+#'
+#' @param coord node by dimension coordinate matrix, name of a network
+#'              attribute containing the coordinate matrix, or vector of
+#'              vertex attribute names containing coordinates (to be used
+#'              jointly); for \code{sphere==TRUE}, first two coordinates
+#'              must be angular units, in lat/lon order
+#' @param metric power to use for the Minkowski metric
+#' @param sphere logical; should great circle distances be used (rather than
+#'               Minkowski distances in free space)?  If \code{TRUE}, the first
+#'               two coordinates are respectively interpreted as latitude and
+#'               longitude coordinates on a sphere of specified radius
+#' @param radius for spherical distances, the radius of the sphere to use;
+#'               defaults to the IUGG mean Earth radius, in km
+#' @param log logical; use log rather than raw distances?
+#' @param mindist for log distances, a lower threshold; for dyads such that
+#'                the observed distance plus any offset is less than the
+#'                minimum, the log of the minimum is used instead.  Ignored for
+#'                raw distances
+#' @param distoff for log distances, an offset to be added to observed
+#'                distances before taking the logarithm
+#' @param scale factor by which raw distances should be rescaled prior to other
+#'              operations (notably, offsetting, logging, and thresholding)
+#' @param pow power to which scaled distances should be raised prior to other
+#'            other operations (notably, offsetting, logging, and thresholding)
+#'
+#' @details Either spherical (\code{sphere=TRUE}) or Minkowski metrics
+#'   (\code{sphere=FALSE}) may be selected.  For the latter, any number of
+#'   dimensions may be supplied, but if \code{coord} is given as a vector the
+#'   space in question is assumed to be one-dimensional.  The choice of
+#'   Minkowski metric is determined by \code{metric}, with the distance being
+#'   given by
+#'   \deqn{
+#'      D(i,j) = (sum_d |x_{id}-x_{jd}|^p)^(1/p)
+#'   }{
+#'      D(i,j) = \left(\sum_d |x_{id}-x_{jd}|^p\right)^{1/p}
+#'   }
+#'   where the sum is over the dimensions of the space, \eqn{x} is the
+#'   coordinate matrix, and \eqn{p} is the metric parameter (i.e.
+#'   \code{metric}).  The default value of \code{metric} is 2, yielding the
+#'   Euclidean distance.  Note that \code{metric==1}, and \code{log==FALSE}
+#'   in the one-dimensional case is identical to the \code{absdiff} term.
+#'
+#'   When \code{sphere=TRUE}, the first two dimensions of \code{coord} must be
+#'   angular coordinates in lat/lon form (i.e., the first column must contain
+#'   units of decimal degrees between -90 and 90, and the second must contain
+#'   units of decimal degrees betweein -180 and 180).  Distance is then
+#'   computed on between the specified angular coordinates on the surface of a
+#'   sphere of radius \code{radius}; by default, this is the IUGG mean Earth
+#'   radius in km, and hence supplying lat/lon coordinates yields geospherical
+#'   distances in kilometers.  As \code{metric} has no meaning here, it is
+#'   ignored.
+#'
+#'   If \code{sphere=TRUE} and more than two coordinates are provided per
+#'   vertex, then the first two are used to compute great circle distances,
+#'   with the remainder being used to compute the selected Minkowski metric.
+#'   These two distances are then added to form the final distance.  While
+#'   this has exotic uses, the most obvious is to add elevation.  It is
+#'   important to ensure that the coordinates are specified such that the
+#'   two distances have the same units!  (By default, spherical distances are
+#'   in kilometers, so unless \code{radius} is changed, the user should ensure
+#'   that additional coordinates are also in kilometers.)  Note that elevation
+#'   effects specified in this way do not take into account factors such as the
+#'   need to go down to the street and then up to a target floor within another
+#'   building (as in the case of artificial elevation in cities), nor the need
+#'   to e.g. get up before you get down in hilly terrain.  Life is full of such
+#'   disappointments.
+#'
+#'   The overall scale of the computed distances can be modified by setting the
+#'   \code{scale} argument.  This is most often useful for numerical or
+#'   interpretational reasons (e.g., to avoid overflow/underflow issues, or to
+#'   change to familiar units).  However, it can also be used in conjunction
+#'   with offsets and log scaling to change the form of the distance effect,
+#'   as described below.  In particular, note that rescaling is performed prior
+#'   to application of offsets or thresholding, where applicable.  Similarly,
+#'   the \code{pow} argument can be used to raise the scaled distance to an
+#'   arbitrary power prior to further calculation.  (Note, however, that power
+#'   law spatial dependence is not realized via this mechanism, but by working
+#'   with log distances, as discussed below.)
+#'
+#'   When \code{log==TRUE}, the logarithm of the distance rather than the
+#'   distance itself is used to form the statistic.  This produces interaction
+#'   functions that roughly approximate the power law spatial interaction
+#'   functions (SIFs) from Butts and Acton (2011).  To prevent divergence of
+#'   the log distance for points at identical positions, \code{distoff} can be
+#'   used to provide an offset to the raw distance (with the computed value
+#'   becoming \code{log(d+distoff)}.  (This can have an effect that is similar
+#'   to an \dQuote{attenuated power law} SIF, per Butts and Acton.)  As an
+#'   additional failsafe, the (offset) distances are thresholded from below by
+#'   \code{mindist}.  With all elements considered, the net effect is thus
+#'   \code{log(max(mindist, distoff + scale*d))}, where \eqn{d} is the raw
+#'   distance.  Setting both \code{mindist} and \code{distoff} to zero is
+#'   allowed, but may produce exciting results if points precisely overlap.
+#'   Note that both \code{distoff} and \code{mindist} are ignored when
+#'   \code{log==FALSE} (but \code{scale} is not).
+#'
+#'   This term can be used for directed or undirected networks.
+#'
+#' @examples
+#' \donttest{
+#' #Create an example network, in a two-dimensional space
+#' #  Effective SIF is 1/(1 + exp(-3 + 2 log(d))); apx inverse square
+#' n <- 300
+#' d <- 2
+#' x <- matrix(runif(d*n,0,50),ncol=d)
+#' net <- simulate(network.initialize(n,directed=FALSE)~edges
+#'    + distance(x), coef=c(3,-2),
+#'    control=control.simulate.formula(MCMC.burnin=n^3))
+#'
+#' #Examine the network
+#' plot(net,coord=x)  #At first, seems not so spatial....
+#' dis <- as.matrix(dist(x))
+#' plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]),
+#'     xlab="Distance", ylab="Edge State")  #But actually quite spatial!
+#' lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF...
+#'     y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2)
+#' lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))),
+#'     lwd=2, lty=3, col=3)                   #Theoretical SIF
+#'
+#' #Recover the parameters
+#' summary(ergm(net~edges+distance(x)))
+#'
+#' #What happens if we use raw instead of log distances?
+#' #  Effective SIF is 1/(1 + exp(-3 + d)); apx exponential decay
+#' net <- simulate(network.initialize(n,directed=FALSE)~edges
+#'    + distance(x, log=FALSE), coef=c(3,-1),
+#'    control=control.simulate.formula(MCMC.burnin=n^3))
+#' plot(net,coord=x)  #Ties are much more local - no long edges
+#'
+#' plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]),
+#'     xlab="Distance", ylab="Edge State")
+#' lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF...
+#'     y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2)
+#' lines((0:100)/100*max(dis), 1/(1+exp(-3+1*((0:100)/100*max(dis)))),
+#'     lwd=2, lty=3, col=3)                   #New SIF (green)
+#' lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))),
+#'     lwd=2, lty=3, col=4)                   #Previous SIF (blue)
+#'
+#' summary(ergm(net~edges+distance(x, log=FALSE))) #Recover parameters
+#'
+#' #A small lat/lon example; begin with some US cities
+#' co<-rbind(
+#'   AustinTX=c(30.26694,-97.74278),
+#'   BaltimoreMD=c(39.29028,-76.6125),
+#'   BatonRougeLA=c(30.45056,-91.15444),
+#'   CharlotteNC=c(35.22694,-80.84333),
+#'   ChicagoIL=c(41.85,-87.65),
+#'   ColumbusOH=c(39.96111,-82.99889),
+#'   DetroitMI=c(42.33139,-83.04583),
+#'   DurhamNC=c(35.99389,-78.89889),
+#'   HoustonTX=c(29.76306,-95.36306),
+#'   IndianapolisIN=c(39.76833,-86.15806),
+#'   IrvineCA=c(33.66944,-117.8222),
+#'   IthacaNY=c(42.44056,-76.49694),
+#'   JacksonvilleFL=c(30.33194,-81.65583),
+#'   KonaHI=c(19.64083,-155.9856),
+#'   LasVegasNV=c(36.175,-115.1364),
+#'   LosAngelesCA=c(34.05222,-118.2428),
+#'   MemphisTN=c(35.14944,-90.04889),
+#'   MiamiFL=c(25.77389,-80.19389),
+#'   MilwaukeeWI=c(43.03889,-87.90639),
+#'   MobileAL=c(30.69417,-88.04306),
+#'   NewOrleansLA=c(29.95444,-90.075),
+#'   NewYorkNY=c(40.71417,-74.00639),
+#'   PhiladelphiaPA=c(39.95222,-75.16417),
+#'   PhoenixAZ=c(33.44833,-112.0733),
+#'   PittsburghPA=c(40.39556,-79.83889),
+#'   SanAntonioTX=c(29.42389,-98.49333),
+#'   SanDiegoCA=c(32.71528,-117.1564),
+#'   SanFranciscoCA=c(37.775,-122.4183),
+#'   SanJoseCA=c(37.33944,-121.8939),
+#'   SeattleWA=c(47.60639,-122.3308)
+#' )
+#'
+#' #By default, spherical distances give us great circle distances
+#' #on the geosphere, in kilometers
+#' net <- network.initialize(30, directed = FALSE)
+#' net[1,8] <- 1  #Create a tie from Austin, TX to Durham, NC
+#' summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 1863 km
+#' net[1,8] <- 0
+#' net[14,22] <- 1  #Now try Kona, HI to NYC, NY
+#' summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 7940 km
+#'
+#' #Model a network among these fine cities
+#' net <- simulate(net ~ edges + distance(co, sphere = TRUE),
+#'     coef=c(5,-1), control=control.simulate.formula(MCMC.burnin=1e5))
+#'
+#' #Plot the cities in lat/lon space
+#' plot(net, coord=co[,2:1], suppress.axes=FALSE, xlab="Longitude",
+#'     ylab="Latitude")
+#'
+#' #Recover the parameters
+#' summary(ergm(net~edges+distance(co, sphere = TRUE)))
+#'}
+#'
+#' @references Butts, Carter T. and Acton, Ryan M.  (2011).  \dQuote{Spatial
+#'   Modeling of Social Networks.}  In Timothy Nyerges, Helen Couclelis, and
+#'   Robert McMaster (Eds.), \emph{The Sage Handbook of GIS and Society
+#'   Research}, 222--250.  SAGE Publications.
+#'
+#' @template ergmTerm-general
+#'
+#' @concept dyad-independent
+#' @concept directed
+#' @concept undirected
+#' @concept quantitative nodal attribute
+InitErgmTerm.distance <- function(nw, arglist, ...) {
+  a <- check.ErgmTerm(nw, arglist,
+      varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale", "pow"),
+      vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric", "numeric"),
+      required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE),
+      defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1,1))
+  #Process arguments
+  if(is.character(a$coord)){  #Some sort of attribute
+    if(length(a$coord)>1){    #Vertex attributes
+      a$coord<-ergm_get_vattr(a$coord,nw,accept="numeric",multiple="matrix")
+    }else{                    #Either vertex or network attribute
+      if(a$coord%in%list.vertex.attributes(nw))
+        a$coord<-ergm_get_vattr(a$coord,nw,accept="numeric")
+      else
+        a$coord<-get.network.attribute(nw,a$coord)
+    }
+    if(is.null(a$coord))
+      stop("Distance term requires either a coordinate matrix, the name of a network attribute containing one, or one or more names of vertex attributes containing coordinates.")
+  }
+  if(length(dim(a$coord))==2){
+    coord<-a$coord
+    if(is.data.frame(coord))
+      coord<-as.matrix(coord)
+  }else if(length(dim(a$coord))==0){
+    coord<-matrix(rep(a$coord,length=network.size(nw)),ncol=1)
+  }else{
+    coord<-as.matrix(coord)
+  }
+  if(NROW(coord)!=network.size(nw))
+    stop("Distance term requires that coordinates be provided for all vertices.")
+  if(any(!apply(coord,1:2,is.numeric)))
+    stop("Distance term requires numeric coordinates.")
+  if(any(is.na(coord)))
+    stop("Missing coordinate values not allowed in distance term.")
+  if(a$sphere){
+    if(NCOL(coord)<2)
+      stop("Spherical coordinates chosen for distance term, but fewer than two dimensions provided; first two coordinates must be lat,lon or equivalent (in that order, and in angular units).")
+    if(any(coord[,1]< -90)||any(coord[,1]>90))
+      stop("Illegal latitude passed to distance term.")
+    if(any(coord[,2]< -180)||any(coord[,2]>180))
+      stop("Illegal longitude passed to distance term.")
+    if(a$radius<=0)
+      stop("Radius of sphere to be used for distance calculations must be positive.")
+  }
+  if(a$log&&(a$mindist<0))
+    stop("mindist cannot be negative in distance term.")
+  if(a$log&&(a$distoff<0))
+    stop("distoff cannot be negative in distance term.")
+  #Set up the call
+  if(a$log)
+    logstr<-".log"
+  else
+    logstr<-""
+  if(a$sphere)
+    basestr<-"dist.S"
+  else
+    basestr<-"dist"
+  if((!a$sphere)||(NCOL(coord)>2))
+    basestr<-paste0(basestr,".L",a$metric)
+  if(a$pow!=1)
+    basestr<-paste0(basestr,".pow",a$pow)
+  list(name = "distance",
+      coef.names = paste0(basestr,logstr),
+      inputs = c(network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, a$scale, a$pow, coord),
+      dependence = FALSE,
+      emptynwstats = 0
+  )
+}
+
+
 ################################################################################
 
 #' @templateVar name dyadcov

R/InitWtErgmTerm.R

@@ -301,6 +301,21 @@ InitWtErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")
   binary_dind_wrap("diff", nw, a, ..., version=version)
 }
 
+#' @templateVar name distance
+#' @template ergmTerm-rdname
+#' @usage
+#' # valued: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714,
+#' #             log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1, form="sum")
+#' @template ergmTerm-form
+InitWtErgmTerm.distance <- function(nw, arglist, ...) {
+  a <- check.ErgmTerm(nw, arglist,
+      varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale", "pow", "form"),
+      vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric", "numeric", "character"),
+      required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE),
+      defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1,1, "sum"))
+
+  binary_dind_wrap("distance", nw, a, ...)
+}
 #' @templateVar name edgecov
 #' @template ergmTerm-rdname
 #' @usage

inst/NEWS.Rd

@@ -79,6 +79,9 @@
 
   \subsection{NEW FEATURES}{
     \itemize{
+      \item{
+        New term \ergmTerm{ergm}{distance}{()} for Minkowski and great circle distances between vertices now available.
+      }
       \item{
         \ergmTerm{ergm}{nodemix}{()} and \ergmTerm{ergm}{mm}{()} are now more user-friendly in how \code{levels2=} arguments are processed.
       }

inst/include/inc/ergm_changestat_common.do_not_include_directly.h

@@ -13,6 +13,9 @@
 double my_choose(double n, int r);
 #define CHOOSE(n,r) ((n)<(r) ? (0) : (my_choose((double)(n),(int)(r)))) 
 
+/* geodesic distance on a sphere */
+double spheredist(double lat0, double lon0, double lat1, double lon1, double r);
+
 /* Comparison macro for doubles: */
 #define EQUAL(a,b) (fabs((a)-(b))<0.0000001)