search

hmsc-r / HMSC

GNU General Public License v3.0
99 stars 36 forks source link

Numerical instability in a very uneven spatial study design #33

Open gtikhonov opened 4 years ago

gtikhonov commented 4 years ago

Dear @gtikhonov ,

I am having a similar problem trying to sampleMcmc a model 2 random effects, one that was the spatial coordinates and the other a regional categorical variable. Since I have >1000 community samples, the spatial method set to GPP, using about 117 knots. It is a hurdle model, so I am first fitting the presence/absence data and then species abundances (individuals ha-1 from 0.001 until ~2000) and relative densities (0-1 values). I log-transformed abundances and relative densities to use distr = "normal".

All my coordinates are in meters and I have jittered coordinates that were too close to each other to avoid problems related to the precision of the spatial coordinates. However, I still have distances between units that vary from 5 meters to 3e+06 meters. Anyways, not sure if this spatial precision is indeed an issue while using GPP.

Notably, the error only appeared when I increased the thin x samples combination and only for the abundance part of the model. The error appears error by the end of chain 1, sometimes of chain 2: Error in chol.default(iU) : the leading minor of order 19 is not positive definite

I run the tests you suggested for keilast here are the results for the computeDataParameters part of the sampleMCMC function. Here are the results: alpha = 4257548 ag = 101 W is positive-definite and has no NAs 'distance' has no NAs but it is not positive-definite. 's' a matrix with 2 columns of numbers and no NAs, with upper off-diagonal = 8916294 and lower off-diagonal = 5766795 9000554.

In conclusion, I think the problem is probably in the small distance between some of my sites. But there may be a possible interaction between these distances and data type (problem not arises with presence/absence data, only abundances).

Any clues? You mentioned some that you had more detailed advice on how you can handle the whole data. I have spent the entire morning truing to solve this issue, so any help would be very welcome.

Thanks in advance, Renato

Originally posted by @LimaRAF in https://github.com/hmsc-r/HMSC/issues/8#issuecomment-562568327

gtikhonov commented 4 years ago

This was originally commented in #8, but I decided to open a new issue, since the nature of the problem is principally different. Here are the first comments. https://github.com/hmsc-r/HMSC/issues/8#issuecomment-562594146 https://github.com/hmsc-r/HMSC/issues/8#issuecomment-562606481

gtikhonov commented 4 years ago

@LimaRAF based on your description I would rather say that the problem is not because of any misspecification of the random effects. Instead, it is due to the very high variation of spatial scales in your study design. While from mathematical side, it still results in a well-posed problem, when we start doing numerics using modern limited float precision, we end up in computational instabilities that most likely cause the error that you have received. We have encountered similarly-flavoured issues in our studies and below I would try to describe the approach that we took. The principal idea is in separating the spatial scales of your study manually. Let's say that we have a study design, where we sampled multiple locations on a number of small islands in a large sea. I would argue that there are 2 types of spatial correlation that can be present in such system:

  1. Nearby islands are generally exhibiting more similar species distributions than those that are far away. This is going to be the "large scale" spatial autocorrelation.
  2. When we look at each island, the samples from nearby locations are more similar than those that are far away from each other, e.g. located on the opposite sides of the island. This is going to be the "small scale" spatial autocorrelation.

So in principle, defining spatial effect only on the samples level should be fine, as then the estimated spatial scale would determine which type of autocorrelation it is capturing. However, if the estimated scale is very large, its correlation between locations within single island is very-very close to 1, which leads to ill-conditioned covariance matrix, which in turn crashes the estimation algorithm. But such covariance matrix can be substituted with another one, where these "very-very close to 1" values are actually replaced with 1. Such matrix can be obtained if you define a spatial effect, defined on islands themselves. Now, this type of spatial effects is well supported by Hmsc - you need to define a spatial effect of islands, e.g. giving the centre mass points of those as coordinates of them, and atop of that define a spatial effect of samples, but for this one you want to restrict the max limit of spatial range to something on the scale of within-island distances.

Of course, the abovedescribed example is somewhat a simplification of the challenges that arise in practice, but I hope it gives you a conceptual understanding. By the way - it would be great if you can post here a 2-D scatterplot (map) of your xy object, preferable coloured by the sites$ecoreg.af- this would give a much better context to figure out what can be the viable strategy.

LimaRAF commented 4 years ago

@gtikhonov Thanks a lot by such a complete explanation. It it indeed very helpful to understand the roots of the problem and the paths of possible solutions. The example of the islands really helped to visualize the idea and I think this may be the case for my study.

As you requested, I am attaching the map of my xy object (colored by ecore.af - left panel) as well as the Knots I am using to build the GPP (just in case - right panel). Some of the colors are repeated, but the are different regions. I have aggregated some of regions to make them more compact but now I see that I may need to simplified it even further.

Thanks again for your help.

Rplot01

LimaRAF commented 4 years ago

@gtikhonov

Some additions that may be important:

So it may be related to the inclusion of more rare species as well, right?

