Lapack exactly singular issue / question

[Ben Padilla]

Hi oSCR crew 

Question for you all about a large dataset my colleagues and I are working on. The basics of the data are that we have multiple regions (management units) that were sampled at varying intervals over the course of 10 years. When a unit was surveyed there were 5-7 grids surveyed using scat dogs. We have done single unit-year analyses and are now trying to combine it all into a comprehensive analysis that includes all units over all years.

I've been working on a subset of the data from a single unit that was sampled over 4 consecutive years. Using each grid-year combination as a session we end up with 28 sessions over 4 years and the model fits fine, albeit slowly. I can even fit this with spatial covariates and sex specific p0 and sigma. 

In all of our data the number of grid-year-wmu sessions would be well more than 100. So, I tried to build a model where all the grids in a given unit-year is the session. I can build the model objects with data2oscr and state space objects. When fitting the NULL ~1 ~1 ~1 model with these objects I run into a singularity error (Lapack routine dgesv: system is exactly singular: u[1,1] = 0). This is the same dataset used to fit the model described above just with the sessions described differently. Is there any reason why multiple grids surveyed over a ~2 week period couldnt be considered a single session? 

If you have any other suggestions on how to efficiently analyze this dataset I would love to hear it. Thanks! 

Ben 

[Daniel Linden]

Hey Ben, conceptually you should be able to define a session however you want.  Your strategy here makes sense.  I assume you don't have individuals that were observed in multiple grids?  If the individual grid-year sessions worked, then something about the combination is causing this error.  I might double check the naming conventions, particularly with traps and individuals.  It could be that something is getting sorted in a strange way.

Can you post a figure of the combined state space and traps?  Just curious about the layout.

--
You received this message because you are subscribed to the Google Groups "oSCR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to oscr_package...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/oscr_package/bbd829b2-ea65-41c7-9273-e229949346c8n%40googlegroups.com.

[Ben Padilla]

Hi Dan, 

We arent entirely sure about individuals across grids. Across all of our data out of more than 15,000 individuals there are 66 that were detected at >1 grid, but the distances between those grids and the timings of surveys has me questioning whether there was a lab error there. To make double sure all the naming conventions and everything were getting sorted properly I added unit-year-grid to each trap and individual ID. Correcting for individuals got me past the Lapack routine error, but instead I ran into "Error in nlm" due to non-finite value supplied to nlm. This is a model fit using starting values derived from using the getStarts=TRUE argument. Below is a plot of one year's surveys in one management unit, where each of those squares is approx a3x3 km2 surveyed by scat dog teams. Then below it shows the spatial scale of our sampling frame. This is really an analysis that needs random effects for the management units, but I was having trouble figuring how to code this model in nimble given the complexity of the state space. oSCR has been working well and I think/hope we can make it work. 

Thanks for the help

Ben

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife
Twitter: @bpdilla ~ Web: www.benpadilla.weebly.com

You received this message because you are subscribed to a topic in the Google Groups "oSCR" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/oscr_package/etRrlDR_W8Y/unsubscribe.
To unsubscribe from this group and all its topics, send an email to oscr_package...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAFWOE77edQsbOPt_Tqd%2B8P%2BobyPdg6BGTsh4XiTFKBeSsyGftA%40mail.gmail.com.

[Daniel Linden]

Have you scaled your distance units?  I think starting values become more important when you have a lot of data like you do here, and with unscaled distance units even more so.  Definitely a tricky problem to troubleshoot though.  I could always give it a go if you send along your data/code.

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAMLk4GX4GzuKP4uEQj6rRBpiveuJPsZcqwxMcXxr3%3DmCzjezYw%40mail.gmail.com.

[Ben Padilla]

We havent scaled the spatial data, that is something I always forget to do! What is your recommendation for best approach? In the past I've generally just reduced the utm coordinates from meters to km or something by moving the decimal just so that the numbers are closer to zero and easier to handle computationally. Or, would you recommend a scale and center transformation? 

Seems like we are making some progress with getting the full dataset together and a null model is spinning without propagating errors yet. Hopefully it doesnt run for 10 hours then run into a singularity! 

Thanks for the help. 

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife

