Laundry list: parallelization, random effects with smooths, prediction intervals, non-convergence...

[dmkaplan2000]

Hi,

I am just starting to familiarize myself with sdmTMB, which I am trying to use for a fishery CPUE standardization exercise (replacing the GAMM models developed here with spatio-temporal GLMM or GAMM models). I have a laundry list of basic questions, some of which are undoubtedly a bit naive, and I was wondering if other sdmTMB users could point my in the right directions. The list is below. Thanks in advance!

1. parallelization: The models I am running involve quite a bit of data, so ideally I would use parallelization to speed them up. I asked ChatGPT how to use parallelization with sdmTMB and it said to execute TMB::openmp(NPROCS) and that it would all work out. I just want to confirm that this is the correct approach?

2. random effects: In my current GAMM models implemented with mgcv, I have a random effect for fishing vessel implemented with s(vessel,bs="re"). I am wondering if this approach is still valid in sdmTMB, or should I use GLM-style random effects +(1|vessel)?

3. model non-convergence: So far, my tests of sdmTMB have not gone too well, either because the models do not pass the sanity checks, suggesting simplifying the model, adjusting the mesh, or adding priors or because the model errors out with singularities. Presumably the models are over determined, but there are lots of different ways I could simplify them, making it hard to know what to do next. Do you have suggestions of how to go about diagnosing these issues? E.g., what to simplify first (spatio-temporal effects, smooth effects in formula, mesh, etc.)?

4. prediction intervals: For my CPUE indices here, I have been using prediction intervals based on the approach developed here to get a better estimate of prediction uncertainty. Am I correct to understand that sdmTMB does not currently include prediction of prediction intervals? Does the function get_index include already this sort of thing (e.g., covariance between prediction uncertainties and unexplained variance in index uncertainty)? If not, is there an equivalent of type="lpmatrix" in the mgcv::predict.gam function?

5. zeroing out intercept for get_index: I notice than in examples using get_index, the intercept is zeroes out (e.g., density ~ 0 + ...). Why is this necessary/appropriate?

Thanks a bunch for the assistance and I apologize in advance if any of these questions are off base...

Cheers,
David

[ericward-noaa]

Here's some responses in-line below:

Hi,

I am just starting to familiarize myself with sdmTMB, which I am trying to use for a fishery CPUE standardization exercise (replacing the GAMM models developed here with spatio-temporal GLMM or GAMM models). I have a laundry list of basic questions, some of which are undoubtedly a bit naive, and I was wondering if other sdmTMB users could point my in the right directions. The list is below. Thanks in advance!

1. parallelization: The models I am running involve quite a bit of data, so ideally I would use parallelization to speed them up. I asked ChatGPT how to use parallelization with sdmTMB and it said to execute TMB::openmp(NPROCS) and that it would all work out. I just want to confirm that this is the correct approach?

Response: please check out the manual on this -- if you look at the sdmTMBcontrol() docs on p 41-42, there's example usage: https://cran.r-project.org/web/packages/sdmTMB/sdmTMB.pdf

2. random effects: In my current GAMM models implemented with mgcv, I have a random effect for fishing vessel implemented with s(vessel,bs="re"). I am wondering if this approach is still valid in sdmTMB, or should I use GLM-style random effects +(1|vessel)?

Response: it's better to use the GLM-style notation here (1|vessel)

3. model non-convergence: So far, my tests of sdmTMB have not gone too well, either because the models do not pass the sanity checks, suggesting simplifying the model, adjusting the mesh, or adding priors or because the model errors out with singularities. Presumably the models are over determined, but there are lots of different ways I could simplify them, making it hard to know what to do next. Do you have suggestions of how to go about diagnosing these issues? E.g., what to simplify first (spatio-temporal effects, smooth effects in formula, mesh, etc.)?

Response: the best approach here would be to completely start from scratch, with a spatial - only model and no main effects. Build the pieces up -- it seems like you've got a complicated model you're trying to work backwards from, and it's hard to identify what might be wrong

4. prediction intervals: For my CPUE indices here, I have been using prediction intervals based on the approach developed here to get a better estimate of prediction uncertainty. Am I correct to understand that sdmTMB does not currently include prediction of prediction intervals? Does the function get_index include already this sort of thing (e.g., covariance between prediction uncertainties and unexplained variance in index uncertainty)? If not, is there an equivalent of type="lpmatrix" in the mgcv::predict.gam function?