Thanks again, Renato

chrishaak commented 4 years ago

@gtikhonov

Hi, I'm working on fitting a spatial model on a large dataset (i.e, 1000s of observations), with distances b/w observations ranging from several m up to >1000000 m. In my fits so far, estimating the scale parameter(s) for my spatial LVs (n=2) has proven to be problematic - with chains that are highly divergent or which often "sit" at zero for the entire duration of the mcmc run. In contrast, a model using the same study design but with unstructured LVs seems to converge just fine. Although I'm not receiving the error described above, I'm wondering if my convergence problem could be related to a similar issue (i.e., the large range of spatial scales in my data). In an attempt to remedy this, I tried to follow your earlier "island" suggestion, assigning my samples to larger "clusters" of samples within a given spatial scale (say, 100km), the centroids of which were then used to specify a "broad scale" spatial random level. Then (also based on your suggestion, as I understand it), I tried to define another spatial random level using the individual observations, but constraining the maximum distance b/w observations to a scale less than my cluster size (< 100km). The only way I can see to accomplish this, however, is by specifying a distance matrix that has been modified to limit the max distance, instead of using the raw coordinates. But then, upon trying to fit, I run into the problem that NNGP or GPP cannot be used with a distance matrix... and I have too many observations to use the "full" estimation for this level.

Its very possible that I am misinterpreting your suggestion - can you offer any more details on this approach?

Alternatively, might there be a way to help address this issue with the priors specified for the spatial random level parameters (Alpha, etc?) I can't find much information in the manual on how to do so.

I'd be grateful for any further guidance you can offer!!

Many thanks in advance for your help... Chris

ovaskain commented 4 years ago

Hi Chris et al,

Gleb can probably give a better answer, but here come my thoughts. I might include two random effects. One for large-scale spatial effects, where you would have clustered the sites e.g. to 100 km scale, and where you would use the default prior for alpha. Then another for small-scale spatial effects, targeted e.g. to variation at 1km to 100km scale. For this I would cluster very nearby sites so that smallest distance among sites is ca. 100 meters, but otherwise keep the sites exactly at locations where they are. Then I would set the prior so that alpha is between 1 and 100km. For this random effects, sites very far from each other (say, 1000 km from each other) would just be independent, and they should not cause numerical problems.

Best,

Otso

On 6.4.2020 17.48, chrishaak wrote:

@gtikhonov https://github.com/gtikhonov

Hi, I'm working on fitting a spatial model on a large dataset (i.e, 1000s of observations), with distances b/w observations ranging from several m up to >1000000 m. In my fits so far, estimating the scale parameter(s) for my spatial LVs (n=2) has proven to be problematic - with chains that are highly divergent or which often "sit" at zero for the entire duration of the mcmc run. In contrast, a model using the same study design but with unstructured LVs seems to converge just fine. Although I'm not receiving the error described above, I'm wondering if my convergence problem could be related to a similar issue (i.e., the large range of spatial scales in my data). In an attempt to remedy this, I tried to follow your earlier "island" suggestion, assigning my samples to larger "clusters" of samples within a given spatial scale (say, 100km), the centroids of which were then used to specify a "broad scale" spatial random level. Then (also based on your suggestion, as I understand it), I tried to define another spatial random level using the individual observations, but constraining the maximum distance b/w observations to a scale less than my cluster size (< 100km). The only way I can see to accomplish this, however, is by specifying a distance matrix that has been modified to limit the max distance, instead of using the raw coordinates. But then, upon trying to fit, I run into the problem that NNGP or GPP cannot be used with a distance matrix... and I have too many observations to use the "full" estimation for this level.

Its very possible that I am misinterpreting your suggestion - can you offer any more details on this approach?

Alternatively, might there be a way to help address this issue with the priors specified for the spatial random level parameters (Alpha, etc?) I can't find much information in the manual on how to do so.

I'd be grateful for any further guidance you can offer!!

Many thanks in advance for your help... Chris

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/hmsc-r/HMSC/issues/33#issuecomment-609840589, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIYMZU5BFEEHQQT3YS4VC3RLHTTHANCNFSM4JW2EO6Q.

chrishaak commented 4 years ago

Hi Otso,

Thanks so much for the quick reply and for sharing your thoughts – It sounds like I am on the right track then.

I’m assuming I use “setpriors” to change “myrandomlevel$alphapw”... which is a matrix of uniformly increasing scale values and their corresponding probabilities (summing to 100). Gonna give that a shot!

thanks again, C

gtikhonov commented 4 years ago

Hi @chrishaak

I generally support the approach that @ovaskain has described. Just a small note on the need to cluster - in certain scenarios it makes good sense to do it, but in some others there may be no natural way to organise your sites in groups. There are multiple potential remedies for that and it is hard to a priori guarantee which will be working and which will be not. E.g., in principle, you can use 2 random levels that are defined on the same set of spatial sites - GPP to capture large-scale spatial autocorrelation and NNGP to capture the short-scale autocorrelation. However, although I would say that such scheme is the most universal for the setting that you have described here, I am not aware of whether anyone has tried it in Hmsc and I have a gut feeling that it could lead to some challenges in MCMC convergence. "Island"-based approach shall be generally easier to fit. To manually set the prior for the spatial range parameter, you may use something like the following lines:

