Problem fitting FEM models

[alice981792]

Hi,

I am working on a project where we are trying to model shops values on a network of streets.

We have isolated the location element of value and are trying to use (very noisy) data to suggest a smooth pattern of changing value with location. We have tried fitting Whittle-Matern models with alpha = 1,2, and the result was not as generalised/smooth as we would have liked. We are therefore trying to explore using higher alpha with the FEM approach in the package.

For reproducibility reasons I want to be able to fix values of range, sigma and sigma_e and fit models under these constraints. I have successfully implemented fixing range, sigma, and sigma_e successfully for alpha = 1 through the model_options parameter in graph_lme. However, when I try to use the same model_options code to fix parameters for higher alpha with the FEM approach, it appears that the model fitting process ignores them.

I’ve included an example simulation below to demonstrate what I mean and have tried to highlight some of the things I have come across that aren’t making sense to me that I need help with.

I’d also be interested in understanding the timelines for whether higher alpha will be implemented in INLA, as we’d be keen to take advantage of things like the formalised priors and inbuilt uncertainty measures.

Cheers in advance for the help!

Kind regards,
Alice

# simulation taken mainly from by https://davidbolin.github.io/MetricGraph/articles/fem_models.html but with different parameters

# Load packages
library(rSPDE)  # rSPDE 2.3.3
library(MetricGraph)  # MetricGraph 1.4.0
library(tidyverse) # tidyverse 2.0.0

# Create example metric graph
edge1 <- rbind(c(0,0),c(1,0))
edge2 <- rbind(c(0,0),c(0,1))
edge3 <- rbind(c(0,1),c(-1,1))
theta <- seq(from=pi,to=3*pi/2,length.out = 20)
edge4 <- cbind(sin(theta),1+ cos(theta))
edges = list(edge1, edge2, edge3, edge4)
graph <- metric_graph$new(edges = edges)
graph$plot()

# Construct mesh
graph$build_mesh(h = 0.01)
graph$plot(mesh=TRUE)

# Specify model for simulation
sigma <- 1.3
range <- 0.5
nu <- 6  # note the high alpha here 
rspde.order <- 2

set.seed(22)
op <- matern.operators(nu = nu, range = range, sigma = sigma, parameterization = "matern", m = rspde.order, graph = graph) 
u <- simulate(op)

graph$plot_function(X = u)

set.seed(2)
# simulate 20 observations of this SPDE on each edge 
obs.per.edge <- 20
obs.loc <- NULL
for(i in 1:graph$nE) {
  obs.loc <- rbind(obs.loc, cbind(rep(i,obs.per.edge), runif(obs.per.edge)))
}
n.obs <- obs.per.edge*graph$nE
A <- graph$fem_basis(obs.loc)

# Add noise and intercept to the observations
sigma.e <- 1.2
intercept <- 10
Y <- intercept + as.vector(A %*% u + sigma.e * rnorm(n.obs))
Y_true <- intercept + as.vector(A %*% u)

# Prepare data for fitting a model
df_data <- data.frame(y = Y, y_true = Y_true, edge_number = obs.loc[,1], distance_on_edge = obs.loc[,2]) %>%
  # sometimes in streets you have gaps in parades of shops, so artificially adding these in to try to match our conditions
  filter(
    !(edge_number == 1 & between(distance_on_edge, 0.5, 0.7))
    & !(edge_number == 3 & between(distance_on_edge, 0.7, 0.9))
    #& !(edge_number == 4 & between(distance_on_edge, 0.5, 0.6))
  )
graph$add_observations(data = df_data, normalized = TRUE, clear_obs = TRUE)

graph$plot(data = "y")

