ngme2.Rmd

---
title: "Ngme2 - A new Flexible R Package for Latent non-Gaussian Models"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ngme2 - A new Flexible R Package for Latent non-Gaussian Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
set.seed(10)
```

## 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
$$

- $m$ is the number of subjects, $n_i$ is the number of observations for each subject,
- $Y$ 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,
- $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:
$$\gamma + \mu V + \sigma \sqrt{V}Z,$$
where $\gamma, \mu, \sigma$ are parameters, $Z\sim N(0,1)$ and is independent of $V$, and $V$ is a positive infinitely divisible random variable.
It results in the following form, where $K$ is the operator part:

$$
KW|V \sim N(\gamma + \mu V, \sigma^2 \, diag(V)),
$$
also, $\mu$ and $\sigma$ can be non-stationary.

One example in `Ngme2` is the normal inverse Gaussian (NIG) noise, in this case, $V$ follows Inverse Gaussian distribution with parameter $\nu$ (IG($\nu$, $\nu$)).
See `?nig` for more details.

### 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](pred_and_est.html) 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:

```{r load-packages, message=FALSE, warning=FALSE}
library(fmesher)
library(splancs)
library(lattice)
library(ggplot2)
library(grid)
library(gridExtra)
library(viridis)
library(ngme2)
# load_all()
```

```{r noise}
noise_normal(sigma = 1)              # normal noise
noise_nig(mu = 1, sigma = 2, nu = 1) # nig noise
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)
```

<!-- to-do, plot some nig noise -->

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:

\[
  - \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`:

```{r f_model}
ngme2::f(1:10, model = "ar1", noise = noise_nig())
```

One useful model would be the SPDE model with Gaussian or non-Gaussian noise, see the vignette
\code{vignette("SPDE-approach", package = "ngme2")} 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

```{r}
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 = \beta_0 + \beta_1 x_1 + W_1(x_2) + W_2(x_3) + \epsilon,
$$
where $W_1$ is an AR(1) process, $W_2$ is a random walk 1 process.
$x_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:

```{r ngme_block}
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
  )
)
```

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

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

```{r simulate_ar1}
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.

```{r fit_ar1, cache=TRUE}
# # 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
  )
)
```

Next we can read the result directly from the object.

```{r out}
ngme_out
```

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.
```{r traceplot, fig.align="center", fig.width=6,fig.height=6,  fig.cap="Parameters of the AR1 model"}
traceplot(ngme_out, "my_ar")
```

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

```{r compare-density}
# 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.

```{r read-parana}
library(INLA)
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:

```{r plot-data, message=FALSE,fig.width=6,fig.height=6,fig.align='center', echo=FALSE, fig.cap="Mean of the rainfall in Octobor 2012 in Paraná"}
ggplot() +
  geom_point(aes(
    x = coords[, 1], y = coords[, 2],
    colour = Y
  ), 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") +
  theme(
    axis.title.x = element_blank(),   # Remove x-axis label
    axis.title.y = element_blank(),   # Remove y-axis label
    axis.line = element_blank(),      # Remove axis lines
    axis.text.x = element_blank(),    # Remove x-axis tick labels
    axis.text.y = element_blank(),    # Remove y-axis tick labels
    legend.title = element_blank()
  )
```

```{r fit-parana, cache=TRUE}
out <- 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 = "nig",
  control_opt = control_opt(
    estimation = T,
    iterations = 1000,
    n_slope_check = 4,
    stop_points = 10,
    std_lim = 0.1,
    n_parallel_chain = 4,
    print_check_info = FALSE,
    seed = 16
  )
)
out

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

## spde model
traceplot(out, "spde")
```

Parameter estimation results:

```{r est-result, echo=FALSE}
res1 <- data.frame(
  intercept    = format(out$replicates[[1]]$feff, digits=3),
  noise_mu     = format(out$replicates[[1]]$noise$theta_mu, digits=3),
  noise_sigma  = format(exp(out$replicates[[1]]$noise$theta_sigma), digits=3),
  noise_nu     = format(exp(out$replicates[[1]]$noise$theta_nu), digits=3),
  rw_sigma     = format(out$replicates[[1]]$models[[1]]$noise$theta_sigma, digits=3),
  ma_kappa     = format(exp(out$replicates[[1]]$models[[2]]$theta_K), digits=3),
  ma_sigma     = format(exp(out$replicates[[1]]$models[[2]]$noise$theta_sigma), digits=3)
)
knitr::kable(res1, caption = "Estimations for the model")
```

### Prediction

```{r predict, cache=TRUE}
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(out, 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.

```{r cv, cache=TRUE}
cross_validation(out, type="k-fold", k=10, print=TRUE)
```