alphapwMatrix = cbind(maxDist*c(0:alphaN)/alphaN, c(0.5,rep(0.5/alphaN,alphaN)))
rL = HmscRandomLevel(xData=data.frame(x1=rep(1,length(TD$X$x1)),x2=TD$X$x2))
rL = setPriors(rL,alphapw=alphapwMatrix)

The first row of alphapwMatrix are the discrete values and the second are the prior weights of the values.

chrishaak commented 4 years ago

Thanks Gleb for weighing in!

My model domain is spatially continuous (i.e., it is not "naturally" divisible into discrete units)... so by defining islands I am imposing somewhat arbitrary boundaries... but In such a case, it sounds like the latter approach you described (simultaneous GPP and NNGP) may be more appropriate... so I will also look into giving that a shot and see what kind of issues I run into.

Wish me luck, c

chrishaak commented 4 years ago

Hi Gleb (and Otso),

I ran both alternative approaches last night and got reasonable results from the former (i.e., "island") approach. The latter (simultaneous GPP & NNGP) ran, but appears to have some issues related to the setup of the randomlevels. The GPP level had a single LV and default priors, while the NNGP level used 2 LVs and an altered scale prior for Alpha. In specifying the model, the two spatial rLs were both associated with the "sample" (i.e., observation) level of the study design.

However, in the fitted object, it appears that the priors from the NNGP level somehow got applied also to the GPP level... because in the fitted object there are 2 LVs for both levels (when GPP should only have 1) and likewise the max scale on the GPP level was similarly constrained to the maximum set in the NNGP prior.

Also of note is that when I run "computepredictedvalues" on the object, I receive the error: "dfPiNew does not contain all the necessary named columns"

I'm guessing these issues might be related to the fact that I "linked" both rLs to the same level of the studydesign? Do I instead need to add an additional (duplicate) level to the study design (i.e., "sample1" and "sample2", because each rL needs to be associated with a discrete level/column?

EDIT: AFTER SOME FIDDLING, I FOUND THAT I NEEDED TO DUPLICATE THE OBSERVATION LEVEL OF studyDesign (i.e., "sampleA" and "sampleB") SO THAT EACH SPATIAL RANDOM LEVEL WAS ASSOCIATED WITH A DISTINCT COLUMN...

cgoetsch commented 3 years ago

@gtikhonov and @ovaskain

I am dealing with a similar issue to chrishaak, in that I am running a spatial model on a large dataset with a large variance in the distances between observations (in fact, knowing who Chris is, it is a subset of the same dataset he is working with, but a different project goal). I have attempted implementing the "island" approach suggested in this thread. I clustered my samples to a 100km grid (for testing the method) and assigned the samples to the appropriate grid centroids to define a "broad" scale spatial random effect. I am using the NNGP method for the spatial random effect. As recommended by Otso, I kept the default prior for alpha for the large-scale random effect. Here is where my problem diverges from chrishaak's. Since I am more interested in the broad scale spatial effect, I did not add a second spatial random effect for the 1km-100km scale.

Instead, following Otso's explanation in issue https://github.com/hmsc-r/HMSC/issues/47#issuecomment-646519920, I added an unstructured random effect at the sample level (tow for my dataset) in order to get the species co-occurence at the sample level instead of the grid level. However, this is where the problem arises. When I run the model with both the unstructured random effect (sample level) and the broad-scale spatial random effect (grid level) - either with or without covariates, the alpha latent factors do not initialize and sit at zero for the duration of the mcmc sampling. If I run the model with only the spatial random effect, the alpha LFs initialize and update. I have read through all the issues on github and checked my model specifications again and again. Either I am missing something, or there is something else going on that is preventing the spatial model from running correctly when I add the unstructured random effect of sample.

Here are my data points classified to the grid: Grid_with_Centroids

Here is my code for specifying the model using a toy dataset of 4 species and 3 years (1086 samples): `X = da %>% select(depth, rugosity)

Y = da %>% select(5:8)

xy_grid = as.matrix(da %>% select(cent_lon, cent_lat) %>% unique()) rownames(xy_grid)<-as.factor(unique(da$grid_id))

studyDesign = data.frame( tow = as.factor(da$UniqueID), grid_id = as.factor(da$grid_id))

unstructured random effect at sample (tow) level

rL.tow = HmscRandomLevel(units = studyDesign$tow)

spatial random effect at grid-level, with method NNGP

rL.spatial = HmscRandomLevel(sData=xy_grid, sMethod = 'NNGP', nNeighbours = 8)

XFormula = ~ depth + rugosity

full model

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula, distr="probit", studyDesign=studyDesign, ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))

spatial only model

m.nngp2 = Hmsc(Y=Y, XData = X, XFormula=~1, distr="probit", studyDesign=studyDesign, ranLevels=list("grid_id" = rL.spatial)) `