# plot each edge on a scatter plot to see the noise distribution around the true value of y
p <- df_data %>%
  ggplot() +
    # black line is the true value of y
    geom_line(
      data = cbind(graph$get_mesh_locations(), intercept + as.vector(u)) %>% 
        as.data.frame() %>% 
        magrittr::set_colnames(c("edge_number", "distance_on_edge", "y_true")),
      aes(distance_on_edge, y_true)
    ) +
    # black dots are the noisy observations
    geom_point(aes(distance_on_edge, y)) +
    facet_wrap(~edge_number)

p

# to make predictions we will use the mesh locations with the correct colnames, so will wrangle this data upfront
mesh_locations <- graph$get_mesh_locations() %>%
  as.data.frame() %>%
  magrittr::set_colnames(c("edge_number", "distance_on_edge"))

# fit model 1: non FEM with alpha = 1 

fit1 <- graph_lme(
  y ~ 1,
  graph = graph,
  model = list(type = "WhittleMatern", alpha = 1, fem = FALSE)
)

# make predictions on the model mesh
# note that predict(fit1, mesh = TRUE) fails with error "In the presence of covariates, you should provide the data"

pred1 <- predict(fit1, newdata = mesh_locations)$mean

# plot model 
p + 
  geom_line(
    data = mesh_locations %>%
      mutate(pred = pred1),
    aes(distance_on_edge, pred),
    colour = "blue"
  )

# Compare the values of the parameters of fit with the true ones:
print(data.frame(
  sigma = c(sigma, fit1$matern_coeff$random_effects[1]), 
  range = c(range, fit1$matern_coeff$random_effects[2]),
  sigma_e = c(sigma.e, fit1$coeff$measurement_error[1]),
  row.names = c("Truth", "Estimates")
))

# fit model 2: non FEM with alpha = 1 and fixed param to the true values for the model (except alpha)

fit2 <- graph_lme(
  y ~ 1,
  graph = graph,
  model = list(type = "WhittleMatern", alpha = 1, fem = FALSE),
  model_options = list(
    fix_sigma = sigma,
    fix_range = range,
    fix_sigma_e = sigma.e
  )
)

# make predictions on the model mesh
pred2 <- predict(fit2, newdata = mesh_locations)$mean

# plot model 
p + 
  geom_line(
    data = mesh_locations %>%
      mutate(pred = pred2),
    aes(distance_on_edge, pred),
    colour = "blue"
  )

# note that for the following models, something is going wrong when predict is applied to a FEM model
# therefore I edited the predict function a bit to try to understand what the problem was 
# see notes later during model fit 3 to explain what I was trying to sort out

