We first load the data and select the species of interest. Here we choose BTBW which is not one of the rare species present in the dataset.
Afterwards we prepare the data in the format that is required by the R-INLA library
Code
data(hbefTrends)revi.data<-hbefTrendssp.names<-dimnames(hbefTrends$y)[[1]]revi.data$y<-revi.data$y[sp.names=='REVI', , , ]revi.data$coords=revi.data$coords/1000#get this in km!## data preparation for inla -----------------------------------------------Y=data.frame(revi.data$y[,1,])Xdet.day=revi.data$det.covs$day[,1,]Xdet.tod=revi.data$det.covs$tod[,1,]for(iin2:9){Y=rbind(Y, data.frame(revi.data$y[,i,]))Xdet.day=rbind(Xdet.day, revi.data$det.covs$day[,i,])Xdet.tod=rbind(Xdet.tod, revi.data$det.covs$tod[,i,])}Xdet.day=(Xdet.day-mean(Xdet.day, na.rm =T))/sd(Xdet.day,na.rm =T)Xdet.tod=(Xdet.tod-mean(Xdet.tod, na.rm =T))/sd(Xdet.tod,na.rm =T)Xdet=cbind(1, Xdet.day[,1], Xdet.tod[,1],1, Xdet.day[,2], Xdet.tod[,2],1, Xdet.day[,3], Xdet.tod[,3])Xocc=data.frame(x =rep(revi.data$coords[,1],9), y =rep(revi.data$coords[,2],9), elev =rep(revi.data$occ.covs$elev,9), scale_elev =scale(rep(revi.data$occ.covs$elev,9)), scale_elev2 =scale(rep(revi.data$occ.covs$elev,9))^2, site =rep(1:373,9), time =rep(1:9,each =373), Int_occ =1, spatialfield =NA)%>%mutate(scale_time =scale(time))
Code
elev_raster=rast(data.frame(x =hbefElev$Easting/1000, y =hbefElev$Northing/1000, elev =hbefElev$val))elev_raster$scale_elev=(elev_raster$elev-mean(revi.data$occ.covs$elev))/sd(revi.data$occ.covs$elev)elev_raster$scale_elev2=elev_raster$scale_elev^2all_elev=inla.group(c(Xocc$elev, values(elev_raster$elev)))val_elev=sort(na.omit(unique(all_elev)))Xocc$group_elev=all_elev[1:dim(Xocc)[1]]elev_group_raster=elev_raster$elevvalues(elev_group_raster)<-all_elev[-c(1:dim(Xocc)[1])]pred_df=data.frame(x =rep(crds(elev_raster, na.rm =F)[,1],9), y =rep(crds(elev_raster, na.rm =F)[,2],9), scale_elev =rep(values(elev_raster$scale_elev),9), scale_elev2 =rep(values(elev_raster$scale_elev2),9), group_elev =rep(values(elev_group_raster),9), time =rep(c(1:9), each=length(values(elev_raster$elev))))%>%dplyr::filter(!is.na(scale_elev))
3 Model fit
3.1 Model 1
This is a separable space time model with linear predictors \[
\eta_{st} = \beta_0 + f_1(\text{elev}) + f_2(t) + f_3(s)
\] where \(f_1(\text{elev})\) is a smooth (RW2) effect of the elevation \(f_2(t)\) is a AR1 effect of time \(f_3(s)\) is an IID effect of location
Warning in normalizePath(path.expand(path), winslash, mustWork):
path[1]="C:/Users/jaf_i/AppData/Local/Temp/RtmpecWRYM/file12d881fcb1e49": The
system cannot find the file specified
3.2 Model 2
Code
# set up spatial model for INLAboundary=inla.nonconvex.hull(points =revi.data$coords, convex =.3)mesh=inla.mesh.2d(boundary =boundary,# loc = cbind(data$X, data$Y), max.edge =c(0.1,0.7), min.angle =20, offset =c(.01, 1), cutoff =0.12,)ggplot()+inlabru::gg(mesh)+coord_equal()+theme_maps
Code
spde<-inla.spde2.pcmatern( mesh =mesh, prior.range =c(5, 0.01), # prior for range prior.sigma =c(1, 0.5))# prior for sd parameter
This is a space-time model with linear predictor as \[
\eta_{st} = \beta_0 + f_1(\text{elev}) + f_2(t) + \omega(s)
\] with \(f_1()\) and \(f_2()\) as before while \(\omega(s)\) is a Gaussian spatial field
Warning in normalizePath(path.expand(path), winslash, mustWork):
path[1]="C:/Users/jaf_i/AppData/Local/Temp/RtmpecWRYM/file12d883d1b74f7": The
system cannot find the file specified
3.3 Model 3
The last model is defined as: \[
\eta_{st} = \beta_0 + f_1(\text{elev}) + \omega(s,t)
\] with \(f_1()\) is as before while \(\omega(s,t)\) is a space-time gaussian spatial field with AR1 time component
Code
h.spec<-list(rho =list(prior ='pc.cor0', param =c(0.5, 0.1)))spde<-inla.spde2.pcmatern( mesh =mesh, prior.range =c(5, 0.7), prior.sigma =c(1, 0.5), constr =T)A_sp2<-inla.spde.make.A(mesh =mesh, group =Xocc$time, loc =cbind(Xocc$x, Xocc$y))formula3_3<-inla.mdata(Y,X)~-1+Int_occ+scale_elev+scale_elev2+f(time, model ="ar1")+f(spatialfield, model=spde, A.local =A_sp2, group =time, control.group =list(model ='ar1',hyper =h.spec))time3=system.time(model3_3<-inla(formula3_3, #the formula data=data_list, #the data stack family='occupancy', #which family the data comes from control.fixed =list(prec =1, prec.intercept =1), control.compute =list(waic =TRUE,config =TRUE), #model diagnostics and config = TRUE gives you the GMRF verbose =F, control.family =list(control.link =list(model ="logit"), link.simple ="logit", hyper =list(beta1 =list(param =c(0,10), initial =0), beta2 =list(param =c(0,10)), beta3 =list(param =c(0,10))))))
df_sf<-Xocc%>%mutate(id =1:nrow(Xocc))%>%st_as_sf(coords =c("x","y"))# create buffer of size 300 (based on estimated range) centred at each sitebuffer_25<-st_buffer(df_sf, dist =0.75)# empty lists to include the indexes of the leave-out-group for each observation iI_i<-list()# loop though each observation and store the leave-out-group based on the bufferfor(iin1:nrow(df_sf)){# Temporal filtering of data within a 2 years of span of observation idf_sf_subset<-df_sf%>%filter(between(time,left =df_sf$time[i]-2, right =df_sf$time[i]+2))# Spatial filtering of the observations that are within the buffer of the ith observationBuffer_i<-df_sf_subset%>%st_intersects(buffer_25[i,],sparse =FALSE)%>%# identify unlist()# obtain the indexes of the leave out groupI_i[[i]]<-df_sf_subset[Buffer_i,]%>%pull(id)}
Doser, Jeffrey W., Andrew O. Finley, Marc Kéry, and Elise F. Zipkin. 2022. “spOccupancy: An R Package for Single-Species, Multi-Species, and Integrated Spatial Occupancy Models.”Methods in Ecology and Evolution 13 (8): 1670–78. https://doi.org/10.1111/2041-210x.13897.