Response: For prediction intervals, get_index is returning confidence intervals. You could do prediction intervals manually -- the difference would be adding on observation / sampling error, and would have to be done at the level of individual grid cells in the prediction dataframe before being aggregated into a total index

5. zeroing out intercept for get_index: I notice than in examples using get_index, the intercept is zeroes out (e.g., density ~ 0 + ...). Why is this necessary/appropriate?

Response: this is one way to get rid of the global intercept (the other being to add -1 to the formula). This isn't always useful, but helpful for example in cases where you want to include fixed year effects and make them more interpretable.

Thanks a bunch for the assistance and I apologize in advance if any of these questions are off base...

Cheers, David

[seananderson]

Adding some additional thoughts:

1. parallelization: The models I am running involve quite a bit of data, so ideally I would use parallelization to speed them up. I asked ChatGPT how to use parallelization with sdmTMB and it said to execute TMB::openmp(NPROCS) and that it would all work out. I just want to confirm that this is the correct approach?

Yes, e.g., TMB::openmp(n = 2, DLL = "sdmTMB") but I haven't tended to see large improvements. Perhaps I should try removing the explicit parallel declarations in the code because I think TMB can try to figure it out on its own now.

Working with lower resolution meshes, especially during testing, will make a large difference. Finer isn't always better even for out of sample prediction.

Also, what will make a massive difference for speed (also for mgcv) is setting up an optimized BLAS.

3. model non-convergence: So far, my tests of sdmTMB have not gone too well, either because the models do not pass the sanity checks, suggesting simplifying the model, adjusting the mesh, or adding priors or because the model errors out with singularities. Presumably the models are over determined, but there are lots of different ways I could simplify them, making it hard to know what to do next. Do you have suggestions of how to go about diagnosing these issues? E.g., what to simplify first (spatio-temporal effects, smooth effects in formula, mesh, etc.)?

Yes, start with a basic model and build up. One thing to note, given you have several penalized smoothers, is that if an SD penalizing the random effect coefficients on the basis functions (ln_smooth_sigma internally) collapses to zero (i.e., leaving you with just the linear fixed effect), that parameter will be on a boundary and you're likely to get a warning about large standard errors on the parameter. That's not necessarily a problem for prediction in this specific case, but you may wish to re-fit the model with that smoother instead fit as a regular linear effect.

4. prediction intervals: For my CPUE indices here, I have been using prediction intervals based on the approach developed here to get a better estimate of prediction uncertainty. Am I correct to understand that sdmTMB does not currently include prediction of prediction intervals? Does the function get_index include already this sort of thing (e.g., covariance between prediction uncertainties and unexplained variance in index uncertainty)? If not, is there an equivalent of type="lpmatrix" in the mgcv::predict.gam function?

Yes, you could basically take the same approach as is in that blog post. The nsim argument in ?predict.sdmTMB() is taking MVN samples from the joint (fixed + random) parameter vector (technically taking samples from the joint precision matrix). You'd then have to add on your own observation error based on the estimated Tweedie parameters. The nsim argument will return a matrix with the number of rows corresponding to the rows of newdata and the columns corresponding to the number of simulations. You can then aggregate them however you would like. There is a get_index_sims() function to help with this but you could also do it yourself. Keep in mind this MVN sampling is an approximation and it can break down. For applied use in assessments it's generally recommended to use TMB's generalized delta method with the generic bias-correction method applied (what get_index(... bias_correct = TRUE) does). You could alternatively sample with tmbstan instead of the MVN approximation, but that is much slower.

I'm still trying to wrap my head around exactly why you'd want to do that here though. The index isn't something you would observe itself (which is where I'd usually expect to make 'prediction intervals' to compare observations to), most assessment models you'd be entering this kind of index into would be assuming any CV or standard error would be on the mean, and because the observation errors are by definition independent, they aren't going to affect the trend, only inflate the uncertainty on the scale and therefore any catchability parameter.

5. zeroing out intercept for get_index: I notice than in examples using get_index, the intercept is zeroes out (e.g., density ~ 0 + ...). Why is this necessary/appropriate?

Yeah, as Eric said, this is just a convenience here. As long as you have year-factors in your formula, this just forces the coefficients on those parameters to represent the mean for a given year as opposed to the intercept representing the first year and the other coefficients representing offsets from that. In the end, the model is identical. Without factor predictors, dropping the intercept will change the model.

[dmkaplan2000]

