Background

Both field seasons are pooled into a single capture history. Individual IDs are suffixed with the year so that the same physical individual caught in both seasons is treated as independent. Year is carried as an individual-level covariate.

Three questions are addressed:

  1. Survival and movement (CJSsecr): Do morphs differ in daily apparent survival \(\phi\) and movement scale \(\sigma\)?

  2. Emergence timing (reverse-time CJSsecr): We reverse the capture history so that apparent survival equals seniority: \(\phi_\text{rev}[t] = \Pr(\text{present at } t-1 \mid \text{present at } t)\). Its complement \(1 - \phi_\text{rev}\) is the per-occasion entry rate. The same spatial kernel as forward CJSsecr is used, so morph differences in \(\sigma\) are explicitly modelled.

  3. Abundance trajectories: \(\hat{N}_{t,m} = n_{t,m}/\hat{p}\) where \(\hat{p}\) comes from the non-spatial CJS (pooled across morphs).

1. Build detector grid and mask

traps_m <- mothTransect |>
  st_as_sf(coords = c("GPS.E", "GPS.N"), crs = 4326) |>
  st_transform(3067) |>
  st_coordinates() |>
  as.data.frame() |>
  rename(x = X, y = Y)

transect_line <- traps_m |>
  st_as_sf(coords = c("x", "y"), crs = 3067) |>
  summarise(geometry = st_combine(geometry)) |>
  st_cast("LINESTRING")

search_area <- st_buffer(transect_line, dist = 200)
bbox        <- st_bbox(search_area)

grid_pts <- expand.grid(
  x = seq(bbox["xmin"], bbox["xmax"], by = 50),
  y = seq(bbox["ymin"], bbox["ymax"], by = 50)
) |>
  st_as_sf(coords = c("x", "y"), crs = 3067)

inside   <- st_within(grid_pts, search_area, sparse = FALSE)[, 1]
grid_pts <- grid_pts[inside, ] |>
  st_coordinates() |>
  as.data.frame() |>
  rename(x = X, y = Y)

traps_obj <- read.traps(data = grid_pts, detector = "proximity")
mask_obj  <- make.mask(traps = traps_obj, buffer = 200,
                       type = "trapbuffer", spacing = 100)

cat("Search area (ha):  ", round(as.numeric(st_area(search_area))/10000, 1), "\n")
## Search area (ha):   42.8
cat("Detector points:   ", nrow(grid_pts), "\n")
## Detector points:    167
cat("Mask points:       ", nrow(mask_obj), "\n")
## Mask points:        101

2. Build combined capture history