I have been running the models for testing purposes with the following: nChains = 2 samples = 250 thin = 1 transient = 50thin verbose = 50thin

I can send an .Rdata file with the relevant data if you need it for testing. I have now spent several days trying to diagnose this and am stuck as I do not see any issues in the code. I'd be grateful for any help to figure out where this is going wrong.

ovaskain commented 3 years ago

Hi,

You mentioned that: ”the alpha latent factors do not initialize and sit at zero for the duration of the mcmc sampling”. This can be a result rather than a bug. The value of alpha=0 is given 0.5 probability in the prior, so having alpha “stuck” to 0 does not necessarily mean that anything goes wrong, it can just mean that there is no spatial signal. As this happens only in the model where you include a sample-level random effect as well, this indicates to me that those two random effects compete partially for the same signal. Maybe modify the prior of the spatial random effect so that zero is excluded because you model non-spatial effects in any case already by the sample-level approach?

Otso

From: cgoetsch @.> Sent: keskiviikko 19. toukokuuta 2021 20:28 To: hmsc-r/HMSC @.> Cc: Ovaskainen, Otso T @.>; Mention @.> Subject: Re: [hmsc-r/HMSC] Numerical instability in a very uneven spatial study design (#33)

@gtikhonovhttps://github.com/gtikhonov and @ovaskainhttps://github.com/ovaskain

I am dealing with a similar issue to chrishaak, in that I am running a spatial model on a large dataset with a large variance in the distances between observations (in fact, knowing who Chris is, it is a subset of the same dataset he is working with, but a different project goal). I have attempted implementing the "island" approach suggested in this thread. I clustered my samples to a 100km grid (for testing the method) and assigned the samples to the appropriate grid centroids to define a "broad" scale spatial random effect. I am using the NNGP method for the spatial random effect. As recommended by Otso, I kept the default prior for alpha for the large-scale random effect. Here is where my problem diverges from chrishaak's. Since I am more interested in the broad scale spatial effect, I did not add a second spatial random effect for the 1km-100km scale.

Instead, following Otso's explanation in issue #47https://github.com/hmsc-r/HMSC/issues/47 (#47 (comment)https://github.com/hmsc-r/HMSC/issues/47#issuecomment-646519920), I added an unstructured random effect at the sample level (tow for my dataset) in order to get the species co-occurence at the sample level instead of the grid level. However, this is where the problem arises. When I run the model with both the unstructured random effect (sample level) and the broad-scale spatial random effect (grid level) - either with or without covariates, the alpha latent factors do not initialize and sit at zero for the duration of the mcmc sampling. If I run the model with only the spatial random effect, the alpha LFs initialize and update. I have read through all the issues on github and checked my model specifications again and again. Either I am missing something, or there is something else going on that is preventing the spatial model from running correctly when I add the unstructured random effect of sample.

Here are my data points classified to the grid: [Grid_with_Centroids]https://user-images.githubusercontent.com/83247406/118856131-9e01d680-b8a4-11eb-866e-ae730c9d4493.PNG

Here is my code for specifying the model using a toy dataset of 4 species and 3 years (1086 samples): `X = da %>% select(depth, rugosity)

Y = da %>% select(5:8)

xy_grid = as.matrix(da %>% select(cent_lon, cent_lat) %>% unique()) rownames(xy_grid)<-as.factor(unique(da$grid_id))

studyDesign = data.frame( tow = as.factor(da$UniqueID), grid_id = as.factor(da$grid_id))

unstructured random effect at sample (tow) level

rL.tow = HmscRandomLevel(units = studyDesign$tow)

spatial random effect at grid-level, with method NNGP

rL.spatial = HmscRandomLevel(sData=xy_grid, sMethod = 'NNGP', nNeighbours = 8)

XFormula = ~ depth + rugosity

full model

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula, distr="probit", studyDesign=studyDesign, ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))

spatial only model

m.nngp2 = Hmsc(Y=Y, XData = X, XFormula=~1, distr="probit", studyDesign=studyDesign, ranLevels=list("grid_id" = rL.spatial)) `

I have been running the models for testing purposes with the following: nChains = 2 samples = 250 thin = 1 transient = 50thin verbose = 50thin

I can send an .Rdata file with the relevant data if you need it for testing. I have now spent several days trying to diagnose this and am stuck as I do not see any issues in the code. I'd be grateful for any help to figure out where this is going wrong.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/hmsc-r/HMSC/issues/33#issuecomment-844314671, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEIYMZWJY4M3HF2BPEA7DX3TOPYKRANCNFSM4JW2EO6Q.

cgoetsch commented 3 years ago

Otso,

Thank you for clarifying that the alpha latent factors at zero does not mean that there is something wrong with the model. This misunderstanding was driving us a bit crazy looking for what we were doing wrong. I will try your suggestion to modify the prior of the spatial random effect so that zero is excluded and see what happens. This may also be resolved when we move to modeling the entire dataset and not just a small subsample.

Once again, I really appreciate your quick feedback.

Chandra

cgoetsch commented 3 years ago

Otso,

I have a follow-up question to the model I detailed above. I am trying to run the HMSC spatial predictions on my model with the spatial random effect at a 100 km grid and unstructured random effect. However, I get the following error when I run the predict function: Error in LFix + Reduce("+", LRan) : non-conformable arrays. I ran the traceback which seemed uninformative:

traceback() 3: predict.Hmsc(m.nngp, Gradient = Gradient, expected = TRUE, predictEtaMean = TRUE) 2: predict(m.nngp, Gradient = Gradient, expected = TRUE, predictEtaMean = TRUE) 1: predict(m.nngp, Gradient = Gradient, expected = TRUE, predictEtaMean = TRUE)

I have included my code below, but am not sure what the problem is.

`# the NNGP model with spatial effect at larger 100km grid scale

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula, distr="probit", studyDesign=studyDesign, ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))`

The prediction grid is at a 4km resolution and appears as so: head(sm_grid_test) pt_id Lon Lat depth slope rugosity seabedforms 1 1 -67.52358 44.46996 -52.81207 0.7794166 1.1219262 2 2 2 -67.47251 44.45993 -67.06577 0.3562204 0.9309982 2 3 3 -67.42146 44.44989 -68.86845 0.4411752 0.8401850 7 4 4 -67.37043 44.43984 -89.23491 0.6646721 1.0424494 7 5 5 -67.31941 44.42976 -103.35296 0.5040390 1.0883795 2 6 6 -67.26841 44.41967 -100.15780 0.2710498 0.7350692 6

I set up the input for the predict function as follows:

xy.pgrid = as.matrix(sm_grid_test %>% dplyr::select(Lon, Lat) )

head(xy.pgrid) Lon Lat [1,] -67.52358 44.46996 [2,] -67.47251 44.45993 [3,] -67.42146 44.44989 [4,] -67.37043 44.43984 [5,] -67.31941 44.42976 [6,] -67.26841 44.41967

XData.grid = data.frame(depth=(sm_grid_test$depth), rugosity=sm_grid_test$rugosity, stringsAsFactors = TRUE)

head(XData.grid) depth rugosity 1 -52.81207 1.1219262 2 -67.06577 0.9309982 3 -68.86845 0.8401850 4 -89.23491 1.0424494 5 -103.35296 1.0883795 6 -100.15780 0.7350692

Gradient = prepareGradient(m.nngp, XDataNew = XData.grid, sDataNew = list(grid_id=xy.pgrid)) image

predY = predict(m.nngp, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

Am I specifying something in the Gradient function or predict function incorrectly? Could this be due to conducting the spatial random effect in the model at the 100km grid level instead of the level of the tow?

I think I am missing where I can link the new spatial coordinates to the grid_ids of the spatial random effect. Each new coordinate of the 4km prediction grid is associated with one of the gird_ids from the model's original study design. original studyDesign: image new prediction grid: image

Thanks again for your help in troubleshooting.

ovaskain commented 3 years ago

Hi,

I cannot spot anything obvious. Jari, can you spot the problem? cgoetsch, maybe you can prepare a small reproducible example that produces the error?

O2

On 14.6.2021 21.30, cgoetsch wrote:

Otso,

I have a follow-up question to the model I detailed above. I am trying to run the HMSC spatial predictions on my model with the spatial random effect at a 100 km grid and unstructured random effect. However, I get the following error when I run the predict function: Error in LFix + Reduce("+", LRan) : non-conformable arrays. I ran the traceback which seemed uninformative:

traceback()
3: predict.Hmsc(m.nngp, Gradient = Gradient, expected = TRUE,
predictEtaMean = TRUE)
2: predict(m.nngp, Gradient = Gradient, expected = TRUE,
predictEtaMean = TRUE)
1: predict(m.nngp, Gradient = Gradient, expected = TRUE,
predictEtaMean = TRUE)

I have included my code below, but am not sure what the problem is.

the NNGP model with spatial effect at larger 100km grid scale

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula, distr="probit", studyDesign=studyDesign, ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))