Thanks so much to both of you for your help. I have a couple of small replies/questions, but in general this was enormously helpful.

Adding some additional thoughts:

1. parallelization: The models I am running involve quite a bit of data, so ideally I would use parallelization to speed them up. I asked ChatGPT how to use parallelization with sdmTMB and it said to execute TMB::openmp(NPROCS) and that it would all work out. I just want to confirm that this is the correct approach?

Yes, e.g., TMB::openmp(n = 2, DLL = "sdmTMB") but I haven't tended to see large improvements. Perhaps I should try removing the explicit parallel declarations in the code because I think TMB can try to figure it out on its own now.

Working with lower resolution meshes, especially during testing, will make a large difference. Finer isn't always better even for out of sample prediction.

Also, what will make a massive difference for speed (also for mgcv) is setting up an optimized BLAS.

For the moment, it is hard for me to tell if this is working for me or not as the models keep failing due to singular matrices (but they do so quickly ;-)). I will try starting small and building up as suggested. I have over 200,000 observations, so one would think that parallelization would help, but only time will tell.

My confusion was mainly due to the fact that sdmTMBcontrol talks about a parallel parameter that is ignored and I wasn't 100% sure when and where I should execute TMB::openmp (but ChatGPT seems to have figured it out however it figures things out...). I get it now in any case.

3. model non-convergence: So far, my tests of sdmTMB have not gone too well, either because the models do not pass the sanity checks, suggesting simplifying the model, adjusting the mesh, or adding priors or because the model errors out with singularities. Presumably the models are over determined, but there are lots of different ways I could simplify them, making it hard to know what to do next. Do you have suggestions of how to go about diagnosing these issues? E.g., what to simplify first (spatio-temporal effects, smooth effects in formula, mesh, etc.)?

Yes, start with a basic model and build up. One thing to note, given you have several penalized smoothers, is that if an SD penalizing the random effect coefficients on the basis functions (ln_smooth_sigma internally) collapses to zero (i.e., leaving you with just the linear fixed effect), that parameter will be on a boundary and you're likely to get a warning about large standard errors on the parameter. That's not necessarily a problem for prediction in this specific case, but you may wish to re-fit the model with that smoother instead fit as a regular linear effect.

I have tried removing all the smooths and just working with parametric linear predictors in the formula, but these models all error out. I will try your approach of building up from very simple models.

4. prediction intervals: For my CPUE indices here, I have been using prediction intervals based on the approach developed here to get a better estimate of prediction uncertainty. Am I correct to understand that sdmTMB does not currently include prediction of prediction intervals? Does the function get_index include already this sort of thing (e.g., covariance between prediction uncertainties and unexplained variance in index uncertainty)? If not, is there an equivalent of type="lpmatrix" in the mgcv::predict.gam function?

Yes, you could basically take the same approach as is in that blog post. The nsim argument in ?predict.sdmTMB() is taking MVN samples from the joint (fixed + random) parameter vector (technically taking samples from the joint precision matrix). You'd then have to add on your own observation error based on the estimated Tweedie parameters. The nsim argument will return a matrix with the number of rows corresponding to the rows of newdata and the columns corresponding to the number of simulations. You can then aggregate them however you would like. There is a get_index_sims() function to help with this but you could also do it yourself. Keep in mind this MVN sampling is an approximation and it can break down. For applied use in assessments it's generally recommended to use TMB's generalized delta method with the generic bias-correction method applied (what get_index(... bias_correct = TRUE) does). You could alternatively sample with tmbstan instead of the MVN approximation, but that is much slower.

I'm still trying to wrap my head around exactly why you'd want to do that here though. The index isn't something you would observe itself (which is where I'd usually expect to make 'prediction intervals' to compare observations to), most assessment models you'd be entering this kind of index into would be assuming any CV or standard error would be on the mean, and because the observation errors are by definition independent, they aren't going to affect the trend, only inflate the uncertainty on the scale and therefore any catchability parameter.

I understand your comment and I go back and forth on this myself, but the reason I do this is that I have noticed that using the standard error approach often produces error bars that are really small, much smaller than the interannual variability in the time series that is undoubtedly partly driven by true uncertainty. By using prediction intervals, I am essentially recreating the fishing data, but this time fishers go out with a spatially and temporally fixed effort distribution. This seems to me a more realistic representation of the uncertainty in standardized CPUE abundance indices that are ultimately the combination of "new" predictions in multiple strata. As to what the assessment model does with those uncertainties, that is another question, but at least visually using prediction interval uncertainties produces a time series with defensible uncertainties.

