spatial correlation

[AlainZuur]

Hello,

It is very nice that gllvm can now do spatial correlation. Thanks!

I have been playing around with it, and it all seems to work. What I would like to do is mimic R-INLA. To be more specific, I want to extract the spatial random effects, interpolate them, and plot them with geom_raster in a ggplot. To do that properly, I need to get the spatial parameters, spatial random effects, and its covariance matrix. Then I can utilise the exponential correlation function to interpolate.

I will start very simple and use the following model:

Y_i ~ Poisson(mu_i)
E[Y_i] = mu_i = exp(Intercept + beta * Covariate + u_i)
u_i ~ N(0, SIGMA)

i is a site index, with more than 200 nicely placed sites. Y_i is univariate, which means that I can compare gllvm results to spaMM, glmmTMB and R-INLA.

To keep things simple, I used a corExp on the row effect in gllvm.

CoordData <- data.frame(Xkm = MI2$Xkm, Ykm = MI2$Ykm)
CoordData <- as.matrix(CoordData)
GroupStruc <- data.frame(Group = 1:nrow(MI2))

TotAbund <- MI2$TotAbund
M1 <- gllvm(y = TotAbund,
X = CovBenthic,
studyDesign = GroupStruc,
row.eff = ~corExp(1 | Group),
dist = list(CoordData),
formula = ~ Temp.s,
num.lv = 0,
family = "poisson")

I have executed the same model in spaMM:

library(spaMM)
M2 <- fitme(TotAbund ~ 1 + Temp.s + Matern(1 | Xkm + Ykm),
family = "poisson",
fixed = list(nu = 0.5),
data = MI2)

where the Matern with fixed nu = 0.5 equals the exponential correlation.

And this is glmmTMB:

MI2$Pos <- numFactor(MI2$Xkm, MI2$Ykm)
MI2$ID <- factor(rep(1, nrow(MI2)))
M3 <- glmmTMB(TotAbund ~ 1 + Temp.s + exp(0 + Pos | ID),
family = "poisson",
data = MI2)

I won't show the INLA code as it is rather long.

The covariance matrix SIGMA is of the form:
sigma^2 if distances between sites are 0
sigma^2 * exp(-d / phi)

where d is the distance between sites. I have also seen sigma^2 * exp(-3 d / phi) in the literature.

Here is the output from gllvm:

M1$params$sigma
Scale|Group Group|Group
0.000209588 0.308296818

The Group|Group is the sigma from the "sigma^2 * exp(-d / alpha)" and that matches the results of spaMM, glmmTMB, INLA.

My question is: What exactly is the Scale|Group parameter? Both INLA and spaMM give a range of about 5. But I can't relate the 000209588 to 5. Or better formulated, how do I go from 000209588 to the phi?

Estimated intercepts by all packages are similar, but the slope obtained by gllvm differs. Random effects obtained by gllvm and spamMM are rather similar.

Kind regards,

Alain

[AlainZuur]

I managed to get a little deeper in this. I changed the data set to something with a stronger spatial correlation. Now I get an estimated intercept and slope from INLA, spaMM, glmmTMB and gllvm, that are all comparable.

I am still struggling with this:

M1$params$sigma
Scale|Group Group|Group
0.03876086 0.23119993

The gllvm model has an exponential correlation on the row effect. All packages give a similar sigma_u.

The covariance matrix is given by: sigma_u^2 * exp(- distance / phi)

I managed to find your corExp function in distrib.h. I am slightly confused why you call this function with Type(0) (that is a 0?). Anyway....I copied your corExp code into ChatGPT (Pro mode). It tells me that you are using alf = 1/exp(0) = 1 in the corExp function. Hence, I guess the s1 in your code is being estimated somewhere else.

It tells me that I have to use:

ϕ = exp(s1/2) where the s1 is the reported Scale|Group. Is this correct?

I realise that this is not the range. INLA's reported range is the distance when the correlation is 0.1. That translates to:

Corr(d) = 0.1 = exp(−d / ϕ).

Which gives a range of -ϕ * ln(0.1). But using ϕ=exp(s1/2) with s1 the reported scale doesn't match the range from spaMM/INLA.

