Ngme2 - A new Flexible R Package for Latent non-Gaussian Models

2025-03-25

Source: vignettes/ngme2.Rmd

Introduction

In this vignette we provide a brief introduction to the ngme2 package.

Ngme2 (https://github.com/davidbolin/ngme2) is the updated version of Ngme, a package for estimating latent non-Gaussian models for repeated measurement data. Ngme2 follows a hierachical structure, differnet components (latent processes, different types of noises) are flexible to change and combine.

1 Features

  1. Support temporal models like AR(1), Ornstein−Uhlenbeck and random walk processes, and spatial models like Matern fields.
  2. Support latent processes constructed by non-Gaussian noises (normal inverse Gaussian(NIG), generalized asymmetric Laplace (GAL)).
  3. Support non-Gaussian, and correlated measurement noises.
  4. Support doing prediction at unknown locations.
  5. Support latent processes and random-effects model for longitudinal data.
  6. Support the bivariate type-G model, which can model 2 non-Gaussian fields jointly (Bolin 2020).
  7. Support the separable space-time model.

2 Model Framework

The package Ngme2 provides methods for mixed effect models in the following form:

$$ {\bf Y}_{ij} = {\bf X}^T_{ij} {\bf \beta} + {\bf D}^T_{ij} {\bf U}_i + W_i(t_{ij}) + \epsilon_{ij}, \qquad j=1 \ldots n_i, i=1,\ldots,m $$

  • mm is the number of subjects, nin_i is the number of observations for each subject,
  • YY is the response variable,
  • ${\bf X}$ is the matrix of fixed effects explanatory variables,
  • ${\bf \beta}$ is the fixed effects,
  • ${\bf D}$ is the matrix of random effects explanatory variables,
  • $\bf U$ is the random effects,
  • Wi(tij)W_i(t_{ij}) is a stochastic process driven by Gaussian or non-Gaussian noise,
  • ϵ\epsilon is measurement error.

Here is a simple template for using the core function ngme to model the single response:

ngme(
  formula=Y ~ x1 + x2 + f(index, model="ar", noise="nig"),
  data=data.frame(Y=Y, x1=x1, x2=x2, index=index),
  noise = noise_normal()
)

Here, function f is for modeling the stochastic process W with Gaussian or non-Gaussian noise, we will discuss this later. noise stands for the measurement noise distribution. In this case, the model will have a Gaussian likelihood.

3 Non-Gaussian Model

Here we assume the non-Gaussian process is a type-G Lévy process, whose increments can be represented as location-scale mixtures: γ+μV+σVZ,\gamma + \mu V + \sigma \sqrt{V}Z, where γ,μ,σ\gamma, \mu, \sigma are parameters, Z∼N(0,1)Z\sim N(0,1) is independent of VV, and VV is a positive infinitely divisible random variable. This results in the following form, where KK is the operator part:

KW|V∼N(γ+μV,σ2diag(V)), KW|V \sim N(\gamma + \mu V, \sigma^2 \, \text{diag}(V)), where μ\mu and σ\sigma can be non-stationary.

One example in ngme2 is the normal inverse Gaussian (NIG) noise, where VV follows an Inverse Gaussian distribution with parameter ν\nu (IG(ν\nu, ν\nu)).

A random variable VV follows an inverse Gaussian distribution with parameters η1\eta_1 and η2\eta_2, denoted by V∼IG(η1,η2)V\sim \text{IG}(\eta_1,\eta_2), if it has probability density function (pdf) given by π(v)=η22πv3exp{−η12v−η22v+η1η2},η1,η2>0.\pi(v) = \frac{\sqrt{\eta_2}}{\sqrt{2\pi v^3}} \exp\left\{-\frac{\eta_1}{2}v - \frac{\eta_2}{2v} + \sqrt{\eta_1\eta_2}\right\},\quad \eta_1,\eta_2>0. We can generate samples from an inverse Gaussian distribution with parameters η1\eta_1 and η2\eta_2 by generating samples from the generalized inverse Gaussian distribution with parameters p=−1/2p=-1/2, a=η1a=\eta_1 and b=η2b=\eta_2. The rGIG function can be used to generate samples from the generalized inverse Gaussian distribution.

If V∼IG(η1,η2)V\sim \text{IG}(\eta_1,\eta_2), and X=γ+μV+σVZX = \gamma +\mu V + \sigma \sqrt{V}Z, with Z∼N(0,1)Z\sim N(0,1) independent of VV, then XX follows a normal inverse Gaussian (NIG) distribution with pdf π(x)=eη1η2+μ(x−γ)/σ2η2μ2/σ2+η1η2πη2σ2+(x−γ)2K1((η2σ2+(x−γ)2)(μ2/σ4+η1/σ2)),\pi(x) = \frac{e^{\sqrt{\eta_1\eta_2}+\mu(x-\gamma)/\sigma^2}\sqrt{\eta_2\mu^2/\sigma^2+\eta_1\eta_2}}{\pi\sqrt{\eta_2\sigma^2+(x-\gamma)^2}} K_1\left(\sqrt{(\eta_2\sigma^2+(x-\gamma)^2)(\mu^2/\sigma^4+\eta_1/\sigma^2)}\right), where K1K_1 is a modified Bessel function of the third kind. In this form, the NIG density is overparameterized, so we set η1=η2=η\eta_1=\eta_2=\eta, which results in E(V)=1E(V)=1. Thus, we have the parameters μ\mu, γ\gamma, and η\eta.

The NIG model assumes that the stochastic variance ViV_i follows an inverse Gaussian distribution with parameters η\eta and ηhi2\eta h_i^2, where hi=∫𝒟φi(𝐬)d𝐬.h_i = \int_{\mathcal{D}} \varphi_i(\mathbf{s}) d\mathbf{s}.

library(fmesher)
library(ngme2)
#> This is ngme2 of version 0.6.0
#> - See our homepage: https://davidbolin.github.io/ngme2 for more details.
library(ggplot2)
library(plyr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:plyr':
#> 
#>     arrange, count, desc, failwith, id, mutate, rename, summarise,
#>     summarize
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(viridis)
#> Loading required package: viridisLite

4 Parameter Estimation

  1. Ngme2 does maximum likelihood estimation through preconditioned stochastic gradient descent.
  2. Multiple chains are run in parallel for better convergence checks.

See Model estimation and prediction for more details.

Ngme Model Structure

Specify the driven noise

There are 2 types of common noise involved in the model, one is the innovation noise of a stochastic process, one is the measurement noise of the observations. They can be both specified by noise_<type> function.

For now we support normal, NIG, and GAL noises.

The R class ngme_noise has the following interface:

noise_normal(sigma = 1)              # normal noise
#> Noise type: NORMAL
#> Noise parameters: 
#>     sigma = 1
noise_nig(mu = 1, sigma = 2, nu = 1) # nig noise
#> Noise type: NIG
#> Noise parameters: 
#>     mu = 1
#>     sigma = 2
#>     nu = 1
noise_nig(            # non-stationary nig noise
  B_mu=matrix(c(1:10), ncol=2),
  theta_mu = c(1, 2),
  B_sigma=matrix(c(1:10), ncol=2),
  theta_sigma = c(1,2),
  nu = 1)
#> Noise type: NIG
#> Noise parameters: 
#>     theta_mu = 1, 2
#>     theta_sigma = 1, 2
#>     nu = 1

The 3rd example is the non-stationary NIG noise, where $\mu = \bf B_{\mu} \bf \theta_{\mu}$, and $\sigma = \exp(\bf B_{\sigma} \bf \theta_{\sigma})$.

ngme_noise(
  type,           # the type of noise
  theta_mu,       # mu parameter
  theta_sigma,    # sigma parameter
  nu,        # nu parameter
  B_mu,           # basis matrix for non-stationary mu
  B_sigma         # basis matrix for non-stationary sigma
)

It will construct the following noise structure:

−𝛍+𝛍V+𝛔VZ - \mathbf{\mu} + \mathbf{\mu} V + \mathbf{\sigma} \sqrt{V} Z

where $\mu = \bf B_{\mu} \bf \theta_{\mu}$, and $\sigma = \exp(\bf B_{\sigma} \bf \theta_{\sigma})$. In this case, we can recover gaussian noise by setting type=“normal and ignoring theta_mu and nu. Or we can simply use helper function noise_normal(sd=1).

Specify stochastic process with f function

The middle layer is the stochastic process, in R interface, it is represented as a f function. The process can be specified by different noise structure. See ?ngme_model_types() for more details.

Some examples of using f function to specify ngme_model:

ngme2::f(1:10, model = "ar1", noise = noise_nig())
#> Model type: AR(1)
#>     rho = 0
#> Noise type: NIG
#> Noise parameters: 
#>     mu = 0
#>     sigma = 1
#>     nu = 1

One useful model would be the SPDE model with Gaussian or non-Gaussian noise, see the vignette for details.

Specifying latent models with formula in ngme

The latent model can be specified additively as a formula argument in ngme function together with fixed effects.

We use R formula to specify the latent model. We can specify the model using f within the formula.

For example, the following formula

formula <- Y ~ x1 + f(
    x2,
    model = "ar1",
    noise = noise_nig(),
    theta_K = 0.5
  ) + f(1:5,
    model = "rw1",
    circular = T,
    noise = noise_normal()
  )

corresponds to the model

Y=β0+β1x1+W1(x2)+W2(x3)+ϵ, Y = \beta_0 + \beta_1 x_1 + W_1(x_2) + W_2(x_3) + \epsilon, where W1W_1 is an AR(1) process, W2W_2 is a random walk 1 process. x2x_2 is random effects.. . By default, we have intercept. The distribution of the measurement error ϵ\epsilon is given in the ngme function.

The entire model can be fitted, along with the specification of the distribution of the measurement error through the ngme function:

ngme(
  formula = formula,
  family = noise_normal(sigma = 0.5),
  data = data.frame(Y = 1:5, x1 = 2:6, x2 = 3:7),
  control_opt = control_opt(
    estimation = FALSE
  )
)
#> *** Ngme object ***
#> 
#> Fixed effects: 
#> (Intercept)          x1 
#>       7.365      -0.869 
#> 
#> Models: 
#> $field1
#>   Model type: AR(1)
#>       rho = 0
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = 0
#>       sigma = 1
#>       nu = 1
#> 
#> $field2
#>   Model type: Random walk (order 1)
#>       No parameter.
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 1
#> 
#> Measurement noise: 
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 0.5
#> 
#> 
#> Number of replicates is  1

It gives the ngme object, which has three parts:

  1. Fixed effects (intercept and x1)
  2. Measurement noise (normal noise)
  3. Latent models (contains 2 models, ar1 and rw1)

We can turn the estimation = TRUE to start estimating the model.

A simple example - AR1 process with nig noise

Now let’s see an example of an AR1 process with nig noise. The process is defined as

Wi=ρWi−1+ϵi, W_i = \rho W_{i-1} + \epsilon_i, Here, ϵ1,..,ϵn\epsilon_1, ..,\epsilon_n is the iid NIG noise. And, it is easy to verify that $$ K{\bf W} = \boldsymbol\epsilon,$$ where K=[1−ρ2−ρ1⋱⋱−ρ1] K = \begin{bmatrix} \sqrt{1-\rho^2} \\ -\rho & 1 \\ & \ddots & \ddots \\ & & -\rho & 1 \end{bmatrix} 

n_obs <- 500
sigma_eps <- 0.5
alpha <- 0.5
mu = 2; delta = -mu
sigma <- 3
nu <- 1

# First we generate V. V_i follows inverse Gaussian distribution
trueV <- ngme2::rig(n_obs, nu, nu, seed = 10)

# Then generate the nig noise
mynoise <- delta + mu*trueV + sigma * sqrt(trueV) * rnorm(n_obs)
trueW <- Reduce(function(x,y){y + alpha*x}, mynoise, accumulate = T)
Y = trueW + rnorm(n_obs, mean=0, sd=sigma_eps)

# Add some fixed effects
x1 = runif(n_obs)
x2 = rexp(n_obs)
beta <- c(-3, -1, 2)
X <- (model.matrix(Y ~ x1 + x2))  # design matrix
Y = as.numeric(Y + X %*% beta)

Now let’s fit the model using ngme. Here we can use control_opt to modify the control variables for the ngme. See ?control_opt for more optioins.

# # Fit the model with the AR1 model
ngme_out <- ngme(
  Y ~ x1 + x2 + f(
    1:n_obs,
    name = "my_ar",
    model = "ar1",
    noise = noise_nig()
  ),
  data=data.frame(x1=x1, x2=x2, Y=Y),
  control_opt = control_opt(
    burnin = 100,
    iterations = 1000,
    std_lim = 0.4,
    n_parallel_chain = 4,
    stop_points = 10,
    print_check_info = FALSE,
    seed = 3,
    sampling_strategy = "ws"
    # verbose = T
  )
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.

Next we can read the result directly from the object.

ngme_out
#> *** Ngme object ***
#> 
#> Fixed effects: 
#> (Intercept)          x1          x2 
#>       -2.88       -1.38        2.03 
#> 
#> Models: 
#> $my_ar
#>   Model type: AR(1)
#>       rho = 0.564
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = 1.85
#>       sigma = 2.99
#>       nu = 0.905
#> 
#> Measurement noise: 
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 0.587
#> 
#> 
#> Number of replicates is  1

As we can see, the model converges in 350 iterations. The estimation results are close to the real parameter.

We can also use the traceplot function to see the estimation traceplot.

traceplot(ngme_out, "my_ar")
Parameters of the AR1 model

Parameters of the AR1 model 

#> Last estimates:
#> $rho
#> [1] 0.5637253
#> 
#> $mu
#> [1] 1.854818
#> 
#> $sigma
#> [1] 2.985931
#> 
#> $nu
#> [1] 0.9143952

We can also do a density comparison with the estimated noise and the true NIG noise:

# ngme_out$replicates[[1]] means for the 1st replicate
plot(
  ngme_out$replicates[[1]]$models[[1]]$noise,
  noise_nig(mu = mu, sigma = sigma, nu = nu)
)

Paraná dataset

The rainfall data from Paraná (Brazil) is collected by the National Water Agency in Brazil (Agencia Nacional de Águas, ANA, in Portuguese). ANA collects data from many locations over Brazil, and all these data are freely available from the ANA website (http://www3.ana.gov.br/portal/ANA).

We will briefly illustrate the command we use, and the result of the estimation.

library(INLA)
#> Loading required package: Matrix
#> This is INLA_25.03.24 built 2025-03-24 16:12:01 UTC.
#>  - See www.r-inla.org/contact-us for how to get help.
#>  - List available models/likelihoods/etc with inla.list.models()
#>  - Use inla.doc(<NAME>) to access documentation
#> 
#> Attaching package: 'INLA'
#> The following object is masked from 'package:ngme2':
#> 
#>     f
data(PRprec)
data(PRborder)

# Create mesh
coords <- as.matrix(PRprec[, 1:2])
prdomain <- fmesher::fm_nonconvex_hull(coords, -0.03, -0.05, resolution = c(100, 100))
prmesh <- fmesher::fm_mesh_2d(boundary = prdomain, max.edge = c(0.45, 1), cutoff = 0.2)

# monthly mean at each location
Y <- rowMeans(PRprec[, 12 + 1:31]) # 2 + Octobor

ind <- !is.na(Y) # non-NA index
Y <- Y_mean <- Y[ind]
coords <- as.matrix(PRprec[ind, 1:2])
seaDist <- apply(spDists(coords, PRborder[1034:1078, ],
  longlat = TRUE
), 1, min)

Plot the data:

Mean of the rainfall in Octobor 2012 in Paraná

Mean of the rainfall in Octobor 2012 in Paraná 

# Define the control options
control = control_opt(
  iterations = 3000,
  n_slope_check = 4,
  stop_points = 10,
  std_lim = 0.1,
  n_parallel_chain = 4,
  print_check_info = FALSE,
  seed = 16
)

m_gauss_nig <- ngme(
  formula = Y ~ 1 +
    f(seaDist, name="rw1", model = "rw1", noise = noise_normal()) +
    f(coords, model = "matern", mesh = prmesh, name="spde", noise = noise_normal()),
  data = data.frame(Y = Y),
  family = noise_nig(),
  control_opt = control
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.
m_gauss_nig
#> *** Ngme object ***
#> 
#> Fixed effects: 
#> (Intercept) 
#>        11.9 
#> 
#> Models: 
#> $rw1
#>   Model type: Random walk (order 1)
#>       No parameter.
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 0.164
#> 
#> $spde
#>   Model type: Matern
#>       kappa = 9.24
#>   Noise type: NORMAL
#>   Noise parameters: 
#>       sigma = 68
#> 
#> Measurement noise: 
#>   Noise type: NIG
#>   Noise parameters: 
#>       mu = 0.54
#>       sigma = 1.97
#>       nu = 1.33
#> 
#> 
#> Number of replicates is  1

# traceplots
## fixed effects and measurement error
traceplot(m_gauss_nig)

#> Last estimates:
#> $mu
#> [1] 0.5375344
#> 
#> $sigma
#> [1] 1.963767
#> 
#> $nu
#> [1] 1.352117
#> 
#> $`fixed effect 1`
#> [1] 11.93805

## spde model
traceplot(m_gauss_nig, "spde")

#> Last estimates:
#> $kappa
#> [1] 9.199731
#> 
#> $sigma
#> [1] 68.39649

Parameter estimation results:

Estimations for the model
intercept noise_mu noise_sigma noise_nu rw_sigma ma_kappa ma_sigma
(Intercept) 11.9 0.54 1.97 1.33 -1.81 9.24 68

Similarily, we can fit some different models:

m_gauss_gauss <- ngme(
  formula = Y ~ 1 +
    f(seaDist, name="rw1", model = "rw1", noise = noise_normal()) +
    f(coords, model = "matern", mesh = prmesh, name="spde", noise = noise_normal()),
  data = data.frame(Y = Y),
  family = noise_normal(),
  control_opt = control
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.

m_nig_gauss <- ngme(
  formula = Y ~ 1 +
    f(seaDist, name="rw1", model = "rw1", noise = noise_nig()) +
    f(coords, model = "matern", mesh = prmesh, name="spde", noise = noise_normal()),
  data = data.frame(Y = Y),
  family = noise_normal(),
  control_opt = control
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.

m_nig_nig <- ngme(
  formula = Y ~ 1 +
    f(seaDist, name="rw1", model = "rw1", noise = noise_nig()) +
    f(coords, model = "matern", mesh = prmesh, name="spde", noise = noise_nig()),
  data = data.frame(Y = Y),
  family = noise_nig(),
  control_opt = control
)
#> Starting estimation... 
#> 
#> Starting posterior sampling... 
#> Posterior sampling done! 
#> Note:
#>       1. Use ngme_post_samples(..) to access the posterior samples.
#>       2. Use ngme_result(..) to access different latent models.

Prediction 

nxy <- c(150, 100)
projgrid <- rSPDE::rspde.mesh.projector(prmesh,
  xlim = range(PRborder[, 1]),
  ylim = range(PRborder[, 2]), dims = nxy
)

xy.in <- inout(projgrid$lattice$loc, cbind(PRborder[, 1], PRborder[, 2]))

coord.prd <- projgrid$lattice$loc[xy.in, ]
plot(coord.prd, type = "p", cex = 0.1)
lines(PRborder)
points(coords[, 1], coords[, 2], pch = 19, cex = 0.5, col = "red")


seaDist.prd <- apply(spDists(coord.prd,
  PRborder[1034:1078, ],
  longlat = TRUE
), 1, min)

# doing prediction by giving the predict location
pds <- predict(m_gauss_nig, map=list(rw1=seaDist.prd, spde=coord.prd))
lp <- pds$mean
ggplot() +
  geom_point(aes(
    x = coord.prd[, 1], y = coord.prd[, 2],
    colour = lp
  ), size = 2, alpha = 1) +
  geom_point(aes(
    x = coords[, 1], y = coords[, 2],
    colour = Y_mean
  ), size = 2, alpha = 1) +
  scale_color_gradientn(colours = viridis(100)) +
  geom_path(aes(x = PRborder[, 1], y = PRborder[, 2])) +
  geom_path(aes(x = PRborder[1034:1078, 1], y = PRborder[
    1034:1078,
    2
  ]), colour = "red")

Cross-validation

We can further validate our model by using cross-validation method.

cv <- cross_validation(
  list(
    gauss_gauss = m_gauss_gauss,
    gauss_nig = m_gauss_nig,
    nig_gauss = m_nig_gauss,
    nig_nig = m_nig_nig
  ), 
  type = "k-fold", 
  k = 10,
  n_gibbs_samples = 3000,
  n_burnin = 200,
  seed = 20
)

# Create a basic table with knitr
cv_table <- knitr::kable(cv$mean.scores, caption = "Cross-validation results")

cv_table
Cross-validation results
MAE MSE neg.CRPS neg.sCRPS
gauss_gauss 1.741812 5.142371 1.249572 1.457003
gauss_nig 1.743243 5.194404 1.249606 1.456385
nig_gauss 1.742106 5.139523 1.249065 1.456654
nig_nig 1.766678 5.461027 1.263476 1.454758

As we can see, the NIG-Gaussian model performs the best in most of the metrics.