Non-euclidean movement cost parameters

[Robby Marrotte]

Hello,

I was wondering whether there is any way to add in some extra parameters that I would like estimated in the userdist function?

As a simple example, using simply the mask I would like to estimate the cost of traversing an NA cell. I would like to do the following in my userdist function:

userdfn2 <- function (xy1, xy2, mask) {

  if (missing(xy1)) return("noneuc","c")

  require(gdistance)

  Sraster <- raster(mask, "noneuc")

  

   # Assume animals can traverse gaps: bridge gaps using an estimated value

  Sraster[is.na(Sraster[])] <- c

  

  # conductance is inverse of friction

  trans <- transition(Sraster, transitionFunction = function(x) 1/mean(x),

                      directions = 16)

  trans <- geoCorrection(trans)

  costDistance(trans, as.matrix(xy1), as.matrix(xy2))

}

I am would like to try to emulate landscape resistance optimization as in the ResistanceGA R package but for SECR models.

The ResistanceGA package uses a gravity model and directly estimates the cost of each pixel using genetic distance between individuals or populations and landscape covariates in raster format.

For a continuous surface such as a DEM, the cost of elevation is estimated from 3 parameters that transform the values of the DEM. 

However, for a binary surface such as the one above, the only parameter to estimate is the cost value of the most costly feature (the NA areas).

Is this at all possible right now in SECR?

Thanks!

Robby

https://github.com/wpeterman/ResistanceGA

https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12984

[Murray Efford]

Interesting, but I'm not sure I have the picture.

In secr-noneuclidean.pdf I read "noneuc’ has its own link scale (default ‘log’) on which it may be modelled as a linear function of anyof the predictors available for density (x, y, x2, y2, xy, session, Session, g, or any mask covariate – see secr-densitysurfaces.pdf)... When the model is fitted by secr.fit, the beta parameters for the ‘noneuc’submodel are estimated along with all the others. To make noneuc available to userdist, ensure that it appears in the ‘model’ argument. Use the formula noneuc ~ 1 if noneuc is constant."

So the coefficients of the linear model for noneuc are 'extra parameters' that might conceivably do the job. I take it you would like to escape the linear model formalism and estimate 3 parameters (coefficients) that are unmodelled (constant across space, time etc.). There is an undocumented feature in secr.fit that allows up to 4 unmodelled parameters, but these are presently not passed to userdist. This could be fixed.

Murray

[Murray Efford]

I was wrong - miscparm is documented and available. See Examples in ?details.

Murray

On Wednesday, July 3, 2019 at 4:23:42 AM UTC+12, Robby Marrotte wrote:

[Robby Marrotte]

Thank you for your quick reply. I'll check it out and see if it does the job!

[Robby Marrotte]

Hello again,

This is great I got it working! One small question though, is it possible to somehow make the user function use a different mask than the one used in the rest of the model?

I would like to do something like:

non_euclidean <- secr.fit(eh,model=list(g0~1, sigma~1), hab_msk, trace=F, CL=T, details = list(miscparm = c(cost = 1), userdist = userdfn2, non_euclid_msk = non_euclid_msk )))

I would like to have the non-euclidean distances calculated dynamically from a more fine-scale mask than the habitat mask.

Thanks

Robby

On Wednesday, July 3, 2019 at 3:39:17 PM UTC-4, Murray Efford wrote:

[Murray Efford]

Obviously, there is no extra mask argument at present. However, I wonder if you can simply write a distance function that uses a named fine mask in the current workspace? Remembering that the distance matrix output needs to be in terms of the coarse mask (from each detector to each coarse-mask point). This would be ad hoc and not easy to generalise or package, if that's your goal.

Murray

[Chris Sutherland]

Robby,

Why not use a density mask that is the resolution you need for the non-euc stuff?

Chris

--
You received this message because you are subscribed to the Google Groups "secr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/cb191f07-f2a4-43dd-94a7-c0cfc830722d%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Robby Marrotte]

That is an easy solution! I'll see if I can get it to work today. Thanks!

[Robby Marrotte]

Hi,

The issue I am having is that cost values outside the mask are given the same cost as obstacles inside the buffer region. I want pixels outside the buffer to set to NA and obstacles within the buffer assigned cost values by ML.

Thanks

Robby

Chris

To unsubscribe from this group and stop receiving emails from it, send an email to secr...@googlegroups.com.

[Chris Sutherland]

So I think the way to do this is the following:

 

-          Make the resolution of the ‘density mask’ the resolution you would like for the cost surface

-          Make the buffer of the ‘density mask’ larger such that all reasonable paths (i.e., least cost paths from activity centers of individuals with non-zero cap prob to the detectors) are contained within it

 

Chris

To unsubscribe from this group and stop receiving emails from it, send an email to secrgroup+...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/secrgroup/44354e44-bd71-420a-9f6e-804e23c27eca%40googlegroups.com.

[Robby Marrotte]

