est_logLik_occu <- function(A_0,thr_0,Y # all settings for theta_0 model
,A,thr,Z){
# Descombes 1999 paper, calculate Z_theta/Z_theta0, using samples Y from theta_0, I will use posterior mean.
A_dif <- as(A_0-A,'dsCMatrix')
rm(A_0)
thr_diff <- thr_0-thr
Ham <- apply(Y,1,function(X,graph,s){(H(graph,X,s))},graph = A_dif,s = thr_diff) # Hamiltonian with parameters as difference and sample of theta_0
robust_cri <- (max(Ham)-min(Ham))/max(Ham)
partition_ratio <- mean(
exp(-Ham)
)# ratio of partitioning function: Z_theta/Z_theta_0
# this is the estimation of a single theta, we need to loop through every sample
logLik_est <- -H(A,Z,thr)-log(partition_ratio)
return(list(logLik_est=logLik_est,robust_cri=robust_cri))
}
logLik_full <- function(theta,envX,detX,distM,link_map,dist_mainland,link_mainland,int_range_intra="nn",int_range_inter="exp",Z, detmat, A_0,thr_0,Y){
## logLik for detection part
nspp <- nrow(theta$spp_mat)
nperiod <- ncol(detmat[[1]])
beta_det <- theta$beta_det
logLdata <- Pdet_Ising_rep(nrep,nperiod,envX,detX,theta$beta_det,theta$spp_mat_det,as.matrix(Z),detmat)
## now estimate the logLik for occu
nsite <- nrow(envX)
beta_occu <- theta$beta_occu
eta_intra <- theta$eta_intra # intra spp, intra island if apply
d_intra <- theta$d_intra
spp_mat <- theta$spp_mat
nspp <- nrow(spp_mat)
nrep <- ncol(Z)
A_in <- getintralayerGraph(distM,link_map$intra,eta_intra,d_intra,int_range = int_range_intra,spp_mat)
eta_inter <- theta$eta_inter # assume there is one
d_inter <- theta$d_inter
A_ex <- getintralayerGraph(distM,link_map$inter,eta_inter,d_inter,int_range = int_range_inter,spp_mat) # graph among islands, if apply, distM should only contain graph
A <- getfullGraph(A_ex,A_in,spp_mat)
rm(A_ex,A_in)
thr <- lapply(1:nspp,
function(i,envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter){
envX %*% beta_occu[1:ncol(envX)+(i-1)*ncol(envX)] +
mainland_thr(dist_mainland,link_mainland,eta_inter[i],d_inter[i],int_range_inter)
},envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter)
thr <- Reduce(rbind,thr)
logLoccu <- est_logLik_occu(A_0,thr_0,Y,A,thr,Z)
# sum them up
return(list(logLik = logLoccu$logLik_est + logLdata,robust = logLoccu$robust_cri))
}
make_list_version_mcmc <- function(mcmc_list,theta_0){ # gonna return a list, with single element is a theta list
n_mcmc <- nrow(mcmc_list[[1]])
nspp <- sqrt( ncol(mcmc_list[[length(mcmc_list)-1]]))
nparas <- length(mcmc_list)
lapply(1:n_mcmc,function(i,mcmc_list,theta_0,nspp,nparas){
temp <- lapply(1:(length(theta_0)-2),function(k,mcmc_list,i){
mcmc_list[[k]][i,]
},mcmc_list,i)
names(temp) <- names(theta_0)[1:(length(theta_0)-2)]
temp$spp_mat <- Matrix(mcmc_list[[nparas-1]][i,],nspp,nspp,sparse = T) # sppmat, need to be formatted
temp$spp_mat_det <- Matrix(mcmc_list[[nparas]][i,],nspp,nspp,sparse = T)
return(temp)
},mcmc_list,theta_0,nspp,nparas)
}
# passed Aug 12 2019
# we calculate the Bayes factor for model selection
logBF <- function(theta_a_mcmc,envX_a,distM,link_map_a,dist_mainland,link_mainland_a,int_range_intra_a="nn",int_range_inter_a="exp",Z_a_mcmc, detX_a, theta_a_point
,theta_mcmc,envX,link_map,link_mainland,int_range_intra="nn",int_range_inter="exp",Z_mcmc, detX, theta_point, detmat, nrep, nY = 3000,nIter = 100,method = "CFTP"){
nsite <- nrow(envX)
beta_occu_point <- theta_point$beta_occu
eta_intra_point <- theta_point$eta_intra # intra spp, intra island if apply
d_intra_point <- theta_point$d_intra
spp_mat_point <- theta_point$spp_mat
nspp <- nrow(spp_mat_point)
A_in_point <- getintralayerGraph(distM,link_map$intra,eta_intra_point,d_intra_point,int_range = int_range_intra,spp_mat_point)
eta_inter_point <- theta_point$eta_inter # assume there is one
d_inter_point <- theta_point$d_inter
A_ex_point <- getintralayerGraph(distM,link_map$inter,eta_inter_point,d_inter_point,int_range = int_range_inter,spp_mat_point) # graph among islands, if apply, distM should only contain graph
A_point <- getfullGraph(A_ex_point,A_in_point,spp_mat_point)
rm(A_ex_point,A_in_point)
thr_point <- lapply(1:nspp,
function(i,envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter){
envX %*% beta_occu[1:ncol(envX)+(i-1)*ncol(envX)] +
mainland_thr(dist_mainland,link_mainland,eta_inter[i],d_inter[i],int_range_inter)
},envX,beta_occu_point,dist_mainland,link_mainland,eta_inter_point,d_inter_point,int_range_inter)
thr_point = Reduce(rbind,thr_point)
beta_occu_0_point <- theta_a_point$beta_occu
eta_intra_0_point <- theta_a_point$eta_intra # intra spp, intra island if apply
d_intra_0_point <- theta_a_point$d_intra
spp_mat_0_point <- theta_a_point$spp_mat
A_in_0_point <- getintralayerGraph(distM,link_map_a$intra,eta_intra_0_point,d_intra_0_point,int_range = int_range_intra_a,spp_mat_0_point)
eta_inter_0_point <- theta_a_point$eta_inter # assume there is a
d_inter_0_point <- theta_a_point$d_inter
A_ex_0_point <- getintralayerGraph(distM,link_map_a$inter,eta_inter_0_point,d_inter_0_point,int_range = int_range_inter_a,spp_mat_0_point) # graph among islands, if apply, distM should only contain graph
A_0_point <- getfullGraph(A_ex_0_point,A_in_0_point,spp_mat_0_point)
rm(A_in_0_point,A_ex_0_point)
thr_0_point <- lapply(1:nspp,
function(i,envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter){
envX %*% beta_occu[1:ncol(envX)+(i-1)*ncol(envX)] +
mainland_thr(dist_mainland,link_mainland,eta_inter[i],d_inter[i],int_range_inter)
},envX_a,beta_occu_0_point,dist_mainland,link_mainland_a,eta_inter_0_point,d_inter_0_point,int_range_inter_a)
thr_0_point <- Reduce(rbind,thr_0_point)
A_mean <- as(0.5*(A_point + A_0_point),'dsCMatrix') # theta in 1999 paper, as the constant, we estimate ratio of partition function of theta_a/theta_mean and theta/theta_mean
thr_mean <- 0.5*(thr_point + thr_0_point)
Y <- IsingSamplerCpp(n=nY,graph = A_mean, thresholds = thr_mean, beta=1, responses = c(-1L, 1L),nIter = nIter,exact = (method=="CFTP"),constrain = NA+thr_mean)
logLik_a <- lapply(1:length(theta_a_mcmc),
function(i,theta,envX,detX
,distM,link_map,dist_mainland
,link_mainland,int_range_intra
,int_range_inter,Z, detmat, A_0,thr_0,Y){
logLik_full(theta[[i]],envX,detX
,distM,link_map,dist_mainland
,link_mainland,int_range_intra,int_range_inter
,matrix(Z[i,],ncol = nrep), detmat, A_0,thr_0,Y)
},theta_a_mcmc,envX_a,detX_a,distM,link_map_a,dist_mainland
,link_mainland_a,int_range_intra_a
,int_range_inter_a,Z_a_mcmc, detmat, A_mean,thr_mean,Y)
logLik_theta <- lapply(1:length(theta_mcmc),
function(i,theta,envX,detX
,distM,link_map,dist_mainland
,link_mainland,int_range_intra
,int_range_inter,Z, detmat, A_0,thr_0,Y){
logLik_full(theta[[i]],envX,detX
,distM,link_map,dist_mainland
,link_mainland,int_range_intra,int_range_inter
,matrix(Z[i,],ncol = nrep), detmat, A_0,thr_0,Y)
},theta_mcmc,envX,detX,distM,link_map,dist_mainland
,link_mainland,int_range_intra
,int_range_inter,Z_mcmc, detmat, A_mean,thr_mean,Y)
rm(Y,A_mean,thr_mean)
pure_logLik_a <- lapply(logLik_a,function(a){a[[1]]})
robust_cri_a <- lapply(logLik_a,function(a){a[[2]]})
pure_logLik_theta <- lapply(logLik_theta,function(a){a[[1]]})
robust_cri_theta <- lapply(logLik_theta,function(a){a[[2]]})
logLik_a <- Reduce(rbind,pure_logLik_a)
logLik_theta <- Reduce(rbind,pure_logLik_theta)
#pD_Gelman04_a = 2 * var(logLik_a)
#pD_Gelman04_t = 2 * var(logLik_theta)
#deltaDIC = - 2 * mean(logLik_a) + 2 * (pD_Gelman04_a) + mean(logLik_theta) - 2 * (pD_Gelman04_t)
deltalogBF <- - log(mean(1/exp(logLik_a),na.rm=T)) + log( mean(1/exp(logLik_theta),na.rm=T) )
robust_cri <- list(robust_a = Reduce(mean,robust_cri_a),robust_theta = Reduce(mean,robust_cri_theta))
return(list(logBF = deltalogBF,robust_criterion = robust_cri))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.