| mcmc.wish.icar {rwc} | R Documentation |
Markov chain Monte Carlo sampler for Generalized Wishart distance matrix data arising from an ICAR spatial model.
Description
Constructs and runs an MCMC sampler to estimate resistance parameters of different landscape features.
Usage
mcmc.wish.icar(Dobs, TL, obs.idx, df=1,
beta.start = rep(0, length(TL)),
beta.prior.mean = rep(0, length(TL)),
beta.prior.cov = diag(10, length(TL)),
tau.start = 0.1, tau.prior.var = 1,
theta.tune = diag(10^-4,length(TL)+1),
n.mcmc=100, adapt.max=10000, adapt.int=100,
print.iter=FALSE, output.trace.plot=FALSE)
Arguments
Dobs |
A square symmetric matrix of observed pairwise distances. For example, a genetic distance matrix. |
TL |
A list of transition matrices for different covariate raster layers, created by get.TL |
obs.idx |
A vector of unique indices of observed nodes in the graph defined by the raster grid. |
df |
An integer specifying the degrees of freedom of Dobs. |
beta.start |
Vector of initial values for conductance parameters beta. Default is a vector of zeros. |
beta.prior.mean |
Vector of prior mean values for conductance parameters beta. Default is a vector of zeros. |
beta.prior.cov |
Matrix of the prior covariance matrix for conductance parameters beta. Default is a diagonal matrix with diagonal entries equal to 10. |
tau.start |
Numeric starting value for the nugget variance tau. Default is 0.1 |
tau.prior.var |
Variance of the half-normal prior for tau. Default is 1. |
theta.tune |
Covariance matrix for the random walk MH sampler for all parameters. Default is a diagonal matrix with variance 10^-4. |
n.mcmc |
Integer number of iterations of the MCMC sampler to run. |
adapt.max |
Integer number (or INF) specifying the last iteration at which the covariance matrix of the proposal distribution will be adapted. Default is 10^5. |
adapt.int |
Interval at which the covariance matrix of the proposal distribution is adapted. Default is every 100 iterations. |
print.iter |
Logical. If TRUE, then the current state of the system will be printed to the console every 100 iterations. |
output.trace.plot |
Logical. If TRUE, then the trace plots of the sampler will be saved to "traceOut.pdf" every 100 iterations. |
Details
Runs an MCMC sampler to draw samples from the posterior distribution of model parameters (beta,tau) from the following model for an observed squared distance matrix Dobs:
-Dobs ~ GenWish(2*Sigma,df)
Sigma = K(Psi)K'+tau*I
where Psi is the covariance matrix of the observed nodes of a graph described by the transition list TL. That is, the total graph has Laplacian (precision) matrix Q, with off-diagonal entries of Q given by
Q_ij = exp( b0 + b1 x_1ij + ... + bp x_pij ), where beta=(b1,b2,...,bp) is a vector of log-conductance values of the covariates x_kij. Each x_kij is equal to (x_ki+x_kj)/2.
The prior on beta is N(beta.prior.mean,beta.prior.cov), and the prior on tau is tau~Half_Normal(0,tau.prior.var).
Value
A list containing output from the MCMC sampler.
beta |
Posterior samples for conductance parameters beta. |
Author(s)
Ephraim M. Hanks
References
Hanks and Hooten 2013. Circuit theory and model-based inference for landscape connectivity. Journal of the American Statistical Association. 108(501), 22-33.
Examples
## Not run:
## The following code conducts a simulation example by
## first simulating a heterogeneous landscape, then
## simulating pairwise distance data on the landscape
## and finally making inference on model parameters.
library(rwc)
library(MASS)
## source("GenWishFunctions_20170901.r")
##
## specify 2-d raster
##
ras=raster(nrow=30,ncol=30)
extent(ras) <- c(0,30,0,30)
values(ras) <- 1
plot(ras,asp=1)
ras
int=ras
cov.ras=ras
## simulate "resistance" raster
TL.int=get.TL(int)
Q.int=get.Q(TL.int,1)
set.seed(1248)
## values(cov.ras) <- as.numeric(rnorm.Q(Q.int
values(cov.ras) <- as.numeric(rnorm.Q(Q.int,zero.constraint=TRUE))
plot(cov.ras)
stack=stack(int,cov.ras)
plot(stack)
TL=get.TL(stack)
## Create "barrier" raster which has no effect, to test model selection
barrier=int
values(barrier) <- 0
barrier[,10:11] <- 1
plot(barrier)
TL.all=get.TL(stack(int,cov.ras,barrier))
##
## sampling locations
##
nsamps=150
set.seed(4567)
samp.xy=cbind(runif(nsamps,0,30),runif(nsamps,0,30))
## samp.xy=samp.xy[rep(1:nsamps,times=5),]
samp.locs=cellFromXY(int,samp.xy)
samp.cells=unique(samp.locs)
nsamps=nrow(samp.xy)
plot(cov.ras)
points(samp.xy)
K=matrix(0,nrow=nsamps,ncol=length(samp.cells))
for(i in 1:nsamps){
K[i,which(samp.cells==samp.locs[i])]=1
}
image(K)
##
## beta values
##
betas=c(-2,-1)
tau=.01
Q=get.Q(TL,betas)
Phi=get.Phi(Q,samp.cells)
## simulate from ibr model
D.rand.ibr=rGenWish(df=20,Sigma=K%*%ginv(as.matrix(Phi))%*%t(K) + diag(nsamps)*tau)
image(D.rand.ibr)
## crude plot of geographic distance vs genetic distance
plot(as.numeric(as.matrix(dist(xyFromCell(ras,samp.locs)))),as.numeric(D.rand.ibr))
##
## fitting using MCMC
##
fit=mcmc.wish.icar(D.rand.ibr,TL,samp.locs,df=20,output.trace.plot=TRUE,
adapt.int=100,adapt.max=100000,n.mcmc=10000)
str(fit)
## End(Not run)