camtrapR (version 3.0.0)

PPC.community: Calculate Community-Level Posterior Predictive Checks for Occupancy Models

Description

A wrapper function that applies Posterior Predictive Checks (PPC) across multiple species in a community occupancy model (from communityModel. It calculates species-specific and community-level Bayesian p-values to assess model fit. The function accepts both predict output format and list format for model parameters.

Usage

PPC.community(
  p,
  psi,
  y,
  input_format = c("predict", "lists"),
  K = NULL,
  model = c("Occupancy", "RN"),
  zhat = NULL,
  z.cond = TRUE,
  type = c("FT", "PearChi2", "Deviance"),
  return.residuals = TRUE,
  show_progress = TRUE,
  ...
)

Value

If return.residuals=TRUE, returns a list containing:

• bp_values - Data frame with species-specific and community-level Bayesian p-values

• species_results - List containing the complete PPC.residuals output for each species

If return.residuals=FALSE, returns only the data frame of Bayesian p-values.

Arguments

p

Either a 4D array [stations, species, iterations, occasions] from predict or a list of 3D arrays [iterations, stations, occasions], one per species

psi

Either a 3D array [stations, species, iterations] from predict or a list of 2D arrays [iterations, stations], one per species

y

List of detection histories, one matrix/vector per species

input_format

Character indicating format of p and psi ("predict" or "lists")

K

Number of occasions as either a scalar or site vector. Calculated automatically if y is a list of matrices.

model

Character indicating model type ("Occupancy" or "RN")

zhat

List of z estimate matrices, one per species (optional). Each matrix should follow the format specified in PPC.residuals.

z.cond

Logical. If TRUE, new data is conditioned on estimated z (testing only detection model fit). If FALSE, generates new z for each posterior sample (testing complete model).

type

Character indicating residual type ("FT", "PearChi2", or "Deviance")

return.residuals

Logical. If TRUE, returns species-specific residuals along with Bayesian p-values. If FALSE, returns only the p-values.

show_progress

Logical; whether to show a progress bar during computation (if package pbapply is available)

...

Additional arguments passed to PPC.residuals.

Author

Rahel Sollmann

Details

This function extends the single-species PPC analysis to the community level by:

• Applying residual calculations to each species in the community

• Aggregating results to assess community-level model fit

• Providing both species-specific and community-level diagnostics

This function provides flexibility in input formats to accommodate different workflows:

• Direct output from camtrapR's predict() function (4D/3D arrays)

• Lists of species-specific arrays (for custom workflows)

See Also

PPC.residuals for details on the underlying single-species calculations

Examples

if (FALSE) {

# Create and fit model 
data("camtraps")

# create camera operation matrix
camop_no_problem <- cameraOperation(CTtable      = camtraps,
                                    stationCol   = "Station",
                                    setupCol     = "Setup_date",
                                    retrievalCol = "Retrieval_date",
                                    hasProblems  = FALSE,
                                    dateFormat   = "dmy"
)

data("recordTableSample")

# make list of detection histories
species_to_include <- unique(recordTableSample$Species)

  DetHist_list <- detectionHistory(
    recordTable         = recordTableSample,
    camOp                = camop_no_problem,
    stationCol           = "Station",
    speciesCol           = "Species",
    recordDateTimeCol    = "DateTimeOriginal",
    species              = species_to_include,
    occasionLength       = 7,
    day1                 = "station",
    datesAsOccasionNames = FALSE,
    includeEffort        = TRUE,
    scaleEffort          = TRUE,
    timeZone             = "Asia/Kuala_Lumpur"
)


# create some fake covariates for demonstration
sitecovs <- camtraps[, c(1:3)]
sitecovs$elevation <- c(300, 500, 600)   

# scale numeric covariates
sitecovs[, c(2:4)] <- scale(sitecovs[,-1])


# bundle input data for communityModel
data_list <- list(ylist = DetHist_list$detection_history,
                  siteCovs = sitecovs,
                  obsCovs = list(effort = DetHist_list$effort))


# create community model for JAGS
modelfile1 <- tempfile(fileext = ".txt")
mod.jags <- communityModel(data_list,
                           occuCovs = list(fixed = "utm_y", ranef = "elevation"),
                           detCovsObservation = list(fixed = "effort"),
                           intercepts = list(det = "ranef", occu = "ranef"),
                           modelFile = modelfile1)

summary(mod.jags)

# fit in JAGS
fit.jags <- fit(mod.jags,
                n.iter = 1000,
                n.burnin = 500,
                chains = 3)   
summary(fit.jags)



# create predictions for p and psi
draws <- 100

p <- predict(object = mod.jags,
             mcmc.list = fit.jags,
             type = "p_array",
             draws = draws)
# output is in order [station, species, draw, occasion]
 
 psi <- predict(object = mod.jags,
                mcmc.list = fit.jags,
                type = "psi_array",
                draws = draws)
 # output is in order [station, species, draw]



ppc_comm <- PPC.community(
  p = p,
  psi = psi,
  y = mod.jags@input$ylist, 
  model = "Occupancy",
  type = "FT")
  
# Bayesian p-values
 ppc_comm$BP
 
   
 str(ppc_comm$residuals)
 
 # get individual species PPC results
 ppc_species <- ppc_comm$residuals[[1]] # first species
 
 plot(apply(ppc_species$res.obs, 2, mean), apply(ppc_species$res.new, 2, mean),
 xlab = "Observed residuals",
 ylab = "Predicted residuals"
 )
 
 abline(0,1)  # diagonal line is not visible due to tiny data set
  
 }