The prediction grid is at a 4km resolution and appears as so: |head(sm_grid_test) pt_id Lon Lat depth slope rugosity seabedforms 1 1 -67.52358 44.46996 -52.81207 0.7794166 1.1219262 2 2 2 -67.47251 44.45993 -67.06577 0.3562204 0.9309982 2 3 3 -67.42146 44.44989 -68.86845 0.4411752 0.8401850 7 4 4 -67.37043 44.43984 -89.23491 0.6646721 1.0424494 7 5 5 -67.31941 44.42976 -103.35296 0.5040390 1.0883795 2 6 6 -67.26841 44.41967 -100.15780 0.2710498 0.7350692 6|

I set up the input for the predict function as follows:

xy.pgrid = as.matrix(sm_grid_test %>% dplyr::select(Lon, Lat) )

head(xy.pgrid)
Lon Lat
[1,] -67.52358 44.46996
[2,] -67.47251 44.45993
[3,] -67.42146 44.44989
[4,] -67.37043 44.43984
[5,] -67.31941 44.42976
[6,] -67.26841 44.41967

XData.grid = data.frame(depth=(sm_grid_test$depth), rugosity=sm_grid_test$rugosity, stringsAsFactors = TRUE)

head(XData.grid)
depth rugosity
1 -52.81207 1.1219262
2 -67.06577 0.9309982
3 -68.86845 0.8401850
4 -89.23491 1.0424494
5 -103.35296 1.0883795
6 -100.15780 0.7350692

