From 54813f53bfb54fefbfd9905bc6c5e78817e09eed Mon Sep 17 00:00:00 2001 From: beth-ross Date: Fri, 19 Jul 2024 12:22:53 -0600 Subject: [PATCH] cleaned up notes --- Analysis/krst_ipm.R | 304 ++------------------------------------------ 1 file changed, 13 insertions(+), 291 deletions(-) diff --git a/Analysis/krst_ipm.R b/Analysis/krst_ipm.R index c0de1af..848d251 100644 --- a/Analysis/krst_ipm.R +++ b/Analysis/krst_ipm.R @@ -8,7 +8,7 @@ library(stringr) cmr = read_xls("NPS CMR data updated.xls",sheet = 1, col_names = TRUE) cmr.df = as.data.frame(cmr) -#need to clean up id numbers +#clean up id numbers #first remove "2 transmitters" and "3 transmitters" cmr.df$id[cmr.df$id=="2 transmitters"] = NA @@ -17,6 +17,7 @@ cmr.df$id[cmr.df$id=="3 transmitters"] = NA #need to add in 114 which got renumbered cmr.df$id[715] = 114 +#correct numbers for individuals identified by genetics idx = str_sort(cmr.df$id,na_last=TRUE,decreasing=TRUE) cmr.df$id[which(cmr.df$id==idx[10])]=813 cmr.df$id[which(cmr.df$id==idx[7])]=814 @@ -77,7 +78,6 @@ cmr.df$id2 = id #fix dates using stringr idx2 = str_split(cmr.df$clutch,"-",simplify = TRUE)[,1] -#cmr.df$date_nest[which(idx2== "MEXICO")] #made up date (kept year) to match other dates cmr.df$date_nest[which(idx2== "MEXICO")][7] = "16 May 2015" @@ -93,50 +93,18 @@ idx2[which(idx2 == idx.tmp[i])] = NA} cmr.df$year = as.numeric(idx2) -#cmr.df$date_nest[11] = "13 June 1997" -#cmr.df$date_nest[21] = "19 April 2000" -#cmr.df$date_nest[22] = "26 April 2002" -#cmr.df$date_nest[23] = "16 May 2002" -#cmr.df$date_nest[24] = "26 April 2006" -#cmr.df$date_nest[25] = "11 May 2006" -#cmr.df$date_nest[26] = "26 April 2008" -#cmr.df$date_nest[27] = "30 April 2008" -#cmr.df$date_nest[28] = "16 May 2008" -#cmr.df$date_nest[29] = "1 June 2008" -#cmr.df$date_nest[30] = "18 June 2012" -#cmr.df$date_nest[137] = "3 June 2008" #was 1008 -#cmr.df$date_nest[215] = "1 May 2015" #no day or month in excel file -#cmr.df$date_nest[468] = "30 April 2007" -#cmr.df$date_nest[1761] = "1 May 2012" #only 'in situ' recorded for 2012 -#cmr.df$date_nest[1818] = "7 May 2007" -#cmr.df$date_nest[1819] = "22 May 2007" -#cmr.df$date_nest[1921] = "1 May 2013" #only 'in situ' recorded for 2013 -#cmr.df$date_nest[1941] = "3 July 2018" -#cmr.df$date_nest[2278] = '1 May 2008' #'in situ' for 2008 -#cmr.df$date_nest[2776] = '8 August 2018' - dat = matrix(0, nrow=822, ncol = (2023-1991)) -#wild = rep(NA, 812) for(i in 1:822){ tmp = subset(cmr.df,id2==i) - #tmp2 = as.Date(tmp$date_nest, format="%d %B %Y") - #years = year(tmp2) - dat[i,(tmp$year-1990)] = 1 - #wild[i] = max(ifelse(tmp$yr_hatch == "W",1,0),na.rm=TRUE) -} - -#wild=wild[-c(114, 227, 291, 506, 556, 588, 604,801)] #renumbered individuals -#wild[!is.finite(wild)] = 1 + dat[i,(tmp$year-1990)] = 1 + } dat = dat[-(which(apply(dat,1,sum)==0)),] #remove individuals that were renumbered in dataset (missing ids: 114, 227, 291, 506, 556, 588, 604,801) tx_nest_dat = read.csv("tx_nest_hatch_data.csv") colnames(tx_nest_dat)[1] = "year" -#nest searching effort -nest.effort = read_xlsx("patrol effort data_NPI_1986 to 2021.xlsx",sheet = 1, col_names = TRUE) -kilo.patrol = nest.effort$`Kilometers patrolled`[nest.effort$Year>1990] #nest protection data #hatch success = prob of hatching #emerge success = prob of being released @@ -165,7 +133,6 @@ krst.ipm <- nimbleCode({ for(t in 1:((Time+Time.fore)-1)){ ln.mu.phi[t] ~ dnorm(mu.phi,var=sd.phi2) - #ln.mu.phi[t] ~ dnorm(mu.phi, sd = 1) } mu.phi ~ dnorm(2.2, sd = 1) sd.phi2 ~ dinvgamma(2,1) @@ -173,17 +140,8 @@ krst.ipm <- nimbleCode({ psiBNB ~ dbeta(psiBNB.mom[1], psiBNB.mom[2]) #transition of Breeders to NonBreeders psiNBB ~ dbeta(psiNBB.mom[1], psiNBB.mom[2]) #transition of NB to B pB ~ dbeta(p.mom[1], p.mom[2]) #det of breeders - #a.0 ~ dnorm(0,sd = 1) - #a.k ~ dnorm(0,sd = 1) - - #model for detection probability by kilo searched - #for(t in 1:(Time+Time.fore-1)){ - #logit(pB[t]) <- a.0 + a.k*kilo.patrol[t] - #} - #priors pre breeding survival - #s.p1 ~ dbeta(1,1) #same table 3, hatch & 1yr surv s.p1 ~ dbeta(s.p.mom[1],s.p.mom[2]) #priors initial abundance @@ -194,11 +152,6 @@ krst.ipm <- nimbleCode({ N.p5[1] ~ dpois(N.p5.ini) N.p6[1] ~ dpois(N.p6.ini) N.p7[1] ~ dpois(N.p7.ini) - #N.p8[1] ~ dpois(N.p8.ini) - #N.p9[1] ~ dpois(N.p9.ini) - #N.p10[1] ~ dpois(N.p10.ini) - #N.b[1] ~ dpois(N.b.ini) - #N.nb[1] ~ dpois(N.nb.ini) N.b[1] <- 1 N.nb[1] <- 1 @@ -226,12 +179,9 @@ krst.ipm <- nimbleCode({ omega[3,1] <- 1 #pr(dead t -> nondetected) omega[3,2] <- 0 #pr(dead t -> detected) - #} - - ## + #####HMM for adult survival and breeding probability#### - ## - + for(i in 1:N){ #latent state at first capture z[i,first[i]] <- 1 @@ -258,8 +208,6 @@ krst.ipm <- nimbleCode({ beta.temp[6] ~ dnorm(insitu.temp,sd=1) beta.temp[7] ~ dnorm(mu.temp, var = var.temp) beta.slr[6] ~ dnorm(insitu.slr,sd=1) - #beta.temp.all ~ dnorm(mu.temp, var = var.temp) - #beta.temp.all2 ~ dnorm(0,sd=2) mu.temp ~ dnorm(0, sd = 1) var.temp ~ dinvgamma(1,1) @@ -294,18 +242,12 @@ krst.ipm <- nimbleCode({ } f[22] ~ dpois(mean(E.ind[1:year.mean.idx[22]])*mean(prob.egg.ind[1:year.mean.idx[22]])) - #f[25] ~ dpois(mean(E.ind[1:year.mean.idx[25]])*mean(prob.egg.ind[1:year.mean.idx[25]])) for(t in 23:Time){ - #for(t in 26:Time){ f[t] ~ dpois(mean(E.ind[(year.mean.idx[t-1]+1):year.mean.idx[t]])*mean(prob.egg.ind[(year.mean.idx[t-1]+1):year.mean.idx[t]])) } for(t in (Time+1):(Time+Time.fore)){ - #logit(prob.egg.fore.tmp[1:6,(t-31)]) <- nest.cov[1:6] + beta.temp[1:6]*clim.cov.fore[(t-31)] + - # beta.slr[1:6]*slr.cov.fore[(t-31)] - #prob.egg.fore[(t-31)] <- sum(prob.egg.fore.tmp[1:6,(t-31)]*nest.mange.fore[1:6])/sum(nest.mange.fore[1:6]) - #f[t] ~ dpois(E.ind.fore*prob.egg.fore[(t-31)]) logit(prob.egg.fore.tmp[1:6,(t-28)]) <- nest.cov[1:6] + beta.temp[1:6]*clim.cov.fore[(t-28)] + beta.slr[1:6]*slr.cov.fore[(t-28)] prob.egg.fore[(t-28)] <- sum(prob.egg.fore.tmp[1:6,(t-28)]*nest.mange.fore[1:6])/sum(nest.mange.fore[1:6]) @@ -316,27 +258,9 @@ krst.ipm <- nimbleCode({ #####number of nests per female per year#### ## - #for(t in 1:Time){ - #n.fem[t] ~ dgamma(mean = mu.n.fem, sd = 1) - #} - mu.n.fem <- mean(n.fem[]) - #mu.n.fem ~ dunif(1,4) - #sd.n.fem ~ dunif(0,10) - - #mu.n.fem <- 2 - #Probability of hatchling being female - #for(t in 1:Time){ - # prob.fem[t] ~ dbeta(r*mean.fem, r*(1-mean.fem)) #Data from Shaver & Calliway - #logit(mu.fem[t]) <- ln.mean.fem - #prob.fem[t] ~ dbeta(prob.fem.mom[1],prob.fem.mom[2]) - #} - - #ln.mean.fem ~ dnorm(0,sd=10) #prior, mean prob - #mean.fem <- exp(ln.mean.fem)/(1+exp(ln.mean.fem)) - #r ~ dgamma(.1,.1) #overdispersion param for beta dist'n mean.fem <- mean(prob.fem[1:18]) #### state space model for abundance/pop matrix #### @@ -347,13 +271,9 @@ krst.ipm <- nimbleCode({ imm[1] <- 1 N.i[1] <- 1 N.tot[1] <- 2 - #mu.imm ~ dnorm(0, sd = 2) - #sd.imm ~ dinvgamma(2,1) - + for(t in 1:(Time+Time.fore)){ eps.imm.tmp[t] ~ dnorm(mu.imm, sd = 1) - #eps.imm[t] <- min(eps.imm.tmp[t],10) - #eps.imm.tmp[t] ~ dnorm(3, sd = 1) eps.imm[t] <- eps.imm.tmp[t] } @@ -362,27 +282,8 @@ krst.ipm <- nimbleCode({ N.tot[t] <- N.b[t] + N.i[t] #s.p1 is transition probability * survival - pre.br1[t] <- f[t-1] * mean.fem * N.tot[t-1] * s.p1 * mu.n.fem #+ - #(N.p1[t-1]*(1-s.p1)) - #pre.br2[t] <- N.p1[t-1]*s.p1 #+ - #N.p2[t-1]*(1-s.p1) - #pre.br3[t] <- N.p2[t-1]*s.p1 #+ - #N.p3[t-1]*(1-s.p1) - #pre.br4[t] <- N.p3[t-1]*s.p1 #+ - #N.p4[t-1]*(1-s.p1) - #pre.br5[t] <- N.p4[t-1]*s.p1 #+ - #N.p5[t-1]*(1-s.p1) - #pre.br6[t] <- N.p5[t-1]*s.p1 #+ - #N.p6[t-1]*(1-s.p1) - #pre.br7[t] <- N.p6[t-1]*s.p1 #+ - #N.p7[t-1]*(1-s.p1) - #pre.br8[t] <- N.p7[t-1]*s.p1 #+ - #N.p8[t-1]*(1-s.p1) - #pre.br9[t] <- N.p8[t-1]*s.p1 #+ - #N.p9[t-1]*(1-s.p1) - #pre.br10[t] <- N.p9[t-1]*s.p1 #+ - #N.p10[t-1]*(1-s.p1) - + pre.br1[t] <- f[t-1] * mean.fem * N.tot[t-1] * s.p1 * mu.n.fem + no.breed[t] <- (N.nb[t-1] * phi[t-1] * (1-psiNBB)) + (N.b[t-1] * phi[t-1] * psiBNB) breed[t] <- (N.p7[t-1] * s.p1) + @@ -390,34 +291,18 @@ krst.ipm <- nimbleCode({ (N.b[t-1] * phi[t-1] * (1-psiBNB)) N.p1[t] ~ dpois(pre.br1[t]) - #N.p2[t] ~ dpois(s.p1*N.p1[t-1]) N.p2[t] ~ dbinom(s.p1,N.p1[t-1]) - #N.p3[t] ~ dpois(N.p2[t-1]*s.p1) N.p3[t] ~ dbinom(s.p1,N.p2[t-1]) - #N.p4[t] ~ dpois(N.p3[t-1]*s.p1) N.p4[t] ~ dbinom(s.p1,N.p3[t-1]) - #N.p5[t] ~ dpois(N.p4[t-1]*s.p1) N.p5[t] ~ dbinom(s.p1,N.p4[t-1]) - #N.p6[t] ~ dpois(N.p5[t-1]*s.p1) N.p6[t] ~ dbinom(s.p1,N.p5[t-1]) - #N.p7[t] ~ dpois(N.p6[t-1]*s.p1) N.p7[t] ~ dbinom(s.p1,N.p6[t-1]) - #N.p8[t] ~ dpois(N.p7[t-1]*s.p1) - #N.p8[t] ~ dbinom(s.p1,N.p7[t-1]) - #N.p9[t] ~ dpois(N.p8[t-1]*s.p1) - #N.p9[t] ~ dbinom(s.p1,N.p8[t-1]) - #N.p10[t] ~ dpois(N.p9[t-1]*s.p1) - #N.p10[t] ~ dbinom(s.p1,N.p9[t-1]) N.nb[t] ~ dpois(no.breed[t]) N.b[t] ~ dpois(breed[t]) log(imm[t]) <- eps.imm[t] N.i[t] ~ dpois(imm[t]) } - #for(t in 1:(Time+Time.fore)){ - # Ntot[t] <- N.b[t] + N.nb[t] - #} - #number of females seen as count data for(t in 1:Time){ @@ -443,50 +328,29 @@ f_ini <- mean(tx_nest_dat$hatch_per_nest) #Nest data starts w/ four years prior to 1991 - 2014, minus total in last row #remove last column (2022) from dat, not complete data -#dat = dat[,1:31] dat = dat[2:814,4:31] y=dat+1 -#z.ini=zinits_dcat+1 first = apply(dat, 1, function(x) min(which(x !=0))) fem.dat = read.csv("hatch_sexes.csv") #forecasting Time = dim(dat)[2] -#Time.fore = 45 Time.fore = 79 #to 2100 -#prob.fem = c(fem.dat$prob_fem[-c(1:4)]/100,rep(NA,8)) prob.fem = c(fem.dat$prob_fem[-c(1:5)]/100,rep(NA,6)) -#prob.fem[(Time+Time.fore)] = NA - -#phi.cov = rnorm(30,0,1) -#fec.cov = rnorm(31,0,1) -#phi.ind = c(rep(1,18), rep(2,3),rep(3,9)) -#phi.ind = c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5),rep(6,5)) n.fem = tx_nest_dat$nests[-c(1:7)]/apply(dat,2,sum) -#n.fem[2:4] <-1 n.fem[1] <- 1 -#n.fem[(Time+Time.fore)] <- NA #Nest success data -#E = (tx_nest_dat$eggs_intact[-c(1:4,25:31)] + tx_nest_dat$eggs_broken[-c(1:4,25:31)]) -#J = tx_nest_dat$hatch_release[-c(1:4,25:31)] E = (tx_nest_dat$eggs_intact[-c(1:7,29:35)] + tx_nest_dat$eggs_broken[-c(1:7,29:35)]) J = tx_nest_dat$hatch_release[-c(1:7,29:35)] - - -#j.e = J/E -#j.e[2:3] = NA #beta data can't include 0 -#E.R = E/tx_nest_dat$nests[-c(1:4,25:31)] -#E.R[2:3] <- 1 E.R = E/tx_nest_dat$nests[-c(1:7,29:35)] incu = nest.pro$incubation.type -#incu = incu[-which(incu=="S")] incu[incu == "I" | incu == "I "] = 1 incu[incu == "C"] = 0 incu[incu == "S"] = 2 @@ -537,26 +401,15 @@ E.ind = E.ind[-1472] E.ind[is.na(E.ind)] = round(mean(E.ind,na.rm=TRUE)) year = nest.pro$Year -#year = year - 1990 year = year - 1993 year = year[-1472] year.mean.idx = as.numeric() -#for(t in 25:31){ for(t in 22:28){ year.mean.idx[t] = max(which(year == t)) } -#nests.fore = 250*Time.fore #250 nests each year in future -#n.fore.fems = (225/250)*nests.fore - -#year[(year.mean.idx[31]+1):(year.mean.idx[31]+nests.fore)] = rep(32:(31+Time.fore),each=250) - -#for(t in 32:(31+Time.fore)){ -# year.mean.idx[t] = max(which(year == t)) -#} - -#inport climate data +#import climate data #historic load("~/KRST/Analysis/KRST climate histical/krst.climate.variables.Rdata") #start with bio5 @@ -571,7 +424,6 @@ bio5.scale = as.numeric(scale(c(bio5$mean,bio5.tmp$mean))) bio5$mean.scale = bio5.scale[1:31] #temp during may (most common nesting month) -#setwd("~/KRST/Analysis/KRST climate histical") load("C:/Users/beross/OneDrive - DOI/Documents/KRST/Analysis/SLR Data/slr.krst.data.Rdata") slr.cov = rep(NA,length(nest.pro$Year)-1) @@ -596,34 +448,19 @@ slr.cov = as.numeric(scale(slr.cov)) #bio5.fore.cov = c(rep(bio5.scale[32],20),rep(bio5.scale[33],20),rep(bio5.scale[34],10)) bio5.fore.cov = bio5.tmp$mean[1:Time.fore] slr.cov.fore = slr.tmp$mean[1:Time.fore] -#100 nests in year T+1 -#E.ind[(year.mean.idx[31]+1):(year.mean.idx[31]+nests.fore)] = round(mean(E.ind)) -#J.ind[(year.mean.idx[31]+1):(year.mean.idx[31]+nests.fore)] = NA #forecast nest management scenarios #%50 in corrals, %50 in situ nest.mange.fore = c(63,62,0,0,0,125) -#nest.mange[(year.mean.idx[31]+1):(year.mean.idx[31]+nests.fore)] = rep(mange.fore,Time.fore) - -#forecast CMR data - 225 new individuals for 250 new nests -#y.fore = rbind(y,matrix(NA,n.fore.fems,31)) -#y.fore = cbind(y.fore,matrix(NA,(814+n.fore.fems),Time.fore)) - -#first[815:12064] = rep(32:81,each=225) #forecast count data y.fem = apply(dat,2,sum) -#y.fem[32:(Time+Time.fore)] = NA y.fem[29:(Time+Time.fore)] = NA -#dat.fore = rbind(dat,matrix(0,100,31)) -#dat.fore = cbind(dat.fore,c(rep(0,795),rep(1,100))) - #need to truncate N.b + N.i to be at least as # large as number of observed individuals. Just # have to use the y.fem data twice, once as data # once as truncation constraint -#trun = c(apply(dat,2,sum),NA) trun = c(apply(dat,2,sum)) #### Data for NIMBLE #### @@ -633,7 +470,6 @@ krst_data <- list(y = y, y.fem = y.fem, psiBNB.mom=psiBNB.mom, p.mom=p.mom, s.p.mom=s.p.mom, - #prob.fem.mom = prob.fem.mom, in.s = in.s, J = J, E.R = E.R, @@ -645,11 +481,6 @@ krst_data <- list(y = y, y.fem = y.fem, N.p5.ini = 33, N.p6.ini = 20, N.p7.ini = 13, - #N.p8.ini = 10, - #N.p9.ini = 5, - #N.p10.ini = 3, - #N.b.ini=2, - #N.nb.ini=2, prob.fem = prob.fem, n.fem = n.fem, J.ind = J.ind, @@ -669,13 +500,6 @@ bio5.585.2 = make.fore.dat(hist.bioc.means = hist.bioc.means, sl.summary = sl.summary, slr.scenario = slr.scenario) -#future patrol effort -k.scale = as.numeric(scale(kilo.patrol)) -kilo.patrol.fore = c(k.scale,rep(k.scale[31],Time.fore)) -kilo.patrol[kilo.patrol > 153250] = 1 -kilo.patrol[kilo.patrol > 1] = 0 -kilo.patrol.fore = c(kilo.patrol,rep(1,Time.fore)) - krst_con <- list(N = nrow(y), Time = ncol(dat), first=first, Time.fore = Time.fore, @@ -683,7 +507,6 @@ krst_con <- list(N = nrow(y), Time = ncol(dat), year = year, year.mean.idx = year.mean.idx, N.nest = length(J.ind), - #kilo.patrol = kilo.patrol.fore, may.tmax.ann = bio5.585.2$bio5.ann[4:31], clim.cov = bio5.585.2$bio5.out, clim.cov.fore = bio5.585.2$bio5.fore.cov, @@ -726,15 +549,10 @@ krst_ini = function() list( N.p5 = rpois((Time+Time.fore),seq(30,400,length.out=(Time+Time.fore))), N.p6 = rpois((Time+Time.fore),seq(20,200,length.out=(Time+Time.fore))), N.p7 = rpois((Time+Time.fore),seq(10,100,length.out=(Time+Time.fore))), - #N.p8 = rpois((Time+Time.fore),seq(8,125,length.out=(Time+Time.fore))), - #N.p9 = rpois((Time+Time.fore),seq(5,60,length.out=(Time+Time.fore))), - #N.p10 = rpois((Time+Time.fore),seq(3,45,length.out=(Time+Time.fore))), - #N.i = rpois((Time+Time.fore),c(trun[1:Time]/2,rep(50,Time.fore))), N.i = c(NA,rpois((Time+Time.fore-1),seq(1,100,length.out=(Time+Time.fore-1)))), r = runif(1,.8,10), ln.mean.fem = runif(1,1,3), mu.n.fem= runif(1,1.5,3), - #sd.n.fem = runif(1,.1,1), mu.phi = runif(1,2,2.2), sd.phi2=runif(1,.1,1.5), mu.egg = runif(1,-1,1), @@ -742,19 +560,12 @@ krst_ini = function() list( insitu.surv = runif(1,.1,.8), mu.temp = runif(1,-1,1), var.temp = runif(1,.1,3), - #beta.temp.all = runif(1,-1,1), beta.temp = runif(7,-1,1), - #mu.imm = runif(1,-1,1), - #sd.imm = runif(1,.01,1), ln.mu.phi = runif((Time+Time.fore-1),1,2.2), eps.imm.tmp = runif((Time+Time.fore),-1,1), J.ind = round(J.ini), prob.fem = prob.fem.ini, - #n.fem = c(rep(NA,Time),rep(2,Time.fore)), - #y.fem = c(rep(NA,Time),rep(150,Time.fore)), beta.slr = c(rep(NA,5),runif(1,-1,1)), - #a.0 = runif(1,1,2), - #a.k = runif(1,-1,1), f = rep(round(f_ini),(Time+Time.fore))) #### nimble run #### @@ -809,25 +620,17 @@ current.5852 <- clusterEvalQ(clus, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) out.mcmc <- as.mcmc.list(current.5852) -#stopCluster(clus) Sys.time() out2.current585 <- clusterEvalQ(clus,{ @@ -837,8 +640,6 @@ out2.current585 <- clusterEvalQ(clus,{ out.current.5852 <- as.mcmc.list(out2.current585) -#stopCluster(clus) - save(out.current.5852, file="results_currentnm_5852_2100.RData") #### current NM, best case global change, insitu beta = -1 #### @@ -892,26 +693,15 @@ out.2455 <- clusterEvalQ(clus2455, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - Sys.time() out2455 <- clusterEvalQ(clus2455,{ @@ -974,28 +764,15 @@ out.2455.is0 <- clusterEvalQ(clus2455.insitu0, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - -Sys.time() - out2455.is0 <- clusterEvalQ(clus2455.insitu0,{ CmodelMCMC$run(100000,thin=10,reset=FALSE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) @@ -1044,8 +821,6 @@ for(j in seq_along(clus5852.insitu25)){ clusterExport(clus5852.insitu25[j],"init") } -Sys.time() - out.5852.is25 <- clusterEvalQ(clus5852.insitu25, { library(nimble) library(coda) @@ -1056,28 +831,15 @@ out.5852.is25 <- clusterEvalQ(clus5852.insitu25, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - -Sys.time() - out5852.is25 <- clusterEvalQ(clus5852.insitu25,{ CmodelMCMC$run(100000,thin=10,reset=FALSE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) @@ -1087,7 +849,6 @@ out.current.5852.is25 <- as.mcmc.list(out5852.is25) save(out.current.5852.is25,file = "results_5852_is25_2100.RData") - #### all in situ NM, best case global change, insitu beta = -1 #### slr.scenario = "0.5 - MED" @@ -1139,28 +900,15 @@ out.2455.is1 <- clusterEvalQ(clus2455.insitu1, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) - CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) + CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - -Sys.time() - out2455.is1 <- clusterEvalQ(clus2455.insitu1,{ CmodelMCMC$run(100000,thin=10,reset=FALSE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) @@ -1226,28 +974,15 @@ out.2455.noincu <- clusterEvalQ(clus2455.noincu, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) - CmodelMCMC$run(15000,reset=TRUE, resetMV=TRUE) + CmodelMCMC$run(15000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - -Sys.time() - out2455.noincu <- clusterEvalQ(clus2455.noincu,{ CmodelMCMC$run(100000,thin=10,reset=FALSE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) @@ -1308,28 +1043,15 @@ out.2455.lowis <- clusterEvalQ(clus2455.lowis, { inits = init) Cmodel <- compileNimble(model) modelConf <- configureMCMC(model) - #modelConf$addMonitors(params) modelConf$addMonitors(c("phi","N.i","N.b","pB","N.nb", "f","eps.imm")) - #configureRJ(modelConf, - # targetNodes = 'beta', - # indicatorNodes = 'indA', - # control = list(mean=0, scale = 2)) modelMCMC <- buildMCMC(modelConf) CmodelMCMC <- compileNimble(modelMCMC, project = model) - #out1 <- runMCMC(CmodelMCMC, niter=15000, nburnin=14000) - #return(as.mcmc(out1)) - CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) + CmodelMCMC$run(100000,reset=TRUE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples))) }) -#out.mcmc <- as.mcmc.list(out) - -#stopCluster(clus) - -Sys.time() - out2455.lowis <- clusterEvalQ(clus2455.lowis,{ CmodelMCMC$run(100000,thin=10,reset=FALSE, resetMV=TRUE) return(as.mcmc(as.matrix(CmodelMCMC$mvSamples)))