# R/par_eCAR.R In eCAR: Eigenvalue CAR Models

#### Documented in par.eCAR.Leroux

```# Wrapper to fit the joint Leroux model
par.eCAR.Leroux <- function(y,x,W,
E=NULL,
C=NULL,
model="Gaussian",
joint_prior_lamx_lamz = FALSE,
lamx.fix.val = NULL, sig2x.fix.val = NULL,
mb=0,s2b=10,
mg=0,s2g=10.0,
alamx=1,blamx=1,
alamz=1,blamz=1,
asig=1,bsig=1,
atau=1,btau=1,
asigx=1,bsigx=1,
mb0=0,s2b0=100,
me=0,s2e=100,
mx=0, s2x=100,
tau_cand_sd=1, sig2_cand_sd=1,
draws=10000,burn=5000,thin=5, verbose=TRUE){
# W - is the neighborhood matrix
# C - is the matrix of confounder covariates that are included not

# Note that y and x are both measure from spatial domain.
# That is, they haven't been projected to spectral domain.
#
# model specifies the data model with options being
# 1 - Gaussian
# 2 - Poisson
# 3 - Binomial
# 4 - Negative Binomial (Poisson with overdispersion)

# HPD calculator (see TeachingDemos)
emp.hpd <- function(x, conf=0.95){
conf <- min(conf, 1-conf)
n <- length(x)
nn <- round( n*conf )
x <- sort(x)
xx <- x[ (n-nn+1):n ] - x[1:nn]
m <- min(xx)
nnn <- which(xx==m)[1]
return( c( x[ nnn ], x[ n-nn+nnn ] ) )
}

cat("A", model, "model is being fit \n")
out <- NULL

nout <- (draws-burn)/thin

nobs <- length(y)

M <- diag(apply(W,1,sum))
R <- M-W

eigendec.R <- eigen(R)
evals <- eigendec.R\$values
evecs <- eigendec.R\$vectors

ystar <- t(evecs) %*% y
xstar <- t(evecs) %*% x

ncov <- ncol(C)
if(is.null(C)){
ncov <- 0
C <- 0
Cstar <- 0
} else {
Cstar <- t(evecs) %*% as.matrix(C)
}

# Assign starting values for sig2x and lamx
# if they are not supplied by the user
sig2x_fix<-sig2x.fix.val; lamx_fix<-lamx.fix.val;
updateXparms <- FALSE;
if(is.null(lamx.fix.val) | is.null(sig2x.fix.val)){
updateXparms <- TRUE
lamx_fix <- 0.5
sig2x_fix <- 1
}

modelPriors = c(mb, s2b, alamx, blamx, alamz, blamz, asig, bsig,
atau, btau, asigx, bsigx, mb0, s2b0, me, s2e, mx, s2x,
mg, s2g)
MHsd <- c(tau_cand_sd, sig2_cand_sd)
beta0 <- beta <- alpha <- tau <- sig2 <- rep(1, nout)
sig2x <- rep(sig2x_fix,nout)
lamx <- lamz <- rep(lamx_fix, nout)
theta <- nb_r <- matrix(0, nrow=nout, ncol=nobs)
eta <- matrix(0, nrow=nout, ncol=ncov)

if(model=="Gaussian"){
# I first center the X and the Y to set the intercept to zero
ystar <- ystar - mean(ystar)
xstar <- xstar - mean(xstar)

C.out <- .C("mcmcloop_leroux_gauss",
as.integer(draws), as.integer(burn), as.integer(thin),
as.integer(nobs), as.double(ystar), as.double(xstar),
as.double(evals), as.double(t(Cstar)), as.integer(ncov),
as.double(modelPriors), as.double(MHsd),
as.integer(verbose), as.integer(joint_prior_lamx_lamz),
as.integer(updateXparms), as.double(lamx_fix),
as.double(sig2x_fix),
beta.out=as.double(beta), alpha.out=as.double(alpha),
tau.out=as.double(tau), sig2x.out=as.double(sig2x),
lamx.out=as.double(lamx), lamz.out=as.double(lamz),
sig2.out=as.double(sig2), eta.out=as.double(eta))

}
if(model!="Gaussian"){
if(model == "Poisson") modelnum = 1
if(model == "Binomial") modelnum = 2
if(model == "Negative Binomial") modelnum = 3

C.out <- .C("mcmcloop_leroux_GLM",
as.integer(draws), as.integer(burn), as.integer(thin),
as.integer(nobs), as.double(y), as.double(x), as.double(E),
as.double(evals), as.double(t(evecs)),
as.double(t(W)), as.double(t(C)), as.integer(ncov),
as.integer(modelnum), as.double(modelPriors),
as.double(MHsd), as.integer(verbose),
as.integer(joint_prior_lamx_lamz),
as.integer(updateXparms),
beta.out=as.double(beta), alpha.out=as.double(alpha),
tau.out=as.double(tau), sig2x.out=as.double(sig2x),
lamx.out=as.double(lamx), lamz.out=as.double(lamz),
theta.out=as.double(theta), beta0.out=as.double(beta0),
eta.out=as.double(eta), r.out=as.double(nb_r))

}

out\$beta <- matrix(C.out\$beta.out, nrow=nout, byrow=TRUE)
out\$gamma <- matrix(C.out\$alpha.out, nrow=nout, byrow=TRUE)
out\$tau <- matrix(C.out\$tau.out, nrow=nout, byrow=TRUE)
out\$sig2x <- matrix(C.out\$sig2x.out, nrow=nout, byrow=TRUE)
out\$lamx <- matrix(C.out\$lamx.out, nrow=nout, byrow=TRUE)
out\$lamz <- matrix(C.out\$lamz.out, nrow=nout, byrow=TRUE)
if(model=="Gaussian") out\$sig2 <- matrix(C.out\$sig2.out, nrow=nout, byrow=TRUE)
out\$rho <- matrix((out\$gamma*sqrt(out\$sig2x))/sqrt(out\$tau + out\$gamma^2*out\$sig2x), nrow=nout, byrow=TRUE)
out\$sig2z <- matrix(out\$tau + out\$gamma^2*out\$sig2x, nrow=nout, byrow=TRUE)
if(model!="Gaussian"){
out\$theta <- matrix(C.out\$theta.out, nrow=nout, byrow=TRUE)
out\$beta0 <- matrix(C.out\$beta0.out, nrow=nout, byrow=TRUE)
}
if(ncov>0) out\$eta <- matrix(C.out\$eta.out, nrow=nout, byrow=TRUE)
if(model=="Negative Binomial") out\$nb_r = matrix(C.out\$r.out, nrow=nout, byrow=TRUE)

# Produce beta as a function of eigen-value
Dseq <- seq(min(evals), max(evals), length=1000)
c.beta <- matrix(NA, nrow=nout, ncol=length(Dseq))
for(t in 1:nout){
c.beta[t,] <- out\$beta[t] + out\$gamma[t]*sqrt((1-out\$lamx[t]+out\$lamx[t]*Dseq)/(1-out\$lamz[t]+out\$lamz[t]*Dseq))
}

if(model=="Gaussian"){
beta.mn <-  matrix(apply(c.beta,2,function(x) mean(x)),nrow=1000,byrow=TRUE)
beta.q025 <-  matrix(apply(c.beta,2,function(x) emp.hpd(x))[1,],nrow=1000,byrow=TRUE)
beta.q975 <-  matrix(apply(c.beta,2,function(x) emp.hpd(x))[2,],nrow=1000,byrow=TRUE)
}
if(model!="Gaussian"){
beta.mn <-  matrix(apply(exp(c.beta),2,function(x) mean(x)),nrow=1000,byrow=TRUE)
beta.q025 <-  matrix(apply(exp(c.beta),2,function(x) emp.hpd(x))[1,],nrow=1000,byrow=TRUE)
beta.q975 <-  matrix(apply(exp(c.beta),2,function(x) emp.hpd(x))[2,],nrow=1000,byrow=TRUE)
}
omega <- Dseq
result <- eCAR.out(data_model = model,
beta_omega = cbind(beta.mn, beta.q025, beta.q975, omega),
posterior_draws = out,
DIC=NULL,
regrcoef = NULL)

return(result)
}
```

## Try the eCAR package in your browser

Any scripts or data that you put into this service are public.

eCAR documentation built on April 5, 2023, 5:09 p.m.