Gradient = prepareGradient(m.nngp, XDataNew = XData.grid, sDataNew = list(grid_id=xy.pgrid)) image https://user-images.githubusercontent.com/83247406/121939302-e208bf80-cd1a-11eb-9cce-a34b94283d0b.png

predY = predict(m.nngp, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

Am I specifying something in the Gradient function or predict function incorrectly? Could this be due to conducting the spatial random effect in the model at the 100km grid level instead of the level of the tow? I didn't think that I needed to assign my prediction points to 100 km grid squares as well as the observation points, since this is already incorporated into the spatial random effect of the fitted model. But maybe I am wrong?

I am using a small test model of 4 species and 3 years. I am concerned with getting my code specified correctly so that it runs, rather than getting accurate model predictions at this point.

Thanks again for your help in troubleshooting.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hmsc-r/HMSC/issues/33#issuecomment-860900573, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIYMZQLTWXLK6GFUCWKY7DTSZDFFANCNFSM4JW2EO6Q.

cgoetsch commented 3 years ago

Otso,

Below is code that produces the error. Attached is an Rdata file with the necessary data objects to run the code. cgoetsch_example_data.zip

# R objects needed for code

#'[da: The presence/absence test data with 1086 observations: 
#'[    3 years (2013, 2015, 2015) 
#'[    4 species (atlher, alewif, bluher, butter)
#'[    3 environmental covariates (depth, slope rugosity)
#'[    UniqueID with Lon and Lat coordinates for each individual tow
#'[    grid_id with centroid Long and Lat coordinates for the 100 km grid squares to define the spatial random effect

#'[sm.grid.test: The test prediction grid with 7,975 new points (only a portion of our final prediction grid):
#'[              pt_id - unique identifier for each prediction point on the 4km prediction grid with associated Lon and Lat coordinates
#'[              4 environmental covariates - depth, slope, rugosity, seabedforms (only using 2 for test model)
#'[              grid_ID to associate each new prediction point with the 100km scale grid we used for the spatial random effect,
#'[                  with the centroid coordinates for that grid square

#Setting up data
load("filepath/cgoetsch_example_data.RData")

#community data
Y = da %>%
  dplyr::select(5:8) 

#environmental data
X = da %>%
  dplyr::select(depth, rugosity)

#spatial coordinates of the centroids of the 100km grid squares

xy_grid = as.matrix(da %>%
                      dplyr::select(cent_lon, cent_lat) %>%
                      unique())

rownames(xy_grid)<-as.factor(unique(da$grid_id))
plot(xy_grid)

##The model

#StudyDesign with an unstructured random effect at the level of the tow and spatial random effect at level of the 100 km grid
studyDesign = data.frame(
  tow = as.factor(da$UniqueID),
  grid_id = as.factor(da$grid_id))

#associates each tow with one of the 42 100km grid squares

#Unstructured random effect of sample level 
rL.tow = HmscRandomLevel(units = studyDesign$tow)

#Spatial random effect at level of the 100 km grid
rL.spatial = HmscRandomLevel(sData=xy_grid, sMethod = 'NNGP', nNeighbours = 8) 

#Test model formula
XFormula = ~ depth + rugosity

#The test model
m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula,
              distr="probit", studyDesign=studyDesign,
              ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))

#Fit the model

set.seed(1)

#testing parameters
nChains = 2
samples = 500
thin = 5
transient = 50*thin
verbose = 50*thin

#fit the model
m.nngp = sampleMcmc(m.nngp, thin = thin, samples = samples, transient = transient,
                    nChains = nChains, verbose = verbose, 
                    updater=list(GammaEta=FALSE))

#Make spatial predictions

#spatial coordinates of the prediction grid points
xy.pgrid = as.matrix(sm_grid_test %>%
                       dplyr::select(Lon, Lat))

#environmental data for the prediction grid points
XData.grid = data.frame(depth=sm_grid_test$depth, rugosity=sm_grid_test$rugosity, stringsAsFactors = TRUE) 

#Format data for spatial predictions
Gradient = prepareGradient(m.nngp, XDataNew = XData.grid, sDataNew = list(grid_id=xy.pgrid))

#compute predictions
predY = predict(m.nngp, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)
cgoetsch commented 3 years ago

