Creating Predicted plots in R from PRESENCE models

[NicolaG]

Hi!
I am trying to create plots from the output of my single-season occupancy model in PRESENCE. I’m using R to create the plots. I am struggling to find somewhere that shows the code for when you haven’t used an R package to create the occupancy model. I wanted to check if my code is correct. I have mostly used linear covariates but have a few models with a quadratic term, and one for month where I used sine and cosine. In these examples, my covariates were standardised in Excel before inputting them into PRESENCE.

For the predicted detection probability plots for a linear covariate, this is the code I used for the model psi(.), p(TempMinST):

Code: Select all

# --- Input values from PRESENCE output ---
# Betas
intercept <- -1.467973   # B1: intercept for detection
slope     <- -0.694926   # B2: slope for TempMinST

# Variance-covariance components
var_intercept <- 0.158806 #VC matrix: B1-B1
var_slope     <- 0.083539 #VC matrix: B2-B2
cov_is        <- 0.042429 #VC matrix: B1-B2

# Raw covariate stats
temp_mean <- 19.56622024
temp_sd   <- 3.11020065
temp_min  <- 13.66666667
temp_max  <- 26

# Sequence of raw temperature values for plotting
temp_vals_raw <- seq(temp_min, temp_max, length.out = 100)

# Standardise these values
temp_vals_std <- (temp_vals_raw - temp_mean) / temp_sd

# Predict detection probability on logit scale
logit_p <- intercept + slope * temp_vals_std

# Compute SE of logit predictions
se_logit <- sqrt(
  var_intercept +
    (temp_vals_std^2 * var_slope) +
    (2 * temp_vals_std * cov_is)
)

# Compute 95% confidence intervals on logit scale
logit_upper <- logit_p + 1.96 * se_logit
logit_lower <- logit_p - 1.96 * se_logit

# Back-transform to probability scale
p_hat   <- plogis(logit_p)
p_upper <- plogis(logit_upper)
p_lower <- plogis(logit_lower)

# --- Plot ---
plot(temp_vals_raw, p_hat, type = "l", lwd = 2, col = "darkblue",
     ylim = c(0, 1),
     xlab = "Minimum Temperature (°C)",
     ylab = "Detection Probability",
     main = "Detection vs. Minimum Temperature")

# Add CI lines
lines(temp_vals_raw, p_upper, lty = 2, col = "gray40")
lines(temp_vals_raw, p_lower, lty = 2, col = "gray40")


I wasn't sure how to deal with the fact standardised values went into the model and how/if I need to standardised the predicted values for the plot. I'd appreciate any advice!

Also, is there a set time when you would report a predicted probability plot vs a plot with the Individual Site Estimates of <P[1]> and their confidence intervals?

Thanks!
Nicola

[jhines]

You are correct in standardizing the temperatures for the plot. When you standardized the input covariates to the model, the beta estimates are based on standardized temperatures. With un-standardized temperature covariates, the beta estimate represents the amount of change in detection per degree of temperature. With standardized temperatures, the beta estimate represents the amount of change in detection per standardized unit of temperature.