graph_lme {MetricGraph}

Metric graph linear mixed effects models

Description

Fitting linear mixed effects model in metric graphs. The random effects can be Gaussian Whittle-Matern fields, discrete Gaussian Markov random fields based on the graph Laplacian, as well as Gaussian random fields with isotropic covariance functions.

Usage

graph_lme(
  formula,
  graph,
  model = list(type = "linearModel"),
  which_repl = NULL,
  optim_method = "L-BFGS-B",
  possible_methods = c("Nelder-Mead", "L-BFGS-B"),
  model_options = list(),
  BC = 0,
  previous_fit = NULL,
  fix_coeff = FALSE,
  parallel = FALSE,
  n_cores = parallel::detectCores() - 1,
  optim_controls = list(),
  improve_hessian = FALSE,
  hessian_args = list(),
  check_euclidean = TRUE
)

Arguments

formula

Formula object describing the relation between the response variables and the fixed effects.

graph

A metric_graph object.

model

The random effects model that will be used (it also includes the option of not having any random effects). It can be either a character, whose options are 'lm', for linear models without random effects; 'WM1' and 'WM2' for Whittle-Matern models with \alpha=1 and 2, with exact precision matrices, respectively; 'WM' for Whittle-Matern models where one also estimates the smoothness parameter via finite-element method; 'isoExp' for a model with isotropic exponential covariance; 'GL1' and 'GL2' for a SPDE model based on graph Laplacian with \alpha = 1 and 2, respectively. 'WMD1' is the directed Whittle-Matern with \alpha=1. There is also the option to provide it as a list containing the elements type, which can be linearModel, WhittleMatern, graphLaplacian or isoCov. linearModel corresponds to a linear model without random effects. For WhittleMatern models, that is, if the list contains type = 'WhittleMatern', one can choose between a finite element approximation of the precision matrix by adding fem = TRUE to the list, or to use the exact precision matrix (by setting fem = FALSE). If fem is FALSE, there is also the parameter alpha, to determine the order of the SPDE, which is either 1 or 2. If fem is TRUE and alpha is not specified, then the default value of alpha=1 will be used. If fem is TRUE and one does not specify alpha, it will be estimated from the data. However, if one wants to have alpha fixed to some value, the user can specify either alpha or nu in the list. See the vignettes for examples. Finally, for type 'WhittleMatern', there is an optional argument, rspde_order, that chooses the order of the rational approximation. By default rspde_order is 2. Finally, if one wants to fit a nonstationary model, then fem necessarily needs to be TRUE, and one needs to also supply the matrices B.tau and B.kappa or B.range and B.sigma. For graph-Laplacian models, the list must also contain a parameter alpha (which is 1 by default). For isoCov models, the list must contain a parameter cov_function, containing the covariance function. The function accepts a string input for the following covariance functions: 'exp_covariance', 'WM1', 'WM2', 'GL1', 'GL2'. For another covariance function, the function itself must be provided as the cov_function argument. The default is 'exp_covariance', the exponential covariance. We also have covariance-based versions of the Whittle-Matern and graph Laplacian models, however they are much slower, they are the following (string) values for 'cov_function': 'alpha1' and 'alpha2' for Whittle-Matern fields, and 'GL1' and 'GL2' for graph Laplacian models. Finally, for Whittle-Matern models, there is an additional parameter version, which can be either 1 or 2, to tell which version of the likelihood should be used. Version is 1 by default.

which_repl

Vector or list containing which replicates to consider in the model. If NULL all replicates will be considered.

optim_method

The method to be used with optim function.

possible_methods

Which methods to try in case the optimization fails or the hessian is not positive definite. The options are 'Nelder-Mead', 'L-BFGS-B', 'BFGS', 'CG' and 'SANN'. By default only 'Nelder-Mead' and 'L-BFGS-B' are considered.

model_options

A list containing additional options to be used in the model. Currently, it is possible to fix parameters during the estimation or change the starting values of the parameters. The general structure of the elements of the list is fix_parname and start_parname, where parname stands for the name of the parameter. If fix_parname is not NULL, then the model with be fitted with the parname being fixed at the value that was passed. If start_parname is not NULL, the model will be fitted using the value passed as starting value for parname. the For 'WM' models, the possible elements of the list are: fix_sigma_e, start_sigma_e, fix_nu, start_nu, fix_sigma, start_sigma, fix_range, start_range. Alternatively, one can use fix_sigma_e, start_sigma_e, fix_nu, start_nu, fix_tau, start_tau, fix_kappa, start_kappa. For 'WM1', 'WM2', 'isoExp', 'GL1' and 'GL2' models, the possible elements of the list are fix_sigma_e, start_sigma_e, fix_sigma, start_sigma, fix_range, start_range. Alternatively, one can use fix_sigma_e, start_sigma_e, fix_tau, start_tau, fix_kappa, start_kappa. For 'isoCov' models, the possible values are fix_sigma_e, start_sigma_e, fix_par_vec, start_par_vec. Observe that contrary to the other models, for 'isoCov' models, both fix_par_vec and start_par_vec should be given as vectors of the size of the dimension of the vector for the input of the covariance function passed to the 'isoCov' model. Furthermore, for 'isoCov' models, fix_par_vec is a logical vector, indicating which parameters to be fixed, and the values will be kept fixed to the values given to start_par_vec, one can also use fix_sigma_e and start_sigma_e for controlling the std. deviation of the measurement error.

BC

For WhittleMatern models, decides which boundary condition to use (0,1). Here, 0 is Neumann boundary conditions and 1 specifies stationary boundary conditions.

previous_fit

An object of class graph_lme. Use the fitted coefficients as starting values.

fix_coeff

If using a previous fit, should all coefficients be fixed at the starting values?

parallel

logical. Indicating whether to use optimParallel() or not.

n_cores

Number of cores to be used if parallel is true.

optim_controls

Additional controls to be passed to optim() or optimParallel().

improve_hessian

Should a more precise estimate of the hessian be obtained? Turning on might increase the overall time.

hessian_args

List of controls to be used if improve_hessian is TRUE. The list can contain the arguments to be passed to the method.args argument in the hessian function. See the help of the hessian function in 'numDeriv' package for details. Observet that it only accepts the "Richardson" method for now, the method "complex" is not supported.

check_euclidean

Check if the graph used to compute the resistance distance has Euclidean edges? The graph used to compute the resistance distance has the observation locations as vertices.

Value

A list containing the fitted model.