Thanks this works with no hiccups.

[Robby Marrotte]

Hi,

I'll give it a try, but I am trying to use the same procedures a colleague previously did. 

I want to see whether adding in the cost of movement produces a better model and see how much this changes the density estimate.

Thanks!

[Robby Marrotte]

Hi,

Any idea what scale these betas are being optimized? I have 1 parameter that can range from 1 to 9 and the optimization tends to jump around by small increments.

I am going to try dividing the value by 10 so the search is within 0.1 and 0.9 to see what happens.

Thanks

Robby

[Murray Efford]

You probably figured this out by now. The key wording in ?details is

"Any transformation (link function) etc. is handled by the user in the userdist function."

So by default the link function is identity, but you can treat the beta value as you like within the userdist function. To maximize on the log scale, exponentiate the value in userdist before using it (and provide the logged starting value in secr.fit(..., details = list(miscparm = ))

[Robby Marrotte]

Hi,

No, unfortunately, I missed that. I'll give it a try. 

I can get the rest of the transformation parameter to work its just the type of function I can't get to work.

Right now I loop through 8 different types of transformation functions (Figure 1 in https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12984) and fit a shape and max_scale parameter to each and compare them with AIC.

It seems to work using the default optimizer. 

Thanks

Robby

[Robby Marrotte]

Hi,

How specifically are the euclidean and non-euclidean distance matrices handled within secr.fit? I am having a difficult time finding this procedure in the code. 

Are the distances scaled to unit variance? I am trying to figure out the impact of fitting a model with accumulated cost distance compared to the euclidean distance in meters.

From my understanding, the cost distance function in gdistance computes the accumulated cost distance of the least-cost path, not the length of the LCP.

Would it be better to use the length, since it is also in meters and not cost?

Although, I haven't noticed any major change in sigma when I use a non-euclidean model.

Thanks again

Robby

[Murray Efford]

'secr' does nothing but apply your function and pass the distances. No scaling. I'm not much use on the LCP stuff, but, really, once you have left the Euclidean world anything goes. Either way, the challenge is to explain to readers what the hell it really means. In this case maybe just an ad hoc way to get a slightly better model of the detection process?

Murray

I find this code in the likelihood function

       ## Apply user-provided distance function or matrix
        ## do not use if detector is one of
        ## polygonX, polygon, transectX, transect, telemetry
        ## The matrix distmat is passed to the C function secrloglik, etc.
        if (is.null(details$userdist))
            distmat <- -1
        else {
           
            if (is.null(covariates(session.mask)))
                covariates(session.mask) <- data.frame(row.names = 1:m)
            if ('noneuc' %in% userdistnames)
                covariates(session.mask)$noneuc <- NE[1:m,,min(dim(NE)[3],sessnum)]
            if ('D' %in% userdistnames)
                covariates(session.mask)$D <- density
            ## pass miscellaneous unmodelled parameter(s)
            if (nmiscparm > 0)
                attr(session.mask, 'miscparm') <- miscparm
           
            distmat <- valid.userdist (details$userdist,
                                       detector(session.traps),
                                       xy1 = session.traps,
                                       xy2 = session.mask,
                                       mask = session.mask,
                                       sessnum = sessnum)
            baddist <- (!is.finite(distmat)) | (distmat<0) | is.na(distmat)
            if (any(baddist)) {
                warning ("replacing infinite, negative and NA userdist values with 1e10")
                distmat[baddist] <- 1e10
            }
        }

which calls

valid.userdist <- function (userdist, detector, xy1, xy2, mask, sessnum) {
    if (is.null(userdist)) {
        ## default to Euclidean distance
        result <- edist(xy1, xy2)
    }
    else {
        if (any(detector %in% .localstuff$polydetectors)) {
            stop ("userdist cannot be used with polygon detector types;")
        }
        if (is.function(userdist))
        {
            OK <- getuserdistnames(userdist) %in% names(covariates(mask))
            if ((length(OK)>0) & !all(OK))
                stop ("covariates required by userdist function not in mask : ",
                      paste(getuserdistnames(userdist)[!OK], collapse=','))
            result <- do.call(userdist, c(list(xy1, xy2, mask)))
        }
        else {
            if (is.character(userdist)) {
                userdist <- get(userdist, pos=-1)
            }

            if (is.list(userdist) & !is.data.frame(userdist)) {
                if (missing(sessnum))
                    stop("This use does not yet allow for session-specific userdist")
                else
                    result <- userdist[[sessnum]]
            }
            else {
                result <- userdist
            }
        }
        if (!all(dim(result) == c(nrow(xy1), nrow(xy2))))
            stop ("invalid distance matrix dim = ", dim(result)[1], ',', dim(result)[2])
        baddist <- (!is.finite(result)) | (result<0) | is.na(result)
        if (any(baddist)) {
            warning ("replacing infinite, negative and NA userdist values with 1e10")
            result[baddist] <- 1e10
        }
    }
    result
}