predict_edit <- function(object, newdata = NULL, loc = NULL, mesh = FALSE, which_repl = NULL, compute_variances = FALSE, posterior_samples = FALSE,
                               n_samples = 100, sample_latent = FALSE, edge_number = "edge_number",
                               distance_on_edge = "distance_on_edge", normalized = FALSE, return_as_list = FALSE, return_original_order = TRUE,
                               alpha_int, # EDIT 1: I added in this argument
                               ...
                               ) {

    
  if (lifecycle::is_present(data)) {
    if (is.null(newdata)) {
      lifecycle::deprecate_warn("2.3.2.9000", "predict(data)", "predict(newdata)",
        details = c("`data` was provided but not `newdata`. Setting `newdata <- data`.")
      )
      newdata <- data
    } else {
      lifecycle::deprecate_warn("2.3.2.9000", "predict(data)", "predict(newdata)",
        details = c("Both `newdata` and `data` were provided. Only `newdata` will be considered.")
      )
    }
    data <- NULL
  }

  data <- newdata

  if(is.null(data)){
    if(!mesh){
      stop("If 'mesh' is false, you should supply data!")
    }
  }

  out <- list()

  coeff_fixed <- object$coeff$fixed_effects %>% as.matrix() 
  coeff_random <- object$coeff$random_effects
  coeff_meas <- object$coeff$measurement_error

  if(inherits(object, "graph_lme")){
    if(object$estimate_nu){
      coeff_random[1] <- coeff_random[1] - 0.5
    }
  }

  if(object$has_graph){
    loc <- cbind(data[[edge_number]], data[[distance_on_edge]])
  } else if(is.character(loc)){
    loc <- data[loc]
  }

  loc <- as.matrix(loc)

  X_cov_initial <- as.matrix(object$model_matrix[,-1])      # EDIT 2: I edited the this bit as it's messing up when you have an intercept only model 
  if(ncol(X_cov_initial) > 0){
    if(mesh){
      stop("In the presence of covariates, you should provide the data, including the covariates at the prediction locations. If you only want predictions for the latent model, set 'only_latent' to TRUE.")
    }
  }

  X_cov_pred <- stats::model.matrix(object$covariates, data)
  if(nrow(X_cov_pred) != nrow(as.matrix(loc))){
    stop("Covariates not found in data.")
  }

  if(sum(duplicated(loc)) > 0){
    warning("There are duplicated locations for prediction, we will try to process the data to extract the unique locations,
    along with the corresponding covariates.")
    if(nrow(X_cov_pred) == nrow(loc)){
      data_tmp <- cbind(loc, X_cov_pred)
    }
    data_tmp <- unique(data_tmp) 
    if(sum(duplicated(cbind(data_tmp[,1:ncol(loc)]))) > 0){
      stop("Data processing failed, please provide a data with unique locations.")
    }
  }

 
  if(!mesh){
    n_prd <- length(loc[,1])
    # Convert data to normalized
    if(object$has_graph && !normalized){
      loc[,2] <- loc[,2] / object$graph$edge_lengths[loc[,1]]
    }
    Aprd <- object$latent_model$make_A(loc)
  } else{
    Aprd <- Matrix::Diagonal(nrow(object$latent_model$C))
  }

  ## 
  repl_vec <- object$repl
  if(is.null(which_repl)){
    u_repl <- unique(repl_vec)
  } else{
    u_repl <- unique(which_repl)
  }

  ##

  if(all(dim(X_cov_pred) == c(0,1))){
    X_cov_pred <- matrix(1, nrow = nrow(loc), ncol=1)
  }
  if(ncol(X_cov_pred) > 0){
    mu_prd <- X_cov_pred %*% coeff_fixed
  } else{
    mu_prd <- matrix(0, nrow = nrow(loc), ncol=1)
  }

  model_matrix_fit <- object$model_matrix

  model_matrix_fit <- as.matrix(model_matrix_fit)

  if(ncol(model_matrix_fit)>1){
    X_cov_fit <- model_matrix_fit[,2:ncol(model_matrix_fit)] %>% as.matrix()
    mu <- X_cov_fit %*% coeff_fixed
  } else{
    mu <- 0
  }

  Y <- model_matrix_fit[,1] - mu

  model_type <- object$latent_model

  sigma.e <- coeff_meas[[1]]
  sigma_e <- sigma.e

  ## construct Q

  if(object$estimate_nu){
      alpha_est <- coeff_random[1]
      tau_est <- coeff_random[2]
      kappa_est <- coeff_random[3]
  } else{
      tau_est <- coeff_random[1]
      kappa_est <- coeff_random[2]
      alpha_est <- NULL
  }

  new_rspde_obj <- update(object$latent_model,
                        user_alpha = alpha_est,
                        user_kappa = kappa_est,
                        user_tau = tau_est,
                        parameterization = "spde")

                   
  Q <- new_rspde_obj$Q                   

  idx_obs_full <- as.vector(!is.na(Y))

  Aprd <- kronecker(matrix(1, 1, object$rspde_order + 1), Aprd)

  for(repl_y in u_repl){
    idx_repl <- repl_vec == repl_y

    idx_obs <- idx_obs_full[idx_repl]

    y_repl <- Y[idx_repl]
    y_repl <- y_repl[idx_obs] 

    A_repl <- object$A_list[[repl_y]]

    if (!alpha_int) {  # EDIT 3: I added this if statement in (as I think it's the bit going wrong?)
      A_repl <- A_repl
    } else {
      A_repl <- kronecker(matrix(1, 1, object$rspde_order + 1), A_repl)     
    }

    Q_xgiveny <- t(A_repl) %*% A_repl/sigma_e^2 + Q

    mu_krig <- solve(Q_xgiveny,as.vector(t(A_repl) %*% y_repl / sigma_e^2)) 

    # mu_krig <- mu_krig[(gap+1):length(mu_krig)]

    mu_krig <- Aprd %*% mu_krig 

    mu_re <- mu_krig

    mu_fe <- mu_prd

    mu_krig <- mu_prd + mu_krig

    mean_tmp <- as.vector(mu_krig)

    mean_fe_tmp <- as.vector(mu_fe)

    mean_re_tmp <- as.vector(mu_re)
        
    if(!return_as_list){
      out$mean <- c(out$mean, mean_tmp)
      out$repl <- c(out$repl, rep(repl_y,n_prd))
      out$fe_mean <- c(out$fe_mean, mean_fe_tmp)
      out$re_mean <- c(out$re_mean, mean_re_tmp)
    } else{
      out$mean[[repl_y]] <- mean_tmp
      out$fe_mean[[repl_y]] <- mean_fe_tmp
      out$re_mean[[repl_y]] <- mean_re_tmp
    }

    if (compute_variances) {
        post_cov <- Aprd%*%solve(Q_xgiveny, t(Aprd))
        var_tmp <- pmax(diag(post_cov),0)

      if(!return_as_list){
        out$variance <- rep(var_tmp, length(u_repl))
      } else {
          for(repl_y in u_repl){
            out$variance[[repl_y]] <- var_tmp
          }
      }
    }

    if(posterior_samples){
      mean_tmp <- as.vector(mu_krig)
      post_cov <- Aprd%*%solve(Q_xgiveny, t(Aprd))
      Z <- rnorm(dim(post_cov)[1] * n_samples)
      dim(Z) <- c(dim(post_cov)[1], n_samples)
      LQ <- tryCatch(chol(post_cov), error =  function(e){chol(post_cov + 1e-8*Matrix::Diagonal(nrow(post_cov)))})
      X <- LQ %*% Z
      X <- X + mean_tmp
      if(!sample_latent){
        X <- X + matrix(rnorm(n_samples * length(mean_tmp), sd = sigma.e), nrow = length(mean_tmp))
      } else{
        X <- X - as.vector(mu_prd)
      }

      if(!return_as_list){
        out$samples <- rbind(out$samples, X)
      } else{
        out$samples[[repl_y]] <- X
      }
    }
  }

  return(out)
}

# fit model 3: FEM with alpha = 1

fit3 <- graph_lme(
  y ~ 1,
  graph = graph,
  model = list(type = "WhittleMatern", alpha = 1, fem = TRUE)
)

# make predictions on the model mesh
# something is going wrong when predict is applied to a FEM model
# predict(fit3, newdata = mesh_locations)$mean errors with Error in if (ncol(X_cov_initial) > 0) { : argument is of length zero
# I identified that this was caused by having an intecept only model, so the function was expecting a matrix where it recieved a vector (EDIT 2)
# EDIT 1 and EDIT 3 are to deal with things erroring when alpha is an integer. Something to do with matrices being of incorrect dimensions. I made a guess at how the function could be corrected to see if it would output anything sensible...
pred3 = predict_edit(fit3, newdata = mesh_locations, alpha_int = TRUE)$mean

# plot model 
p + 
  geom_line(
    data = mesh_locations %>%
      mutate(pred = pred3),
    aes(distance_on_edge, pred),
    colour = "blue"
  )

# I'm a bit confused as to why this looks so smooth in comparison to the other alpha = 1 models I've seen fitted the non FEM way...
# It looks like this might be using a higher alpha??

# Compare the values of the parameters of fit with the true ones:
print(data.frame(
  sigma = c(sigma, fit3$matern_coeff$random_effects[1]), 
  range = c(range, fit3$matern_coeff$random_effects[2]),
  sigma_e = c(sigma.e, fit3$coeff$measurement_error[1]),
  row.names = c("Truth", "Estimates")
))

# fit model 4: FEM with alpha = 1 and attempting to fix params to true values

fit4 <- graph_lme(
  y ~ 1,
  graph = graph,
  model = list(type = "WhittleMatern", alpha = 1, fem = TRUE),
  model_options = list(
    fix_sigma = sigma,
    fix_range = range,
    fix_sigma_e = sigma.e
  )
)

# make predictions on model mesh
pred4 = predict_edit(fit4, newdata = mesh_locations, alpha_int = TRUE)$mean

# plot model 
p + 
  geom_line(
    data = mesh_locations %>%
      mutate(pred = pred4),
    aes(distance_on_edge, pred),
    colour = "blue"
  )

# This comes out looking the same?

# Compare the values of the parameters of fit with the true ones:
print(data.frame(
  sigma = c(sigma, fit4$matern_coeff$random_effects[1]), 
  range = c(range, fit4$matern_coeff$random_effects[2]),
  sigma_e = c(sigma.e, fit4$coeff$measurement_error[1]),
  row.names = c("Truth", "Estimates")
))

# model_options argument has been ignored?

# Please help :)

[AlexandreBSimas]

Hi Alice,

You are correct, it is ignoring the model_options when using fem. Basically, internally it calls rspde_lme, and these options are not yet available for rspde_lme. However, I believe this week we will update rSPDE package by adding such functionalities to rspde_lme, and then, I will let you know here, so you can use it.

Thank you for calling attention for this.

Best,
Alexandre

[alice981792]

Hi Alexandre,

Cheers for getting back - that sounds fantastic. I'll keep an eye on this thread for when it's ready to test out.

Kind regards,
Alice

[alice981792]

Hi Alexandre,

Any update on timings for this? Super excited to be able to implement this in my modelling.

Kind regards,
Alice

[AlexandreBSimas]

Hi Alice,

I am sorry, I did not receive any notification for this message and only saw now.

But, in any case, I just finished all the corrections now. I messaged you as soon as the updates were available.

Hope you enjoy :)