In other words....I am missing the link between the reported scale and the phi in the exponential correlation.

Alain

PS....Matern correlation is next...:-)

[BertvanderVeen]

The function is called with Type(0), where Type(0) is just a 0 of TMB "type", because the coordinates are already scaled by the parameter in the correlation function before hand (admittedly, slightly confusing and not well documented in the code). As log_s1 = 0, and exp(log_s1) = 1, we thus ignore s1 in the corExp function.

Since the matrix of coordinates is first scaled by 1/Scale|Group before calculating euclidean distances, this is equivalent to calculating the euclidean distances of the unscaled coordinates and dividing by Scale|Group in the correlation function. This makes for a parameterization where Scale|Group is the range that you are looking for.

To summarize, the covariance matrix is given by: 1|Group^2 * exp(- distance / Scale|Group)

[AlainZuur]

Thanks. That is what I needed to know. I will simulate some univariate data based on a covariance matrix 1|Group^2 * exp(- distance / Scale|Group), apply gllvm/spaMM/Friends, and see how similar the correlation is (quantified via Scale|Group and 1|Group).

Alain

[AlainZuur]

Hello Bert, Are you sure that your comment above from 19 December is correct?

I think that the reported Scale|Group (range) is with respect to the scaled distances. My coordinates are in km..... your code scales them via:

}else if((cstruc(re) == 4) || (cstruc(re) == 2)){ // corMatern, corExp
// Distance matrix calculated from the coordinates for rows
matrix DiSc(dc(dccounter).cols(),dc(dccounter).cols()); DiSc.fill(0.0);
matrix dc_scaled(dc(dccounter).rows(),dc(dccounter).cols()); dc_scaled.fill(0.0);
DiSc.setZero();
DiSc.diagonal().array() += 1/sigma(sigmacounter);
sigmacounter++;
dc_scaled = dc(dccounter)*DiSc;
if(cstruc(re)==2){ // corExp
Sr = gllvm::corExp(sigma(sigmacounter), Type(0), nr(re), dc_scaled);
sigmacounter++;
} else if(cstruc(re)==4) { // corMatern

If I understand this code well, then it uses:

DiSc = 1/sigma;
dc_scaled = dc * DiSc;

So...the coordinates dc_scaled are scaled by dividing by sigma (which makes them bigger when sigma<1), where sigma is $params$sigma["Group|Group"].

The dc_scaled and sigma then go into corExp where you calculate this:

S(d,j)=s0*exp(-sqrt( (((dc.row(d)-dc.row(j))alf)(dc.row(d)-dc.row(j)).transpose()).sum() ) )*s0;//

where the dc.row are the scaled coordinates.

So..you scale by ["Group|Group"].

I am now trying to express the range in the original units.

With covariance = (Group|Group)^2 * exp(- scaled distance / (Scale|Group) ) I would have thought that

covariance = (Group|Group)^2 * exp(- ((1 / sigma) * distance) / (Scale|Group) ) =
(Group|Group)^2 * exp(- distance / ( Group|Group * Scale|Group ))

In other words: corrected range (or phi in the expression above) = reported Group|Group * reported Scale|Group

But that value does not make sense for the distances I am working with. It is about a factor 1000 or so too small.

What am I missing? Are you dividing the coordinates somewhere by 1000?
Alain

[BertvanderVeen]

dc is the coordinates object on the C++ side. These are not divided by anything elsewhere. You can verify this by model$TMBfn$env$data$dc
dc_scaled = dc(dccounter)*DiSc; scales the coordinates by the inverse of sigma, same as dividing by it on the R-side.
In corExp this is again scaled (on the outside of exp) by the square of group|group

i.e., the same as I wrote in my previous post

[AlainZuur]

Bert....it must have been the data set that I used. I have simulated 1 response variable Y with corExp correlated noise. And then I applied gllvm with a row-level spatial correlated term (also corExp).

The estimated range and sigma obtained by gllvm are comparable to what I used to simulate the data. Hence, the software does what it is supposed to do.

Thanks.
Alain