5. zeroing out intercept for get_index: I notice than in examples using get_index, the intercept is zeroes out (e.g., density ~ 0 + ...). Why is this necessary/appropriate?

Yeah, as Eric said, this is just a convenience here. As long as you have year-factors in your formula, this just forces the coefficients on those parameters to represent the mean for a given year as opposed to the intercept representing the first year and the other coefficients representing offsets from that. In the end, the model is identical. Without factor predictors, dropping the intercept will change the model.

OK, I see. I was initially using a smooth for the non-factor year effect, so in this case I understand that I would definitely want to include an intercept. On the other hand, the simpler models I am currently trying to run are using factor time variables, so in this case, I should zero out the intercept to get an error estimate for year 0.

Thanks again!

[dmkaplan2000]

Also, what will make a massive difference for speed (also for mgcv) is setting up an optimized BLAS.

Regarding optimized BLAS, I think I already have an optimized version. I tried running the suggested code on that page and it gave an elapsed time of ~0.08s on my oldish desktop computer and a similar time for the oldish cluster I am using. This is more than indicated on the web page, but still considerably less than a second.

[dmkaplan2000]

I have a few more questions. I have started looking into why my models die and as soon as I add a spatio-temporal effect, I have problems. This led me to look at the meshes I am using and I am confused by the results.

1. When I plotted the mesh, I found that it doesn't cover the actual data itself. Is this perhaps because a subset of the full dataset is used for large datasets? The mesh was created with the following code:

mymesh = make_mesh(data |> as.data.frame(),c("lon","lat"),cutoff=10,type="cutoff")

Note that I am explicitly using a low resolution mesh to start with and I am using lon,lat as spatial coordinates as the zone of the data is too large for UTM, but relatively close to the equator so distortion from true distance is small. This produced the following mesh:

I "fixed" this by identifying points on the convex hull of the dataset and moving them to the top of the dataset (i.e., the first few rows). After this change, the mesh looks better, though some of the points are still very close to the mesh edge:

Is this the correct approach or is there a better way to assure that the mesh covers the data?

2. One reason that I can imagine that the models are having problems is that there may be little or no data for certain combinations of year-triangle (or year-node; I am not sure if the fit is by center of triangles or by triangle vertex). Is there a convenient way to look a this, for example extracting the number of data points per triangle-year? Alternatively, is there a way to convert the mesh into a set of triangular polygons that I can intersect with the original dataset to do these counts?

3. A minor query is that I noticed in the help of the plot.sdmTMBmesh that this function is supposed to return a ggplot2 object if that package is installed, but looking at the actual code of getAnywhere(plot.sdmTMBmesh), it makes no mention of ggplot2 and when I execute plot(mymesh) it does not seem to return a ggplot2 object. Is this normal? Am I misunderstanding something?

[dmkaplan2000]

After several tests, it seems like the mesh not covering all the data was a large part of my problem. After fixing that, the model convergences with a fairly complex set of predictors. Based on these working models, I have a few additional questions (sorry to ask so many, but I am in the thick of it...):

1. One of the more complex models I have run fails one aspect of the sanity check: `b_j` gradient > 0.001. I tried running run_extra_optimization(model,nlminb_loops = 1, newton_loops = 1), but this didn't resolve the issue. I will try standardizing covariates, but do you have any comments regarding this particular sanity check that might help me decide how serious a problem it is?

2. There is one important factor that definitely causes the models to fail - a boolean variable that is true for all years after the imposition of a quota in 2017. It isn't surprising to me that this causes issues as this variable is redundant of the year variable, but what I wanted to do is include a step function for the period before and after imposing the quota that should help standardize for this management change when predicting. Do you have any suggestions for this type of predictor? I have looked at various options, such as breakpoint models or thresholds, but as far as I can tell none of it seems to fit this use case, though I am not sure to have captured all the subtleties.

3. I have been using parallelization and it definitely seems to be occupying several CPUs. The strange thing, however, is that sometimes it seems to go overboard, using essentially all of the CPUs on the machine. For example, if I give it 4 processors, sometimes the linux top command tells me that the R process is using 700% of CPU. I am not sure how top calculates CPU usage, so this might not be accurate, but it seemed high to me. Does this make sense to you?

As always, thanks for the help!

[dmkaplan2000]

