Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error in hM$rL[[r]]$s[levels(hM$dfPi[, r]), ] : subscript out of bounds #47

Closed
rhettdharrison opened this issue Jun 17, 2020 · 8 comments
Labels
help wanted Extra attention is needed

Comments

@rhettdharrison
Copy link

Hello again,

I am trying to run a model with spatial latent effects following Vignette 3 and the Manual, but getting this error. If I remove the spatial effect then the error disappears. I've checked the xycoords and it is fine - as far as I can see (I plotted it and checked the values). Below the error is the model code in full (and attached data). The xycoords are for plot locations. There are 2 samples per plot (at 1 m and at 18 m height).

m.pa.spatial = sampleMcmc(m.pa.spatial, thin = thin, samples = samples, transient = transient,

  •               nChains = nChains, verbose = verbose)
    

Error in hM$rL[[r]]$s[levels(hM$dfPi[, r]), ] : subscript out of bounds

traceback()
2: computeDataParameters(hM)
1: sampleMcmc(m.pa.spatial, thin = thin, samples = samples, transient = transient,
nChains = nChains, verbose = verbose)

xycoords <- matrix(c(plot.char$UTM_E,plot.char$UTM_N),ncol=2)
colnames(xycoords) <- c("UTM_E","UTM_N")
rownames(xycoords) <- plot.id
xycoords <- apply(xycoords,2, function(x) sapply(x, function(m) m - min(x)))

studyDesign = data.frame(plot = plot.id, sample = sample.id)
rL.spatial = HmscRandomLevel(sData = xycoords)
rL.spatial = setPriors(rL.spatial,nfMin=1,nfMax=1)
rL2 = HmscRandomLevel(units = studyDesign$sample)

#set presence-absence model
m.pa.spatial <- Hmsc(Y = sign(lh_sp), XData = plot.char, XFormula = ~ HEIGHT + ELEV + SLOPE + LAI, distr = "probit", studyDesign = studyDesign, ranLevels = list("plot" = rL.spatial))

plot.char.txt
lh_sp.txt

@ovaskain
Copy link
Collaborator

ovaskain commented Jun 17, 2020 via email

@rhettdharrison
Copy link
Author

rhettdharrison commented Jun 18, 2020

That is correct. However, I have tried re-specifying the model and cannot get it to work.

If I correct the xycoords by removing the repeats, I get the same error. If, alternatively, I assign the sample.id to the rownames of xycoords and remove the plot level from the design, I get the following error (I rechecked the data).

Error in chol.default(W) : 
  the leading minor of order 2 is not positive definite
> traceback()
4: chol.default(W)
3: chol(W)
2: computeDataParameters(hM)
1: sampleMcmc(m.ab.spatial, thin = thin, samples = samples, transient = transient, 
       nChains = nChains, verbose = verbose)

xycoords <- matrix(c(plot.char$UTM_E,plot.char$UTM_N),ncol=2)
colnames(xycoords) <- c("UTM_E","UTM_N")
rownames(xycoords) <- sample.id
xycoords <- apply(xycoords,2, function(x) sapply(x, function(m) m - min(x)))

studyDesign = data.frame(sample = sample.id)
rL.spatial = HmscRandomLevel(sData = xycoords)
rL.spatial = setPriors(rL.spatial,nfMin=1,nfMax=1)

#set presence-absence model
m.pa.spatial <- Hmsc(Y = sign(lh_sp), XData = plot.char, XFormula = ~ HEIGHT + ELEV + SLOPE + LAI,
             distr = "probit", studyDesign = studyDesign,
             ranLevels = list("sample" = rL.spatial))

@jarioksa
Copy link
Collaborator

I have tried to reproduce this, but failed. In other words: the commands work when I try them and I cannot trigger an error. So this is guessing: the error message you get is the same as earlier with chol (leading minor of order 2 is not positive definite) when you had non-numeric data. In this case the error comes from chol() of spatial weight matrix W which you find from your spatial coordinates. So check these first. In this case, I think you should first see what you get from dist(xycoords) and then how data are read into Hmsc. With your model the spatial coordinates should be in m.pa.spatial$rL["sample"]$s, and the spatial weights are W = dist(m.pa.spatial$rL["sample"]$s).

This is blind guessing though, as we do not have a reproducible example. Can you work out such an example with, e.g., TD data in Hmsc?

@rhettdharrison
Copy link
Author

The data in xycoords was fine (plotted and checked the sum), but obviously coordinates of each of the pair of samples taken from the same plot were identical, which produced zeros in the distance matrix. I am not sure why that creates a problem, but adding 5 m to the second x of each pair coodinates (a trivial change when the plots are 100s to kms apart) corrected the error.

So it seems that it is not possible to specify a spatial model where multiple samples are taken from the sample sites, at least not without jittering the points.

@ovaskain
Copy link
Collaborator

ovaskain commented Jun 19, 2020 via email

@jarioksa jarioksa added the help wanted Extra attention is needed label Jun 19, 2020
@ewmorr
Copy link

ewmorr commented Jun 19, 2020

Thanks for this explanation. This answers my questions in #49. I'll update that issue by way of comment and close.

@rhettdharrison
Copy link
Author

Thanks from here too. Now very clear.

@smithja16
Copy link

Sorry about posting on a closed thread, but I wanted to post my own experience with this error (using Hmsc v3.0-13). It arose when I was doing my own form of cross-validation (i.e. fitting a model to a data subset).

The solution was to use droplevels() to ensure that the data subset did not reference levels of a random effect (here, 'Vessel') missing in that subset:
studyDesign = data.frame(sample=as.factor(1:nrow(train_data)), vessel=droplevels(train_data$Vessel))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
help wanted Extra attention is needed
Projects
None yet
Development

No branches or pull requests

5 participants