Array for multi-season occupancy model

21 views
Skip to first unread message

Joaquín Aldabe

unread,
Oct 22, 2025, 9:49:44 AMOct 22
to spOccupancy and spAbundance users
This is Joaquín from Uruguay and I´m trying to fit a multi-season detection-non detection model using spOccupancy to a dataset with a set of J sites are sampled across a set of 3 seasons/years. Within each season t, each site j is sampled K j,t times. As far as I understand, I should use  stPGOcc function and a 3D Array with Site, Year, Visit.  Is it correct?

I´m having trouble with the 3d array creation. I appreciate it if you could check my dataset (attached csv file) and suggest a code to make my array. 

Thanks in advance. 

Joaquín.    

Marc Kéry

unread,
Oct 22, 2025, 10:29:04 AMOct 22
to Joaquín Aldabe, spOccupancy and spAbundance users
Hola Joaquin,

If you send the data set and I can have a try. 

In general, formatting detection/nondetection data for more complex occupancy models can be a major stumbling block for a user (I know this very well from when I started with such models). In the AHM2 book, we have chapter section 5.4. which is entitled "DATA FORMATTING FOR THE DCM: FUN WITH MULTIDIMENSIONAL ARRAYS"; this may be helpful. If you don't have the book then you can download the code from here.

Best regards  --- Marc

--
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/7967f59d-7dd0-4ee4-8d56-a6f78167478fn%40googlegroups.com.


--
______________________________________________________________
 
Marc Kéry
Tel. ++41 41 462 97 93
marc...@vogelwarte.ch
www.vogelwarte.ch
 
Swiss Ornithological Institute | Seerose 1 | CH-6204 Sempach | Switzerland
______________________________________________________________

*** Hierarchical modeling in ecology ***

José Wagner Ribeiro Júnior

unread,
Oct 22, 2025, 10:30:54 AMOct 22
to Joaquín Aldabe, spOccupancy and spAbundance users
Hi Joaquín, 

For multi-season occupancy models in the spOccupancy package, the detection/non-detection data y should actually be a four-dimensional array, with dimensions corresponding to species, sites, primary time periods (e.g., years or seasons), and replicates (visits). You can find more details in vignette("spaceTimeModels") from the function.  

Also, I wasn’t able to find the CSV file you mentioned. Could you please check if it was correctly attached?  

Best regards,
José

--
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/7967f59d-7dd0-4ee4-8d56-a6f78167478fn%40googlegroups.com.


--
José Wagner Ribeiro Jr (Xuleta)PhD         
Quantitative Ecologist

Joaquín Aldabe

unread,
Oct 22, 2025, 10:40:57 AMOct 22
to spOccupancy and spAbundance users
Thanks Marc and José! I´m attaching the dataset. Please let me know if you can get it. 

My model is for one species, so i guess the array is 3D instead of 4D?

All my best,
Joaquín. 

dataset.csv

Marc Kéry

unread,
Oct 22, 2025, 12:10:53 PMOct 22
to Joaquín Aldabe, spOccupancy and spAbundance users
Dear Joaquin,

this seems to work and produces a Site x Year x Visit array for your data. I note that your data are extremely sparse and would warn that you may perhaps not have much luck fitting such a complicated model.

Best regards  ---- Marc

Formatting data into a 3d array.R

Joaquín Aldabe

unread,
Oct 22, 2025, 1:57:40 PMOct 22
to Marc Kéry, spOccupancy and spAbundance users
Thank you so much Marc! I will fit the model tomorrow and see if it works. 

All my best,
Joaquín. 
--
Joaquín Aldabe, PhD. 
MSc. Ecología y Evolución
Lic. Ciencias Biológicas

Departamento de Sistemas Agrarios y Paisajes Culturales
Centro Universitario de la Región Este, Universidad de la República, Uruguay
Ruta 15 (y Ruta 9), Km 28.500, Departamento de Rocha  

Joaquín Aldabe

unread,
Oct 29, 2025, 11:44:03 AM (8 days ago) Oct 29
to Marc Kéry, spOccupancy and spAbundance users
Hi all, I´m trying to fit a multi-season single species spatial occupancy model. My dataset is attached, and below is the code I tried. Marc kindly created the array.