Hi, I imagine you have moved on to other things, but the issue with meshes not covering the data seems important to me. I have tried various ways of reordering the data and they do significantly change the resulting mesh, which seems strange to me. Is this a known issue documented somewhere? If not, then perhaps a reordering of datasets to place the convex hull points at the beginning of the dataset before running make_mesh would be a useful addition to the function? I can share the code I used to do this if you would like.

[seananderson]

Sorry I missed replying to this. Answers below:

When I plotted the mesh, I found that it doesn't cover the actual data itself. Is this perhaps because a subset of the full dataset is used for large datasets? The mesh was created with the following code:

I "fixed" this by identifying points on the convex hull of the dataset and moving them to the top of the dataset (i.e., the first few rows). After this change, the mesh looks better, though some of the points are still very close to the mesh edge:

Is this the correct approach or is there a better way to assure that the mesh covers the data?

... the issue with meshes not covering the data seems important to me. I have tried various ways of reordering the data and they do significantly change the resulting mesh, which seems strange to me. Is this a known issue documented somewhere? If not, then perhaps a reordering of datasets to place the convex hull points at the beginning of the dataset before running make_mesh would be a useful addition to the function? I can share the code I used to do this if you would like.

This is all an issue with INLA, now fmesher, unfortunately, and out of my hands. The not covering the data part is not something I've commonly seen. The default sdmTMB::make_mesh() is pretty basic and calls fmesher::fm_rcdt_2d_inla() with either the specified cutoff or knots specified via a kmeans algorithm. In retrospect, I'd design this differently. A "better" option is likely to develop your own mesh with fmesher::fm_mesh_2d() or fmesher::fm_mesh_2d_inla() and pass that in via the mesh argument. You can also specify that within make_mesh() by specifying the fmesher_func argument and passing arguments. There are some examples in the ?make_mesh and the last preprint version. https://doi.org/10.1101/2022.03.24.485545 (I'll update that shortly too)

One reason that I can imagine that the models are having problems is that there may be little or no data for certain combinations of year-triangle (or year-node; I am not sure if the fit is by center of triangles or by triangle vertex).

In general, that shouldn't be a major issue. The vertices become the random effects. It's OK if they don't all have data associated with them.

Is there a convenient way to look a this, for example extracting the number of data points per triangle-year? Alternatively, is there a way to convert the mesh into a set of triangular polygons that I can intersect with the original dataset to do these counts?

See this example:

library(sdmTMB)
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 5, type = "cutoff")

# vertices:
head(mesh$mesh$loc[,1:2])
#>          [,1]     [,2]
#> [1,] 446.4752 5793.426
#> [2,] 446.4594 5800.136
#> [3,] 436.9157 5802.305
#> [4,] 420.6101 5771.055
#> [5,] 408.2088 5771.287
#> [6,] 414.3656 5760.894

loc <- mesh$mesh$loc

# every data point gets a weighted average of the nearest 3 vertices
# based on this matrix:
A <- mesh$A_st

# so, for row of data 3, bilinear interpolation is the combination
# of these vertices:
vert <- which(A[3, ] != 0)
vert
#> [1]   2 535 563

And these are the weights:

A[3, vert]
#> [1] 0.6333086 0.2156787 0.1510126

plot(mesh$mesh)
points(loc[vert,], col = "blue", pch = 20)
points(pcod[3, c("X", "Y")], col = "red", pch = 20, cex = 0.8)

^{Created on 2024-07-03 with reprex v2.1.0}

A minor query is that I noticed in the help of the plot.sdmTMBmesh that this function is supposed to return a ggplot2 object if that package is installed, but looking at the actual code of getAnywhere(plot.sdmTMBmesh), it makes no mention of ggplot2 and when I execute plot(mymesh) it does not seem to return a ggplot2 object. Is this normal? Am I misunderstanding something?

Thanks, I just fixed the docs. I had to pull out that ggplot code when I removed INLA as an imported package because it relies on inlabru. There's now a note about how you can pass your_mesh$mesh (or any INLA mesh) to inlabru::gg().

One of the more complex models I have run fails one aspect of the sanity check: b_j gradient > 0.001. I tried running run_extra_optimization(model,nlminb_loops = 1, newton_loops = 1), but this didn't resolve the issue. I will try standardizing covariates, but do you have any comments regarding this particular sanity check that might help me decide how serious a problem it is?

There's no magic number. If that's the only issue and the gradient is still < 0.01, I probably wouldn't be too worried. The default when fitting is now newton_loops = 1, which in general should be sufficient.