Web: www.benpadilla.weebly.com

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAFWOE75oZXg%2BbN%3DFcah0c1Ryo0f-%2BKZcxSW41Sr7zwjhtb3EqQ%40mail.gmail.com.

[Daniel Linden]

It mostly depends on what sigma might be, but usually having distance units in km instead of m is enough. When the log-scale model coefficients are small single digits, the optimization tends to work better.

So yeah, making sure your coordinates for the state space and trap locations are in km is all you probably need to do.

Glad you’ve got something working, or at least not breaking quickly!

On Apr 17, 2025, at 5:21 PM, Ben Padilla <benjamin...@gmail.com> wrote:



We havent scaled the spatial data, that is something I always forget to do! What is your recommendation for best approach? In the past I've generally just reduced the utm coordinates from meters to km or something by moving the decimal just so that the numbers are closer to zero and easier to handle computationally. Or, would you recommend a scale and center transformation? 

Seems like we are making some progress with getting the full dataset together and a null model is spinning without propagating errors yet. Hopefully it doesnt run for 10 hours then run into a singularity! 

Thanks for the help. 

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife
Web: www.benpadilla.weebly.com

On Thu, Apr 17, 2025 at 12:25 PM Daniel Linden <danl...@gmail.com> wrote:

Have you scaled your distance units?  I think starting values become more important when you have a lot of data like you do here, and with unscaled distance units even more so.  Definitely a tricky problem to troubleshoot though.  I could always give it a go if you send along your data/code.

On Wed, Apr 16, 2025 at 2:34 PM Ben Padilla <benjamin...@gmail.com> wrote:

Hi Dan, 

We arent entirely sure about individuals across grids. Across all of our data out of more than 15,000 individuals there are 66 that were detected at >1 grid, but the distances between those grids and the timings of surveys has me questioning whether there was a lab error there. To make double sure all the naming conventions and everything were getting sorted properly I added unit-year-grid to each trap and individual ID. Correcting for individuals got me past the Lapack routine error, but instead I ran into "Error in nlm" due to non-finite value supplied to nlm. This is a model fit using starting values derived from using the getStarts=TRUE argument. Below is a plot of one year's surveys in one management unit, where each of those squares is approx a3x3 km2 surveyed by scat dog teams. Then below it shows the spatial scale of our sampling frame. This is really an analysis that needs random effects for the management units, but I was having trouble figuring how to code this model in nimble given the complexity of the state space. oSCR has been working well and I think/hope we can make it work. 

Thanks for the help

Ben

<image.png>

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAMLk4GXb1_nCW9%3Djq7BNxtZfkRnjx9ug-UGKna7rf4pD0VBBBw%40mail.gmail.com.

[Chris Sutherland]

Sorry to be late to the party here.

With such a large data set, there are a few things I would do:

1. Fit a bunch of single session models to get a sense of the session-specific variability in parameters and identify any quirks.

2. Do this with a relatively coarse resolution state space and collapsed occasions (with Poisson encounters) at first for computational efficiency.

3. Then, once satisfied all is good, start assembling the MS model using informed start values from steps 1 & 2 which will speed up model fitting.

Hope that's of some use. By the way, it's an extremely decent looking data set. Once you have the model up and running, I can see loadsa potential for cool questions! Good luck with it all.

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/7F78C471-A191-464D-A310-B44C27BC2C96%40gmail.com.

[Ben Padilla]

Hi Chris, 

No need to apologize, welcome to the party! 

These are all great ideas. We have fit all the single year-unit models with and without covariates, great idea to use those as starting values for the full model, I hadnt thought of that. We will give these suggestions a try. And yes, this is a great dataset with a lot of potential. It took a lot of work to get everything sorted out and in working order but it will be well worth the efforts in the end! 

Ben

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CANjB3vSUSV8cbqT1XmxRegtUhFVPbsoz14XS8ixHXOvKk2sxEw%40mail.gmail.com.

[Ben Padilla]

Hi Chris and Dan, 

A follow up question here about running in parallel. We have some code to fit our models in parallel using foreach and doParallel. I'm pretty sure it is essentially what is going on under the hood in the oscr.parfit function. But, based on this thread this effectively runs each model in the list on a different core, but doesnt speed up runs of a given model by spread it out across multiple cores. 

