pkern_GLS: Generalized least squares (GLS) estimator

In deankoch/pkern: Fast Kriging with Product Kernels

Generalized least squares (GLS) estimator

Description

Computes coefficients b of the linear predictor E(Z) = Xb using the GLS equation

Usage

pkern_GLS(g_obs, pars, X = NA, fac = NULL, fac_method = "eigen", out = "b")

Arguments

g_obs	list of form returned by pkern_grid (with entries 'gdim', 'gres', 'gval')
pars	list of form returned by pkern_pars (with entries 'y', 'x', 'eps', 'psill')
X	matrix or NA, the linear predictors (in columns) excluding intercept
fac	matrix or list, (optional) pre-computed covariance matrix factorization
fac_method	character, factorization method: 'eigen' (default) or 'chol' (see pkern_var)
out	character, either 'b' (coefficients) or 'z' (linear predictor)

Details

The GLS equation is: b = ( X^T V^-1 X )^-1 X^T V^-1 z,

where V is the covariance matrix for data z, and X is a matrix of covariates. V is generated from the covariance model specified in pars on the grid g_obs, where the non-NA values of g_obs$gval form the vector z. Operations with V^-1 are computed using the factorization fac, or else as specified in fac_method.

When out='z', the function returns the product Xb instead of b.

Argument X should provide matrix X without the intercept column (a vector of 1's). DO NOT include an intercept column in argument X or you will get collinearity errors (the function appends it without checking if its there already). The columns of X must be independent, and its rows should match the entries of g_obs$gval, or its non-NA data values, in order. Use X=NA to specify an intercept-only model; ie to fit a spatially constant mean. This replaces X in the GLS equation by a vector of 1's.

g_obs$gval can be a matrix whose columns are multiple repetitions (layers) of the same spatial process (see pkern_grid), in which case the covariates in X are recycled for each layer. Layers are assumed mutually independent and the GLS equation is evaluated using the corresponding block-diagonal V. Note that this is equivalent to (but faster than) calling pkern_GLS separately on each layer with the same X and averaging the estimated b's.

Value

numeric vector, either b (default) or Xb (if out='z')

Examples

# set up example grid, and covariance parameters
gdim = c(45, 31)
n = prod(gdim)
g = pkern_grid(gdim)
pars = modifyList(pkern_pars(g, 'gau'), list(psill=2))

# generate spatial noise
z = pkern_sim(g, pars)
pkern_plot(modifyList(g, list(gval=z)))

# generate some covariates and data
n_betas = 3
betas = rnorm(n_betas, 0, 10)
X_all = cbind(1, matrix(rnorm(n*(n_betas-1)), n))
lm_actual = as.vector(X_all %*% betas)
g_obs = modifyList(g, list(gval=z+lm_actual))

# exclude intercept column in calls to pkern_GLS
X_pass = X_all[,-1]

# find the GLS coefficients
betas_est = pkern_GLS(g_obs, pars, X_pass)
print(betas_est)
print(betas)

# compute trend as product of betas with matrix X_all, or by setting out='z'
lm_est = X_all %*% betas_est
max( abs( pkern_GLS(g_obs, pars, X_pass, out='z') - lm_est ) )

# repeat with pre-computed eigen factorization
fac_eigen = pkern_var(g_obs, pars, fac_method='eigen', sep=TRUE)
betas_est_compare_eigen = pkern_GLS(g_obs, pars, X_pass, fac=fac_eigen)
max( abs( betas_est_compare_eigen - betas_est ) )

# missing data example
n_obs = 10
idx_rem = sort(sample.int(n, n-n_obs))
g_miss = g_obs
g_miss$gval[idx_rem] = NA
pkern_plot(g_miss)
betas_est = pkern_GLS(g_miss, pars, X_pass)
print(betas_est)
print(betas)

# set X to NA to estimate the a spatially constant trend (the adjusted mean)
pkern_GLS(g_miss, pars, X=NA)
mean(g_miss$gval, na.rm=TRUE)

# generate some extra layers
z_extra = lapply(seq(9), function(x) pkern_sim(g, pars))
z_multi = lm_actual + do.call(cbind, c(list(z), z_extra))

# multi-layer example with missing data
is_obs = !is.na(g_miss$gval)
map_sparse = match(seq(n), which(is_obs))
g_sparse = modifyList(g_miss, list(gval=z_multi[is_obs,], idx_grid=map_sparse))
betas_sparse = pkern_GLS(g_obs=g_sparse, pars, X=X_pass)
print(betas_sparse)
print(betas)

pkern_GLS(g_sparse, pars, NA)
mean(g_sparse$gval, na.rm=TRUE)

deankoch/pkern documentation built on Oct. 26, 2023, 8:54 p.m.