sanity(fit, gradient_thresh = 0.01)

There is one important factor that definitely causes the models to fail - a boolean variable that is true for all years after the imposition of a quota in 2017. It isn't surprising to me that this causes issues as this variable is redundant of the year variable

Yes...

but what I wanted to do is include a step function for the period before and after imposing the quota that should help standardize for this management change when predicting. Do you have any suggestions for this type of predictor? I have looked at various options, such as breakpoint models or thresholds, but as far as I can tell none of it seems to fit this use case, though I am not sure to have captured all the subtleties.

You'll need some kind of time process if you want to estimate another parameter that changes as a factor variable between years. E.g. time_varying = ~ 1, time_varying_type = "ar1", which is like formula = as.factor(year) assuming time = "year" but with the intercepts by year following an AR1 process. A random walk but be a slightly less restrictive option that wouldn't mean-revert. That or assume the only difference between those years is explained by your management change, i.e., manipulate your factors. Or put a strong prior on the step change.

I have been using parallelization and it definitely seems to be occupying several CPUs. The strange thing, however, is that sometimes it seems to go overboard, using essentially all of the CPUs on the machine. For example, if I give it 4 processors, sometimes the linux top command tells me that the R process is using 700% of CPU. I am not sure how top calculates CPU usage, so this might not be accurate, but it seemed high to me. Does this make sense to you?

If you're using OpenBLAS, I believe the default is to use all cores. You may be able to control this with:

RhpcBLASctl::blas_set_num_threads(1)
RhpcBLASctl::omp_set_num_threads(1)

and put that in your .Rprofile.

[dmkaplan2000]

Thanks. I will have to chew on it a bit to understand the comments about the quota step function. Regarding the mesh, one potential reason that I am seeing this issue, but not commonly in other cases, is that I am using lon,lat instead of UTM coordinates, which means that my scale is different from typical scales by about a factor of 100 (I guess I could try multiplying lon,lat by 111 km to test, but I haven't done that). Perhaps this interacts with some of the defaults in the mesh generation algorithm. At some point, I tried changing the defaults regarding mesh-data distance, but didn't get it to behave the way I wanted (a somewhat larger mesh envelope). In any case, putting some of the points near the edge of the convex hull at the top of the data frame does seem to fix the issue. In my latest incarnations, I don't just use the convex hull points, but also all positions within a certain distance of the convex hull, and this produces a mesh with no points outside or on the edge of the mesh.

[seananderson]

I am using lon,lat instead of UTM coordinates, which means that my scale is different from typical scales by about a factor of 100

This can affect model convergence if it means the Matérn range parameter is very small (or in other cases, very big). The range parameter is distance in units of your coordinates.

If the data span much wider than a UTM zone, you could use an equal area projection for your area of interest projection such as Albers in km. There's an example with Snowy Owls across North America in the preprint. There's an example in the Atlantic here. My understanding is you want lat_1 and lat_2 to be roughly 1/6 and 5/6 of the latitude range. lat_0 and lon_0 just control what is zero in the new coordinates.

[dmkaplan2000]

Thanks for the info. In the case of our zone, pure lon,lat only has a maximum distance error of <5%, so it is acceptable in my view, but if I get a chance, I will test the albers projection. To test just how scale impacts results, the easiest thing is to center the lon,lat (i.e., remove the means) and then multiple by 111 to convert to km. I will give it a try...

More generally, is there some way to standardize the Matérn range so that it is unit-less? It seems strange that the convergence should depend on the choice of units (e.g., I imagine using meters instead of km could also cause differences even if this would just be moving the decimal place).

[seananderson]

The Matérn range is defined as the distance at which correlation has decayed to ~0.13, so it will always be in the units of the coordinates. That said, the model should be identical except for the scale of the Matérn range parameter. The only issue is computational at times if the units are far too big or far too small.

[dmkaplan2000]

In that case, I wouldn't think it should make much difference as the units are smaller when using lon,lat, but nothing dramatic, about a factor of 100.

[dmkaplan2000]

If scale is a real issue, would it be logical to internally in sdmTMB rescale the spatial units so that they are always centered and of order 1?

[dmkaplan2000]

And by the way, I recently tried generating the mesh on a machine running a newer version of Linux and it covers all the data without having to reorganize that data frame to put the convex hull at the beginning. Perhaps this is an issue with an older version of one of the C libraries used to generate the mesh...