And if you run into any problem, please let me know.

Best,
Alexandre

[AlexandreBSimas]

Hi Alice,

Everything has been fixed now! You can now use it, and everything should work.

But you will need to install the latest development version of rSPDE:

devtools::install_github("davidbolin/rSPDE", ref = "devel")

and the latest development version of MetricGraph:

 devtools::install_github("davidbolin/MetricGraph", ref = "devel")

Then, your code should directly work without the need for your predict_edit, but by directly using the predict. The problem was that it was an integer case, in which there is no rational approximation, so the dimensions did not match as it was assuming a rational approximation.

Please, let me know if everything is ok now.

Best,
Alexandre.

[AlexandreBSimas]

Ahhh, another thing. I have updated the vignettes for rSPDE:
https://davidbolin.github.io/rSPDE/articles/rspde_cov.html#further-details-on-rspde_lme

and for fem models in metric graphs:
https://davidbolin.github.io/MetricGraph/articles/fem_models.html (this one is the online version is not updated yet, but if not you will probably just need to wait a little bit, as it is currently building).

So you can check a bit of the details how to choose the options, but since you were already using them correctly, maybe it is not super relevant to you, but just wanted to let you know about this.

[alice981792]

Hi Alexandre,

Cheers for getting in touch and thanks so much for the updated functionality!

I've been having a play around today and all looks to be working great - can't wait to spend more time seeing what is possible with different types of alpha for my project.

Kind regards,
Alice

[AlexandreBSimas]

Hi Alice,

Great to hear that everything is working smoothly! I hope you have a good time with the package while doing your project.

Feel free to reach out if you run into anything else.

Best,
Alexandre