Plotting spPGOcc predictions over spatial extent

18 views
Skip to first unread message

Kristin Engebretsen

unread,
Mar 24, 2026, 3:32:07 PM (7 days ago) Mar 24
to spOccupancy and spAbundance users
Hi all,

I am trying to use PGOcc and spPGOcc to model marten occupancy over a section of National Forest land in Idaho. I have successfully fitted short simple runs of both models with a small number of occ.covs and det.covs to practice using the functions, and I can back transform to look at the occupancy estimates and detection estimates for my covariates. However, I am having trouble getting a prediction surface for occupancy using the covariates in my model. 

I'm familiar with using terra::predict(raster, fitted.mod) to predict out across a raster surface. However, it seems that the model output from (sp)PGOcc require that you specify the coordinates to predict over as a df or matrix. I've gone back to try to get this to work on the non-spatial model first, shown here.

I can successfully generate an occupancy prediction by just giving the covariate values at the camera locations, which I'll attached the output figure for as well. 

The problem seems to be the memory size generating the prediction at all coordinates of my intended raster. I'm interested if there is a workaround for this? Is there a way to predict the fitted model with a raster of interest as one can do with other fitted model objects?
 
Thanks!
Kristin


# package data. occ.covs contains several covariates extracted at camera locations and then scaled, including DEM.sc, which I am going to use in the model. 

marten_data <- list( y=station_det_wk_bin,
                     occ.covs = occ.covs,
                     det.covs = list(effort=station_det_wk_eff,
                                     Lure = station_metadata.df$Lure,
                                     ScentFamily = station_metadata.df$ScentFamily,
                                     Bait.Type = station_metadata.df$Bait.Type),
                     coords = station.coords3)


n.samples <- 3000
n.burn <- 1000
n.thin <- 2
n.chains <- 1

mod1 <- PGOcc(occ.formula = ~DEM.sc,  
             det.formula = ~Lure,
             data = marten_data,
             # inits = oven.inits,
             n.samples = n.samples,
             # priors = oven.priors,
             n.omp.threads = 1,
             verbose = TRUE,
             n.report = 1000,
             n.burn = n.burn,
             n.thin = n.thin,
             n.chains = n.chains)

summary(mod1)
summary(mod1) Call: PGOcc(occ.formula = ~DEM.sc, det.formula = ~Lure, data = marten_data, n.samples = n.samples, n.omp.threads = 1, verbose = TRUE, n.report = 1000, n.burn = n.burn, n.thin = n.thin, n.chains = n.chains) Samples per Chain: 3000 Burn-in: 1000 Thinning Rate: 2 Number of Chains: 1 Total Posterior Samples: 1000 Run Time (min): 0.0392 Occurrence (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) 0.7201 0.2149 0.3134 0.7187 1.1447 NA 661 DEM.sc 0.3549 0.2092 -0.0469 0.3542 0.7530 NA 861 Detection (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) -1.3776 0.1048 -1.5932 -1.3725 -1.1739 NA 1000 LureNone -0.0746 0.2981 -0.6727 -0.0535 0.4797 NA 604 LureSingleDose -0.0139 0.1409 -0.2761 -0.0233 0.2604 NA 887


# predict and plot data at camera locations
pred1.MOD2 <- predict(mod1, X.0=cbind(1,marten_data$occ.covs$DEM.sc))

plot_df.2 <- data.frame( x = marten_data$coords$long,
                       y = marten_data$coords$lat,
                       mean.psi = apply(pred1.MOD2$psi.0.samples,2,mean),
                       sd.psi = apply(pred1.MOD2$psi.0.samples, 2, sd))

dat.stars <- st_as_stars(plot_df.2, dims = c('x', 'y'))
head(dat.stars)

ggplot() +
  geom_stars(data=dat.stars, aes(x, y, fill = mean.psi)) +
  scale_fill_viridis_c(na.value = "transparent") +
  labs(x = "Longitude", y = "Latitude", fill = "Predicted\nOccupancy") 

* see attached fig*


# try to plot over area of interest; modify the original DEM layer into a df to use predict() and then scale it
coords.df <- crds(DEM, df=TRUE)
DEM.df <- as.data.frame(DEM)
elev.as.df <- data.frame(coords.df,DEM.df)

elev.pred <- (elev.as.df$Elev - mean(station.coords2$DEM$Elev) / sd(station.coords2$DEM$Elev))

X.0 <- cbind(1, elev.pred)
head(X.0)
head(X.0) elev.pred [1,] 1 1300.863 [2,] 1 1306.006 [3,] 1 1309.755 [4,] 1 1315.756 [5,] 1 1322.356 [6,] 1 1328.449

pred.MOD <- predict(mod1, X.0)
Error: cannot allocate vector of size 194.1 Gb




Kristin Engebretsen

unread,
Mar 24, 2026, 3:33:44 PM (7 days ago) Mar 24
to spOccupancy and spAbundance users
Would help if I attached the figure too :D 
predict_PGOcc.png

Jeffrey Doser

unread,
Mar 27, 2026, 3:53:21 AM (5 days ago) Mar 27
to Kristin Engebretsen, spOccupancy and spAbundance users
Hi Kristin, 

Thanks for the message, and hope all is well. There is no functionality in spOccupancy that allows you to generate predictions directly from a raster object, so you need to convert the raster to a matrix/df of coordinates as you are doing. In order to get around memory issues, you will likely need to split up the coordinates into smaller chunks, predict individually at the chunks, and then combine them all together at the end. There are other ways you can reduce the memory required for the prediction objects as well. See this thread on the listserv for some more discussion on this, as well as an example script that shows how to go about splitting the coordinates up into groups to aid in prediction. Let me know if that thread doesn't solve your problems. 

Kind regards,

Jeff

--
You received this message because you are subscribed to the Google Groups "spOccupancy and spAbundance users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to spocc-spabund-u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/spocc-spabund-users/95ef4ad8-128b-4c91-8186-ee9f9d7a1c3bn%40googlegroups.com.


--
Jeffrey W. Doser, Ph.D.
Assistant Professor
Department of Forestry and Environmental Resources
North Carolina State University
Pronouns: he/him/his
Reply all
Reply to author
Forward
0 new messages