1 Introduction
One of the key features ofubms is the ability toinclude random effects in model formulas. This is not possible inunmarked, but it is possible with custom models usingJAGS. To explore howestimates of parameters and uncertainty intervals compared betweenubms andJAGS, I simulated a basic occupancydataset with known parameter values and fit models with both. In thissimulation, individual sites were nested within groups, and I estimatedgroup-level random intercepts. This is a common scenario for so-calledstackeddata from multiple years, where the “groups” are unique sites and“sites” represent unique site-year combinations.
2 Simulating thedata
First I set the true parameter values for the vector of fixed effects\(\boldsymbol{\beta}\) and the randomeffects standard deviation\(\sigma\).
beta<-c(0.4,-0.5,0.3,0.5)sigma<-1.2
Then, I simulated covariate data for the occupancy and detectionsubmodels. This included a random group assignment (26 possible groups)for each site.
dat_occ<-data.frame(x1=rnorm(500),group=factor(sample(letters[1:26],500,replace=T)))dat_p<-data.frame(x2=rnorm(500*5))
Next I simulated the random effect value for each group (with aneffects parameterization centered on 0).
Finally, I simulated an occupancy datasety using thedata and parameter values obtained above.
y<-matrix(NA,500,5)z<-rep(NA,500)group_idx<-as.numeric(dat_occ$group)#convert group assignment to numeric indexidx<-1for (iin1:500){#True latent state of site i in group group_idx[i] z[i]<-rbinom(1,1,plogis(beta[1]+ beta[2]*dat_occ$x1[i]+ re[group_idx[i]]))#Observation processfor (jin1:5){#Observed data at site i for observation j y[i,j]<- z[i]*rbinom(1,1,plogis(beta[3]+ beta[4]*dat_p$x2[idx])) idx<- idx+1 }}
3 Fitting the model inubms
I combined the covariates and simulated data into anunmarkedFrame.
library(ubms)umf<-unmarkedFrameOccu(y=y,siteCovs=dat_occ,obsCovs=dat_p)
Usingstan_occu, I fit a model with a random effect ofgroup on occupancy ((1|group)). I probably did not useenough iterations but it’s fine for this example.
ubms_fit<-stan_occu(~x2~x1+(1|group), umf,chains=3,iter=300,cores=3,seed=123)
4 Fitting the model inJAGS
A bit more work is required to fit this model in JAGS. First Ireorganized the data into a list so it could be sent to JAGS.
nsites<-nrow(umf@y)jags_data<-list(y=umf@y,x1=dat_occ$x1,group=as.numeric(dat_occ$group),ngroups=nlevels(dat_occ$group),x2=matrix(dat_p$x2,nrow=nsites,byrow=TRUE),nsites=nsites,nobs=ncol(umf@y))
Next I wrote a simple occupancy model with group-level interceptsgint in the BUGS language, and saved it to a temporaryfile.
modfile<-tempfile()writeLines("model{#Likelihoodfor (i in 1:ngroups){ gint[i] ~ dnorm(beta[1], tau.group)}for (i in 1:nsites){ z[i] ~ dbern(psi[i]) logit(psi[i]) <- gint[group[i]] + beta[2]*x1[i] for (j in 1:nobs){ y[i,j] ~ dbern(p[i,j]*z[i]) logit(p[i,j]) <- beta[3] + beta[4]*x2[i,j] }}#Priorsfor (i in 1:4){ beta[i] ~ dnorm(0,0.01)}sd.group ~ dunif(0,10)tau.group <- pow(sd.group, -2)}",con=modfile)
JAGS also requires a list of parameters you want to save, and afunction to generate initial values. Generally you need to providereasonable initial values for the latent state\(z\). Note that I’m saving both theparameter values ("beta" and"sd.group"),along with the actual random effects estimates("gint").
params<-c("beta","sd.group","gint")inits<-function(){list(z=apply(y,1, max,na.rm=TRUE))}
Finally, I fit the model using the R packagejagsUI.
set.seed(123)library(jagsUI)jags_fit<-jags(jags_data, inits, params, modfile,n.chains=3,n.adapt=100,n.iter=2000,n.burnin=1000,n.thin=2,parallel=TRUE,verbose=FALSE)
Running the model in JAGS is significantly slower, although you couldspeed up by marginalizing the likelihood.
5 Compare fixed effectestimates
With the model fit in bothubms andJAGS, Icompared the fixed effect parameter point estimates and uncertaintyintervals visually.
#Data frame of JAGS parameter estimates and UIsjags_sum<- jags_fit$summarybeta_jags<-as.data.frame(jags_sum[1:5,c(1,3,7)])names(beta_jags)<-c("mean","lower","upper")beta_jags$source<-"JAGS"#Data frame of ubms parameter estimates and UIsubms_sum<-summary(ubms_fit,"state")beta_ubms<-rbind(summary(ubms_fit,"state")[,c(1,4,8)],summary(ubms_fit,"det")[,c(1,4,8)])beta_ubms<- beta_ubms[c(1:2,4:5,3),]names(beta_ubms)<-c("mean","lower","upper")beta_ubms$source<-"ubms"#Data frame of true parameter valuesbeta_truth<-data.frame(mean=c(beta, sigma),lower=NA,upper=NA,source="Truth")#Bind togetherbeta_plot<-rbind(beta_jags, beta_ubms, beta_truth)beta_plot$Parameter<-rep(rownames(beta_jags),3)#Plotlibrary(ggplot2)dodge<-position_dodge(0.4)ggplot(data=beta_plot,aes(x=Parameter,y=mean,col=source))+geom_errorbar(aes(ymin=lower,ymax=upper),position=dodge)+geom_point(position=dodge,size=3)+labs(x='Parameter',y='Value',col='Source')+theme_bw()+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.text=element_text(size=12),axis.title=element_text(size=14),legend.text=element_text(size=12),legend.title=element_text(size=14))
Fixed-effects parameter estimates are very similar betweenubms andJAGS, as are the uncertaintyintervals. In both cases the intervals overlap the true parametervalues.
6 Compare random effectestimates
Finally, I did the same for the random effect estimates. Theranef method extracts random effects terms from a fittedubms model. Although we specified the random effect inubms as an “effects” parameterization,ranefautomatically adds the random effect to the intercept term to get thecomplete random intercept at the group level. The group-level randomintercepts in JAGS are parameter"gint".
#Get random effect estimates from ubmsre_ubms<-ranef(ubms_fit,"state",summary=TRUE)[[1]][[1]][,c(1,3,4)]re_ubms$source<-"ubms"#Get random effects estimates from JAGSre_jags<-as.data.frame(jags_sum[grepl("gint",rownames(jags_sum)),c("mean","2.5%","97.5%")])re_jags$source<-"JAGS"names(re_jags)<-names(re_ubms)#Get truthtruth<-data.frame(Estimate=re,`2.5%`=NA,`97.5%`=NA,source="Truth",check.names=FALSE)#Combineplot_dat<-rbind(re_ubms, re_jags, truth)plot_dat$group<-rep(letters[1:26],3)levels(plot_dat$source)<-c("Truth","JAGS","ubms")library(ggplot2)dodge<-position_dodge(0.4)#Plot a subset of random effectsplot_dat_sub<- plot_dat[plot_dat$group%in% letters[1:6],]ggplot(data=plot_dat_sub,aes(x=group,y=Estimate,col=source))+geom_errorbar(aes(ymin=`2.5%`,ymax=`97.5%`),position=dodge)+geom_point(position=dodge,size=3)+labs(x='Group',y='Random effect value',col='Source')+theme_bw()+theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(),axis.text=element_text(size=12),axis.title=element_text(size=14),legend.text=element_text(size=12),legend.title=element_text(size=14))
As with the fixed effects, the estimated random intercepts are verysimilar between the two methods. Success!