diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index 5d927d602..23784ebdd 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -25,6 +25,7 @@ # C: = # D: # +# # E: # G: # @@ -2412,6 +2413,296 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) ) } +################################################################################ + +#' @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 diff --git a/R/InitWtErgmTerm.R b/R/InitWtErgmTerm.R index df2666d70..3878b7f05 100644 --- a/R/InitWtErgmTerm.R +++ b/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 diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 7210c5bba..c3f711e37 100644 --- a/inst/NEWS.Rd +++ b/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. } diff --git a/inst/include/inc/ergm_changestat_common.do_not_include_directly.h b/inst/include/inc/ergm_changestat_common.do_not_include_directly.h index 83442a489..74ee3db2e 100644 --- a/inst/include/inc/ergm_changestat_common.do_not_include_directly.h +++ b/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) diff --git a/man/distance-ergmTerm-d5141986.Rd b/man/distance-ergmTerm-d5141986.Rd new file mode 100644 index 000000000..2fc18b7d0 --- /dev/null +++ b/man/distance-ergmTerm-d5141986.Rd @@ -0,0 +1,250 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/InitErgmTerm.R, R/InitWtErgmTerm.R +\name{distance-ergmTerm} +\alias{distance-ergmTerm} +\alias{InitErgmTerm.distance} +\alias{InitWtErgmTerm.distance} +\title{Inter-point distances} +\usage{ +# binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +# log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1) + +# valued: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +# log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1, form="sum") +} +\arguments{ +\item{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} + +\item{metric}{power to use for the Minkowski metric} + +\item{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} + +\item{radius}{for spherical distances, the radius of the sphere to use; +defaults to the IUGG mean Earth radius, in km} + +\item{log}{logical; use log rather than raw distances?} + +\item{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} + +\item{distoff}{for log distances, an offset to be added to observed +distances before taking the logarithm} + +\item{scale}{factor by which raw distances should be rescaled prior to other +operations (notably, offsetting, logging, and thresholding)} + +\item{pow}{power to which scaled distances should be raised prior to other +other operations (notably, offsetting, logging, and thresholding)} + +\item{form}{how to aggregate tie values in a valued ERGM: \code{"sum"} +(the default) for a statistic of the form \eqn{\sum_{i,j} x_{i,j} +y_{i,j}}{sum[i,j] x[i,j]*y[i,j]}, where \eqn{y_{i,j}}{y[i,j]} is +the value of dyad \eqn{(i,j)} and \eqn{x_{i,j}}{x[i,j]} is the +term's covariate associated with it; and \code{"nonzero"} with the +edge considered to be present if its value is not 0. See +\code{\link{ergmTerm}} for more information.} +} +\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. +} +\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. +} +\seealso{ +\code{\link{ergmTerm}} for index of model terms currently visible to the package. + +\Sexpr[results=rd,stage=render]{ergm:::.formatTermKeywords("ergmTerm", "distance", "subsection")} +} +\concept{directed} +\concept{dyad-independent} +\concept{quantitative nodal attribute} +\concept{undirected} diff --git a/src/changestat.c b/src/changestat.c index daaae9db7..e85867028 100644 --- a/src/changestat.c +++ b/src/changestat.c @@ -31,3 +31,26 @@ double my_choose(double n, int r) { ans*=(n--); return ans; } + +/* +Utility function to calculate distances on a sphere of radius r, between two points +in spherical coordinates. Our notation betrays our assumption that these are +lat/lon points on the geosphere, but technically you could use any sphere you +wanted, by passing an alternative radius (so long as your coordinates were given +in angular units). +*/ +double spheredist(double lat0, double lon0, double lat1, double lon1, double r) { + double rlat0,rlat1,rlon0,rlon1,cosrl0,cosrl1,cosdl,sinrl0,sinrl1,sindl; + + rlat0=lat0/180.0*M_PI; + rlat1=lat1/180.0*M_PI; + rlon0=lon0/180.0*M_PI; + rlon1=lon1/180.0*M_PI; + cosrl0=cos(rlat0); + cosrl1=cos(rlat1); + cosdl=cos(rlon0-rlon1); + sinrl0=sin(rlat0); + sinrl1=sin(rlat1); + sindl=sin(rlon0-rlon1); + return r*atan2(sqrt((cosrl0*sindl) * (cosrl0*sindl) + (cosrl1*sinrl0-sinrl1*cosrl0*cosdl) * (cosrl1*sinrl0-sinrl1*cosrl0*cosdl)), sinrl0*sinrl1+cosrl0*cosrl1*cosdl); +} diff --git a/src/changestats.c b/src/changestats.c index ee0ba70d8..869884ca5 100644 --- a/src/changestats.c +++ b/src/changestats.c @@ -1463,6 +1463,7 @@ C_CHANGESTAT_FN(c_degree_w_homophily) { } } + /***************** changestat: d_dyadcov *****************/ diff --git a/src/changestats_dyad_ind.c.template.do_not_include_directly.h b/src/changestats_dyad_ind.c.template.do_not_include_directly.h index 8abbdc218..078991b85 100644 --- a/src/changestats_dyad_ind.c.template.do_not_include_directly.h +++ b/src/changestats_dyad_ind.c.template.do_not_include_directly.h @@ -435,3 +435,66 @@ ETYPE(C_CHANGESTAT_FN)(SVARIANT(c_sociality)) { } } + +/***************** + changestat: c_distance +*****************/ +ETYPE(C_CHANGESTAT_FN)(SVARIANT(c_distance)) { + Vertex t,h; + int j; + int nv,dim,logd,sphd; + double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis,pw; + + /*Set things up*/ + nv = (int)INPUT_PARAM[0]; /*Number of vertices*/ + dim = (int)INPUT_PARAM[1]; /*Dimension of coordinate space*/ + dexp = INPUT_PARAM[2]; /*Minkowski exponent*/ + logd = (int)INPUT_PARAM[3]; /*Should we use the log distance?*/ + sphd = (int)INPUT_PARAM[4]; /*Should we use lat/lon distances on the geosphere?*/ + sphr = INPUT_PARAM[5]; /*Radius to use for spherical coordinates*/ + logmind = INPUT_PARAM[6]; /*Distance minimum for the log case*/ + logdoff = INPUT_PARAM[7]; /*Distance offset for the log case*/ + scale = INPUT_PARAM[8]; /*Scaling adjustment for raw distance*/ + pw = INPUT_PARAM[9]; /*Power to which distance should be raised*/ + coord = INPUT_PARAM+10; /*Pointer to the coordinate matrix*/ + + /*Compute the changescore*/ + t=tail-1; /*0-indexed values, for C convenience*/ + h=head-1; + if(sphd){ /*Compute great circle distances on the sphere*/ + dis=spheredist(coord[t],coord[t+nv],coord[h],coord[h+nv],sphr); + }else{ + dis=0.0; + } + /*Compute Minkowski distances in free space (adding extra dimensions to + great circle distance if sphd and we were given dim>2)*/ + mdis=0.0; + for(j=0+2*sphd;j0]), ignore_attr=TRUE) +}) +test_that("L1 works", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE)), sum(d1[g>0]), ignore_attr=TRUE) +}) +test_that("L2 works", { + expect_equal(summary(net~distance(co)), sum(log(d2[g>0])), ignore_attr=TRUE) +}) +test_that("log L1 works", { + expect_equal(summary(net~distance(co,metric=1)), sum(log(d1[g>0])), ignore_attr=TRUE) +}) +test_that("scaling works", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE,scale=3)), 3*sum(d1[g>0]), ignore_attr=TRUE) +}) +test_that("powers work", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE,pow=3)), sum(d1[g>0]^3), ignore_attr=TRUE) +}) +test_that("thresholding works", { + expect_equal(summary(net~distance(co,metric=1,log=TRUE,scale=1e-10,mindist=5)), log(5)*sum(g>0), ignore_attr=TRUE) +}) + +#Verify that alternative specifications work +test_that("character matches matrix", { + expect_equal(summary(net~distance(co)), summary(net~distance("co")), ignore_attr=TRUE) +}) +test_that("character matches data.frame", { + expect_equal(summary(net~distance(as.data.frame(co))), summary(net~distance("co")), ignore_attr=TRUE) +}) +test_that("vertex attributes matches matrix", { + expect_equal(summary(net~distance(co)), summary(net~distance(c("x","y","z"))), ignore_attr=TRUE) +}) +test_that("single attribute defaults to vertex", { + expect_equal(summary(net~distance(co[,1])), summary(net~distance("x")), ignore_attr=TRUE) +}) + + +#Test for lat/lon calculations +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 +test_that("first great circle distance is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 1862.879, ignore_attr=TRUE) +}) +net[1,8] <- 0 +net[14,22] <- 1 #Now try Kona, HI to NYC, NY +test_that("second great circle distance is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 7940.019, ignore_attr=TRUE) +}) +net[1,8] <- 1 +test_that("great circle distance sum is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 9802.898, ignore_attr=TRUE) +}) +test_that("log great circle distance sum is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = TRUE)),3), 16.51, ignore_attr=TRUE) +}) diff --git a/tests/testthat/test-ergm-term-doc.R b/tests/testthat/test-ergm-term-doc.R index b62c354a1..5677a5618 100644 --- a/tests/testthat/test-ergm-term-doc.R +++ b/tests/testthat/test-ergm-term-doc.R @@ -44,7 +44,7 @@ test_that("test search ergm term", { expect_equal(length(search.ergmTerms(keywords = 'bipartite', packages='ergm')), 36) ## expect_gt(length(search.ergmTerms(keywords = 'valued')), 44) - expect_equal(length(search.ergmTerms(keywords = 'valued', packages='ergm')), 50) + expect_equal(length(search.ergmTerms(keywords = 'valued', packages='ergm')), 51) ## expect_gt(length(search.ergmTerms(keywords = 'valued', packages=c('ergm', 'ergm.count'))), 44) })