I have 90 sites sampled for three different seasons (Years). In each year, I revisited each site. Sites over the years are the same. When I run the model, an error says that I need the same nrow of y as for coordinates: Error in stPGOcc(occ.formula = ~vor + esp_1m + perc.fact, det.formula = ~1, :
data$grid.index must be specified if nrow(data$coords) != nrow(data$y)

This is because in my dataset each row corresponds to a site and year, and each site was sampled three times (one per season). How can I fix this?

The other issue is that my coordinates are in Grades, Decimals (e.g. -53.83907 -33.32663). Does this work for soOcuppancy or I need UTM? If UTM is the only option, I have another issue. My sites cover two different UTM zones (21s and 22s). In the case I have to use UTM format, do you know how I can manage this issue?

Thanks a lot for your helps.

Joaquín.  


# Read in data set
dat <- read.table("dataset.csv", header = TRUE, sep = ',')
str(dat)

# Renumber sites and years consecutively and get dimensions
dat$Site                     # Look at data
site <- as.numeric(as.factor(dat$Site))
nsite <- max(site)
dat$Year                     # Look at data (2013, 2016, 2018)
year <- as.numeric(as.factor(dat$Year))
nyear <- max(year)
nvisit <- 5

# Create empty array of right dimensionality
y3d <- array(NA, dim = c(nsite, nyear, nvisit))
dimnames(y3d) <- list(unique(dat$Site), unique(dat$Year), paste('visit', 1:5, sep = ''))

# Fill it
for(i in 1:nrow(dat)){
  y3d[site[i], year[i], 1:5] <- as.numeric(dat[i, c('Visit.1', 'Visit.2', 'Visit.3', 'Visit.4', 'Visit.5')])
}
y3d                          # Look at array

# Quick sum test
sum(y3d, na.rm = TRUE)  ;  sum(dat[, c('Visit.1', 'Visit.2', 'Visit.3', 'Visit.4', 'Visit.5')], na.rm = TRUE)

#Build spOccupancy objects with appropiate names

y<-y3d
str(y)

det.covs <- list()

# Covariables de ocupación
occ.covs <- dat %>%
  dplyr::select(Vor, Esp_1m, perc.fact)

# Coordenadas (para modelos espaciales)
coords <- dat%>%
  dplyr::select(x, y) %>%
  as.matrix()
str(coords)

# Package all data into list object
data.list <- list(y = y,
                  occ.covs = occ.covs,
                  det.covs = det.covs,
                  coords = coords)

#Multiseason-occupancy spatial model

out <- stPGOcc(occ.formula = ~ vor + esp_1m + perc.fact,
               det.formula = ~ 1,
               data = data.list,
               NNGP = TRUE,
               n.burn = 10000, n.thin = 3, n.chains = 3)
dataset.csv

Susana Requena

unread,
Nov 1, 2025, 1:38:49 PM (5 days ago) Nov 1
to Joaquín Aldabe, Marc Kéry, spOccupancy and spAbundance users
Hi Joaquín.
Yes, your coordinates must be projected, in meters. Not lon lat. And obviously, every point has the same coordinates over the years; the coordinates must have the same length as the points. This looks to me like the source of the error (I've been there before :).

In my experience, working in UTM can be very cumbersome, especially if your data is spread across different grids.

My approach in these cases is to project the coordinates from degrees into a more 'general' projection. Which to choose depends on your latitude and the extent of the area.  An equidistant Conic projection works well for median latitudes; another, more "local" one will work just as well. But you should centre the projection on the sample to minimise distortion, setting standard parallels 1 and 2 to the northmost and southmost coordinates of your study area. Then, if you need it, you could reproject your maps and results to UTM. This works for me.

See my code below (I sacrificed coding style for clarity). I tested it on your dataset, and it worked; it should be fine. Happy to have suggestions, 

Cheers!
Susana

P.D. 
(tidyverse)

dataset <- read_csv("C:/Users/SusanaRequena/Downloads/dataset.csv")
str(dataset)

dataset$lon<- as.numeric((dataset$x), rm.NA=TRUE)# change names to avoid confusions later
dataset$lat<- as.numeric((dataset$y), rm.NA=TRUE)

range(dataset$lon)
range(dataset$lat) # check there's nothing extrange

#Projecting and packing coordinates-----

library(terra)
library(sp)

centralpoint<-data.frame(cbind(mean(dataset$lon), mean(dataset$lat))) # points centroid (could be median to avoid effect of extremes)
centralpoint

# Northmost and southmost latitudes from the dataset range
lat_north <- -33.03455 # Northmost latitude
lat_south <- -33.90712 # Southmost latitude

# Convert original points to a SpatialPoints object in WGS84 (assuming it's your datum)
coords.lonlat <- SpatialPoints(cbind(dataset$lon, dataset$lat), proj4string = CRS("+proj=longlat +datum=WGS84") )

# Create your own plocal projection
DgProj <- CRS(paste("+proj=eqdc +lon_0=", centralpoint$X1, " +lat_0=", centralpoint$X2, " +lat_1=", lat_south, " +lat_2=", lat_north, sep="")) # Centrering the general projection in your study area

# Project them into your local  system, derived from the Equidistant Conic in this case
coords_proj <- spTransform(coords.lonlat, DgProj)

# Add projected coordinates back into the data
dataset$X <- coordinates(coords_proj)[, 1]
dataset$Y <- coordinates(coords_proj)[, 2]

# --- Coordinates reduced to unique X,Y by point---
coords <- dataset %>%
  group_by(Site) %>%
  summarize(
    X = mean(X), # getting rid of small variations if the GPS coords were registered in every replica
    Y = mean(Y),
    .groups = "drop") %>%
  arrange(Site)%>%
  as.data.frame() # convert tibble to data.frame

str(coords) #same lenght as points (ie. 90)

# Convert to a numeric matrix with X, Y columns

rownames(coords) <- coords$Site ## Set row names to point_id

# Now convert to matrix
coords.matrix <- as.matrix(coords[, c("X", "Y")])

# Check structure
str(coords.matrix)

class(coords.matrix)

dim(coords.matrix) # as many as points (90): 2 (X,Y)

coords.matrix ["1" ,1 ] #correct, test some of them

#6.4 Package all data into list object----

data.model <- list(y = y.dat,
                  occ.covs = occ.covs,
                  det.covs = det.covs,
                  coords = coords.matrix) # directly


Reply all
Reply to author
Forward
0 new messages