build_combined_capthist <- function() {

  dat <- mothCaptures |>
    filter(Year %in% c(2014, 2015), Morph %in% c("W", "Y"),
           !is.na(GPS.E), !is.na(GPS.N)) |>
    st_as_sf(coords = c("GPS.E", "GPS.N"), crs = 4326) |>
    st_transform(3067) |>
    mutate(cx = st_coordinates(geometry)[, 1],
           cy = st_coordinates(geometry)[, 2]) |>
    st_drop_geometry() |>
    group_by(Year) |>
    mutate(wday = as.integer(censusday - min(censusday))) |>
    ungroup() |>
    transmute(
      ID       = paste0(as.character(GenID), "_", Year),
      occasion = as.integer(censusday),
      cx       = cx,
      cy       = cy,
      Morph    = factor(Morph, levels = c("W", "Y")),
      Year     = as.integer(Year),
      wday     = wday
    ) |>
    group_by(ID, occasion) |>
    slice(1) |>
    ungroup() |>
    mutate(occasion = as.integer(match(occasion, sort(unique(occasion)))))

  n_occ <- n_distinct(dat$occasion)

  cap_sf     <- st_as_sf(dat, coords = c("cx", "cy"), crs = 3067)
  trap_sf    <- st_as_sf(grid_pts, coords = c("x", "y"), crs = 3067)
  dat$trapID <- st_nearest_feature(cap_sf, trap_sf)

  cat(sprintf("Combined: %d captures, %d individuals, %d occasions\n",
              nrow(dat), n_distinct(dat$ID), n_occ))

  captures_file <- tempfile(fileext = ".txt")
  traps_file    <- tempfile(fileext = ".txt")

  write.table(
    data.frame(session=1L, ID=dat$ID, occasion=dat$occasion, trapID=dat$trapID),
    captures_file, row.names=FALSE, col.names=FALSE, sep="\t", quote=FALSE)

  write.table(
    data.frame(trapID=seq_len(nrow(grid_pts)), x=grid_pts$x, y=grid_pts$y),
    traps_file, row.names=FALSE, col.names=FALSE, sep="\t", quote=FALSE)

  ch <- read.capthist(captfile=captures_file, trapfile=traps_file,
                      detector="proximity", noccasions=n_occ)

  cov_df <- dat |>
    group_by(ID) |> slice(1) |> ungroup() |>
    mutate(Year = factor(Year)) |>
    select(ID, Morph, Year, wday) |>
    as.data.frame()
  rownames(cov_df) <- cov_df$ID
  covariates(ch) <- cov_df[, c("Morph", "Year", "wday")]

  # Assign each occasion to a year based on which year's individuals
  # were captured there. Since IDs are GenID_Year, extract year from ID.
  occ_year_map <- dat |>
    mutate(ind_year = as.integer(sub(".*_", "", ID))) |>
    group_by(occasion) |>
    summarise(
      year_2014 = sum(ind_year == 2014),
      year_2015 = sum(ind_year == 2015),
      Year      = as.character(ifelse(year_2014 >= year_2015, 2014, 2015)),
      wday      = first(wday),
      .groups   = "drop"
    ) |>
    arrange(occasion)

  cat("Occasion-year assignment:\n"); print(occ_year_map)

  season_starts <- c("2014" = as.Date("2014-06-27"),
                     "2015" = as.Date("2015-06-22"))

  occ_meta <- occ_year_map |>
    select(occasion, Year, wday) |>
    mutate(
      date    = season_starts[Year] + wday,
      wday_sc = as.numeric(scale(wday))
    ) |>
    group_by(Year) |>
    mutate(occ_yr = seq_len(n())) |>
    ungroup()

  counts_by_occ <- dat |>
    group_by(occasion, Morph) |>
    summarise(n_caught = n(), .groups = "drop")

  list(ch=ch, occ_meta=occ_meta, counts_by_occ=counts_by_occ)
}

result        <- build_combined_capthist()
## Combined: 225 captures, 151 individuals, 18 occasions
## Occasion-year assignment:
## # A tibble: 18 × 5
##    occasion year_2014 year_2015 Year   wday
##       <int>     <int>     <int> <chr> <int>
##  1        1         6         1 2014      0
##  2        2         6         0 2014      1
##  3        3         5         0 2014      2
##  4        4         2         0 2014      3
##  5        5         6         3 2014      4
##  6        6         4         1 2014      5
##  7        7         6        10 2015      6
##  8        8         5         6 2015      7
##  9        9         4        22 2015      8
## 10       10         4        12 2015      9
## 11       11         4        15 2015     10
## 12       12         1         5 2015     11
## 13       13         9        14 2015     12
## 14       14         1        25 2015     13
## 15       15         5        10 2015     14
## 16       16         3         7 2015     15
## 17       17         2        10 2015     16
## 18       18         3         8 2015     17
ch_both       <- result$ch
occ_meta      <- result$occ_meta
counts_by_occ <- result$counts_by_occ

