# mcmc.wish.icar: Markov chain Monte Carlo sampler for Generalized Wishart... In rwc: Random Walk Covariance Models

## Description

Constructs and runs an MCMC sampler to estimate resistance parameters of different landscape features.

## Usage

 ```1 2 3 4 5 6 7 8``` ```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.

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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107``` ```## 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) ```

rwc documentation built on May 2, 2019, 3:34 p.m.