# R/network.R In geonet: Intensity Estimation on Geometric Networks with Penalized Splines

#### Defines functions network_locationnetwork_ISEnetwork_integralnetwork_intensityinternalbin_datanetwork_binsnetwork_knotsdelta_h_global

```#' Computes a Global Knot Distance from the Input
#'
#' @param G An object of class \code{gn}) or a point pattern on a geometric
#' network (object of class \code{gnpp}).
#' @param delta A numeric vector of length one which already defines
#' the global knot distance. Alternatively, delta can be supplied in terms of a
#' quantile of the curve lengths of the network, i.e. a number in the unit
#' interval. In the latter case, delta must be supplied as a character vector of
#' length one, see the examples. By default, delta is chosen to be half of the
#' minimal curve length.
#' @param h A numeric vector of length one which already defines the
#' global knot distance h. Alternatively, h can be supplied in terms of a
#' fraction of delta, i.e. a number in the interval (0,1]. In the latter case,
#' h must be supplied as character vector of length one, see the examples.
#' By default, h is chosen to be half of the global knot distance.
#' @return The global knot distance delta.
#' @author Marc Schneble \email{marc.schneble@@stat.uni-muenchen.de}
#' @export
#' @examples
#' G <- as_gn(montgomery)
#' # use default arguments
#' setup <- delta_h_global(G)
#' setup
#' # set numeric value for delta and fraction for h
#' setup <- delta_h_global(G, delta = 0.1, h = "0.25")
#' setup
#' # set quantile for delta
#' setup <- delta_h_global(G, delta = "0.05")
#' setup

delta_h_global <- function(G, delta = NULL, h = NULL) {
#global knot distance
stopifnot(inherits(G, "gn"))
if (is.null(delta)) {
delta <- min(G\$d)/2
} else if (is.character(delta) & !is.na(as.numeric(delta))) {
if (is.character(delta) & !is.na(as.numeric(delta))) {
quant <- as.numeric(delta)
if (quant >= 0 & quant <= 1) {
delta <- as.numeric(quantile(G\$d, quant)/2)
if (delta > quantile(G\$d, 0.5)/2) {
warning(paste0("Many network segments have no curve specific ",
"B-splines.\nChoose a lower value for the quantile."))
}
} else {
stop("The quantile supplied for delta must be in the interval [0,1]")
}
} else {
stop("The character must be convertible to class numeric.")
}
} else if (is.numeric(delta)) {
if (delta <= 0) {
stop("delta must be positive!")
}
if (delta > quantile(G\$d, 0.5)/2) {
warning(paste0("Many network segments have no curve specific ",
"B-splines.\nChoose a lower value for delta."))
}
} else {
stop("Supply delta in correct format, see the documentation.")
}
if (is.null(h)) {
h <- delta/2
} else if (is.numeric(h)) {
if (h <= 0) stop("h must be positive!")
if (h > delta) stop("h must be at most as large as delta!")
} else if (is.character(h) & !is.na(as.numeric(h))) {
mult <- as.numeric(h)
h <- delta*mult
} else {
stop("Supply h in correct format, see the documentation.")
}
list(delta = delta, h = h)
}

#' Defining knots on a Geometric Network
#'
#' \code{network_knots} defines knots on a geometric network (object of class
#' \code{gn}) which can be used to construct linear B-splines on it.
#'
#' @param G An object of class \code{gn}) or a point pattern on a geometric
#' network (object of class \code{gnpp}).
#' @param delta The global knot distance delta.
#' @return A list which contains the knot sequence of every curve of the
#' geometric network.
#' @author Marc Schneble \email{marc.schneble@@stat.uni-muenchen.de}

network_knots <- function(G, delta){
if (inherits(G, "gnpp")) G <- as_gn(G)
if (!inherits(G, "gn")) stop("G muss be of class 'gn'")
# line specific knot distances
delta <- G\$d*(delta > G\$d) + delta*(delta <= G\$d)
delta <- pmin(G\$d/floor(G\$d/delta)*(G\$d/delta - floor(G\$d/delta) < 0.5) +
G\$d/ceiling(G\$d/delta)*(G\$d/delta - floor(G\$d/delta) >= 0.5),
G\$d/2)
# initializing...
tau <- vector("list", G\$M)
J <- rep(0, G\$M)

# do for every line segment
for (m in 1:G\$M) {
# knot sequences tau
tau[[m]] <- seq(0, G\$d[m], delta[m])
# count of linear B-splines on line segment
J[m] <- length(tau[[m]]) - 2
}
knots <- list(delta = delta, tau = tau, J = J)
knots
}

#' Defining bins on a Geometric Network
#'
#' \code{network_bins} subdivides each curve segment into several bins.
#'
#' @param G An object of class \code{gn}) or a point pattern on a geometric
#' network (object of class \code{gnpp}).
#' @param h The global bin width h.
#' @return A list which contains the bins of every curve of the geometric
#' network.
#' @author Marc Schneble \email{marc.schneble@@stat.uni-muenchen.de}

network_bins <- function(G, h = NULL){
if (inherits(G, "gnpp")) G <- as_gn(G)
if (!inherits(G, "gn")) stop("G muss be of class 'gn'")
# line specific bin widths
h <- pmin(G\$d/floor(G\$d/h)*(G\$d/h - floor(G\$d/h) < 0.5) +
G\$d/ceiling(G\$d/h)*(G\$d/h - floor(G\$d/h) >= 0.5), G\$d/2)
# initializing...
b <- z <- vector("list", G\$M)
N <- rep(0, G\$M)

# do for every line segment
for (m in 1:G\$M) {
# bin boundaries b
b[[m]] <- seq(0, G\$d[m], h[m])
# characterization of bins by midpoints z
z[[m]] <- (b[[m]][1:(length(b[[m]])-1)] + b[[m]][2:length(b[[m]])])/2
# total count of bins in the geometric network
N[m] <- length(z[[m]])
}
bins <- list(h = h, b = b, z = z, N = N)
bins
}

#' Bin Point Pattern on a Geometric Network
#'
#' \code{bin_data} bins the data on the supplied point pattern according to
#' all possible combination of covariates.
#'
#' @param X Point pattern on a geometric network (object of class \code{gnpp})
#' @param bins A list containing the bins of the geometric network.
#' @param vars A character vector containing the name of all covariates in the
#' model.
#' @param vars_internal A character vector containing the name of all interval
#' covariates in the model.
#' @param scale A named list which specifies the rescaling of network related
#' covariates. Currently, internal covariates "x", "y", and "dist2V" can be
#' scaled.
#' @return The binned data.
#' @import dplyr
#' @importFrom stats setNames
#' @author Marc Schneble \email{marc.schneble@@stat.uni-muenchen.de}

bin_data <- function(X, bins = NULL, vars = NULL, vars_internal = NULL,
scale = NULL){

# get covariates from every point on the network (if applicable)
if (is.null(bins)) bins <- network_bins(X\$network)
vars_external <- setdiff(vars, vars_internal)
if (length(vars_external) > 0) {
covariates <- as_tibble(X\$data) %>% select(all_of(vars_external))
} else {
covariates <- tibble(a = rep(1, nrow(X\$data)))
}

# get all combinations of covariates and calculate the number of rows of the data matrix
vars_comb <- covariates %>% distinct() %>% expand.grid() %>%
distinct() %>% as_tibble()

# sort data frame
for (a in 1:length(vars_comb)) {
vars_comb <- arrange(vars_comb, !!sym(names(vars_comb[a])))
}

# initializing
data <- tibble(id = 1:sum(bins\$N), count = NA, h = NA)
if (length(vars_internal) > 0) {
data <- bind_cols(data, as_tibble(internal(vars_internal, X, bins, scale)))
}
data <- data %>% slice(rep(1:n(), nrow(vars_comb))) %>%
bind_cols(vars_comb %>% slice(rep(1:n(), each = sum(bins\$N))))

ind <- 1
for (j in 1:nrow(vars_comb)) {
ind_comb <- which(do.call(paste, covariates) == do.call(paste, vars_comb[j, ]))
data_sub <- X\$data[ind_comb, ]
for (m in 1:X\$network\$M) {
# positions of data on curve e
ind_e <- which(data_sub\$e == m)
y_e <- sort(as.numeric(data_sub[ind_e, ]\$tp_e))*X\$network\$d[m]

# bin data
y_b <- rep(0, length(bins\$z[[m]]))
for (k in 1:length(bins\$z[[m]])) {
y_b[k] <- length(which(y_e < bins\$b[[m]][k+1] & y_e > bins\$b[[m]][k]))
}

# stack into one vector
data\$count[ind:(ind + length(y_b) - 1)] <- y_b
data\$h[ind:(ind + length(y_b) - 1)] <- bins\$h[m]
ind <- ind + length(y_b)
}
# add bin id for every row
data\$id[((j-1)*sum(bins\$N) + 1):(j*sum(bins\$N))] <- 1:sum(bins\$N)
}
data\$offset <- 1
data
}

#' Internal Covariates
#'
#' \code{internal} computes the values of internal covariates at the midpoints
#' of the bins of the network. Internal covariates can either be supplied via
#' the point pattern or they are a function of the network. Currently, x- and y-
#' coordinates are supported for the latter.
#'
#' @param vars The name of the covariates which should go into the model
#' as linear internal covariates.
#' @param X Point pattern on a geometric network (object of class \code{gnpp})
#' @param bins A list containing the bins of the geometric network.
#' @param scale A named list which specifies the rescaling of network related
#' covariates. Currently, only x- and y-coordinates can be scaled.
#' @return A data frame with the number of rows equal to the number of bins of the
#' geometric network (\code{sum(bins\$N)}) and the number of columns equal to
#' the length of \code{vars}.
#' @author Marc Schneble \email{marc.schneble@@stat.uni-muenchen.de}

internal <- function(vars, X, bins, scale){
e <- NULL
out <- list()
if (is.element("x", vars)) {
x <- rep(NA, sum(bins\$N))
ind <- 1
test <- rep(NA, X\$network\$M)
for (m in 1:X\$network\$M) {
lins_m <- filter(X\$network\$lins, e == m)
tp <- bins\$z[[m]]/X\$network\$d[m]
dx <- lins_m\$a2_x - lins_m\$a1_x
for (i in 1:nrow(lins_m)) {
ind_id <- which(tp >= lins_m\$frac1[i] & tp < sum(lins_m\$frac2[1:i]))
x[ind:(ind + length(ind_id) - 1)] <- lins_m\$a1_x[i] + tp[ind_id]*dx[i]
ind <- ind + length(ind_id)
}
test[m] <- ind - 1
}
if (!is.null(scale\$x)) out\$x <- x*scale\$x
else out\$x <- x
}
if (is.element("y", vars)) {
y <- rep(NA, sum(bins\$N))
ind <- 1
for (m in 1:X\$network\$M) {
lins_m <- filter(X\$network\$lins, e == m)
tp <- bins\$z[[m]]/X\$network\$d[m]
dy <- lins_m\$a2_y - lins_m\$a1_y
for (i in 1:nrow(lins_m)) {
ind_id <- which(tp >= lins_m\$frac1[i] & tp < sum(lins_m\$frac2[1:i]))
y[ind:(ind + length(ind_id) - 1)] <- lins_m\$a1_y[i] + tp[ind_id]*dy[i]
ind <- ind + length(ind_id)
}
}
if (!is.null(scale\$y)) out\$y <- y*scale\$y
else out\$y <- y
}
if (is.element("dist2V", vars)) {
d <- rep(NA, sum(bins\$N))
ind <- 1
for (m in 1:X\$network\$M) {
d[ind:(ind + bins\$N[m] - 1)] <- pmin(bins\$z[[m]], X\$network\$d[m] - bins\$z[[m]])
ind <- ind + bins\$N[m]
}
if (!is.null(scale\$dist2V)) out\$dist2V <- d*scale\$dist2V
else out\$dist2V <- d
}
X_internal <- colnames(X\$network\$lins[-(1:11)])
if (length(X_internal) > 0) {
for (k in 1:length(X_internal)) {
if (is.element(X_internal[k], vars)) {
if (inherits(X\$network\$lins[[X_internal[k]]], "factor")) {
val <- factor(rep(NA, sum(bins\$N)),
levels = levels(X\$network\$lins[[X_internal[k]]]))
} else {
val <- rep(NA, sum(bins\$N))
}
ind <- 1
for (m in 1:X\$network\$M) {
lins_m <- filter(X\$network\$lins, e == m)
if (length(unique(lins_m[[paste(X_internal[k])]])) > 1) {
#stop("Internal covariates must be unique on each curve.")
}
val[ind:(ind + bins\$N[m] - 1)] <- lins_m[[paste(X_internal[k])]][1]
ind <- ind + bins\$N[m]
}
out[[paste(X_internal[k])]] <- val
}
}
}
out
}

#' Fitted Intensity on a Geometric Network
#'
#' @param z The shortest path distance from the beginning of the network
#' segment.
#' @param m The network segment index.
#' @param fit1 A fitted geometric network.
#' @param fit2 A second fitted geometric network. If specified, the function
#' returns the squared difference of the intensity fits at the specified point
#' of the network.
#' @param scale A numeric vector of length two which determines the scaling
#' of the two intensity functions.
#'
#' @return A numeric vector of length one, indicating the intensity (or the
#' squared difference of two intensities) at the specified point.

network_intensity <- function(z, m, fit1, fit2 = NULL, scale = NULL){

ind_v1 <- which(fit1\$network\$incidence[, m] == -1)
ind_v2 <- which(fit1\$network\$incidence[, m] == 1)

B1 <- matrix(0, length(z), fit1\$knots\$J[m] + 2)
B1[, 1] <- (1-z/fit1\$knots\$delta[m])*(1-z/fit1\$knots\$delta[m] > 0)
B1[, 2:(ncol(B1)-1)] <- splineDesign(knots = fit1\$knots\$tau[[m]], x = z, ord = 2, outer.ok = TRUE)
B1[, ncol(B1)] <- (1-(fit1\$network\$d[m]-z)/fit1\$knots\$delta[m])*(1-(fit1\$network\$d[m]-z)/fit1\$knots\$delta[m] > 0)

gamma1 <- fit1\$coefficients[fit1\$ind[[1]]]

gamma1_m <- c(gamma1[sum(fit1\$knots\$J) + ind_v1],
gamma1[((cumsum(fit1\$knots\$J)-fit1\$knots\$J)[m]+1):cumsum(fit1\$knots\$J)[m]],
gamma1[sum(fit1\$knots\$J) + ind_v2])

if (!is.null(fit2)) {
B2 <- matrix(0, length(z), fit2\$knots\$J[m] + 2)
B2[, 1] <- (1-z/fit2\$knots\$delta[m])*(1-z/fit2\$knots\$delta[m] > 0)
B2[, 2:(ncol(B2)-1)] <- splineDesign(knots = fit2\$knots\$tau[[m]], x = z, ord = 2, outer.ok = TRUE)
B2[, ncol(B2)] <- (1-(fit2\$network\$d[m]-z)/fit2\$knots\$delta[m])*(1-(fit2\$network\$d[m]-z)/fit2\$knots\$delta[m] > 0)

gamma2 <- fit2\$coefficients[fit2\$ind[[1]]]

gamma2_m <- c(gamma2[sum(fit2\$knots\$J) + ind_v1],
gamma2[((cumsum(fit2\$knots\$J)-fit2\$knots\$J)[m]+1):cumsum(fit2\$knots\$J)[m]],
gamma2[sum(fit2\$knots\$J) + ind_v2])
}

if (min(z) >= 0 & max(z) <= fit1\$network\$d[m]) {
if (!is.null(fit2)) {
if (is.null(scale)) scale = c(1, 1)
return(as.vector((exp(B1%*%gamma1_m)/scale[1] - exp(B2%*%gamma2_m)/scale[2])^2))
} else {
return(as.vector(exp(B1%*%gamma1_m)))
}
} else {
stop("Wrong domain of z on this edge!")
}
}

#' Integral of a fitted intensity
#'
#' @param fit A fitted point process on a geometric network.
#'
#' @return A numeric vector of the length one, the integral.

network_integral <- function(fit) {
int <- rep(NA, fit\$network\$M)
for (m in 1:fit\$network\$M) {
int[m] <- integrate(network_intensity, lower = 0, upper = fit\$network\$d[m],
m = m, fit1 = fit)\$value
}
sum(int)
}

#' Computation of the Integrated Squared Error
#'
#' Computes the integrated squared error between a true between a true
#' intensity \code{fit1} and an estimate \code{fit2}.
#'
#' @param fit1 The true intensity.
#' @param fit2 The estimated intensity.
#'
#' @return The (normalized) integrated squared error.
#' @export

network_ISE <- function(fit1, fit2) {

scale <- c(nrow(fit1\$data), network_integral(fit2))
int <- rep(NA, fit1\$network\$M)
for (m in 1:fit1\$network\$M) {
int[m] <- integrate(network_intensity, lower = 0, upper = fit1\$network\$d[m],
m = m, fit1 = fit1, fit2 = fit2, scale = scale)\$value
}
sum(int)
}

#' Find Location of a Point on a Geometric Network
#'
#' @param G A geometric network.
#' @param m The segment index.
#' @param z The shortest path distance from the beginning of the network
#' segment.
#'
#' @return A list with names l, tp_l, x and y.

network_location <- function(G, m, z){
e <- NULL
lins_m <- filter(G\$lins, e == m)
frac <- c(lins_m\$frac1, 1)
tp_e <- z/G\$d[m]
lin <- findInterval(tp_e, frac, rightmost.closed = TRUE, all.inside = TRUE)
l <- lins_m\$l[lin]
dx <- lins_m\$a2_x[lin] - lins_m\$a1_x[lin]
dy <- lins_m\$a2_y[lin] - lins_m\$a1_y[lin]
tp_l <- (z - lins_m\$frac1[lin]*G\$d[m])/(lins_m\$length[lin])
x <- lins_m\$a1_x[lin] + tp_l*dx
y <- lins_m\$a1_y[lin] + tp_l*dy
list(l = l, tp_l = tp_l, x = x, y = y)
}
```

## Try the geonet package in your browser

Any scripts or data that you put into this service are public.

geonet documentation built on July 11, 2022, 9:08 a.m.