Otso,

Update: Today, I tried to make predictions on a version of the same model without the spatial random effect. All of the above code is the same except the following:

#model without the spatial random effect; only includes the unstructured random effect and covariates

m.nngp4 = Hmsc(Y=Y, XData = X, XFormula=XFormula,
               distr="probit", studyDesign=studyDesign,
               ranLevels=list("tow"= rL.tow))

#fit model with unstructured random effect with covariates
m.nngp4 = sampleMcmc(m.nngp4, thin = thin, samples = samples, transient = transient,
                     nChains = nChains, verbose = verbose, 
                     updater=list(GammaEta=FALSE))

#for the model with no spatial random effect
Gradient = prepareGradient(m.nngp4, XDataNew = XData.grid)

#here I left out the SDataNew argument, since there is no spatial random effect in this version of the model

#make predictions of occupancy to the new data
predY = predict(m.nngp4, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

Again, I get the same error as the full model: Error in LFix + Reduce("+", LRan) : non-conformable arrays.

Am I not specifying the code correctly to get predictions on a model with an unstructured random effect and covariates, but no spatial random effect? It should be possible to make spatial predictions for a model even if there is no explicit spatial random effect, correct?

Cheers, Chandra

ovaskain commented 3 years ago

Hi,

Below one version of code that goes through without error message.

Some issues that I noticed:

-Concerning your update, prepareGradient is meant to be used only with spatial prediction. For non-spatial prediction, use constructGradient.

-The reason behind the error message is that your XData.grid has NA for rugosity for row 4610.

Best,

Otso

setwd("U:/all stuff/hmsc/hmsc tmp") load("cgoetsch_example_data.RData") library(Hmsc) library(dplyr)

community data

Y = da %>% dplyr::select(5:8)

environmental data

X = da %>% dplyr::select(depth, rugosity)

spatial coordinates of the centroids of the 100km grid squares

xy_grid = as.matrix(da %>% dplyr::select(cent_lon, cent_lat) %>% unique()) rownames(xy_grid)<-as.factor(unique(da$grid_id)) plot(xy_grid)

The model

StudyDesign with an unstructured random effect at the level of the tow and spatial random effect at level of the 100 km grid

studyDesign = data.frame( tow = as.factor(da$UniqueID), grid_id = as.factor(da$grid_id))

associates each tow with one of the 42 100km grid squares

Unstructured random effect of sample level

rL.tow = HmscRandomLevel(units = studyDesign$tow)

Spatial random effect at level of the 100 km grid

rL.spatial = HmscRandomLevel(sData=xy_grid, sMethod = 'NNGP', nNeighbours = 8)

Test model formula

XFormula = ~ depth + rugosity

The test model

studyDesign = data.frame(grid_id = as.factor(da$grid_id))

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula,

distr="probit", studyDesign=studyDesign,

ranLevels=list("tow"= rL.tow, "grid_id" = rL.spatial))

m.nngp = Hmsc(Y=Y, XData = X, XFormula=XFormula, distr="probit", studyDesign=studyDesign, ranLevels=list("grid_id" = rL.spatial))

Fit the model

set.seed(1)

testing parameters

nChains = 1 samples = 5 thin = 1 transient = 1 verbose = 1

fit the model

m.nngp = sampleMcmc(m.nngp, thin = thin, samples = samples, transient = transient, nChains = nChains, verbose = verbose, updater=list(GammaEta=FALSE))

Make spatial predictions

spatial coordinates of the prediction grid points

xy.pgrid = as.matrix(sm_grid_test %>% dplyr::select(Lon, Lat))

environmental data for the prediction grid points

XData.grid = data.frame(depth=sm_grid_test$depth, rugosity=sm_grid_test$rugosity, stringsAsFactors = TRUE)

Format data for spatial predictions

nn = 4609 xy.pgrid2 = xy.pgrid[1:nn,] XData.grid2 = XData.grid[1:nn, ] Gradient = prepareGradient(m.nngp, XDataNew = XData.grid2, sDataNew = list(grid_id=xy.pgrid2))

compute predictions