We have a long list of models we'd like to fit, so spreading them out across multiple cores is great, but each individual model will still take a week or more at best. We are using scaled state space, covariates, and starting values (from getStarts). Do you know if there is a way to spread a single model across multiple cores with the way the oSCR.fit function works? 

Thanks for the help! 

Ben

[Daniel Linden]

Hi Ben, we use "nlm" for the likelihood minimization process and I do not know of any parallel versions.  Technically, we could use "optim" instead and there is a parallelized version of that ("optimParallel").  That would require getting into the guts of the oSCR.fit function, which is not impossible.  The code is explicit and designed to allow folks to poke into once you understand the steps.  But this would not be an insignificant venture...

The other option is to make sure that your state space resolution is as coarse as possible, to reduce the computation time.  This kind of tuning can be tricky and you'd probably have to settle for an approximation to make it worthwhile (i.e., having a state space that is a little too coarse).

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/e2054bba-b201-4b51-9313-ebd03838a8a2n%40googlegroups.com.

[Víctor Masías Hinojosa]

Hi,

Its not tested but its close to replace "nlm" by optimParallel.

Best.

oSCR.fit <-
function (model = list(D ~ 1, p0 ~ 1, sig ~ 1, asu ~1), scrFrame, ssDF,
          encmod = c("B", "P", "CLOG","M")[1], multicatch = FALSE, theta = 2,
          trimS = NULL, DorN = c("D", "N")[1], sexmod = c("constant", "session")[1],
          costDF = NULL, distmet = c("euc", "user", "ecol")[1], directions = 8,
          PROJ = NULL, rsfDF = NULL, RSF = FALSE, telemetry.type = c("none","ind","dep")[1],
          se = TRUE, predict = FALSE, start.vals = NULL, getStarts = FALSE, pxArea = 1,
          plotit = F, mycex = 1, nlmgradtol = 1e-06, nlmstepmax = 10, nlmiterlim = 200,
          smallslow = FALSE, print.level = 0){
 
  # Load optimParallel package
  if (!require(optimParallel, quietly = TRUE))
    stop("need to install package 'optimParallel'")
 
  # Rest of the original code remains the same until the optimization part
 
  # ... [previous code library(parallel)
library(optimParallel)

oSCR.fit <- function(model = list(D ~ 1, p0 ~ 1, sig ~ 1, asu ~ 1), scrFrame, ssDF,
          encmod = c("B", "P", "CLOG","M")[1], multicatch = FALSE, theta = 2,
          trimS = NULL, DorN = c("D", "N")[1], sexmod = c("constant", "session")[1],
          costDF = NULL, distmet = c("euc", "user", "ecol")[1], directions = 8,
          PROJ = NULL, rsfDF = NULL, RSF = FALSE, telemetry.type = c("none","ind","dep")[1],
          se = TRUE, predict = FALSE, start.vals = NULL, getStarts = FALSE, pxArea = 1,
          plotit = FALSE, mycex = 1, nlmgradtol = 1e-06, nlmstepmax = 10, nlmiterlim = 200,
          smallslow = FALSE, print.level = 0) {

    if (!requireNamespace("optimParallel", quietly = TRUE))
        stop("Package 'optimParallel' is required. Install it using install.packages('optimParallel')")

    # If optimization is required
    if (!getStarts && !predict) {
        message("Fitting model: D", paste(model)[1],
                            ", p0", paste(model)[2],
                         ", sigma", paste(model)[3],
                           ", asu", paste(model)[4], sep = " ")
        if (!anySex) {
            if (telem) {
                message("Telemetry: ", telemetry.type)
            }
            message("Using likelihood function 'msLL.nosex' \nHold on tight!")
            message(Sys.time())
            message(paste(pn, " ", sep = " | "))
            message(" ")

            # OptimParallel Implementation
            cl <- makeCluster(detectCores())
            setDefaultCluster(cl = cl)
            myfit <- suppressWarnings(
                optimParallel(par = pv, fn = msLL.nosex, pn = pn, D = D, nG = nG, nK = nK,
                    hiK = hiK, dm.den = dm.den, dm.trap = dm.trap,
                    hessian = hessian, control = list(maxit = nlmiterlim, trace = TRUE),
                    method = "BFGS"))
            stopCluster(cl)
        } else {
            if (telem) {
                message("Telemetry: ", telemetry.type)
            }
            message("Using likelihood function 'msLL.sex' \nHold on tight!")
            message(Sys.time())
            message(paste(pn, " ", sep = " | "))
            message(" ")

            # OptimParallel Implementation
            cl <- makeCluster(detectCores())
            setDefaultCluster(cl = cl)
            myfit <- suppressWarnings(
                optimParallel(par = pv, fn = msLL.sex, scrFrame = scrFrame, pn = pn, D = D, nG = nG, nK = nK,
                    hiK = hiK, dm.den = dm.den, dm.trap = dm.trap,
                    hessian = hessian, control = list(maxit = nlmiterlim, trace = TRUE),
                    method = "BFGS"))
            stopCluster(cl)
        }

        # Adjust output to match nlm structure
        myfit$estimate <- myfit$par
        myfit$minimum <- myfit$value
    }
}

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAFWOE76yYeaQ_0h2%3DBAHf_N2Ng%2BDcmwZ7Reh_%3D0D54dQa2W0Ww%40mail.gmail.com.

--

Víctor Masías Hinojosa (Cargo) - (Unidad) vma...@fen.uchile.cl  Teléfono: +56 2 297x xxxx | Móvil +56 9xxxx xxxx Universidad de Chile, Facultad de Economía y Negocios Diagonal Paraguay 257 - Torre 26 fen.uchile.cl Por un Campus Sustentable, no imprima este correo si no es realmente necesario. En FEN vivimos la sustentabilidad.

[Ben Padilla]

Hi Victor, 

Thanks for sharing. I'll have to have a look at this and test it out. Have you tested it out at all? Did it decrease run times? 

Ben

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CABsnA6io4-pKPpuYtcupQey-Fxho-sasTtX5uKWjh5HeSZXB7g%40mail.gmail.com.

[Ben Padilla]

Hi all, 

Quick follow up question regarding run time and when we should think about checking on things. We have 10 models that have been running simultaneously on 10 cores using parallel and lapply for 4 weeks now. As noted above, it is a large dataset, but we've scaled coordinates, added starting values, etc. to try and speed things up. We've also reduced the model complexity as much as we can. Our most complex model for density includes session (17 sessions), and an interaction between canopy cover and land ownership (3 categories). No model has more than 3 covariates for density (screenshot of our model list attached). Is >4 weeks a reasonable run time for this set of models? Presumably each one is running at the same time on a unique core so we are just limited by the slowest of the bunch. 

If we stop the process will all progress be lost? I understand that this is more a question of how the parallelization is working rather than oSCR itself, but curious if you have any experience? 

Thanks. 

Ben

[Daniel Linden]

Hey Ben, it can be really hard to tell how long the models will take to fit.  It depends on many things, from the spatial resolution to how well your starting values are able to reduce the iterations needed for optimization.

My best advice is to try a version with a very coarse state space, much coarser than is ideal.  Or fitting a model for half the survey area, etc.  Such exploratory attempts can give you some idea of how computational time increases with finer resolutions.  A model with <500 pixels in the state space is usually pretty fast but you would need to test.  I would keep reducing the resolution until you know what you're working with.

Unfortunately, stopping the model fitting will lose everything as far as I know.  Similar to MCMC, you would need some data writing mechanism embedded in your function to get "current values" saved somewhere.  Otherwise, the progress is lost.

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/e3e49619-fdbf-4849-b0c4-47dae757ef96n%40googlegroups.com.

[Ben Padilla]

Hi Dan, 

Thanks for the idea, we'll try that. We had already coarsened the resolution, but clearly not as much as we should have. We are testing out some coarse resolution models. 

Ben

Benjamin Padilla Ph.D
Wildlife Research Supervisor
Oregon Department of Fish and Wildlife

Twitter: @bpdilla ~ Web: www.benpadilla.weebly.com

To view this discussion visit https://groups.google.com/d/msgid/oscr_package/CAFWOE777WychgqOY5F2xe4Y-CBHQj9rV0YR2zbdFSp99F_WJRw%40mail.gmail.com.