Nothing
#==================
# define functions
#==================
tnorm <- function(n, lo, hi, mu, sig) {
# generates truncated normal variates based on cumulative normal distribution normal truncated lo and hi
if (length(lo) == 1 & length(mu) > 1)
lo <- rep(lo, length(mu))
if (length(hi) == 1 & length(mu) > 1)
hi <- rep(hi, length(mu))
q1 <- pnorm(lo, mu, sig) #cumulative distribution
q2 <- pnorm(hi, mu, sig) #cumulative distribution
z <- runif(n, q1, q2)
z <- qnorm(z, mu, sig)
z[z == -Inf] <- lo[z == -Inf]
z[z == Inf] <- hi[z == Inf]
z
} # tnorm
fix.MH <- function(lo, hi, old1, new1, jump) {
jold <- pnorm(hi, mean = old1, sd = jump) - pnorm(lo, mean = old1, sd = jump)
jnew <- pnorm(hi, mean = new1, sd = jump) - pnorm(lo, mean = new1, sd = jump)
log(jold) - log(jnew) #add this to pnew
} # fix.MH
get.logl <- function(theta, phi, pobs.mat, y, nmat) {
prob <- theta %*% phi * pobs.mat
cond <- prob < 1e-05
prob[cond] <- 1e-05
cond <- prob > 0.99999
prob[cond] <- 0.99999
dbinom(y, size = nmat, prob = prob, log = T)
} # get.logl
acceptMH <- function(p0, p1, x0, x1, BLOCK) {
# accept for M, M-H if BLOCK, then accept as a block, otherwise, accept individually
nz <- length(x0) #no. to accept
if (BLOCK)
nz <- 1
a <- exp(p1 - p0) #acceptance PR
z <- runif(nz, 0, 1)
keep <- which(z < a)
if (BLOCK & length(keep) > 0)
x0 <- x1
if (!BLOCK)
x0[keep] <- x1[keep]
accept <- length(keep)
list(x = x0, accept = accept)
} # acceptMH
update.theta <- function(param, jump, nl, nc, ns, y, x, nmat, id.binomial) {
phi <- param$phi
beta <- param$beta
mu <- x %*% beta
pobs <- param$pobs
pobs.mat <- matrix(pobs, nl, ns, byrow=T)
pobs.mat[id.binomial,] <- 1
vmat.ori <- vmat.old <- param$vmat
vmat.tmp <- exp(rnorm(nl*(nc-1), mean=log(vmat.old[,-nc]), sd=jump[,-nc]))
proposed <- cbind(matrix(vmat.tmp, nl, nc-1), 1)
for (j in 1:(nc-1)) {
# last column has to be 1
vmat.new <- vmat.old
vmat.new[,j] <- proposed[,j]
theta.old <- vmat.old / rowSums(vmat.old)
theta.new <- vmat.new / rowSums(vmat.new)
prob.old <- get.logl(theta=theta.old, phi=phi, pobs.mat=pobs.mat, y=y, nmat=nmat)
prob.new <- get.logl(theta=theta.new, phi=phi, pobs.mat=pobs.mat, y=y, nmat=nmat)
pold <- rowSums(prob.old) + dnorm(log(vmat.old[,j]), mu[,j], sigma, log = T)
pnew <- rowSums(prob.new) + dnorm(log(vmat.new[,j]), mu[,j], sigma, log = T)
k <- acceptMH(p0=pold, p1=pnew, x0=vmat.old[,j], x1=vmat.new[,j], BLOCK=F)
vmat.old[,j] <- k$x
}
vmat <- vmat.old
theta <- vmat / rowSums(vmat)
list(theta=theta, vmat=vmat, accept=vmat.ori!=vmat.old)
} # update.theta
update.phi <- function(param, jump, nl, nc, ns, y, nmat, id.binomial, a.phi, b.phi) {
theta <- param$theta
pobs <- param$pobs
pobs.mat <- matrix(pobs, nl, ns, byrow=T)
pobs.mat[id.binomial,] <- 1
phi.ori <- phi.old <- param$phi
proposed <- matrix(tnorm(nc*ns, lo=0, hi=1, mu=phi.old, sig=jump), nc, ns)
adj <- fix.MH(lo=0, hi=1, old1=phi.old, new1=proposed, jump=jump)
for (j in 1:nc) {
phi.new <- phi.old
phi.new[j,] <- proposed[j,]
prob.old <- get.logl(theta=theta, phi=phi.old, pobs.mat=pobs.mat, y=y, nmat=nmat)
prob.new <- get.logl(theta=theta, phi=phi.new, pobs.mat=pobs.mat, y=y, nmat=nmat)
pold <- colSums(prob.old) + dbeta(phi.old[j,], a.phi, b.phi, log=T)
pnew <- colSums(prob.new) + dbeta(phi.new[j,], a.phi, b.phi, log=T)
k <- acceptMH(p0=pold, p1=pnew + adj[j,], x0=phi.old[j,], x1=phi.new[j,], BLOCK=F)
phi.old[j,] <- k$x
}
phi <- phi.old
list(phi=phi, accept=phi.ori!=phi.old)
} # update.phi
update.lpmu <- function(param, jump, ns) {
pobs <- param$pobs
lpmu.old <- param$lpmu
lpsigma.old <- param$lpsigma
lpmu.new <- rnorm(1, lpmu.old, jump)
lpsigma.new <- exp(rnorm(1, log(lpsigma.old), jump))
pold <- sum(dnorm(logit(pobs), lpmu.old, lpsigma.old, log=T)) +
dnorm(lpmu.old, 0, 10, log=T) + dnorm(log(lpsigma.old), 0, 10, log=T)
pnew <- sum(dnorm(logit(pobs), lpmu.new, lpsigma.new, log=T)) +
dnorm(lpmu.new, 0, 10, log=T) + dnorm(log(lpsigma.new), 0, 10, log=T)
k <- acceptMH(p0=pold, p1=pnew, x0=c(lpmu.old,lpsigma.old), x1=c(lpmu.new,lpsigma.new), BLOCK=F)
lpmu <- k$x[1]
lpsigma <- k$x[2]
list(lpmu=lpmu, lpsigma=lpsigma)
} # update.lpmu
update.pobs <- function(param, jump, nl, ns, y, nmat, id.binomial) {
theta <- param$theta
phi <- param$phi
lpmu <- param$lpmu
lpsigma <- param$lpsigma
pobs.old <- param$pobs
pobs.new <- inv.logit(rnorm(ns, logit(pobs.old), jump))
pobs.new <- ifelse(pobs.new > .99, .99, pobs.new)
pobs.new <- ifelse(pobs.new < .01, .01, pobs.new)
pobs.old.mat <- matrix(pobs.old, nl, ns, byrow=T)
pobs.old.mat[id.binomial,] <- 1
pobs.new.mat <- matrix(pobs.new, nl, ns, byrow=T)
pobs.new.mat[id.binomial,] <- 1
prob.old <- get.logl(theta=theta, phi=phi, pobs.mat=pobs.old.mat, y=y, nmat=nmat)
prob.new <- get.logl(theta=theta, phi=phi, pobs.mat=pobs.new.mat, y=y, nmat=nmat)
pold <- sum(prob.old) + sum(dnorm(logit(pobs.old), lpmu, lpsigma, log=T))
pnew <- sum(prob.new) + sum(dnorm(logit(pobs.new), lpmu, lpsigma, log=T))
k <- acceptMH(p0=pold, p1=pnew, x0=pobs.old, x1=pobs.new, BLOCK=F)
pobs <- k$x
list(pobs=pobs, accept=pobs!=pobs.old)
} # update.pobs
update.beta <- function(param, jump, nl, nc, x) {
nx <- dim(x)[2]
vmat <- param$vmat
beta.old <- param$beta
beta.new <- matrix(rnorm(length(beta.old), beta.old, jump), nx, nc)
mu.old <- x %*% beta.old
mu.new <- x %*% beta.new
prob.old <- dnorm(log(vmat[,-nc]), mu.old[,-nc], sigma, log=T)
prob.new <- dnorm(log(vmat[,-nc]), mu.new[,-nc], sigma, log=T)
pold <- sum(prob.old) + sum(dnorm(beta.old[-nc], 0, 5, log=T))
pnew <- sum(prob.new) + sum(dnorm(beta.new[-nc], 0, 5, log=T))
k <- acceptMH(p0=pold, p1=pnew, x0=beta.old, x1=beta.new, BLOCK=F)
beta <- k$x
list(beta=beta, accept=beta!=beta.old)
} # update.beta
thetaphi.jumpTune <- function(accept, jump, ni, adapt=2000, low=.3, high=.8) {
nstart <- ifelse(ni>=100, ni-99, 1)
if (ni > adapt) {
jump.new <- jump
} else if (ni %% 10 > 0) {
jump.new <- jump
} else {
accept.rate <- apply(accept[,,nstart:ni], 1:2, mean)
jump.new <- jump
jump.new[which(accept.rate < low)] <-
jump[which(accept.rate < low)] * rnorm(length(which(accept.rate < low)),.5,.01)
jump.new[which(accept.rate > high)] <-
jump[which(accept.rate > high)] * rnorm(length(which(accept.rate > high)),2,.01)
}
jump <- jump.new
list(jump=jump)
} # jumpTune
pobsbeta.jumpTune <- function(accept, jump, ni, adapt=2000, low=.3, high=.8) {
nstart <- ifelse(ni>=100, ni-99, 1)
if (ni > adapt) {
jump.new <- jump
} else if (ni %% 10 > 0) {
jump.new <- jump
} else {
accept.rate <- mean(accept[nstart:ni], na.rm=T)
if (accept.rate < low) {
jump.new <- jump * rnorm(1,.5,.05)
} else if (accept.rate > high) {
jump.new <- jump * rnorm(1,2,.05)
} else {
jump.new <- jump
}
}
jump <- jump.new
list(jump=jump)
} # pobsbeta.jumpTune
#============================
# define regression function
#============================
rlda.regression <- function(y=y, x=x, nmat=nmat, id.binomial=id.binomial, nc=nc,
niter=niter, nburn=nburn, adapt=adapt,
vmat.init=vmat.init, phi.init=phi.init, beta.init=beta.init,
logit_pobs_mu.init=logit_pobs_mu.init,
logit_pobs_sigma.init=logit_pobs_sigma.init,
pobs.init=pobs.init) {
#===========
# run model
#===========
nl <- dim(y)[1] # number of locations
ns <- dim(y)[2] # number of species
nx <- dim(x)[2] # number of covariates
theta.jump <- matrix(.05, nrow=nl, ncol=nc)
phi.jump <- matrix(.05, nrow=nc, ncol=ns)
pobs.jump <- .05
beta.jump <- .05
param <- list()
param$vmat <- vmat.init
param$theta <- vmat.init / rowSums(vmat.init)
param$phi <- phi.init
param$beta <- beta.init
param$lpmu <- logit_pobs_mu.init
param$lpsigma <- logit_pobs_sigma.init
param$pobs <- pobs.init
vmat.post <- theta.post <- theta.jump.post <- theta.accept <- array(, dim=c(nl, nc, niter))
vmat.post[,,1] <- param$vmat
theta.post[,,1] <- param$theta
theta.jump.post[,,1] <- theta.jump
theta.accept[,,1] <- FALSE
phi.post <- phi.jump.post <- phi.accept <- array(, dim=c(nc, ns, niter))
phi.post[,,1] <- param$phi
phi.jump.post[,,1] <- phi.jump
phi.accept[,,1] <- FALSE
lpmu.post <- lpsigma.post <- numeric(niter)
lpmu.post[1] <- param$lpmu
lpsigma.post[1] <- param$lpsigma
pobs.post <- matrix(, ns, niter)
pobs.post[,1] <- param$pobs
pobs.jump.post <- numeric(niter)
pobs.jump.post[1] <- pobs.jump
pobs.accept <- logical(niter)
beta.post <- array(, dim=c(nx,nc,niter))
beta.post[,,1] <- param$beta
beta.jump.post <- numeric(niter)
beta.jump.post[1] <- beta.jump
beta.accept <- logical(niter)
for (i in 2:niter) {
theta.up <- update.theta(param, jump=theta.jump, nl, nc, ns, y, x, nmat, id.binomial)
vmat.post[,,i] <- param$vmat <- theta.up$vmat
theta.post[,,i] <- param$theta <- theta.up$theta
theta.accept[,,i] <- theta.up$accept
phi.up <- update.phi(param, jump=phi.jump, nl, nc, ns, y, nmat, id.binomial, a.phi=1, b.phi=1)
phi.post[,,i] <- param$phi <- phi.up$phi
phi.accept[,,i] <- phi.up$accept
theta.jump.up <- thetaphi.jumpTune(accept=theta.accept, jump=theta.jump, ni=i, adapt=adapt, low=.3, high=.8)
theta.jump.post[,,i] <- theta.jump <- theta.jump.up$jump
phi.jump.up <- thetaphi.jumpTune(accept=phi.accept, jump=phi.jump, ni=i, adapt=adapt, low=.3, high=.8)
phi.jump.post[,,i] <- phi.jump <- phi.jump.up$jump
lpmu.up <- update.lpmu(param, jump=.1, ns)
lpmu.post[i] <- param$lpmu <- lpmu.up$lpmu
lpsigma.post[i] <- param$lpsigma <- lpmu.up$lpsigma
pobs.up <- update.pobs(param, jump=pobs.jump, nl, ns, y, nmat, id.binomial)
pobs.post[,i] <- param$pobs <- pobs.up$pobs
pobs.accept[i] <- pobs.up$accept
pobs.jump.up <- pobsbeta.jumpTune(accept=pobs.accept, jump=pobs.jump, ni=i, adapt=adapt, low=.3, high=.8)
pobs.jump.post[i] <- pobs.jump <- pobs.jump.up$jump
beta.up <- update.beta(param, jump=beta.jump, nl, nc, x)
beta.post[,,i] <- param$beta <- beta.up$beta
beta.accept[i] <- beta.up$accept
beta.jump.up <- pobsbeta.jumpTune(accept=beta.accept, jump=beta.jump, ni=i, adapt=adapt, low=.3, high=.8)
beta.jump.post[i] <- beta.jump <- beta.jump.up$jump
}
#==============
# save results
#==============
theta.est <- apply(theta.post[,,nburn:niter], 1:2, median)
phi.est <- apply(phi.post[,,nburn:niter], 1:2, median)
beta.est <- apply(beta.post[,,nburn:niter], 1:2, median)
lpmu.est <- median(lpmu.post[nburn:niter])
lpsigma.est <- median(lpsigma.post[nburn:niter])
pobs.est <- apply(pobs.post[,nburn:niter], 1, median)
res <- list()
res$theta <- theta.est
res$phi <- phi.est
res$beta <- beta.est
res$lpmu <- lpmu.est
res$lpsigma <- lpsigma.est
res$pobs <- pobs.est
class(res) <- c("rlda", "list")
res
} # rlda.regression
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.