svcFit: Fit a spatially varying coefficient model

In jmhewitt/telefit: Estimation and Prediction for Remote Effects Spatial Process Models

Description

Fit a spatially varying coefficient model

Usage

svcFit(
  y,
  X,
  z,
  coords,
  miles = T,
  priors,
  nSamples,
  thin = 1,
  rw.initsd = 0.1,
  inits = list(),
  C = 1,
  alpha = 0.44
)

Arguments

y	vector containing responses for each timepoint. vector is blocked by timepoint.
X	matrix containing local covariates for each timepoint. each row are the covariates for one location and timepoint. matrix is blocked by timepoint.
z	matrix containing remote covariates. each column has remote covariates for one timepoint.
coords	n x 2 matrix containing lon-lat coordinates for locations.
miles	T/F for whether to compute great circle distances in miles (T) or km (F)
priors	A list containing parameters for the prior distributions. The list needs to contain the following values T list(Psi=matrix, nu=double) specifying the Inverse wishart prior distribution for the spatially varying coefficient process. beta list(Linv=matrix) specifying the prior precision matrix for the fixed local covariates. sigmasq list(a=double, b=double) specifying the prior shape and scale parameters for the covariance scale and nugget parameters. rho list(L=double, U=double) specifying the lower and upper bounds for the spatial range parameter. cov list(nu=double) specifying the smoothness for the matern covariance.
nSamples	number of MCMC iterations to run
thin	MCMC thinning; defaults to no thinning (thin=1)
rw.initsd	Initial proposal standard deviation for RW samplers
inits	optional list containing starting parameters for MCMC sampler
C	scaling factor used in adapting random walk proposal variances.
alpha	target acceptance rate for random walk proposals.

Examples

library(fields)
library(mvtnorm)

set.seed(2018)

# set key parameters
dims = list(N=100, nt=3, k=2, p=2)
params = list(sigmasq=.2, rho=.3, eps=.5, nu=.5)

# generate parameters and data
coords = matrix( runif(2 * dims$N), ncol = 2 )
X = matrix( rnorm(dims$p * dims$N * dims$nt), ncol = dims$p )
beta = c(-1, .5)
z = matrix( rnorm(dims$k * dims$nt), ncol = dims$nt)
H = maternCov(rdist.earth(coords), scale = params$sigmasq, range = params$rho,
              smoothness = params$nu, nugget = params$sigmasq * params$eps)
Hinv = solve(H)
Tm = matrix(c(.5,.2, .2, .5), ncol=2)/2
theta = kronSamp(Hinv, Tm)


# generate response
xb = X %*% beta
zt = as.numeric(apply(z, 2, function(d) { 
  kronecker(diag(dims$N), t(d)) %*% theta }))
w = kronSamp(diag(dims$nt), H)
y =  xb + zt + w


# fit model
it = 100
priors = list(
  T = list(Psi = .1*diag(dims$k), nu = dims$k),
  beta = list(Linv = diag(dims$p) * 1e-2),
  sigmasq = list(a=2, b=1),
  rho = list(L=0, U=1),
  cov = list(nu=.5)
)

fit = svcFit(y=y, X=X, z=z, coords=coords, priors=priors, nSamples=it)


#
# predict at new timepoints
#

# generate parameters and data
nt0 = 3
Xn = matrix( rnorm(dims$p * dims$N * nt0), ncol = dims$p )
zn = matrix( rnorm(dims$k * nt0), ncol = nt0)

# generate response
xbn = Xn %*% beta
ztn = as.numeric(apply(zn, 2, function(d) { 
  kronecker(diag(dims$N), t(d)) %*% theta }))
wn = kronSamp(diag(nt0), H)
yn =  xbn + ztn + wn

# predict responses
pred = svcPredict(fit, Xn, zn, burn = 50)

jmhewitt/telefit documentation built on Feb. 9, 2020, 7:15 p.m.