Heteroskedastic noise ( again ! )

[Sébastien Coube]

Hello, 

I am coming back to you in order to ask (again !) how to model spatial heteroskedastic noise using INLA SPDE. The R code below shows what I mean; the objective is to estimate noise_variance and spatial_field. 

A few years ago, I have been in this forum to ask the same question. The answer that was kindly suggested to me was to run an INLA model with two Gaussian SPDEs : one for the Matérn field, one for the Gaussian noise field. The two SPDEs were differring in the size of their A-matrices, the noise SPDE having a diagonal A-matrix. The solution did work, but encountered difficulties when the number of observations and the number of variables explaining the noise variance augmented. 

I would like to try again, and to be sure that I am not missing a new method that has been introduced in INLA since the last time. 

Again, thank you all for the terrific work you are doing ! 

Kind regards

Sébastien Coube-Sisqueille

# spatial locations
locs = as.matrix(expand.grid(seq(1e2), seq(1e2)))

# Gaussian field explaining the mean
spatial_field = matrix(rnorm(1e4), 1e2)
sqexp_cov_chol = exp(-as.matrix((dist(seq(1e2), diag = T)/10)^2))
diag(sqexp_cov_chol) = 1.001
sqexp_cov_chol = chol(sqexp_cov_chol)
spatial_field = t(sqexp_cov_chol) %*% spatial_field %*% sqexp_cov_chol
plot(
  locs,
  col = heat.colors(100)[cut(spatial_field, seq(min(spatial_field), max(spatial_field), length.out = 101))],
  pch = 15,
  main = "spatial field"
     )

# Gaussian field explaining the mean
noise_variance = matrix(rnorm(1e4), 1e2)
noise_variance = - 2 + t(sqexp_cov_chol) %*% noise_variance %*% sqexp_cov_chol
noise_variance = exp(noise_variance)
noise_realization = sqrt(noise_variance) * rnorm(1e4)
plot(
  locs,
  col = heat.colors(100)[cut(noise_variance, seq(min(noise_variance), max(noise_variance), length.out = 101))],
  pch = 15,
  main = "variance of the noise"
     )
plot(
  locs,
  col = heat.colors(100)[cut(noise_realization, seq(min(noise_realization), max(noise_realization), length.out = 101))],
  pch = 15,
  main = "realization of the noise"
     )

# combining spatial field and noise
observed_field = spatial_field + noise_realization
plot(
  locs,
  col = heat.colors(100)[cut(observed_field, seq(min(observed_field), max(observed_field), length.out = 101))],
  pch = 15,
  main = "observed field"
     )

[Helpdesk (Haavard Rue)]

if one of the models (mean and log variance) is simple, then there is a new
likelihood

inla.doc("ggaussian")

you might be able to use ?

> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/dc9d7fac-0924-4054-bf36-6c626be6251bn%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Sébastien Coube]

Hello Haavard, 

I will definitely look into this. THank you very much !

Sébastien

[Sébastien Coube]

Dear Haavard, and rest of the INLA community, 

I have had a look into the ggaussian family of INLA, and it does what I have in mind. However, I believe that the number of predictors for the noise variance is limited to 10. Is there a way to loosen that requirement ? Ans subsidiarily, is it possible to pass matrices to inla.mdata() instead of individual vectors ?

Again, as usual : thank you all for your amazing work and your help with the software. 

Kind regards

Seb

n <- 1000
x <- rnorm(n)
xx <- rnorm(n)
eta <- 0.1 + 1.1 * x + 2.2 * xx
s <- rep(1, n)
z <- rnorm(n)
zz <- rnorm(n)

# try with 9 +1 predictors for the noise
X_noise = matrix(rnorm(n * 9), n)
beta_noise = rnorm(ncol(X_noise))
eta.prec <- -3 + X_noise %*% beta_noise
y <- eta + 1/sqrt(s * exp(eta.prec)) * rnorm(n)
Y <- inla.mdata(y, s, 1,
                X_noise[,1],
                X_noise[,2],
                X_noise[,3],
                X_noise[,4],
                X_noise[,5],
                X_noise[,6],
                X_noise[,7],
                X_noise[,8],
                X_noise[,9]
                )
r <- inla(Y ~ 1 + x + xx,
          data = list(Y = Y, x = x, xx = xx, z = z, zz = zz, s = s),
          family = "ggaussian")
summa = summary(r)
plot(
summa$hyperpar[-1,1],
beta_noise
)
abline(a=0, b=1)


# try with 10 +1 predictors for the noise
X_noise = matrix(rnorm(n * 10), n)
beta_noise = rnorm(ncol(X_noise))
eta.prec <- -3 + X_noise %*% beta_noise
y <- eta + 1/sqrt(s * exp(eta.prec)) * rnorm(n)
Y <- inla.mdata(y, s, 1,
                X_noise[,1],
                X_noise[,2],
                X_noise[,3],
                X_noise[,4],
                X_noise[,5],
                X_noise[,6],
                X_noise[,7],
                X_noise[,8],
                X_noise[,9],
                X_noise[,10])

r <- inla(Y ~ 1 + x + xx,
          data = list(Y = Y, x = x, xx = xx, z = z, zz = zz, s = s),
          family = "ggaussian")
summa = summary(r)
plot(
summa$hyperpar[-1,1],
beta_noise
)
abline(a=0, b=1)

[Helpdesk (Haavard Rue)]

yes, you can use

Y <- inla.mdata(y, s, 1, X_noise)

but it will show individuals names in the output, one for each column of X_noise

yes, there is currently a restriction of 10. it seems rare that more is
required... ?

of'course this could be increased, but just adding more and more will just add
more problems elsewhere, so...

[Sébastien Coube]

Hello Haavard, 

I was on holidays so I did not answer sooner. Thank you for your answer. 

I was thinking about including a few dozen of spatial basis functions to model the precision of the Gaussian observations. I have some evidence that this can really help with spatial data sets. 

Do you think that this is doable in the current framework ?

Kind regards !

Sébastien

[Havard Rue]

You can.  Just keep them centered and scaled 

-- 

Håvard Rue 

hr...@r-inla.org

You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/27b2bed7-ab64-481d-81e0-e140d7ea0336n%40googlegroups.com.