minor speed-up to logistic code when distance is a covariate

NEWS

@@ -3,6 +3,7 @@ mrds 2.1.14
 
  * updated initialvalues calculation for hazard-rate -- now uses Beavers & Ramsay method to scale parameters for hazard-rate
  * automatic parameter rescaling for covariate models when covariates are poorly scaled. Now default for nlminb method
+ * minor speed-up to logistic code when distance is a covariate
 
 
 mrds 2.1.13

R/g0.R

@@ -6,5 +6,5 @@
 #' @return vector of p(0) values
 #' @author Jeff Laake
 g0 <- function(beta, z){
-  exp(as.matrix(z) %*% beta)/(1 + exp(as.matrix(z) %*% beta))
+  exp(z %*% beta)/(1 + exp(z %*% beta))
 }

R/integratelogisticdup.R

@@ -1,7 +1,27 @@
 integratelogisticdup <- function(x1, x2, models, beta, lower=0, width, point){
   # numerical integral of product of logistic detection functions
 
-  integrate(logisticdupbyx, lower=lower, upper=width,
-            subdivisions=10, rel.tol=0.01, abs.tol=0.01,
-            x1=x1, x2=x2, models=models, beta=beta, point=point)$value
+  # computation speed-up when there is distance in the formula
+  # but only if there are no interactions with distance
+  if(sum(grepl("distance",names(beta)))==1){
+
+    # set parameter to be zero for distance
+    beta_distance <- beta[grepl("distance",names(beta))]
+    beta[grepl("distance",names(beta))] <- 0
+
+    # calculate the rest of the linear predictor
+    x1 <- setcov(x1,models$g0model)%*%beta
+    x2 <- setcov(x2,models$g0model)%*%beta
+
+    # do some integration
+    integrate(logisticdupbyx_fast, lower=lower, upper=width,
+              subdivisions=10, rel.tol=0.01, abs.tol=0.01,
+              x1=x1, x2=x2, models=models, beta=beta, point=point,
+              beta_distance=beta_distance)$value
+  }else{
+  # Otherwise just go ahead and do the numerical integration
+    integrate(logisticdupbyx, lower=lower, upper=width,
+              subdivisions=10, rel.tol=0.01, abs.tol=0.01,
+              x1=x1, x2=x2, models=models, beta=beta, point=point)$value
+  }
 }

R/logisticdupbyx.R

@@ -14,20 +14,24 @@
 #' @return vector of probabilities
 #' @author Jeff Laake
 logisticdupbyx <- function(distance, x1, x2, models, beta, point){
+
+  # avoid using g0 which calls exp and matrix multiplication twice
+  ologit <- function(p) p/(1+p)
+
   #  Functions used: g0, setcov
   xlist <- as.list(x1)
   xlist$distance <- distance
   xmat <- expand.grid(xlist)
 
-  gx1 <- g0(beta, setcov(xmat, models$g0model))
+  gx1 <- ologit(exp(setcov(xmat, models$g0model) %*% beta))
 
   xlist <- as.list(x2)
   xlist$distance <- distance
   xmat <- expand.grid(xlist)
 
   if(!point){
-    return(gx1 * g0(beta, setcov(xmat, models$g0model)))
+    return(gx1 * ologit(exp(setcov(xmat, models$g0model) %*% beta)))
   }else{
-    return(gx1 * g0(beta, setcov(xmat, models$g0model))*2*distance)
+    return(gx1 * ologit(exp(setcov(xmat, models$g0model) %*% beta))*2*distance)
   }
 }

R/logisticdupbyx_fast.R

@@ -0,0 +1,24 @@
+#' Logistic for duplicates as a function of covariates (fast)
+#'
+#' As \code{\link{logisticdupbyx}}, but faster when distance is a covariate (but no interactions with distance occur.
+#'
+#' @inheritParams logisticdupbyx
+#' @param beta_distance parameter for distance
+#' @param x1 linear predictor for 1, without distance
+#' @param x2 linear predictor for 2, without distance
+#' @author David L Miller
+logisticdupbyx_fast <- function(distance, x1, x2, models, beta, point, beta_distance){
+
+  # function to calculate p/(1+p)
+  ologit <- function(p) p/(1+p)
+
+  # first part of the function
+  gx1 <- ologit(exp(x1 + distance*beta_distance))
+
+  # calculate second and return
+  if(!point){
+    return(gx1 * ologit(exp(x2 + distance*beta_distance)))
+  }else{
+    return(gx1 * ologit(exp(x2 + distance*beta_distance))*2*distance)
+  }
+}

R/predict.io.fi.R

@@ -58,19 +58,19 @@ predict.io.fi <- function(object,newdata=NULL,compute=FALSE, int.range=NULL,
 
     # now int.range is a vector with lower and upper bounds
     if(is.null(int.range)){
-      pdot.list <- pdot.dsr.integrate.logistic(width,width, model$mr$coef,
-                     newdata,integral.numeric, FALSE, models,GAM, point=point)
+      pdot.list <- pdot.dsr.integrate.logistic(width, width, model$mr$coef,
+                     newdata, integral.numeric, FALSE, models,GAM, point=point)
     }else{
       pdot.list <- pdot.dsr.integrate.logistic(int.range,width, model$mr$coef,
-                       newdata,integral.numeric, FALSE, models,GAM, point=point)
+                      newdata, integral.numeric, FALSE, models,GAM, point=point)
     }
 
     # if there is left truncation, take that off the integral
     if(left !=0){
       pdot.list$pdot <- pdot.list$pdot -
                     pdot.dsr.integrate.logistic(left, width, model$mr$coef,
-                                   newdata,integral.numeric, FALSE, models,GAM,
-                                   point=point)$pdot
+                                   newdata, integral.numeric, FALSE, models,
+                                   GAM, point=point)$pdot
     }
 
     fitted <- pdot.list$pdot

man/logisticdupbyx_fast.Rd

@@ -0,0 +1,30 @@
+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/logisticdupbyx_fast.R
+\name{logisticdupbyx_fast}
+\alias{logisticdupbyx_fast}
+\title{Logistic for duplicates as a function of covariates (fast)}
+\usage{
+logisticdupbyx_fast(distance, x1, x2, models, beta, point, beta_distance)
+}
+\arguments{
+\item{distance}{vector of distance values}
+
+\item{x1}{linear predictor for 1, without distance}
+
+\item{x2}{linear predictor for 2, without distance}
+
+\item{models}{model list}
+
+\item{beta}{logistic parameters}
+
+\item{point}{\code{TRUE} for point transect data}
+
+\item{beta_distance}{parameter for distance}
+}
+\description{
+As \code{\link{logisticdupbyx}}, but faster when distance is a covariate (but no interactions with distance occur.
+}
+\author{
+David L Miller
+}
+