predY = predict(m.nngp, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

From: cgoetsch @.> Sent: tiistai 15. kesäkuuta 2021 21:55 To: hmsc-r/HMSC @.> Cc: Ovaskainen, Otso T @.>; Mention @.> Subject: Re: [hmsc-r/HMSC] Numerical instability in a very uneven spatial study design (#33)

Otso,

Update: Today, I tried to make predictions on a version of the same model without the spatial random effect. All of the above code is the same except the following:

model without the spatial random effect; only includes the unstructured random effect and covariates

m.nngp4 = Hmsc(Y=Y, XData = X, XFormula=XFormula,

           distr="probit", studyDesign=studyDesign,

           ranLevels=list("tow"= rL.tow))

fit model with unstructured random effect with covariates

m.nngp4 = sampleMcmc(m.nngp4, thin = thin, samples = samples, transient = transient,

                 nChains = nChains, verbose = verbose,

                 updater=list(GammaEta=FALSE))

for the model with no spatial random effect

Gradient = prepareGradient(m.nngp4, XDataNew = XData.grid)

here I left out the SDataNew argument, since there is no spatial random effect in this version of the model

make predictions of occupancy to the new data

predY = predict(m.nngp4, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

Again, I get the same error as the full model: Error in LFix + Reduce("+", LRan) : non-conformable arrays.

Am I not specifying the code correctly to get predictions on a model with an unstructured random effect and covariates, but no spatial random effect? It should be possible to make spatial predictions for a model even if there is no explicit spatial random effect, correct?

Cheers, Chandra

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/hmsc-r/HMSC/issues/33#issuecomment-861752525, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEIYMZQRVKEE4LTN2IYEEUTTS6OYVANCNFSM4JW2EO6Q.

jarioksa commented 3 years ago

@cgoetsch This peculiar error is caused by one missing value in variable rugosity in sm_grid_test (row 4610). This file is used to build the Gradient in prepareGradient function, and later in predict this case with missing value is removed from LFix but kept in LRan and – voilà we have "non-conformable arrays".

The following seems to fix this case:

kk <- complete.cases(Gradient$XDataNew)
Gradient$XDataNew <- Gradient$XDataNew[kk,]
Gradient$studyDesignNew <- Gradient$studyDesignNew[kk,]
predY = predict(m.nngp, Gradient=Gradient, expected = TRUE, predictEtaMean = TRUE)

Naturally, this should be handled more smoothly by the software. We do check against missing values in XData in setting the Hmsc model, but not in constructing gradients.

cgoetsch commented 3 years ago

Otso,

Thanks so much for your help. I should have caught that myself before bothering you with it, but I thought that all the NAs had already been expunged. Everything is working great now and I have some fantastic prediction maps.

As for making spatial predictions on a model without an explicit spatial random effect, I am a somewhat confused, so maybe you can clarify.

I am interested in making predictions for new spatial points - getting the occurrence probabilities for each species across our spatial prediction grid based on the same model - just without the spatial random effect. With constructGradient, from the book example for Finnish birds (Ch 11), it appears as if I have to select a focal variable and that the output of constructGradient is meant to feed into plotGradient to see patterns in the data by covariate independent of space. This is not what I am interested in for this question. We want to compare the predicted species distributions of model A (full model with covariates and spatial random effect) to model B (model with covariates but no spatial random effect). But getting the predicted species distribution for both versions of the model are both spatial predictions, correct? Am I misunderstanding the application of constructGradient or is it not possible to get spatial predictions in a model without a spatial random effect?

Thanks again for your time and patience, Chandra

ovaskain commented 3 years ago

Hi,

With Hmsc, predictions are always done with the predict function. The functions constructGradient & prepareGradient are just auxiliary functions that help you to prepare the predictors for special cases. If you have a non-spatial model and you have the XDataNew for the new prediction units, the easiest is probably to use directly the predict function.

Best,

Otso

From: cgoetsch @.> Sent: keskiviikko 16. kesäkuuta 2021 19:00 To: hmsc-r/HMSC @.> Cc: Ovaskainen, Otso T @.>; Mention @.> Subject: Re: [hmsc-r/HMSC] Numerical instability in a very uneven spatial study design (#33)

Otso,

Thanks so much for your help. I should have caught that myself before bothering you with it, but I thought that all the NAs had already been expunged. Everything is working great now and I have some fantastic prediction maps.

As for making spatial predictions on a model without an explicit spatial random effect, I am a somewhat confused, so maybe you can clarify.

I am interested in making predictions for new spatial points - getting the occurrence probabilities for each species across our spatial prediction grid based on the same model - just without the spatial random effect. With constructGradient, from the book example for Finnish birds (Ch 11), it appears as if I have to select a focal variable and that the output of constructGradient is meant to feed into plotGradient to see patterns in the data by covariate independent of space. This is not what I am interested in for this question. We want to compare the predicted species distributions of model A (full model with covariates and spatial random effect) to model B (model with covariates but no spatial random effect). But getting the predicted species distribution for both versions of the model are both spatial predictions, correct? Am I misunderstanding the application of constructGradient or is it not possible to get spatial predictions in a model without a spatial random effect?

Thanks again for your time and patience, Chandra

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/hmsc-r/HMSC/issues/33#issuecomment-862503444, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEIYMZQGMNGGMJ36UQ7V7BTTTDDBNANCNFSM4JW2EO6Q.

jarioksa commented 3 years ago

@ovaskain The canonical name (in R) for new data in predict is newdata: currently the documentation is not very clear for the meaning of XData. It seems that predict suffers from the same NA-in-XData problem with the cryptic error message of "non-conformable arrays". I'll add tests against NA in predict.Hmsc, too.

jarioksa commented 3 years ago

This issue was started long time ago and has rambled far away from the original message. Please post here only comments that are directly related to the first message, and report all other problems as new issues so that they won't be buried deep in sediments.