constructGradient error and conditional estimates of species associations

[jiaangou]

Hello,

I am experimenting with the current Hmsc R package (from CRAN). I have a time-series (6 time points) abundance data set of terrestrial arthropods (100spp), 14 plots (w/ spatial structure), 102 sampling units, and 1 categorical trait variable (9 levels). I'm examining how species respond to 2 covariates: (1) time/julian date (continuous) and (2) Farm type (categorical with 2 levels: O & C).

My model structure is as follows:

hmsc_mod = Hmsc(Y = sp_comp, XData = XData, XFormula = XFormula,
                TrData = trait_data, TrFormula = TrFormula,
                studyDesign = studyDesign, 
                ranLevels = list(farm = randPlot, sample = randSample, space = randSpace),
                distr = 'poisson')

At the moment, I have 2 questions:

1) Error message from constructGradient()

When I try to use the constructGradient() function to create a new dataframe for prediction, I get the follow error message:

 Gradient <- constructGradient(output, focalVariable = "julian_date",
        non.focalVariables = list('farm_type'=list(3,"O")))
Error in `[.data.frame`(hM$XData, , vars[i]) : undefined columns selected
> traceback()
5: stop("undefined columns selected")
4: `[.data.frame`(hM$XData, , vars[i])
3: hM$XData[, vars[i]]
2: is.factor(hM$XData[, vars[i]])
1: constructGradient(output, focalVariable = "julian_date", non.focalVariables = list(farm_type = list(3, "O")))

I tried running the source code directly but still could not find where the issue is.

2) Covariate-dependent species-to-species associations

I am aware that the Matlab version has this functionality (Tikhonov et al 2017). I'm wondering if there is a way to do this in the R package.

Any help on either of the two questions would be much appreciated. Thanks!

[ovaskain]

Hi,

Concerning 1, see if all.vars(XFormula) gives extra variables on top of those covariates that you actually have. For example, if you have second order polynomial for which you say "raw=T", the extra variable "T" may appear. This can be fixed by "raw=TRUE". This was the source of that same error message in one previous case.

Concerning 2, this is implemented, with the xData option for defining random effects.

Otso

[jarioksa]

Also check the spelling of your variables – including case. For instance, farm_type and Farm_type would be different variables and only one of these will be found in the data frame. We do not have your data and cannot reconstruct your problem locally. Without a reproducible example, all we can do is to guess. The variable names (especially for the focalVariable) you use in constructGradient() should be found in names(output$XData) or colnames(output$XData) where output is your private result data that you used in your constructGradient() command.

[jiaangou]

Thanks for the responses Otso and Jari.

Hi, Concerning 1, see if all.vars(XFormula) gives extra variables on top of those covariates that you actually have. For example, if you have second order polynomial for which you say "raw=T", the extra variable "T" may appear. This can be fixed by "raw=TRUE".

Changing "raw=T" to "raw=TRUE" fixed the problem!

[jiaangou]

Hi,

I tried constructing a covariate-dependent random effect using the xData formulation. The model constructs fine with no error but once I run MCMC, I keep getting this error message:

Error in (Eta[[r]][Pi[, r], ] rL[[r]]$x[as.character(dfPi[, r]), k]) %% : non-conformable arguments In addition: Warning message: In Ops.factor(Eta[[r]][Pi[, r], ], rL[[r]]$x[as.character(dfPi[, : ‘*’ not meaningful for factors

This isn't just happening on my dataset, but also on the simulated "TD" dataset within the Hmsc package.

rL = HmscRandomLevel(xData=data.frame(x1=rep(1,length(TD$X$x1)),x2=TD$X$x2))

test_hmsc = Hmsc(Y = as.matrix(TD$Y),XData = TD$X, XFormula = ~x1+x2, TrData = TD$Tr, TrFormula = ~T1+T2, studyDesign = TD$studyDesign, ranLevels = list('sample' = rL), distr = 'poisson')

test_output = sampleMcmc(test_hmsc,samples=samples,transient=transient,thin=thin, verbose=verbose,nChains = 2)

Any idea why this might be? Thanks!

[gtikhonov]

Hi @jiaangou,

Please have a look at the other issue #31, where we point out some current limitations of the package regarding the covariate-dependent species associations.

[aquila06]

Hi, Concerning 1, see if all.vars(XFormula) gives extra variables on top of those covariates that you actually have. For example, if you have second order polynomial for which you say "raw=T", the extra variable "T" may appear. This can be fixed by "raw=TRUE". This was the source of that same error message in one previous case. Concerning 2, this is implemented, with the xData option for defining random effects. Otso

Hi,

I encountered the same issue and I suspect this means that the model has to be fitted again with the raw=TRUE, instead of raw=T, command in the formula ? Does this only affect variable names or other aspects of the model fit as well?

Thanks in advance for any hint.

[ovaskain]

Examine if m$X has the first and second order terms of the predictor as they should be in the raw order (e.g. first order = 5 and second order = 25). If yes, the model has been fitted correctly, and you can just re-define the XFormula by m$XFormula = ...

On 7.5.2020 17.13, aquila06 wrote:

Hi, Concerning 1, see if all.vars(XFormula) gives extra variables
on top of those covariates that you actually have. For example, if
you have second order polynomial for which you say "raw=T", the
extra variable "T" may appear. This can be fixed by "raw=TRUE".
This was the source of that same error message in one previous
case. Concerning 2, this is implemented, with the xData option for
defining random effects. Otso

Hi,

I encountered the same issue and I suspect this means that the model has to be fitted again with the raw=TRUE, instead of raw=T, command in the formula ? Does this only affect variable names or other aspects of the model fit as well?

Thanks in advance for any hint.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/hmsc-r/HMSC/issues/29#issuecomment-625281030, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEIYMZXJXGCVAZPLN5E3IWDRQK62BANCNFSM4JI7M35A.

[aquila06]

Thanks for your answer