cat("\nCapture history summary:\n"); print(summary(ch_both))
## 
## Capture history summary:
## Object class       capthist 
## Detector type      proximity 
## Detector number    167 
## Average spacing    50 m 
## x-range            431943.8 432843.8 m 
## y-range            6896024 6896774 m 
## 
## Counts by occasion 
##                     1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
## n                   7   6   5   2   9   5  16  11  26  16  19   6  23  26  15
## u                   7   5   5   0   5   1  13   7  22  11  11   5  17  15   7
## f                 104  27  14   5   1   0   0   0   0   0   0   0   0   0   0
## M(t+1)              7  12  17  17  22  23  36  43  65  76  87  92 109 124 131
## losses              0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
## detections          7   6   5   2   9   5  16  11  26  16  19   6  23  26  15
## detectors visited   4   5   4   2   9   5   8   8  11  10  10   5  15  11   9
## detectors used    167 167 167 167 167 167 167 167 167 167 167 167 167 167 167
##                    16  17  18 Total
## n                  10  12  11   225
## u                   6   7   7   151
## f                   0   0   0   151
## M(t+1)            137 144 151   151
## losses              0   0   0     0
## detections         10  12  11   225
## detectors visited   5   9   8   138
## detectors used    167 167 167  3006
## 
## Individual covariates
##  Morph     Year         wday       
##  W:102   2014:55   Min.   : 0.000  
##  Y: 49   2015:96   1st Qu.: 7.000  
##                    Median : 9.000  
##                    Mean   : 9.437  
##                    3rd Qu.:13.000  
##                    Max.   :17.000
cat("\nOccasion metadata:\n"); print(occ_meta)
## 
## Occasion metadata:
## # A tibble: 18 × 6
##    occasion Year   wday date       wday_sc occ_yr
##       <int> <chr> <int> <date>       <dbl>  <int>
##  1        1 2014      0 2014-06-27 -1.59        1
##  2        2 2014      1 2014-06-28 -1.40        2
##  3        3 2014      2 2014-06-29 -1.22        3
##  4        4 2014      3 2014-06-30 -1.03        4
##  5        5 2014      4 2014-07-01 -0.843       5
##  6        6 2014      5 2014-07-02 -0.656       6
##  7        7 2015      6 2015-06-28 -0.468       1
##  8        8 2015      7 2015-06-29 -0.281       2
##  9        9 2015      8 2015-06-30 -0.0937      3
## 10       10 2015      9 2015-07-01  0.0937      4
## 11       11 2015     10 2015-07-02  0.281       5
## 12       12 2015     11 2015-07-03  0.468       6
## 13       13 2015     12 2015-07-04  0.656       7
## 14       14 2015     13 2015-07-05  0.843       8
## 15       15 2015     14 2015-07-06  1.03        9
## 16       16 2015     15 2015-07-07  1.22       10
## 17       17 2015     16 2015-07-08  1.40       11
## 18       18 2015     17 2015-07-09  1.59       12

3. Sanity check — non-spatial CJS

fit_cjs_test <- openCR.fit(ch_both, type="CJS",
                            model=list(phi ~ Year, p ~ Year), trace=FALSE)
cat("Non-spatial CJS:\n"); print(predict(fit_cjs_test))
## Non-spatial CJS:
## $p
##    session Year  estimate SE.estimate       lcl       ucl
## 1        1 2014        NA          NA        NA        NA
## 2        2 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 3        3 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 4        4 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 5        5 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 6        6 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 7        7 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 8        8 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 9        9 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 10      10 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 11      11 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 12      12 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 13      13 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 14      14 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 15      15 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 16      16 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 17      17 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 18      18 2014 0.2177676  0.06305558 0.1187593 0.3651189
## 
## $phi
##    session Year estimate SE.estimate       lcl       ucl
## 1        1 2014 0.665787  0.06805774 0.5224155 0.7839204
## 2        2 2014 0.665787  0.06805774 0.5224155 0.7839204
## 3        3 2014 0.665787  0.06805774 0.5224155 0.7839204
## 4        4 2014 0.665787  0.06805774 0.5224155 0.7839204
## 5        5 2014 0.665787  0.06805774 0.5224155 0.7839204
## 6        6 2014 0.665787  0.06805774 0.5224155 0.7839204
## 7        7 2014 0.665787  0.06805774 0.5224155 0.7839204
## 8        8 2014 0.665787  0.06805774 0.5224155 0.7839204
## 9        9 2014 0.665787  0.06805774 0.5224155 0.7839204
## 10      10 2014 0.665787  0.06805774 0.5224155 0.7839204
## 11      11 2014 0.665787  0.06805774 0.5224155 0.7839204
## 12      12 2014 0.665787  0.06805774 0.5224155 0.7839204
## 13      13 2014 0.665787  0.06805774 0.5224155 0.7839204
## 14      14 2014 0.665787  0.06805774 0.5224155 0.7839204
## 15      15 2014 0.665787  0.06805774 0.5224155 0.7839204
## 16      16 2014 0.665787  0.06805774 0.5224155 0.7839204
## 17      17 2014 0.665787  0.06805774 0.5224155 0.7839204
## 18      18 2014       NA          NA        NA        NA
p_hat_cjs <- predict(fit_cjs_test)$p |>
  filter(is.finite(estimate)) |>
  pull(estimate) |>
  mean()
cat(sprintf("\nPooled per-occasion p = %.4f\n", p_hat_cjs))
## 
## Pooled per-occasion p = 0.2178

4. Part 1 — CJSsecr: survival and movement scale

fit_or_load <- function(path, expr) {
  if (!refit && file.exists(path)) {
    cat("Loading cached:", path, "\n"); return(readRDS(path))
  }
  fit <- eval(expr); saveRDS(fit, path); fit
}

fit_null <- fit_or_load("fitD_null.rds", quote(
  openCR.fit(ch_both, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year, p ~ Year, sigma ~ 1),
             ncores=ncores, trace=FALSE)
))
## Loading cached: fitD_null.rds
fit_sigma_morph <- fit_or_load("fitD_sigma.rds", quote(
  openCR.fit(ch_both, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year, p ~ Year, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitD_sigma.rds
fit_phi_sigma <- fit_or_load("fitD_phi_sigma.rds", quote(
  openCR.fit(ch_both, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year + h2, p ~ Year, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitD_phi_sigma.rds
fit_all_morph <- fit_or_load("fitD_all.rds", quote(
  openCR.fit(ch_both, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year + h2, p ~ Year + h2, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitD_all.rds

5. CJSsecr model comparison

aic_raw <- AIC(fit_null, fit_sigma_morph, fit_phi_sigma, fit_all_morph) |>
  as.data.frame()

aic_col <- names(aic_raw)[sapply(aic_raw, is.numeric) &
                           grepl("AIC", names(aic_raw), ignore.case=TRUE)][1]

aic_tbl <- aic_raw |>
  mutate(
    model     = c("year only", "sigma ~ morph",
                  "phi + sigma ~ morph", "phi + p + sigma ~ morph"),
    delta_AIC = round(.data[[aic_col]] - min(.data[[aic_col]]), 2),
    AIC       = round(.data[[aic_col]], 2)
  ) |>
  select(model, AIC, delta_AIC) |>
  arrange(delta_AIC)

kable(aic_tbl, col.names=c("Model","AIC","ΔAIC"),
      caption="CJSsecr model comparison (both years pooled)") |>
  kableExtra::kable_styling(full_width=FALSE) |>
  kableExtra::row_spec(1, bold=TRUE, background="#f0f7e6")
CJSsecr model comparison (both years pooled)
Model AIC ΔAIC
fit_sigma_morph year only 1200.63 0.00
fit_phi_sigma sigma ~ morph 1200.66 0.02
fit_all_morph phi + sigma ~ morph 1200.66 0.02
fit_null phi + p + sigma ~ morph 1205.02 4.39

6. CJSsecr parameter estimates

best_fit <- fit_sigma_morph
preds    <- predict(best_fit)

sigma_tbl <- preds$sigma |>
  filter(is.finite(estimate), is.finite(SE.estimate)) |>
  group_by(h2) |> slice(1) |> ungroup() |>
  mutate(morph = ifelse(as.character(h2) == "1", "W", "Y")) |>
  select(morph, estimate, SE.estimate, lcl, ucl)

phi_tbl <- predict(fit_phi_sigma)$phi |>
  filter(is.finite(estimate)) |>
  group_by(h2) |>
  summarise(estimate=mean(estimate), SE.estimate=mean(SE.estimate),
            lcl=mean(lcl), ucl=mean(ucl), .groups="drop") |>
  mutate(morph = ifelse(as.character(h2) == "1", "W", "Y"))

kable(sigma_tbl, digits=1, caption="Movement scale sigma (m) by morph") |>
  kableExtra::kable_styling(full_width=FALSE)
Movement scale sigma (m) by morph
morph estimate SE.estimate lcl ucl
W 92.6 19.5 61.3 139.8
Y 47.0 5.1 37.9 58.3 
kable(phi_tbl |> select(morph, estimate, SE.estimate, lcl, ucl),
      digits=3, caption="Daily apparent survival phi by morph") |>
  kableExtra::kable_styling(full_width=FALSE)
Daily apparent survival phi by morph
morph estimate SE.estimate lcl ucl
W 0.536 0.179 0.221 0.825
Y 0.803 0.063 0.652 0.898 
data_dir <- system.file("data", package="plantaginisCMR")
saveRDS(list(sigma_tbl=sigma_tbl, phi_tbl=phi_tbl,
             aic_tbl=aic_tbl, best_fit=best_fit, preds=preds),
        file.path(data_dir, "modD_objects.rds"))
cat("Saved modD_objects.rds\n")
## Saved modD_objects.rds

7. Part 2 — Emergence timing (reverse-time CJSsecr)

extract_aic_val <- function(fit) {
  tryCatch({
    a   <- AIC(fit)
    col <- names(a)[sapply(a, is.numeric) &
                    grepl("AIC", names(a), ignore.case=TRUE)][1]
    if (is.na(col)) return(NA_real_)
    val <- a[[col]]
    if (!is.finite(val)) NA_real_ else val
  }, error=function(e) NA_real_)
}

reverse_capthist <- function(ch) {
  n_occ   <- ncol(ch)
  arr     <- unclass(ch)
  arr_rev <- arr[, n_occ:1, , drop=FALSE]
  class(arr_rev)      <- class(ch)
  traps(arr_rev)      <- traps(ch)
  covariates(arr_rev) <- covariates(ch)
  arr_rev
}

ch_rev <- reverse_capthist(ch_both)
cat(sprintf("Reversed capthist: %d individuals, %d occasions\n",
            nrow(ch_rev), ncol(ch_rev)))
## Reversed capthist: 151 individuals, 18 occasions
fit_rev_null <- fit_or_load("fitR_null.rds", quote(
  openCR.fit(ch_rev, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year, p ~ Year, sigma ~ 1),
             ncores=ncores, trace=FALSE)
))
## Loading cached: fitR_null.rds
fit_rev_sigma <- fit_or_load("fitR_sigma.rds", quote(
  openCR.fit(ch_rev, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year, p ~ Year, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitR_sigma.rds
fit_rev_phi_sigma <- fit_or_load("fitR_phi_sigma.rds", quote(
  openCR.fit(ch_rev, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year + h2, p ~ Year, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitR_phi_sigma.rds
fit_rev_all <- fit_or_load("fitR_all.rds", quote(
  openCR.fit(ch_rev, mask=mask_obj, type="CJSsecr",
             model=list(phi ~ Year + h2, p ~ Year + h2, sigma ~ h2),
             groups="Morph", ncores=ncores, trace=FALSE)
))
## Loading cached: fitR_all.rds
aic_rev_raw <- AIC(fit_rev_null, fit_rev_sigma,
                   fit_rev_phi_sigma, fit_rev_all) |> as.data.frame()
aic_rev_col <- names(aic_rev_raw)[sapply(aic_rev_raw, is.numeric) &
                grepl("AIC", names(aic_rev_raw), ignore.case=TRUE)][1]

aic_rev_tbl <- aic_rev_raw |>
  mutate(
    model     = c("year only", "sigma ~ morph",
                  "seniority + sigma ~ morph",
                  "seniority + p + sigma ~ morph"),
    delta_AIC = round(.data[[aic_rev_col]] - min(.data[[aic_rev_col]]), 2),
    AIC       = round(.data[[aic_rev_col]], 2)
  ) |>
  select(model, AIC, delta_AIC) |>
  arrange(delta_AIC)

kable(aic_rev_tbl, col.names=c("Model","AIC","ΔAIC"),
      caption="Reverse-time CJSsecr: seniority model comparison") |>
  kableExtra::kable_styling(full_width=FALSE) |>
  kableExtra::row_spec(1, bold=TRUE, background="#f0f7e6")
Reverse-time CJSsecr: seniority model comparison
Model AIC ΔAIC
fit_rev_sigma year only 1196.57 0.00
fit_rev_phi_sigma sigma ~ morph 1198.00 1.43
fit_rev_all seniority + sigma ~ morph 1198.00 1.43
fit_rev_null seniority + p + sigma ~ morph 1201.45 4.89 
# Always use morph-specific seniority model for the emergence comparison
preds_rev <- predict(fit_rev_phi_sigma)
n_occ_fwd <- ncol(ch_both)

phi_rev <- preds_rev$phi |>
  filter(is.finite(estimate)) |>
  mutate(
    entry_rate   = 1 - estimate,
    entry_lcl    = 1 - ucl,
    entry_ucl    = 1 - lcl,
    fwd_occasion = (n_occ_fwd + 1) - as.integer(session),
    Morph        = factor(ifelse(as.character(h2)=="1","W","Y"),
                          levels=c("W","Y"))
  )

cat("\nEntry rates by morph:\n")
## 
## Entry rates by morph:
print(phi_rev |>
      select(fwd_occasion, Morph, Year, estimate, entry_rate,
             entry_lcl, entry_ucl) |>
      head(12))
##    fwd_occasion Morph Year  estimate entry_rate  entry_lcl entry_ucl
## 1            18     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 2            17     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 3            16     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 4            15     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 5            14     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 6            13     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 7            12     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 8            11     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 9            10     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 10            9     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 11            8     W 2014 0.7091831  0.2908169 0.09490887 0.6159229
## 12            7     W 2014 0.7091831  0.2908169 0.09490887 0.6159229

8. Part 3 — Abundance trajectories

ht_df <- counts_by_occ |>
  mutate(
    N_hat    = n_caught / p_hat_cjs,
    N_hat_se = sqrt(n_caught * (1 - p_hat_cjs)) / p_hat_cjs,
    N_hat_lo = pmax(n_caught, N_hat - 1.96 * N_hat_se),
    N_hat_hi = N_hat + 1.96 * N_hat_se
  ) |>
  left_join(occ_meta |> select(occasion, date, Year, wday, occ_yr),
            by="occasion") |>
  mutate(Morph = factor(Morph, levels=c("W","Y")))

cat(sprintf("p_hat = %.4f\n", p_hat_cjs))
## p_hat = 0.2178
cat("Year split:\n")
## Year split:
print(ht_df |> group_by(Year) |>
      summarise(n=n(), dates=paste(min(date), max(date), sep=" to "),
                .groups="drop"))
## # A tibble: 2 × 3
##   Year      n dates                   
##   <chr> <int> <chr>                   
## 1 2014     11 2014-06-27 to 2014-07-02
## 2 2015     24 2015-06-28 to 2015-07-09

9. Figure 1 — Survival and emergence timing

p_phi <- ggplot(phi_tbl, aes(x=morph, y=estimate, colour=morph)) +
  geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.15, linewidth=0.8) +
  geom_point(size=4) +
  scale_colour_manual(values=morph_colours, labels=c(W="White", Y="Yellow")) +
  scale_y_continuous(limits=c(0, 1),
                     labels=scales::number_format(accuracy=0.01)) +
  labs(x=NULL, y=expression("Daily apparent survival ("*phi*")")) +
  theme_moth() +
  theme(legend.position="none")

entry_morph <- phi_rev |>
  filter(is.finite(entry_rate), fwd_occasion > 1) |>
  group_by(Morph) |>
  summarise(estimate = mean(entry_rate),
            lcl      = mean(entry_lcl),
            ucl      = mean(entry_ucl),
            .groups  = "drop")

cat("Mean entry rates:\n"); print(entry_morph)
## Mean entry rates:
## # A tibble: 2 × 4
##   Morph estimate    lcl   ucl
##   <fct>    <dbl>  <dbl> <dbl>
## 1 W        0.291 0.0949 0.616
## 2 Y        0.189 0.0996 0.331
p_entry <- ggplot(entry_morph, aes(x=Morph, y=estimate, colour=Morph)) +
  geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.15, linewidth=0.8) +
  geom_point(size=4) +
  scale_colour_manual(values=morph_colours, labels=c(W="White", Y="Yellow")) +
  scale_y_continuous(limits=c(0, 1),
                     labels=scales::number_format(accuracy=0.01)) +
  labs(x=NULL, y="Mean per-capita entry rate (1 \u2212 seniority)") +
  theme_moth() +
  theme(legend.position="none")

p_phi + p_entry +
  plot_annotation(tag_levels="A") &
  theme(plot.tag=element_text(size=14, face="bold"))

ggsave("figD1_survival_emergence.pdf", width=10, height=5)

10. Figure 2 — Movement kernel and female encounter rate

hr_radius_W <- 2.45 * sigma_tbl$estimate[sigma_tbl$morph=="W"]
hr_radius_Y <- 2.45 * sigma_tbl$estimate[sigma_tbl$morph=="Y"]

build_kernel <- function(sigma_est, sigma_se, morph_label) {
  if (!is.finite(sigma_est) || !is.finite(sigma_se) || sigma_est <= 0)
    return(tibble())
  dist_seq <- seq(0, (sigma_est + sigma_se) * 5, length.out=400)
  tibble(dist=dist_seq, morph=morph_label,
         density=dnorm(dist_seq, 0, sigma_est),
         dens_lo=dnorm(dist_seq, 0, max(sigma_est - sigma_se, 0.1)),
         dens_hi=dnorm(dist_seq, 0, sigma_est + sigma_se))
}

kernel_df <- sigma_tbl |>
  filter(is.finite(estimate), is.finite(SE.estimate)) |>
  rowwise() |>
  do(build_kernel(.$estimate, .$SE.estimate, .$morph)) |>
  ungroup()

hr_lines <- tibble(
  morph  = c("W","Y"),
  radius = c(hr_radius_W, hr_radius_Y),
  label  = sprintf("95%% HR\n%.0f m", c(hr_radius_W, hr_radius_Y))
)

p_kernel <- ggplot(kernel_df, aes(x=dist, colour=morph, fill=morph)) +
  geom_ribbon(aes(ymin=dens_lo, ymax=dens_hi), alpha=0.2, colour=NA) +
  geom_line(aes(y=density), linewidth=1.1) +
  geom_vline(data=hr_lines,
             aes(xintercept=radius, colour=morph),
             linetype="dashed", linewidth=0.7) +
  geom_text(data=hr_lines,
            aes(x=radius, y=Inf, label=label, colour=morph),
            vjust=1.4, hjust=-0.1, size=3, show.legend=FALSE) +
  scale_colour_manual(values=morph_colours, labels=c(W="White", Y="Yellow")) +
  scale_fill_manual(values=morph_fills,     labels=c(W="White", Y="Yellow")) +
  labs(x="Distance moved between occasions (m)",
       y="Probability density", colour="Morph", fill="Morph") +
  theme_moth()

modF_path <- local({
  # Check working directory first (where F vignette saves it)
  wd_path <- "modF_objects.rds"
  if (file.exists(wd_path)) return(wd_path)
  # Fall back to package data directory
  pkg_path <- system.file("data", "modF_objects.rds", package="plantaginisCMR")
  if (nchar(pkg_path) > 0 && file.exists(pkg_path)) return(pkg_path)
  return("")
})
cat("modF path:", modF_path, "\n")
## modF path:
if (nchar(modF_path) > 0 && file.exists(modF_path)) {
  modF_data        <- readRDS(modF_path)
  encounter_scores <- modF_data$encounter_scores
  modF1            <- modF_data$modF1

  # Extract marginal mean encounter score per morph from the mixed model
  # on the log scale, then back-transform to the original scale.
  # W is the reference level; Y effect is the MorphY fixed effect.
  fe      <- fixef(modF1)
  se_fe   <- sqrt(diag(vcov(modF1)))

  # Mean log(n_females+1) across the dataset (to marginalise over covariate)
  log_nfem_mean <- mean(log(encounter_scores$n_females + 1), na.rm=TRUE)

  lp_W  <- fe["(Intercept)"] + fe["log(n_females + 1)"] * log_nfem_mean
  lp_Y  <- lp_W + fe["MorphY"]

  # SE via delta method (covariance between intercept and MorphY)
  V       <- vcov(modF1)
  se_lp_W <- sqrt(V["(Intercept)","(Intercept)"] +
                  log_nfem_mean^2 * V["log(n_females + 1)","log(n_females + 1)"] +
                  2 * log_nfem_mean * V["(Intercept)","log(n_females + 1)"])
  se_lp_Y <- sqrt(se_lp_W^2 + V["MorphY","MorphY"] +
                  2 * V["(Intercept)","MorphY"] +
                  2 * log_nfem_mean * V["log(n_females + 1)","MorphY"])

  encounter_est <- tibble(
    Morph    = factor(c("W","Y"), levels=c("W","Y")),
    lp       = c(lp_W, lp_Y),
    lp_se    = c(se_lp_W, se_lp_Y),
    estimate = exp(lp) - 0.001,          # back-transform offset
    lcl      = exp(lp - 1.96*lp_se) - 0.001,
    ucl      = exp(lp + 1.96*lp_se) - 0.001
  )

  cat("Model-estimated encounter scores:\n"); print(encounter_est)

  p_encounter <- ggplot(encounter_est,
                        aes(x=Morph, y=estimate, colour=Morph)) +
    geom_errorbar(aes(ymin=lcl, ymax=ucl), width=0.15, linewidth=0.8) +
    geom_point(size=4) +
    scale_colour_manual(values=morph_colours, labels=c(W="White", Y="Yellow")) +
    scale_y_continuous(labels=scales::number_format(accuracy=0.001)) +
    labs(x=NULL,
         y=expression(sigma*"-corrected female encounter score")) +
    theme_moth() +
    theme(legend.position="none")

  p_kernel + p_encounter +
    plot_annotation(tag_levels="A") &
    theme(plot.tag=element_text(size=14, face="bold"))

} else {
  cat("Module F output not found — run F_female_proximity.Rmd first\n")
  print(p_kernel)
}
## Module F output not found — run F_female_proximity.Rmd first

ggsave("figD2_movement_encounter.pdf", width=10, height=5)

11. Figure 3 — Estimated population size through time

# Two panels: one per year, x-axis = actual date within season
ht_plot <- ht_df |>
  select(-any_of(c("occ_yr","date","Year","wday"))) |>
  left_join(occ_meta |> select(occasion, Year, date, occ_yr), by="occasion")

cat("Year split in ht_plot:\n")
## Year split in ht_plot:
print(ht_plot |> group_by(Year) |>
      summarise(n=n(), from=min(date), to=max(date), .groups="drop"))
## # A tibble: 2 × 4
##   Year      n from       to        
##   <chr> <int> <date>     <date>    
## 1 2014     11 2014-06-27 2014-07-02
## 2 2015     24 2015-06-28 2015-07-09
make_abund_panel <- function(df, yr) {
  ggplot(df, aes(x=date, y=N_hat, colour=Morph, fill=Morph)) +
    geom_ribbon(aes(ymin=N_hat_lo, ymax=N_hat_hi), alpha=0.2, colour=NA) +
    geom_line(linewidth=1) +
    geom_point(size=2.5) +
    scale_colour_manual(values=morph_colours, labels=c(W="White", Y="Yellow")) +
    scale_fill_manual(values=morph_fills,     labels=c(W="White", Y="Yellow")) +
    scale_x_date(date_labels="%d %b", date_breaks="2 days") +
    labs(x=NULL, y=expression(hat(N)[t]),
         colour="Morph", fill="Morph", title=yr) +
    theme_moth() +
    theme(axis.text.x=element_text(angle=45, hjust=1))
}

years_present <- sort(unique(ht_plot$Year))
cat("Years:", years_present, "\n")
## Years: 2014 2015
plot_list <- lapply(years_present, function(yr)
  make_abund_panel(filter(ht_plot, Year==yr), yr))

wrap_plots(plot_list, nrow=1, guides="collect") +
  plot_annotation(
    tag_levels="A",
    caption=sprintf("N_hat = n/p,  p = %.3f (non-spatial CJS, pooled morphs)",
                    p_hat_cjs)) &
  theme(plot.tag=element_text(size=14, face="bold"),
        legend.position="right")

ggsave("figD3_abundance.pdf", width=12, height=5)

12. Save all outputs

saveRDS(
  list(sigma_tbl   = sigma_tbl,
       phi_tbl     = phi_tbl,
       aic_tbl     = aic_tbl,
       aic_rev_tbl = aic_rev_tbl,
       phi_rev     = phi_rev,
       ht_df       = ht_df,
       p_hat_cjs   = p_hat_cjs),
  "modD_summary.rds"
)
cat("Saved modD_summary.rds\n")
## Saved modD_summary.rds

13. Notes for paper

CJSsecr (forward): \(\sigma_W = 93\) m vs \(\sigma_Y = 47\) m — white males range ~twice as far. Their lower recapture rate reflects wider ranging, not lower survival. Apparent survival: \(\phi_W = 0.54\) vs \(\phi_Y = 0.80\), CIs overlap.

Reverse-time CJSsecr: \(\phi_\text{rev}[t] = \Pr(\text{present at } t-1 \mid \text{present at } t)\). Entry rate \(= 1 - \phi_\text{rev}\). The spatial kernel accounts for morph differences in \(\sigma\), avoiding the detection confound of non-spatial Pradel models. Always use fit_rev_phi_sigma for morph comparison regardless of AIC ranking.

Abundance (HT): \(\hat{N}_{t,m} = n_{t,m}/\hat{p}\), \(\hat{p} \approx 0.22\) (non-spatial CJS, pooled morphs). Since \(\hat{p}\) is the same for both morphs, the W/Y ratio is unbiased. Treat absolute values with caution.

Home range: 95% HR radius \(= 2.45\sigma\): \(\approx 228\) m (white) vs \(\approx 115\) m (yellow).

Set refit <- FALSE after first successful run.