Applied analyses: justice, health, environment and rescue

Courts and crime: flow, time, outcomes, place and case mix

Court and police statistics are published as counts: cases received and disposed each month, shares disposed within so many months, incidents with a date and a place. This vignette works through the questions those counts are asked to answer, on the Ontario Court of Justice’s criminal statistics and the Toronto Police Service’s open data (see ?justice-data for every file and its source). Each section uses the excerpt shipped with the package; the full tables, every edition and every court, load from the data hub:

clearance <- dload("ocj_statistics/clearance_rate_criminal")
offences  <- dload("ocj_statistics/offence_based_criminal")
mci       <- dload("tps_open_data/major_crime_indicators_open_data_0")

Is the court keeping up? Case flow and clearance

A court keeps up when it disposes of as many cases as it receives. caseflow() turns the monthly counts into the clearance rate (disposed per case received) with an exact interval, its rolling twelve-month value, and, from the pending count, how many months the inventory would take to clear at the recent pace.

f <- ext("ocj_flow_ontario.csv")
cf <- caseflow(f$received, f$disposed, pend = f$pending, prd = f$month)
cf[
  c(1, 12, 24),
  c("period", "received", "disposed", "clearance", "rolclr", "pertoclear")
]
#>    period received disposed clearance    rolclr pertoclear
#> 1 2024-04    21265    22053  1.037056        NA   6.752097
#> 2 2025-03    20688    20996  1.014888 0.9742875   7.466076
#> 3 2026-03    20602    23136  1.122998 0.9937306   7.196190
plot(cf)

Over these 24 months the court disposed of 489,448 cases against 497,393 received, a clearance rate of 0.984. The monthly rate swings well beyond its interval – summer months clear less – which is why the rolling rate is the one to read.

How long do cases take? A life table from banded times

The Court publishes the share of cases disposed within 10 months, 10 to 15, 15 to 18 and over 18. lifetbl() reads those bands directly and gives the share still open at each band edge, with a Greenwood interval, and the median.

d <- ext("ocj_time_to_disposition.csv")
last <- nrow(d)
shares <- unlist(d[last, c("under_10", "m10_15", "m15_18", "over_18")])
lt <- lifetbl(
  c(
    "Under 10 Months",
    "Between 10-15 Months",
    "Between 15-18 Months",
    "Over 18 Months"
  ),
  events = shares * f$disposed[last]
)
lt
#> ╭── Life table (actuarial) ────────────────────────────────────────────────╮
#> │   Bands     4                                                            │
#> │   Cases     23,136                                                       │
#> │   Events    23,136                                                       │
#> │   Censored  0                                                            │
#> │   Interval  95% log-log Greenwood                                        │
#> │                                                                          │
#> │   lower  upper  entering  events  censored  atrisk      q  survival      │
#> │   -----  -----  --------  ------  --------  ------  -----  --------      │
#> │       0     10    23,136  18,185         0  23,136  0.786   0.21400      │
#> │      10     15     4,951   2,869         0   4,951  0.579   0.09000      │
#> │      15     18     2,082     879         0   2,082  0.422   0.05200      │
#> │      18    Inf     1,203   1,203         0   1,203  1.000   0.00000      │
#> │                                                                          │
#> │   lower  95% CI              hazard                                      │
#> │   -----  ------------------  ------                                      │
#> │       0  [0.20874, 0.21931]  0.1295                                      │
#> │      10  [0.08636, 0.09373]  0.1632                                      │
#> │      15  [0.04919, 0.05491]  0.1784                                      │
#> │      18                                                                  │
#> │                                                                          │
#> │   Median time: 6.36 (linear within its band)                             │
#> │   survival is the share still without the event at the end of each band. │
#> │                                                                          │
#> │   A censored case counts as half a case at risk in its band.             │
#> ╰──────────────────────────────────────────────────────────────────────────╯
plot(lt)

The median of 6.4 months assumes cases are spread evenly within the first band; the published average for the month is 185 days. Only disposed cases are in these statistics, so the table describes the cases that ended, not every case received.

How do cases end? The mix of outcomes

x <- ext("ocj_dispositions.csv")
mix <- dispmix(x[, -(1:2)], prd = x$year, scores = x$year)
mix$trend[, c("outcome", "statistic", "pholm", "change")]
#>                         outcome    statistic        pholm        change
#> 1 withdrawn_stayed_before_trial 30150.915355 0.000000e+00  0.1292036827
#> 2      guilty_plea_before_trial 10825.914695 0.000000e+00 -0.0805960172
#> 3            other_before_trial     7.296441 6.909133e-03 -0.0002635828
#> 4            withdrawn_at_trial  2890.841438 0.000000e+00 -0.0127761314
#> 5          guilty_plea_at_trial  6782.412494 0.000000e+00 -0.0218969033
#> 6                other_at_trial   175.420561 9.691647e-40 -0.0006064793
#> 7               following_trial  2044.569303 0.000000e+00 -0.0130645687
plot(mix)

Every share comes with Goodman’s simultaneous interval (in mix$shares). With a few hundred thousand cases a year every trend is “significant”; the change column says how large each shift is. The years are the scores, so the gaps at 2020 and 2023 (see ?justice-data) are spaced as they should be.

Is a change in the clearance rate real? Simpson’s paradox and case mix

An aggregate clearance rate is a weighted average of the rates of its parts – offence groups, courts – weighted by their share of the work. It can fall while most parts improve, if the work shifts toward the parts that clear more slowly. Demers and Rossmo (2015) found 210 such reversals in fifteen years of clearance rates for Canada’s fifty largest police services. clrdcmp() splits a change exactly into the change in the parts’ rates and the change in the mix (Kitagawa, 1955), and tests the within-part change with a Mantel-Haenszel rate ratio.

g <- ext("ocj_offence_groups.csv")
g <- g[g$year %in% c(2014, 2022), ]
cd <- clrdcmp(g$disposed, g$received, g$offence_group, g$year)
cd
#> ╭── Clearance rate, 2022 against 2014 (Kitagawa) ─────────────────────────────╮
#> │   Rate 2014                0.976                                            │
#> │   Rate 2022                0.9749                                           │
#> │   Change                   -0.0011 = rate effect +0.0096 + mix effect       │
#> │                            -0.0107                                          │
#> │   2022 at the 2014 mix     0.9912                                           │
#> │   2022 at the average mix  0.9830                                           │
#> │   Parts                    6                                                │
#> │                                                                             │
#> │   part                       rate 2014  rate 2022  share 2014  share 2022   │
#> │   -------------------------  ---------  ---------  ----------  ----------   │
#> │   Administration of Justice     0.9788      1.038      0.2689      0.2818   │
#> │   Crimes Against Property       0.9851      0.969      0.2279      0.1973   │
#> │   Crimes Against the Person     0.9180      0.858      0.2620      0.3187   │
#> │   Criminal Code Traffic         1.0105      1.021      0.0808      0.0948   │
#> │   Federal Statute               1.0232      1.173      0.1148      0.0534   │
#> │   Other Criminal Code           1.0675      1.082      0.0455      0.0539   │
#> │                                                                             │
#> │   part                       ratio  95% CI            rateeff     mixeff    │
#> │   -------------------------  -----  --------------  ---------  ---------    │
#> │   Administration of Justice  1.060  [1.048, 1.072]   0.016263   0.000424    │
#> │   Crimes Against Property    0.984  [0.971, 0.997]  -0.003375  -0.000052    │
#> │   Crimes Against the Person  0.934  [0.923, 0.945]  -0.017526  -0.004966    │
#> │   Criminal Code Traffic      1.011  [0.990, 1.031]   0.000939   0.000566    │
#> │   Federal Statute            1.147  [1.123, 1.171]   0.012631  -0.007535    │
#> │   Other Criminal Code        1.013  [0.987, 1.040]   0.000713   0.000834    │
#> │                                                                             │
#> │   4 of 6 parts rose, 2 fell.                                                │
#> │   Simpson's paradox through the mix: the aggregate fell, but at a fixed     │
#> │   mix it rose                                                               │
#> │   Rate effect: the change at a fixed case mix; mix effect: the change       │
#> │   from the mix of cases alone.                                              │
#> │                                                                             │
#> │   Rate ratio 2022/2014: crude 0.9989 (0.9929 to 1.0050); within parts       │
#> │   (Mantel-Haenszel) 1.0085 (1.0024 to 1.0146), p 0.006                      │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
plot(cd)

# a group's share of the cases received in a year, as a percentage
caseshare <- function(group, yr) {
  sprintf(
    "%.0f%%",
    100 *
      sum(g$received[g$offence_group == group & g$year == yr]) /
      sum(g$received[g$year == yr])
  )
}

Between 2014 and 2022 the province-wide clearance rate moved by -0.0011. Held at a fixed mix it moved by +0.0096; the shift of work contributed -0.0107. Crimes against the person, the group with the lowest clearance rate in both years, grew from 26% to 32% of the cases received, and federal-statute cases, the group with the highest rate in 2022, fell from 11% to 5%. 4 of the 6 groups cleared a larger share in 2022 and 2 a smaller one. The crude ratio’s interval includes 1; the Mantel-Haenszel ratio within offence groups is 1.0085 (1.0024 to 1.0146). Read alone, the aggregate would say nothing changed; within the groups, taken together, clearance improved, though not in every group.

The grouping matters, and so does the pair of years. By court, with Toronto’s courthouses taken together (they were consolidated in 2023), comparing 2012 with 2024:

c <- ext("ocj_courts.csv")
c <- aggregate(cbind(received, disposed) ~ year + court_series, c, sum)
c <- c[c$year %in% c(2012, 2024), ]
c <- c[c$court_series %in% names(which(table(c$court_series) == 2)), ]
cc <- clrdcmp(c$disposed, c$received, c$court_series, c$year)
cc$overall
#>       ratea     rateb       change     rateeff      mixeff staonbase  staonave
#> 1 0.9990086 0.9891387 -0.009869871 -0.01247238 0.002602504 0.9852411 0.9871899
cc$simpson$type
#> [1] "none"

Here the aggregate (-0.0099) and the rate at a fixed mix (-0.0125) move the same way: no reversal, and 27 of 54 courts cleared a smaller share in 2024. A decomposition does not always find a paradox; it says which way the parts moved and how much the mix did.

Where do incidents cluster? Kernel density

crmdens() estimates the density of incidents in incidents per square kilometre and outlines the smallest areas holding half and nine-tenths of them.

b <- ext("tps_break_enter_d14_2024.csv")
kd <- crmdens(b$lon, b$lat, ngrid = 150)
kd
#> ╭── Kernel density of incidents ───────────────────────────────────────╮
#> │   Incidents     508                                                  │
#> │   Coordinates   longitude and latitude, centred on -79.4114, 43.6508 │
#> │   Grid          150 x 150 cells                                      │
#> │   Bandwidth     0.337 km east-west, 0.412 km north-south             │
#> │   Peak density  77.92 incidents per km2                              │
#> │                                                                      │
#> │   The smallest areas holding each share of incidents                 │
#> │   incidents  areakm2  shrofmaparea  densatedge                       │
#> │   ---------  -------  ------------  ----------                       │
#> │   50%           5.50  11.9%              31.88                       │
#> │   90%          14.84  32.1%              11.23                       │
#> │                                                                      │
#> │   508 incidents; bandwidth 0.337 km east-west, 0.412 km north-south  │
#> │   50% of the incidents lie in 11.9% of the mapped area.              │
#> ╰──────────────────────────────────────────────────────────────────────╯
plot(kd)

Do incidents repeat nearby, soon after? The Knox test

Break-ins close in space and time to an earlier one are the near-repeat pattern. The Knox test compares the pairs close in both with the number expected if the dates were shuffled among the places; the expectation is exact, and a permutation test gives each cell’s p-value.

b <- b[!is.na(b$lon), ]
kn <- knoxtst(b$lon, b$lat, as.Date(b$occ_date), spcbands = c(100, 200, 400),
                tmbands = c(7, 14, 28), nsim = 199, seed = 1)
kn
#> ╭── Knox test of near-repeat incidents ───────────────────────────────────────╮
#> │   Incidents     508                                                         │
#> │   Permutations  199                                                         │
#> │   Close pairs   5,585                                                       │
#> │   Space bands   100, 200, 400                                               │
#> │   Time bands    7, 14, 28                                                   │
#> │                                                                             │
#> │   Knox ratios, observed / expected pairs                                    │
#> │   distance   0-7 d   7-14 d  14-28 d                                        │
#> │   ---------  ------  ------  -------                                        │
#> │   0-100 m    1.66**  1.16    1.09                                           │
#> │   100-200 m  1.25    0.91    1.17                                           │
#> │   200-400 m  1.07    1.03    1.01                                           │
#> │                                                                             │
#> │   A ratio above 1: more pairs this close in space and in time than          │
#> │   chance alone gives, the near-repeat pattern.                              │
#> │                                                                             │
#> │   * p < 0.05, ** p < 0.01, one-sided per cell (not adjusted for the number  │
#> │   of cells); - nothing expected                                             │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
plot(kn)

Do the hot spots last?

A hot spot that stays put is a place problem; one that moves is not. hotpers() counts the incidents in grid cells each period, reports the share of places that holds half of them, which cells are hot (the busiest 5% of the cells with any incident), and how much each period’s hot spots overlap the previous period’s:

h <- hotpers(b$lon, b$lat, quarters(as.Date(b$occ_date)), cell = 0.3)
h
#> ╭── Hot spots ──────────────────────────────────────────────────────────────╮
#> │   Active cells             139                                            │
#> │   Cell side                0.3 km                                         │
#> │   Periods                  4                                              │
#> │   Hot cells                the busiest 7 (5%)                             │
#> │   Hot in half the periods  7                                              │
#> │                                                                           │
#> │   By period                                                               │
#> │   period  incidents  actcells  cellsforhalf  shrforhalf  hotcells  tiedin │
#> │   ------  ---------  --------  ------------  ----------  --------  ------ │
#> │   Q1            128        77            22  15.8%             11       4 │
#> │   Q2            133        74            22  15.8%              7       0 │
#> │   Q3            112        70            19  13.7%              7       0 │
#> │   Q4            135        72            23  16.5%              4       0 │
#> │                                                                           │
#> │   period  tiedout  hotshr  jacpre                                         │
#> │   ------  -------  ------  ------                                         │
#> │   Q1            0  32.8%                                                  │
#> │   Q2            0  22.6%   0.2857                                         │
#> │   Q3            0  26.8%   0.1667                                         │
#> │   Q4           12  13.3%   0.1000                                         │
#> │                                                                           │
#> │   The busiest 7 cells hold 13.3% to 32.8% of a period's incidents.        │
#> │   Their overlap with the period before averages 0.184: the hot set is     │
#> │   mostly new each period, so these spots do not persist. 7 cells are      │
#> │   hot in at least half the periods.                                       │
#> │                                                                           │
#> │   shrforhalf is the share of active cells that hold half the incidents,   │
#> │   hotshr the share of incidents in the hot cells, and jacpre the          │
#> │   overlap of the hot cells with the period before (1 = the same cells).   │
#> │   Cells tied with the busiest 7 were left out of the hot set: the         │
#> │   counts are too even (or the cells too small) to single out hot spots.   │
#> │   period Q4: 12 cells tied                                                │
#> ╰───────────────────────────────────────────────────────────────────────────╯
plot(h)

The hot set is the busiest 7 cells and every cell tied with the last of them. When a tie would make more than twice that many cells hot, the counts do not single those cells out, so the tied group is left out and counted under tiedout: a quarter can then show fewer hot cells than 7, a smaller hot share, and an overlap with its neighbours that cannot be computed (left blank), as in Q4 here.

Pairs of break-ins within 100 metres and a week of each other occur 1.66 times as often as chance would give.

Indigenous life expectancy

lifetbl() works from counts of events in age intervals. Where only published summaries exist, as for Indigenous life expectancy, the comparison is made on the published estimates and their intervals instead:

le <- indiex("lifeexp")
le[le$period == "2016/2021", c("group", "estimate", "lower", "upper")]
#>                                                group estimate lower upper
#> 15                                             Métis    80.47 79.98 80.95
#> 16                                             Inuit    71.91 70.80 73.03
#> 17              First Nations (Registered or Treaty)    73.98 73.63 74.33
#> 18  First Nations (Registered or Treaty), on reserve    71.42 71.02 71.81
#> 19 First Nations (Registered or Treaty), off reserve    76.04 75.37 76.71
#> 20          First Nations, without registered status    79.23 78.42 80.04
#> 21                                    Non-Indigenous    84.99 84.93 85.04

References

Demers, S. and Rossmo, D. K. (2015). Simpson’s paradox in Canadian police clearance rates. Canadian Journal of Criminology and Criminal Justice, 57(3), 424-434.

Kitagawa, E. M. (1955). Components of a difference between two rates. Journal of the American Statistical Association, 50(272), 1168-1194.

Knox, E. G. (1964). The detection of space-time interactions. Journal of the Royal Statistical Society, Series C, 13(1), 25-30.

Parsing SIU director’s reports

Ontario’s Special Investigations Unit publishes a director’s report for every investigation as an HTML page. The reports share a structure but not a template: fields appear under varying headings, counts are written in words or figures, and the text is long. The package parses them into a fixed, panel-reviewed schema with a parser written for hostile input.

The schema

sch <- siusch()
head(sch$name, 20)
#>  [1] "police_service"                "date_of_incident_iso"         
#>  [3] "date_siu_notified_iso"         "date_of_director_decision_iso"
#>  [5] "siu_investigators"             "siu_forensics_investigators"  
#>  [7] "number_of_witness_officials"   "number_of_civilian_witnesses" 
#>  [9] "number_of_subject_officials"   "age_affected"                 
#> [11] "sex_gender_affected"           "charges_recommended"          
#> [13] "directors_name"                "location_of_call"             
#> [15] "specific_injuries"             "relevant_legislation"
nrow(sch)
#> [1] 16

Parsing a report

The package ships a synthetic report with the structure of a real one and no real person in it:

f <- siuprs(system.file("extdata", "siu_synthetic_report.html",
                                      package = "rmoriebricklayer"))
f[["number_of_subject_officials"]]
#> [1] "2"
f[["date_of_incident_iso"]]
#> [1] "2023-01-05"

The pieces

HTML to plain text, with the layout tags turned into line structure and entities decoded:

siutxt("<p>Number of SIU Investigators assigned: 3</p>")
#> [1] " Number of SIU Investigators assigned: 3\n"

The subject-official count, resolved from the enumerated list rather than a headline number that reports sometimes omit:

siuso("Subject Officials\nSO #1 Interviewed\nSO #2 Declined interview")
#> $count
#> [1] 2
#> 
#> $reason
#> [1] "section: max(ordinal 2, entries 2)"

Dates in the report’s prose, normalised to ISO:

siudate("March 4, 2024")
#> [1] "2024-03-04"
siudate("not a date")
#> [1] ""

Written for hostile input

A parser that reads 5,000 public web pages will meet broken HTML, megabyte-long lines, unclosed tags and entities of every shape. Three review rounds on this parser established and tested the following, each pinned by a test over the sources:

Fetching live reports

siufetch() fetches a report by its identifier through the package’s checked transport (address validation, redirect policy, resource limits) with an Internet Archive fallback, and siufprs() fetches and parses in one call:

f <- siufprs("24-OCI-123")

The corpus of parsed reports, reviewed round by round, is distributed as a dataset by the companion package rmoriedata; this package provides the parser so the corpus can be rebuilt and audited.

Indigenous police services

Ontario’s First Nations police services, among them the Nishnawbe-Aski Police Service, the Anishinabek Police Service and the Treaty Three Police Service, are investigated by the SIU like any other. The parsed police_service field selects their cases (this fetches live reports, so it is not run here):

fnservices <- c(
  "Nishnawbe-Aski Police Service",
  "Anishinabek Police Service",
  "Treaty Three Police Service",
  "UCCM Anishnaabe Police Service"
)
reports <- lapply(reporturls, siufprs)
mine <- Filter(
  function(r) any(grepl(paste(fnservices, collapse = "|"), r$police_service)),
  reports
)
length(mine)

The report’s timeline in a picture

siuprs() returns each date in the report in ISO form. Drawn on one axis, they give the investigation’s timeline at a glance.

rpt <- siuprs(siutxt(paste(
  readLines(
    system.file(
      "extdata",
      "siu_synthetic_report.html",
      package = "rmoriebricklayer"
    ),
    warn = FALSE
  ),
  collapse = "\n"
)))
iso <- rpt[grepl("_iso$", names(rpt)) & grepl("^\\d{4}-\\d{2}-\\d{2}$", rpt)]
when <- as.Date(iso)
o <- order(when)
lab <- gsub("_", " ", sub("_iso$", "", names(iso)[o]), fixed = TRUE)
span <- diff(range(when))
op <- par(mar = c(4, 1, 3, 1))
plot(
  when[o],
  seq_along(o),
  pch = 19,
  col = blpal()[1],
  yaxt = "n",
  xaxt = "n",
  xlab = "",
  ylab = "",
  xlim = range(when) + c(-0.1, 0.1) * as.numeric(span),
  ylim = c(0.5, length(o) + 0.5),
  main = "Dates in the synthetic report"
)
axis.Date(1, when, format = "%Y-%m-%d")
text(
  when[o],
  seq_along(o),
  paste0(lab, "  ", format(when[o])),
  cex = 0.75,
  pos = ifelse(seq_along(o) < length(o), 4, 2)
)

par(op)

Air pollution: dispersion, health burden, equity, and the footprint of the computation

Three questions about air pollution recur in administrative and public-health work: where does a plume go, how many deaths does a concentration imply, and who bears it. A fourth has become part of reproducible research: what did the computation itself cost in electricity and carbon. This vignette walks the four tools, each with the formula it implements and the paper it comes from, on small offline examples.

1. Where the plume goes

Dispersion coefficients

The Gaussian plume model needs the horizontal and vertical spread of the plume as a function of downwind distance and atmospheric stability. Pasquill’s stability classes A (very unstable) to F (very stable) with Briggs’s (1973) curve fits are the standard parameterisation; rural and urban coefficients differ because cities are rougher and warmer.

pgsig(1000, "D")
#> $sigmay
#> [1] 76.27701
#> 
#> $sigmaz
#> [1] 37.94733
pgsig(1000, "D", set = "urban")
#> $sigmay
#> [1] 135.2247
#> 
#> $sigmaz
#> [1] 122.7881

Plume rise

A hot, fast stack gas rises before it levels off. Briggs’s plume-rise equations give the buoyancy flux, the final rise and the distance at which it is reached; before that distance the plume is still rising:

r <- briggs(c(100, 500, 2000), u = 4, diam = 1, exitvel = 12,
                       stktmp = 420, ambtmp = 285)
r$flux
#> [1] 9.45594
r$finalrise
#> [1] 28.88281
r$rise
#> [1] 18.22335 28.88281 28.88281

Ground-level concentration

With the effective height (stack height plus rise), the plume equation gives the concentration at any receptor; ground reflection is included, and a mixing lid adds image terms. Concentration first rises with distance as the plume reaches the ground, then falls as it dilutes:

x <- seq(200, 5000, by = 100)
h <- 30 + r$finalrise
cd <- gplume(q = 100, u = 4, h = h, rcpt = cbind(x, 0, 0), stab = "D")
cb <- gplume(q = 100, u = 4, h = h, rcpt = cbind(x, 0, 0), stab = "B")
plot(
  x,
  cd,
  type = "l",
  xlab = "downwind distance (m)",
  ylab = "concentration (g/m^3)",
  main = "Ground-level centreline concentration"
)
lines(x, cb, col = 2)
legend(
  "topright",
  c("class D (neutral)", "class B (unstable)"),
  col = 1:2,
  lty = 1,
  bty = "n"
)

gspuff() is the instantaneous-release counterpart; advdif2d() and lagrprt() are grid and particle models for cases the Gaussian assumptions (steady wind, flat terrain) do not fit. The particle model is deterministic from its seed and takes any normal generator through normals=, so another package’s random stream can be reproduced exactly.

p <- lagrprt(
  n = 500,
  x0 = 0,
  y0 = 0,
  u = 1,
  v = 0,
  kx = 0.5,
  ky = 0.5,
  dt = 1,
  nsteps = 50,
  seed = 1,
  grid = c(0, 100, -30, 30, 20, 12)
)
c(mean_x = p$mean_x, meany = p$meany)
#>       mean_x        meany 
#> 49.981881150  0.009562424

2. From concentration to deaths

Concentration-response

The relative risk of death at an ambient concentration, relative to a counterfactual concentration below which no excess risk is assumed. The all-cause PM2.5 function is log-linear with the pooled cohort estimate of the WHO 2021 guideline review, RR 1.08 per 10 µg/m³ (Chen and Hoek 2020); NO2 uses RR 1.02 per 10 µg/m³ (Huangfu and Atkinson 2020). The cause-specific PM2.5 functions are the Integrated Exposure-Response curves of Burnett et al. (2014), which saturate at high concentrations:

crfpm25(12)$rr
#> [1] 1.048873
exp(log(1.08) * (12 - 5.8) / 10)          # the same, by hand
#> [1] 1.048873
crfpm25(12, y = "ihd")$rr
#> [1] 1.553634
crfno2(25)$rr
#> [1] 1.03015

A vector of exposures gives the mean RR and keeps the per-unit values:

v <- crfpm25(c(6, 10, 15, 22))
v$rr
#> [1] 1.060136
v$extra$rrperunit
#> [1] 1.001540 1.032852 1.073371 1.132782

Attributable fraction and avoided deaths

Levin’s population attributable fraction turns a relative risk and an exposure prevalence into the share of deaths the exposure accounts for; the BenMAP health-impact function gives the deaths avoided by a counterfactual reduction:

attrfrac(crfno2(25)$rr, expoprev = 1)
#> ╭── Population attributable fraction ─────────────────────────────────────────╮
#> │   Attributable fraction  2.93%                                              │
#> │   Relative risk          1.03                                               │
#> │   Exposed                100.00%                                            │
#> │   Formula                Levin: p (RR - 1) / (1 + p (RR - 1))               │
#> │                                                                             │
#> │   About 2.9% of the cases in this population would not occur without the    │
#> │   exposure, if the relative risk is causal.                                 │
#> │                                                                             │
#> │   It holds only if the relative risk is causal and unconfounded, and        │
#> │   only for a population exposed as this one is (Rothman et al. 2008,        │
#> │   chapter 5).                                                               │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
mortdisp(expdel = 10, pop = 1e6, brate = 0.008,
                    betaperunit = log(1.02) / 10)
#> [1] 156.8627

The burden chain

polbrdn() chains concentration-response, attributable fraction and the baseline case count, as the Global Burden of Disease risk-factor estimates do:

b <- polbrdn(expomean = 25, expoprev = 1, brate = 0.008,
                      pop = 1e6, poll = "NO2")
c(rr = b$extra$rr, paf = b$paf, attributable = b$attrcases)
#>           rr          paf attributable 
#>   1.03014950   0.02926711 234.13691783
b$citation
#> [1] "Huangfu & Atkinson (2020) Environ Int 144:105998; WHO (2021) Global AQ Guidelines; Rothman et al. (2008) section 5"

Per area, sorted worst first:

areas <- data.frame(
  fsa = c("M6H", "M5V", "M4W", "M1B"),
  exposure = c(28, 18, 15, 24),
  population = c(40000, 60000, 25000, 70000),
  brate = 0.008,
  stringsAsFactors = FALSE
)
polbrdarea(areas)
#> ╭── Pollution burden by area ─────────────────────────────────────────────╮
#> │   Areas               4                                                 │
#> │   Attributable cases  36.03                                             │
#> │   Baseline cases      1,560                                             │
#> │   Share attributable  2.31%                                             │
#> │                                                                         │
#> │   fsa  exposure  population     brate     rr      paf  attrcases        │
#> │   ---  --------  ----------  --------  -----  -------  ---------        │
#> │   M1B        24      70,000  0.008000  1.028  0.02734     15.312        │
#> │   M6H        28      40,000  0.008000  1.036  0.03502     11.205        │
#> │   M5V        18      60,000  0.008000  1.016  0.01572      7.544        │
#> │   M4W        15      25,000  0.008000  1.010  0.00985      1.970        │
#> │                                                                         │
#> │   fsa  basecases                                                        │
#> │   ---  ---------                                                        │
#> │   M1B        560                                                        │
#> │   M6H        320                                                        │
#> │   M5V        480                                                        │
#> │   M4W        200                                                        │
#> │                                                                         │
#> │   The areas are sorted by the cases attributable to the exposure,       │
#> │   largest first: rr is the relative risk at the area's exposure and paf │
#> │   the share of its cases the exposure accounts for.                     │
#> ╰─────────────────────────────────────────────────────────────────────────╯

3. Who bears it

The concentration index of Wagstaff, Paci and van Doorslaer (1991) measures whether exposure concentrates among lower- or higher-income people: it is twice the covariance of exposure with the fractional income rank, divided by mean exposure. Negative means the poor bear more.

d <- data.frame(exposure = c(30, 27, 22, 18, 14), income = 1:5)
e <- expconcidx(d, "exposure", "income")
e$concidx
#> [1] -0.1477477
e$interpretation
#> [1] "CI = -0.1477. Pro-poor exposure burden: lower-income individuals bear disproportionately higher pollution."

4. Estimating the exposure effect from your own data

When individual or area data with confounders are available, the partially linear model of Chernozhukov et al. (2018) estimates the exposure effect with cross-fitted nuisance functions, so the confounder adjustment cannot overfit the same rows it is evaluated on. exprspplr() adds a percentile bootstrap for a sensitivity interval. The default estimator uses ordinary-least-squares nuisances; any estimator with the same interface can be passed.

set.seed(1)
n <- 300
sim <- data.frame(income = rnorm(n), density = rnorm(n))
sim$no2 <- 20 - 2 * sim$income + 3 * sim$density + rnorm(n)
sim$asthmarate <- 5 + 0.3 * sim$no2 - 0.5 * sim$income + rnorm(n)
fit <- exprspplr(sim, y = "asthmarate", expo = "no2",
                             confs = c("income", "density"), nboot = 30)
c(ate = fit$ate, lower = fit$cilobs, upper = fit$ciupbs)
#>       ate     lower     upper 
#> 0.2464713 0.1437978 0.3132657

The true coefficient is 0.3 and the interval covers it.

5. The whole pipeline, with an assumption log

vrfpol() runs the chain end to end and records every assumption it checked, so a report says not only the number but the conditions under which it holds. Synthetic demo data stand in here:

r <- vrfpol("no2", demo = TRUE, reg = "demo", years = "2024")
r$status
#> [1] "ok"
cat(polrpttxt(r))
#> ==================================================================
#>   verify-pollution -- NO2 -> allcause
#>   region: demo
#>   years:  2024
#>   data:   demo (synthetic)
#> ==================================================================
#> 
#> Inputs
#>   exposure mean:      24.261
#>   exposure prevalence:0.981
#>   baseline rate/100k: 500.00
#>   population:         1,000,000
#>   reference conc:     10
#> 
#> Assumption log
#>   [PASS] exposure finite and non-negative -- mean 24.26113743 vs ref 10
#>   [PASS] exposure within the concentration-response range (at most 10,000 ug/m3) -- exposure_mean=24.26114
#>   [PASS] prevalence in [0,1] -- exposure_prevalence=0.981
#>   [PASS] brate finite: zero, or between 1 and 10,000 per 100k -- brate=500 per 100k per year
#>   [PASS] population finite, positive and at most 10 billion -- population=1,000,000
#>   [PASS] reference finite and non-negative -- reference=10
#>   [PASS] pollutant supported by the CRFs -- Current CRFs: NO2 (log-linear), PM2.5 (log-linear all-cause; Burnett IER for IHD and stroke). Other pollutants reject.
#> 
#> Concentration-response
#>   RR:       1.0286
#>   source:   Huangfu & Atkinson (2020) Environ Int 144:105998; WHO (2021) Global AQ Guidelines
#> 
#> Attributable fraction (PAF): 0.0273  (95% CI 0.0138 to 0.0534, from the risk's limits)
#> 
#> Mortality displaced
#>   expected avoided deaths: 136.6
#> 
#> Burden of pollution
#>   attributable deaths:   136.7
#> 
#> Equity analysis
#>   concentration index: -0.0647
#> 
#> STATUS: ok

An input that violates an assumption stops the pipeline and says why, with a non-zero exit status for scripts (a negative exposure here); an exposure at or below the reference is a valid run with zero attributable cases, and a call with no exposure data at all is an error, not a clean zero:

bad <- vrfpol("pm25", expomean = -3, expoprev = 0.5)
bad$status
#> [1] "assumption_failure"
attr(bad, "exitstat")
#> [1] 1
none <- vrfpol("pm25")
none$status
#> [1] "error"
attr(none, "exitstat")
#> [1] 2

A CSV with an exposure column (and optionally income) replaces the demo data; a NAPS pull with value in ppb is converted to µg/m³ at 1.88 for NO2.

6. What did this computation cost?

Computation uses electricity, and electricity has a carbon intensity. co2fp() reports both for any expression, two ways. On Linux with readable RAPL counters it measures the CPU package energy, as codecarbon does. Everywhere else it models energy with the Green Algorithms formula (Lannelongue, Grealey and Inouye 2021), \(E = t\,(P_{cpu} u + M \cdot 0.3725\,\mathrm{W/GB})\,PUE\). Every default that stood in for a measurement is listed:

raplok()
#> [1] FALSE
fp <- co2fp(sum(sqrt(seq_len(3e5))), loc = "CA-ON", cpupwrw = 45,
                        memgb = 16)
print(fp)
#> ╭── Computation footprint (modelled) ─────────────────────────────────────────╮
#> │   Time         0.004 s wall, 0.004 s CPU                                    │
#> │   Utilisation  utilisation 0.12 (process) on 8 cores                        │
#> │   Memory       16 GB                                                        │
#> │   Energy       1.29e-08 kWh (cpu 6.25e-09, memory 6.62e-09)                 │
#> │   CO2e         1.17e-06 g at 90.97 gCO2e/kWh (Canada, Ontario, 2024)        │
#> │   Source       Electricity Maps (2025). Canada 2024 Yearly Carbon           │
#> │                Intensity Data (Version January 27, 2025). Electricity Maps  │
#> │                (via Green Algorithms data v3.1)                             │
#> │                                                                             │
#> │   An estimate of the energy this computation used and the emissions of      │
#> │   making it on the local grid; fpeqv() puts the grams in everyday terms.    │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
fp$assumptions
#> character(0)

Utilisation: measured, per process or per machine

The Green Algorithms calculator asks the user for the utilisation \(u\) and defaults to 1; codecarbon measures it in one of two tracking modes. co2fp() measures it too, in the same two ways. usage = "process" (the default) is this R process’s own CPU time, children included, over wall time and cores, so an idle wait is not charged as computation. usage = "machine" is the whole machine’s busy share over the run, read from /proc/stat on Linux, for the case where the computation is all the machine is doing or the question is what the machine drew; where there is no /proc/stat the process figure is used and the assumptions say so. A number fixes it, as the calculator does. loadcur chooses how power follows utilisation: linearly (Green Algorithms) or along codecarbon’s curves, a 10% floor for a process and \(0.1 + 0.9u^3\) for a machine, so a figure can be made comparable with either tool:

busy <- function() sum(sqrt(seq_len(2e5)))
a <- co2fp(busy(), loc = "CA-ON", cpupwrw = 45, memgb = 16)
b <- co2fp(busy(), loc = "CA-ON", cpupwrw = 45, memgb = 16,
                       usage = "machine", loadcur = "codecarbon")
rbind(prolin = c(usage = a$usage, cpukwh = a$enekwh$cpu),
      maccod = c(usage = b$usage, cpukwh = b$enekwh$cpu))
#>        usage     cpukwh
#> prolin 0.125 4.6875e-09
#> maccod 0.000 5.0000e-09
b$usagemode
#> [1] "machine"

The grid: the latest year, every country, offline

The carbon intensity comes from a table bundled with the package, so nothing is fetched and it works in any country: every country at the latest year Our World in Data publishes from Ember’s yearly electricity data (2024 or 2025 at this release for all but a few grids whose latest row is older, Ukraine’s 2022 among them; every row carries its year), sub-national zones from Electricity Maps’ 2024 yearly data as redistributed by Green Algorithms (Canadian provinces, US balancing authorities, Australian states, Indian and Japanese regions, …), and the world average. Each row carries its year and source, so a report can cite them. A location is matched by code, alpha-3 code or name; an unknown zone falls back to its country and says so; a number of your own is accepted when you have a better figure:

co2int("CA-ON")[c("g_per_kwh", "year", "source")]
#> $g_per_kwh
#> [1] 90.97
#> 
#> $year
#> [1] 2024
#> 
#> $source
#> [1] "Electricity Maps (2025). Canada 2024 Yearly Carbon Intensity Data (Version January 27, 2025). Electricity Maps (via Green Algorithms data v3.1)"
co2int("Kenya")[c("location", "g_per_kwh", "year")]
#> $location
#> [1] "KE"
#> 
#> $g_per_kwh
#> [1] 95.44
#> 
#> $year
#> [1] 2025
co2int("CA-XX")$note
#> [1] "no row for zone CA-XX; its country CA used"
tab <- co2inttab()
nrow(tab)
#> [1] 378
table(nchar(tab$location) == 2, tab$year)[, c("2024", "2025")]
#>        
#>         2024 2025
#>   FALSE  155    3
#>   TRUE   125   91

With location = NULL (the default) the location is detected offline, in order, from the RMBL_LOCATION environment variable, codecarbon’s CODECARBON_COUNTRY_ISO_CODE, the system time zone (America/Toronto is Canada, Asia/Kolkata India) and the locale, and the detection is written into the assumptions; when nothing is found the world average is used and the result says so, so one country’s grid is never silently applied to another’s computation:

dtctloc()
#> $location
#> [1] "US"
#> 
#> $method
#> [1] "system time zone America/New_York"
dtctloc(envloc = "", codloc = "", tz = "Asia/Kolkata")
#> $location
#> [1] "IN"
#> 
#> $method
#> [1] "system time zone Asia/Kolkata"
co2fp(1 + 1, cpupwrw = 45, memgb = 16)$assumptions[1]
#> [1] "location US via system time zone America/New_York"

For a long analysis the figures become meaningful; the equivalents are the two the Green Algorithms calculator prints:

fpeqv(fp$co2eg * 1e4)   # if this ran ten thousand times
#> ╭── Footprint in everyday terms ────────────────────────────────────────╮
#> │   Car travel   6.691e-05 km                                           │
#> │   Tree uptake  1.277e-05 tree-months                                  │
#> │                                                                       │
#> │   The same emissions as driving a European passenger car this far, or │
#> │   what this many trees take up in a month.                            │
#> │                                                                       │
#> │   car: 175 gCO2/km, European passenger car (Green Algorithms          │
#> │   context.csv). tree: 917 gCO2 per tree-month (Green Algorithms       │
#> │   context.csv).                                                       │
#> ╰───────────────────────────────────────────────────────────────────────╯

Cross-checked against codecarbon. On the same 8-second CPU-bound workload with the same forced inputs (65 W CPU, 16 GB at 0.3725 W/GB), codecarbon 3.3.1 and co2fp() charged the memory identically (5.96 W) and the CPU within 8%; loadcur = "codecarbon" reproduces codecarbon’s power-from-load curves, so what remains of the difference is the load measurement itself (codecarbon samples it during the run, this function integrates the CPU time). The carbon intensity they applied to Canada differed (170 g/kWh from codecarbon’s 2023 Our World in Data figure; 190.7 g/kWh from the 2025 figure bundled here), which is the point of carrying the year with every row: grid intensities move, and a report should cite the year of the one it used. CodeCarbonR on CRAN wraps the same Python tracker through reticulate; this function needs no Python and no network.

What this does and does not cover. The measured mode reads the CPU packages only; memory, storage, network and GPU are not in the RAPL counters (GPU power needs the vendor library; this package has no GPU dependency). The modelled mode inherits Green Algorithms’ assumptions: memory draws 0.3725 W/GB regardless of how it is populated (codecarbon v3 moved to a per-DIMM heuristic for the same reason), and the grid intensity is an annual average, not the hour’s. Embodied emissions of the hardware are excluded, as in both reference tools. The place to improve a figure is the inputs: a measured cpupwrw for your CPU, the actual memgb your job used, the data centre’s pue, and, if you have it, the hour’s carbon intensity passed as a number.

First Nations, Métis and Inuit populations

The 2021 Census gives the populations a burden estimate needs. The bundled Indigenous Population Profile excerpt has them by province and territory; the concentrations below are scenarios (the WHO 2021 guideline of 5 µg/m³ against 8 and 12 µg/m³), not measurements of any community.

cp <- indiex("censuspop")
tot <- cp[
  cp$GEO_LEVEL != "Country" &
    cp$AGE_GROUP == "Total - Age groups" &
    cp$RESIDENCE == "Total - Residence on or off reserve" &
    grepl("^Total - Average age", cp$CHARACTERISTIC_NAME),
]
pop <- stats::setNames(tot$FIRST_NATIONS, tot$GEO_NAME)
head(pop[order(-pop)], 5)
#>          Ontario British Columbia          Alberta         Manitoba 
#>           251030           180085           145640           134890 
#>     Saskatchewan 
#>           121175
# relative risk of all-cause mortality from PM2.5, and the share of deaths
# attributable to the excess above the guideline if the whole population
# breathes the scenario level
vapply(c(8, 12), function(conc) {
  rr <- crfpm25(conc, refconc = 5)$rr
  c(pm25 = conc, rr = rr, attributable = attrfrac(rr, 1))
}, numeric(3))
#>                    [,1]       [,2]
#> pm25         8.00000000 12.0000000
#> rr           1.02335691  1.0553503
#> attributable 0.02282382  0.0524473

References

Briggs (1973). Diffusion estimation for small emissions. ATDL 79. Briggs (1975). Plume rise predictions. AMS. Burnett et al. (2014). An integrated risk function for estimating the global burden of disease attributable to ambient fine particulate matter exposure. EHP 122(4). Chen, Hoek (2020). Long-term exposure to PM and all-cause and cause-specific mortality. Environment International 143. Chernozhukov et al. (2018). Double/debiased machine learning. Econometrics Journal 21(1). Huangfu, Atkinson (2020). Long-term exposure to NO2 and O3 and all-cause and respiratory mortality. Environment International 144. Lannelongue, Grealey, Inouye (2021). Green Algorithms: quantifying the carbon footprint of computation. Advanced Science 8(12), 2100707. Ember (2026). Yearly Electricity Data, via Our World in Data, ourworldindata.org/grapher/carbon-intensity-electricity. Electricity Maps (2025). 2024 Yearly Carbon Intensity Data, via Green Algorithms data v3.1. Levin (1953). The occurrence of lung cancer in man. Acta Unio Internationalis Contra Cancrum 9. Wagstaff, Paci, van Doorslaer (1991). On the measurement of inequalities in health. Social Science and Medicine 33(5). WHO (2021). Global Air Quality Guidelines.

First Nations, Métis and Inuit health: published tables, distinctions and governance

This vignette covers the tools for First Nations, Métis and Inuit health, mental health and wellness data, as Statistics Canada publishes them. It also covers the rules those tools follow, because with these data, how you work matters as much as what you compute.

Every example runs on published values. indiex() bundles excerpts of eight Statistics Canada tables, reproduced under the Statistics Canada Open Licence:

Excerpt Table Survey and period
"health" 13-10-0924-01 Canadian Community Health Survey, 2015/2018 and 2019/2022
"mentalhealth" 41-10-0067-01 Indigenous Peoples Survey, 2022
"wellbeing" 41-10-0082-01 2024
"suiide" 41-10-0021-01 Aboriginal Peoples Survey, 2017
"lifeexp" 17-10-0160-01 Projected life expectancy at birth
"infmort" 17-10-0161-01 Infant mortality
"censuspop" 98-510-X2021001 2021 Census Indigenous Population Profile, population
"censushh" 98-510-X2021016 2021 Census Indigenous Population Profile, households

The 13-10-0924 excerpt is the table behind The Daily of 2025-02-13. The full tables, and the other Indigenous-titled tables, are at the data service: see indtabs().

Some of this vignette concerns suicide. If you or someone you know needs support, the Hope for Wellness Help Line (1-855-242-3310, or the online chat at hopeforwellness.ca) offers counselling and crisis support to Indigenous people across Canada, in English and French, and in Cree, Ojibway and Inuktitut on request. You can also call or text 9-8-8.

Whose data these are

Two sets of principles govern Indigenous data in Canada.

Walter and Andersen (2013) set out what a quantitative methodology looks like when it starts from Indigenous standpoints. Walter (2016) names the pattern it pushes against: statistics about Indigenous peoples clustered around five Ds, namely disparity, deprivation, disadvantage, dysfunction and difference.

The package follows these principles in six ways.

Reading a table as published

A Statistics Canada table puts each estimate and its limits on separate rows. scests() puts them side by side. It reads every layout the Indigenous tables use:

d <- data.frame(
  REF_DATE = "2022",
  GEO = "Canada",
  `Indigenous identity` = rep(c("First Nations", "M\u00e9tis"), each = 4),
  Statistics = rep(
    c(
      "Number of persons",
      "Percent",
      "Low 95% confidence interval",
      "High 95% confidence interval"
    ),
    2
  ),
  UOM = rep(c("Persons", "Percent", "Percent", "Percent"), 2),
  VALUE = c(415000, 47.1, 45.2, 49.0, 320000, 50.3, 48.1, 52.5),
  STATUS = "",
  check.names = FALSE,
  stringsAsFactors = FALSE
)
scests(d)
#>   REF_DATE    GEO Indigenous.identity           measure estimate lower upper
#> 1     2022 Canada       First Nations Number of persons 415000.0    NA    NA
#> 2     2022 Canada       First Nations           Percent     47.1  45.2  49.0
#> 3     2022 Canada               Métis Number of persons 320000.0    NA    NA
#> 4     2022 Canada               Métis           Percent     50.3  48.1  52.5
#>       cisou quality
#> 1      <NA>      ok
#> 2 published      ok
#> 3      <NA>      ok
#> 4 published      ok

A mirrored table loads by its id. This needs the data service key, so it is not run here:

indtabs("mental health")
e <- scests("41100067")

Perceived mental health, against non-Indigenous adults

The share of adults who rate their mental health very good or excellent, 2019/2022, from table 13-10-0924-01:

h <- indiex("health")
mh <- h[h$period == "2019/2022" & h$age == "Total, 18 years and over" &
          h$indicator == "Perceived mental health, very good or excellent", ]
g <- indgap(mh, "group", "estimate", "lower", "upper", qual = "quality")
g
#> ╭── Gap against a reference group ──────────────────────────────────────────╮
#> │   Reference  Non-Indigenous                                               │
#> │   Groups     3                                                            │
#> │   Intervals  95%                                                          │
#> │                                                                           │
#> │   group          estimate  refest   ratio  ratio CI        difference     │
#> │   -------------  --------  ------  ------  --------------  ----------     │
#> │   First Nations     45.80   61.50  0.7447  [0.700, 0.792]      -15.70     │
#> │   Métis             52.80   61.50  0.8585  [0.813, 0.907]       -8.70     │
#> │   Inuit             47.80   61.50  0.7772  [0.679, 0.890]      -13.70     │
#> │                                                                           │
#> │   group          difference CI                                            │
#> │   -------------  --------------                                           │
#> │   First Nations  [-18.5, -12.9]                                           │
#> │   Métis          [-11.6, -5.8]                                            │
#> │   Inuit          [-20.2, -7.2]                                            │
#> │                                                                           │
#> │   Below the reference (interval below 1): First Nations, Métis and Inuit. │
#> │                                                                           │
#> │   ratio: the group's estimate over the reference's; difference: the       │
#> │   group's less the reference's, in the estimate's own units.              │
#> │   Report a gap beside each group's own level, not alone.                  │
#> ╰───────────────────────────────────────────────────────────────────────────╯
plot(g)

Before pooling the three groups into one “Indigenous” figure, test whether they differ:

dt <- dstntst(
  mh[mh$group != "Non-Indigenous", ],
  "group",
  "estimate",
  "lower",
  "upper"
)
dt
#> ╭── Distinctions test ───────────────────────────────────────────────────╮
#> │   Groups       3                                                       │
#> │   Cochran's Q  11.97 on 2 df                                           │
#> │   p-value      0.003                                                   │
#> │   I-squared    83% of the variation is beyond chance                   │
#> │   Pooled       49.11 (se 0.973, inverse-variance mean)                 │
#> │                                                                        │
#> │   group          estimate     se  95% CI                               │
#> │   -------------  --------  -----  --------------                       │
#> │   First Nations     45.80  1.429  [43.00, 48.60]                       │
#> │   Métis             52.80  1.454  [49.95, 55.65]                       │
#> │   Inuit             47.80  3.291  [41.35, 54.25]                       │
#> │                                                                        │
#> │   Q says the groups differ: report each one, not the pooled figure.    │
#> │   Q = 11.97, p = 0.00252.                                              │
#> │                                                                        │
#> │   Distinctions-based: each group is reported on its own; a             │
#> │   pan-Indigenous total is left out of the test, never used as a group. │
#> │   Higgins and Thompson (2002) for Q and I-squared.                     │
#> ╰────────────────────────────────────────────────────────────────────────╯
plot(dt)

The whole scale, not one cut of it

“Very good or excellent” cuts a five-point scale at one place (the table publishes it in four categories, with fair and poor together). ordcmp() compares the whole distribution. It gives the probability that a person in one group rates their mental health higher than a person in the reference group, with ties counted half (Agresti 1980). The data are from the 2022 Indigenous Peoples Survey, table 41-10-0067-01:

r <- indiex("mentalhealth")
r <- r[r$region == "Canada", ]
oc <- ordcmp(
  r,
  "group",
  "rating",
  "estimate",
  "lower",
  "upper",
  lvls = c("Fair or poor", "Good", "Very good", "Excellent"),
  ref = "M\u00e9tis",
  qual = "quality"
)
oc
#> ╭── Ordinal comparison against a reference group ────────────────────────╮
#> │   Reference       Métis                                                │
#> │   Reference neff  3,884                                                │
#> │   Levels          4                                                    │
#> │   Measure         chance a person in the group rates higher, ties half │
#> │                                                                        │
#> │   group          probbet  95% CI             delta   neff              │
#> │   -------------  -------  --------------  --------  -----              │
#> │   First Nations   0.4653  [0.453, 0.477]  -0.06940  4,021              │
#> │   Inuit           0.5062  [0.493, 0.520]   0.01247  2,792              │
#> │                                                                        │
#> │   Rate lower than Métis (interval below 0.5): First Nations.           │
#> │   Not distinguishable from 0.5 (interval includes it): Inuit.          │
#> │                                                                        │
#> │   probbet 0.5 is no difference; above 0.5 the group rates higher.      │
#> │   delta (Cliff's) = share rating higher less share rating lower.       │
#> │   neff: the effective sample size the published limits imply (it       │
#> │   carries the survey's design effect).                                 │
#> ╰────────────────────────────────────────────────────────────────────────╯
plot(oc)

Age standardisation

First Nations, Métis and Inuit populations are younger than the non-Indigenous population, and many outcomes vary steeply with age. stdprev() standardises survey percentages to one age structure. Here it is applied to serious thoughts of suicide in the past year, from the 2017 Aboriginal Peoples Survey (table 41-10-0021-01):

si <- indiex("suicide")
si <- si[si$measure == "Seriously contemplated suicide in the past 12 months" &
           si$age != "Total, 15 years and over", ]
si <- si[order(si$group, si$age), ]
sp <- stdprev(si$estimate, si$lower, si$upper, std = c(20, 18, 17, 17, 28),
                             grp = si$group, qual = si$quality)
sp
#> ╭── Age-standardised prevalence ──────────────────────────────────────────────╮
#> │   Groups  3                                                                 │
#> │   Method  direct standardisation to the standard population given           │
#> │                                                                             │
#> │   group            estimate  95% CI              se                         │
#> │   ---------------  --------  --------------  ------                         │
#> │   First Nations *     5.792  [5.141, 6.443]  0.3321                         │
#> │   Inuit *             5.940  [5.036, 6.844]  0.4612                         │
#> │   Métis *             4.742  [4.125, 5.359]  0.3148                         │
#> │                                                                             │
#> │   Highest prevalence: Inuit * (5.94); lowest: Métis * (4.742).              │
#> │   Their intervals overlap, which alone does not mean they are the same:     │
#> │   the difference has its own, narrower interval (svycmp() for survey        │
#> │   domains).                                                                 │
#> │                                                                             │
#> │   Each group's prevalence as if it had the standard population's ages,      │
#> │   so the groups compare without their age structures.                       │
#> │   Compare only rates standardised to the same population.                   │
#> │   * use with caution: the source flags the estimate (E, high sampling       │
#> │   variability).                                                             │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
plot(sp)

For counts of events (deaths, hospitalisations), agestd() gives directly standardised rates with the Fay-Feuer gamma interval, and rategap() compares two of them.

Wellness beside distress

Indigenous scholars of mental health describe it as wellness, and as a process of healing, grounded in community, cultural identity, a holistic view and interdependence (Stewart 2008). They ask that strengths be reported together with disparities (Stewart 2020).

The estimators in this section and the next put that into how published tables are analysed. They are our translation, inspired by that work, not methods its authors state, and they say what they cannot show.

welprof() puts wellness indicators first and distress after them. This example uses the 2024 well-being table, 41-10-0082-01:

wb <- indiex("wellbeing")
w <- welprof(
  wb,
  "group",
  "indicator",
  "estimate",
  "lower",
  "upper",
  well = c(
    "Life satisfaction rating of 8, 9 or 10",
    "Sense of meaning and purpose rating of 8, 9 or 10",
    "Always or often has a hopeful outlook"
  ),
  distress = c(
    "Life satisfaction rating between 0 and 5",
    "Rarely or never has a hopeful outlook"
  ),
  qual = "quality"
)
plot(w)

frmaud(w)
#> ╭── Framing audit ────────────────────────────────────────────────────╮
#> │   Groups               3                                            │
#> │   Wellness indicators  9                                            │
#> │   Distress indicators  6                                            │
#> │   Deficit-framed       40%                                          │
#> │                                                                     │
#> │   By group (shrdef: the deficit-framed share)                       │
#> │   group          wellness  distress  shrdef                         │
#> │   -------------  --------  --------  ------                         │
#> │   First Nations         3         2  0.4000                         │
#> │   Métis                 3         2  0.4000                         │
#> │   Inuit                 3         2  0.4000                         │
#> │                                                                     │
#> │   40% of the indicators reported are deficit-framed.                │
#> │                                                                     │
#> │   Walter (2016): statistics about Indigenous peoples cluster around │
#> │   disparity, deprivation, disadvantage, dysfunction and difference. │
#> │   This audits the report, not the people it describes.              │
#> ╰─────────────────────────────────────────────────────────────────────╯

dualcont() places each group on two axes, wellness and distress. They are related but distinct (Keyes 2002, 2005):

dc <- dualcont(
  w,
  "Always or often has a hopeful outlook",
  "Life satisfaction rating between 0 and 5"
)
dc
#> ╭── Dual continuum: wellness beside distress ────────────────────────────────╮
#> │   Wellness  Always or often has a hopeful outlook                          │
#> │   Distress  Life satisfaction rating between 0 and 5                       │
#> │   Groups    3                                                              │
#> │   Ratio     wellness over distress, log-normal interval                    │
#> │                                                                            │
#> │   group          wellness  interval        distress  interval        ratio │
#> │   -------------  --------  --------------  --------  --------------  ----- │
#> │   First Nations     50.20  [45.90, 54.50]     31.40  [27.40, 35.70]  1.599 │
#> │   Métis             52.30  [47.60, 57.00]     28.40  [24.00, 33.10]  1.842 │
#> │   Inuit             54.40  [45.70, 62.80]     25.20  [19.40, 31.80]  2.159 │
#> │                                                                            │
#> │   group          ratio interval                                            │
#> │   -------------  --------------                                            │
#> │   First Nations  [1.37, 1.87]                                              │
#> │   Métis          [1.53, 2.21]                                              │
#> │   Inuit          [1.61, 2.90]                                              │
#> │                                                                            │
#> │   A ratio above 1: more people report the wellness than the distress.      │
#> │   3 of 3 groups have the whole interval above 1.                           │
#> │                                                                            │
#> │   Both estimates come from the same respondents and their covariance       │
#> │   is not published: the interval treats them as independent.               │
#> ╰────────────────────────────────────────────────────────────────────────────╯
plot(dc)

Holistic balance

Wellness is a balance among the domains of a person; illness is a loss of that balance (Stewart 2008; Gabriel, Mirza and Stewart 2022). An average of domains cannot show balance. holbal() measures it as the ratio of the geometric mean of the domain levels to their arithmetic mean. That ratio is 1 when every domain is level, and falls as one domain falls behind the others.

Which indicator stands for each domain is a choice for the people whose data these are to make. The choice below is one reading of it:

dom <- c(
  "Perceived health, very good or excellent" = "physical",
  "Perceived mental health, very good or excellent" = "mental",
  "Life satisfaction, satisfied or very satisfied" = "emotional",
  "community"
)
names(dom)[4] <- paste(
  "Sense of belonging to local community,",
  "somewhat strong or very strong"
)
d <- h[
  h$period == "2019/2022" &
    h$age == "Total, 18 years and over" &
    h$indicator %in% names(dom),
]
d$domain <- dom[d$indicator]
hb <- holbal(
  d,
  "group",
  "domain",
  "estimate",
  "lower",
  "upper",
  qual = "quality"
)
hb
#> ╭── Holistic balance across the domains ──────────────────────────────────────╮
#> │   Groups   4                                                                │
#> │   Weights  equal, over 4 domains                                            │
#> │   Balance  weighted geometric over arithmetic mean of the domain levels     │
#> │                                                                             │
#> │   group           level  95% level CI    balance  95% balance CI  weakest   │
#> │   --------------  -----  --------------  -------  --------------  --------  │
#> │   First Nations   59.52  [58.25, 60.80]   0.9641  [0.958, 0.970]  physical  │
#> │   Métis           63.73  [62.44, 65.01]   0.9752  [0.971, 0.979]  physical  │
#> │   Inuit           63.80  [60.74, 66.86]   0.9700  [0.956, 0.981]  mental    │
#> │   Non-Indigenous  69.80  [69.60, 70.00]   0.9847  [0.984, 0.985]  physical  │
#> │                                                                             │
#> │   Domain levels (percent)                                                   │
#> │   group           physical  mental  community  emotional                    │
#> │   --------------  --------  ------  ---------  ---------                    │
#> │   First Nations      43.90   45.80      63.70      84.70                    │
#> │   Métis              49.70   52.80      65.10      87.30                    │
#> │   Inuit              49.40   47.80      73.50      84.50                    │
#> │   Non-Indigenous     58.90   61.50      67.40      91.40                    │
#> │                                                                             │
#> │   Highest balance: Non-Indigenous (0.9847); lowest: First Nations (0.9641). │
#> │   Their intervals do not overlap: the gap is more than sampling error.      │
#> │   The weakest domain is physical for 3 of 4 groups.                         │
#> │                                                                             │
#> │   balance is 1 when every domain is at the same level and falls toward      │
#> │   0 as one falls behind; level is the weighted mean, reported beside it.    │
#> │   Our method, inspired by Stewart (2008), not taken from it.                │
#> │   The weights say how much each domain counts; a community may set          │
#> │   its own (CARE: authority to control, Carroll et al. 2020).                │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
plot(hb)

A community that weighs the domains differently sets weights. The result says which weights it used.

The four directions

The First Nations Mental Wellness Continuum Framework (Health Canada and Assembly of First Nations 2015) describes mental wellness in four directions: hope, belonging, meaning and purpose. mhcont() lays out the published proxies for each direction. It says which directions the tables do not measure, and does not fill them in:

mw <- mhcont(
  wb,
  "group",
  "indicator",
  "estimate",
  "lower",
  "upper",
  hope = "Always or often has a hopeful outlook",
  meaning = "Sense of meaning and purpose rating of 8, 9 or 10",
  purp = "Sense of meaning and purpose rating of 8, 9 or 10",
  qual = "quality"
)
mw[, c("group", "indicator", "domain", "estimate", "shared")]
#>            group indicator                                            domain
#> 1  First Nations      hope             Always or often has a hopeful outlook
#> 2  First Nations belonging                                no published proxy
#> 3  First Nations   meaning Sense of meaning and purpose rating of 8, 9 or 10
#> 4  First Nations   purpose Sense of meaning and purpose rating of 8, 9 or 10
#> 5          Métis      hope             Always or often has a hopeful outlook
#> 6          Métis belonging                                no published proxy
#> 7          Métis   meaning Sense of meaning and purpose rating of 8, 9 or 10
#> 8          Métis   purpose Sense of meaning and purpose rating of 8, 9 or 10
#> 9          Inuit      hope             Always or often has a hopeful outlook
#> 10         Inuit belonging                                no published proxy
#> 11         Inuit   meaning Sense of meaning and purpose rating of 8, 9 or 10
#> 12         Inuit   purpose Sense of meaning and purpose rating of 8, 9 or 10
#>    estimate shared
#> 1      50.2  FALSE
#> 2        NA  FALSE
#> 3      47.5   TRUE
#> 4      47.5   TRUE
#> 5      52.3  FALSE
#> 6        NA  FALSE
#> 7      51.6   TRUE
#> 8      51.6   TRUE
#> 9      54.4  FALSE
#> 10       NA  FALSE
#> 11     50.5   TRUE
#> 12     50.5   TRUE
plot(mw)

Healing as a process

A process shows in change. healchg() gives the change between survey periods, wellness first. It uses the two periods of table 13-10-0924-01:

hc <- healchg(
  h[h$age == "Total, 18 years and over", ],
  "group",
  "indicator",
  "period",
  "estimate",
  "lower",
  "upper",
  from = "2015/2018",
  to = "2019/2022",
  qual = "quality",
  well = c(
    "Perceived mental health, very good or excellent",
    "Sense of belonging to local community, somewhat strong or very strong"
  )
)
hc <- hc[
  hc$indicator %in%
    c(
      "Perceived mental health, very good or excellent",
      "Sense of belonging to local community, somewhat strong or very strong",
      "Mood disorder",
      "Anxiety disorder"
    ),
]
plot(hc)

Between the two periods, which span the COVID-19 pandemic, the share rating their mental health very good or excellent fell for First Nations, Métis and non-Indigenous adults. The Inuit interval includes no change. Because the reference group fell too, a group’s change is read against the reference group’s change, not on its own.

Connection across places

The land is “an instrumental part of all aspects of Indigenous life including health and well-being” (Mashford-Pringle and Stewart 2019). Health and education are interconnected (Mashford-Pringle 2016). Across places, connections of this kind leave a trace: do two indicators rise together where people live?

connassoc() measures this with Spearman’s rho. Its p-value is exact for up to eight places, and a Monte Carlo estimate beyond. Its band is a sensitivity band, not a confidence interval: it shows how far the published sampling error alone moves rho for these places. Here, by region, it asks whether First Nations and Métis people’s mental health rise together, using table 41-10-0067-01:

r <- indiex("mentalhealth")
r <- r[r$rating == "Excellent or very good" & r$region != "Canada", ]
pair <- merge(
  r[r$group == "First Nations", ],
  r[r$group == "M\u00e9tis", ],
  by = "region",
  suffixes = c("fn", "me")
)
ca <- connassoc(
  pair,
  "region",
  "estimatefn",
  "estimateme",
  "lowerfn",
  "upperfn",
  "lowerme",
  "upperme"
)
ca
#> ╭── Connection across places ──────────────────────────────────────────────╮
#> │   Places            11                                                   │
#> │   Indicators        estimatefn and estimateme                            │
#> │   Rank correlation  rho = 0.52                                           │
#> │   p-value           0.106 (permutation)                                  │
#> │   Interval          95% [-0.19, 0.77], with the published sampling error │
#> │                                                                          │
#> │   Rank correlation across 11 places: rho = 0.52, permutation p = 0.106   │
#> │   Across these places the two rank moderately together.                  │
#> │                                                                          │
#> │   An association across places only:                                     │
#> │   it does not show that one causes the other for any person.             │
#> ╰──────────────────────────────────────────────────────────────────────────╯
plot(
  ca,
  xlab = "First Nations: excellent or very good (%)",
  ylab = "M\u00e9tis: excellent or very good (%)"
)

This is an association across places. It cannot say that one thing makes a person well.

Community-set weights, and the Community Well-Being index

commidx() combines protective and risk factors with weights the community sets. It follows the structure of Inuit Tapiriit Kanatami’s National Inuit Suicide Prevention Strategy (2016). The function never chooses the weights. A factor the tables do not publish leaves the index NA; it is not filled in.

h2 <- h[h$period == "2019/2022" & h$age == "Total, 18 years and over", ]
ci <- commidx(
  h2,
  "group",
  "indicator",
  "estimate",
  "lower",
  "upper",
  qual = "quality",
  w = c(
    "Sense of belonging to local community, somewhat strong or very strong" = 2,
    "Has a regular healthcare provider" = 1,
    "Perceived life stress, most days quite a bit or extremely stressful" = -2,
    "Heavy drinking" = -1
  )
)
plot(ci)

cwbscr() computes Indigenous Services Canada’s Community Well-Being index from census profile figures: income, education, housing and labour force activity. The index measures socio-economic conditions only. It is not a measure of wellness as First Nations and Inuit define it.

The 2021 Census

The Indigenous Population Profile of the 2021 Census gives counts for every geography, by gender, age group and about two thousand characteristics, with First Nations, Métis, Inuit and non-Indigenous people in their own columns. In the provincial table, residence on and off reserve is split out too. cnsprof() reads it with the hierarchy of the characteristics. profshr() takes a share with an interval for the 25% long-form sample and the random rounding of the counts.

Foster children among children under 15, Canada, 2021:

p <- cnsprof(indiex("censuspop"))
fc <- profshr(p, "Foster children", "Total - Children in census families")
fc <- fc[
  fc$geo == "Canada" &
    fc$age == "0 to 14 years" &
    fc$group %in% c("First Nations", "M\u00e9tis", "Inuit", "Non-Indigenous"),
]
fc[, c(
  "residence",
  "group",
  "numerator",
  "denominator",
  "estimate",
  "lower",
  "upper"
)]
#>                              residence          group numerator denominator
#> 1  Total - Residence on or off reserve  First Nations     11735      265405
#> 2                           On reserve  First Nations      2350       90270
#> 3                          Off reserve  First Nations      9385      175135
#> 4  Total - Residence on or off reserve          Métis      1975      134670
#> 5                           On reserve          Métis        15         610
#> 6                          Off reserve          Métis      1960      134060
#> 7  Total - Residence on or off reserve          Inuit       770       21525
#> 8                           On reserve          Inuit        35         100
#> 9                          Off reserve          Inuit       735       21425
#> 10 Total - Residence on or off reserve Non-Indigenous     12705     5502310
#> 11                          On reserve Non-Indigenous        20        4070
#> 12                         Off reserve Non-Indigenous     12690     5498235
#>      estimate      lower      upper
#> 1   4.4215444  4.2880409  4.5590064
#> 2   2.6033012  2.5989608  2.6076487
#> 3   5.3587233  5.1789488  5.5443733
#> 4   1.4665479  1.3593905  1.5820167
#> 5   2.4590164  1.8933392  3.1882099
#> 6   1.4620319  1.3548094  1.5776046
#> 7   3.5772358  3.1713963  4.0328470
#> 8  35.0000000 30.9738945 39.2518422
#> 9   3.4305718  3.0324927  3.8788168
#> 10  0.2309030  0.2240599  0.2379546
#> 11  0.4914005  0.4040293  0.5975523
#> 12  0.2308013  0.2239573  0.2378541
plot(fc[fc$residence == "Total - Residence on or off reserve", ])

The overrepresentation of First Nations, Métis and Inuit children in care is the subject of the first five Calls to Action of the Truth and Reconciliation Commission.

Knowledge of an Indigenous language (all ages):

know <- paste(
  "Total - Knowledge of Indigenous languages for the population in private",
  "households - 25% sample data"
)
la <- profshr(p, "Indigenous languages", know)
la[
  la$geo == "Canada" &
    la$age == "Total - Age groups" &
    la$residence == "Total - Residence on or off reserve" &
    la$group %in% c("First Nations", "M\u00e9tis", "Inuit"),
  c("group", "estimate", "lower", "upper")
]
#>           group  estimate    lower     upper
#> 1 First Nations 17.530439 17.40473 17.656859
#> 2         Métis  1.332062  1.28369  1.382232
#> 3         Inuit 58.130139 57.49824 58.759377

The census age structure of the Indigenous population is the standard to use when age-standardising First Nations, Métis and Inuit rates (Kumar and Tjepkema 2019):

cnsstd(p)
#>                  age population     weight
#> 1       0 to 4 years     137510 0.07608819
#> 2       5 to 9 years     155565 0.08607853
#> 3     10 to 14 years     166135 0.09192722
#> 4     15 to 19 years     148820 0.08234633
#> 5     20 to 24 years     136070 0.07529140
#> 6     25 to 34 years     263655 0.14588780
#> 7     35 to 44 years     217990 0.12062006
#> 8     45 to 54 years     206360 0.11418485
#> 9     55 to 64 years     202785 0.11220670
#> 10    65 to 74 years     119410 0.06607295
#> 11 75 years and over      52945 0.02929597

A Community Well-Being score from the public profile (income, education, housing, labour force), by distinction. The household table splits Indigenous from non-Indigenous households only, so First Nations, Métis and Inuit people share the Indigenous households’ housing component. ?cnscwb explains where it departs from the Department’s own index:

h <- cnsprof(indiex("censushh"))
cw <- cnscwb(p, h)
cw[
  cw$geo %in% c("Canada", "Nunavut"),
  c("geo", "group", "income", "education", "housing", "labour", "cwb")
]
#> ╭── Community Well-Being index, from the census ──────────────────────────────╮
#> │   Rows         8                                                            │
#> │   Geographies  Canada; Nunavut                                              │
#> │   Groups       First Nations; Métis; Inuit; Non-Indigenous                  │
#> │   Not scored   2                                                            │
#> │   Scale        0 to 100 for each component and the score                    │
#> │                                                                             │
#> │   geo      group           income  education  housing  labour    cwb        │
#> │   -------  --------------  ------  ---------  -------  ------  -----        │
#> │   Canada   First Nations    71.03      51.62    91.88   75.85  72.59        │
#> │   Canada   Métis            78.73      59.56    91.88   83.15  78.33        │
#> │   Canada   Inuit            72.13      38.51    91.88   74.75  69.32        │
#> │   Canada   Non-Indigenous   84.38      68.01    95.99   86.30  83.67        │
#> │   Nunavut  First Nations                                                    │
#> │   Nunavut  Métis                                                            │
#> │   Nunavut  Inuit            70.22      28.49    70.82   69.85  59.85        │
#> │   Nunavut  Non-Indigenous  100.00      79.89    93.30   94.85  92.01        │
#> │                                                                             │
#> │   Highest score: Nunavut, Non-Indigenous (92.01); lowest: Nunavut, Inuit    │
#> │   (59.85).                                                                  │
#> │   Within a geography the groups' scores are furthest apart in Nunavut       │
#> │   (32.2 points) and closest in Canada (14.4).                               │
#> │                                                                             │
#> │   cwb is the mean of income, education, housing and labour force            │
#> │   activity (Indigenous Services Canada's Community Well-Being index).       │
#> │   CWB-style, close to but not the Department's index: the public            │
#> │   profile's inputs differ (ages 25 and over; housing from households,       │
#> │   pooled for First Nations, Métis and Inuit people).                        │
#> │   Never a measure of wellness as First Nations, Métis or Inuit              │
#> │   people define it. A group below minpop people gets no score.              │
#> │   Housing needs the household table (dwell).                                │
#> ╰─────────────────────────────────────────────────────────────────────────────╯
plot(cw)

Intergenerational exposure

Studies of the residential schools ask whether harm passes down families. Among them are Bombay, Matheson and Anisman (2014), Kaspar (2014), Hackett, Feeny and Tompa (2016), and McQuaid et al. (2017).

The counts below are made up, to show the computation:

genexpo(
  parent = c(TRUE, FALSE, TRUE, NA),
  gparent = c(TRUE, FALSE, FALSE, TRUE)
)
#> [1]  2  0  1 NA
et <- exptrnd(cases = c(60, 95, 70), n = c(500, 520, 300))
plot(et)

Individual-level data of this kind belong to the communities who collected them. These functions are for analyses done under their authority, or on tables they have published.

The Calls to Action

Call to Action 19 asks for “measurable goals to identify and close the gaps in health outcomes between Aboriginal and non-Aboriginal communities”, for indicators that include infant mortality, suicide, mental health, addictions, life expectancy and chronic diseases. trcind() lists which published table measures each of them. trchlthgap() measures every gap against non-Indigenous people on one scale, the disadvantage ratio: the ratio of the two estimates, turned so that above 1 means the group fares worse, with its interval. It is a ratio of published estimates, not a relative risk.

g <- trchlthgap()
g[
  g$measure %in% c("Life expectancy at birth", "Infant mortality rate"),
  c(
    "measure",
    "period",
    "group",
    "estimate",
    "refest",
    "disratio",
    "lower",
    "upper"
  )
]
#>                     measure    period
#> 1  Life expectancy at birth 2016/2021
#> 2  Life expectancy at birth 2016/2021
#> 3  Life expectancy at birth 2016/2021
#> 4  Life expectancy at birth 2016/2021
#> 5  Life expectancy at birth 2016/2021
#> 6  Life expectancy at birth 2016/2021
#> 7     Infant mortality rate 2014/2016
#> 8     Infant mortality rate 2014/2016
#> 9     Infant mortality rate 2014/2016
#> 10    Infant mortality rate 2014/2016
#> 11    Infant mortality rate 2014/2016
#>                                                group estimate refest disratio
#> 1                                              Métis    80.47  84.99 1.056170
#> 2                                              Inuit    71.91  84.99 1.181894
#> 3               First Nations (Registered or Treaty)    73.98  84.99 1.148824
#> 4   First Nations (Registered or Treaty), on reserve    71.42  84.99 1.190003
#> 5  First Nations (Registered or Treaty), off reserve    76.04  84.99 1.117701
#> 6           First Nations, without registered status    79.23  84.99 1.072700
#> 7                                              Métis     5.20   4.20 1.238095
#> 8                                              Inuit    13.50   4.20 3.214286
#> 9               First Nations (Registered or Treaty)     9.60   4.20 2.285714
#> 10  First Nations (Registered or Treaty), on reserve    10.00   4.20 2.380952
#> 11 First Nations (Registered or Treaty), off reserve     8.10   4.20 1.928571
#>        lower    upper
#> 1  1.0497867 1.062592
#> 2  1.1636936 1.200379
#> 3  1.1433513 1.154323
#> 4  1.1833943 1.196648
#> 5  1.1078697 1.127620
#> 6  1.0617669 1.083745
#> 7  0.8366176 1.832235
#> 8  2.0315414 5.085613
#> 9  1.7776500 2.938987
#> 10 1.4796282 3.831323
#> 11 1.3632810 2.728262
plot(g)

Life expectancy at birth in 2016/2021 was 71.4 years for First Nations people living on reserve and 71.9 for Inuit, against 85.0 for non-Indigenous people. Infant mortality among Inuit was 13.5 per 1,000 live births in 2014/2016, against 4.2.

How fast the Calls themselves are being completed:

The Truth and Reconciliation Commission made 94 Calls to Action in 2015. trccmpl() holds the number completed each year, as the Yellowhead Institute counts them (Jewell and Mosby 2021, 2022). trcprog() takes that pace as a Poisson rate:

p <- trcprog()
p
#> ╭── Truth and Reconciliation Calls to Action: progress ────────────────────╮
#> │   Completed        13 of 94                                              │
#> │   Remaining        81                                                    │
#> │   Years counted    7.04                                                  │
#> │   Calls a year     1.85 [0.983, 3.160]                                   │
#> │   All complete by  2066 (2048 to 2105)                                   │
#> │                                                                          │
#> │   At this pace the remaining 81 Calls are complete around 2066.          │
#> │                                                                          │
#> │   Completions taken as a Poisson process at a constant rate; the range   │
#> │   comes from the rate's interval. Trackers count completion differently. │
#> ╰──────────────────────────────────────────────────────────────────────────╯
plot(p)

The projection only extends the pace so far. It is not a forecast.

Cultural continuity

Chandler and Lalonde (1998) found that British Columbia First Nations communities with more markers of cultural continuity had far fewer youth suicides. Hallett, Chandler and Lalonde (2007) found the same for knowledge of an Indigenous language. Their unit was the community. Across provinces and territories, the published tables allow a coarser question: where more First Nations people know an Indigenous language, do more rate their mental health highly?

prof <- cnsprof(indiex("censuspop"))
la <- profshr(prof, "Indigenous languages", know)
la <- la[
  la$age == "Total - Age groups" &
    la$group == "First Nations" &
    la$geo != "Canada" &
    la$residence == "Total - Residence on or off reserve",
]
mh2 <- indiex("mentalhealth")
mh2 <- mh2[
  mh2$rating == "Excellent or very good" & mh2$group == "First Nations",
]
both <- merge(
  la[c("geo", "estimate", "lower", "upper")],
  mh2[c("region", "estimate", "lower", "upper")],
  by.x = "geo",
  by.y = "region",
  suffixes = c("lang", "mh")
)
cc <- connassoc(
  both,
  "geo",
  "estimatelang",
  "estimatemh",
  "lowerlang",
  "upperlang",
  "lowermh",
  "uppermh"
)
cc
#> ╭── Connection across places ──────────────────────────────────────────────╮
#> │   Places            8                                                    │
#> │   Indicators        estimatelang and estimatemh                          │
#> │   Rank correlation  rho = 0.47                                           │
#> │   p-value           0.247 (exact, every permutation)                     │
#> │   Interval          95% [-0.31, 0.79], with the published sampling error │
#> │                                                                          │
#> │   Rank correlation across 8 places: rho = 0.47, exact p = 0.247          │
#> │   Across these places the two rank moderately together.                  │
#> │                                                                          │
#> │   An association across places only:                                     │
#> │   it does not show that one causes the other for any person.             │
#> ╰──────────────────────────────────────────────────────────────────────────╯
plot(
  cc,
  xlab = "First Nations who know an Indigenous language (%)",
  ylab = "First Nations: mental health excellent or very good (%)"
)

Across these eight provinces and territories, the rank correlation is positive but far from settled: its exact p-value is large and its sensitivity band includes zero. That is not evidence for or against the communities’ findings. Eight places are far too few, a province mixes many communities, and the two measures come from different surveys and years. The community-level design needs the census profile of each First Nations community (98-510-X2021002 at the data service), and outcomes measured for the same communities.

Your own data, step by step

Nothing here is tied to Canada or to these tables. Each estimator takes a data frame with a group column, an estimate and its limits, and the arguments name those columns: grp, est, lower, upper, with ind (or cat, domain) for the indicator and qual for quality flags.

1. A published table. A table from another agency needs two steps first: limits from whatever it publishes in their place (limfrom() takes a standard error, a coefficient of variation or a margin of error), and its quality flags mapped onto suppressed and use-with-caution (qltyrcd()). Here, made-up numbers in an Australian layout, with relative standard errors and that agency’s flags:

own <- data.frame(
  group = c(
    "Aboriginal",
    "Torres Strait Islander",
    "Non-Indigenous Australians"
  ),
  pct = c(31.2, 27.5, 18.9),
  rse = c(4.1, 22.0, 1.2),
  flag = c("", "*", ""),
  stringsAsFactors = FALSE
)
own <- cbind(own, limfrom(own$pct, cv = own$rse)[c("lower", "upper")])
own$quality <- qltyrcd(own$flag, suppr = "**", caution = "*")
indgap(
  own,
  "group",
  "pct",
  "lower",
  "upper",
  ref = "Non-Indigenous Australians",
  qual = "quality"
)
#> ╭── Gap against a reference group ────────────────────────────────────────────╮
#> │   Reference  Non-Indigenous Australians                                     │
#> │   Groups     2                                                              │
#> │   Intervals  95%                                                            │
#> │                                                                             │
#> │   group                     estimate  refest  ratio  ratio CI               │
#> │   ------------------------  --------  ------  -----  ------------           │
#> │   Aboriginal                   31.20   18.90  1.651  [1.52, 1.80]           │
#> │   Torres Strait Islander *     27.50   18.90  1.455  [0.92, 2.31]           │
#> │                                                                             │
#> │   group                     difference  difference CI                       │
#> │   ------------------------  ----------  -------------                       │
#> │   Aboriginal                     12.30  [9.8, 14.8]                         │
#> │   Torres Strait Islander *        8.60  [-3.3, 20.5]                        │
#> │                                                                             │
#> │   Above the reference (interval above 1): Aboriginal.                       │
#> │   Not distinguishable from 1 (interval includes it): Torres Strait          │
#> │   Islander *.                                                               │
#> │                                                                             │
#> │   ratio: the group's estimate over the reference's; difference: the         │
#> │   group's less the reference's, in the estimate's own units.                │
#> │   Report a gap beside each group's own level, not alone.                    │
#> │   * use with caution: the source flags the estimate (E, high sampling       │
#> │   variability).                                                             │
#> ╰─────────────────────────────────────────────────────────────────────────────╯

Name the groups as the peoples concerned name themselves, and give pooled any row that pools several of them, so it is never set beside its own parts.

2. Microdata. With one row per respondent, estimate first, then compare. svydes() describes the design (here only a weight; strata, clusters and replicate weights are its other arguments), and svyest() gives each group’s estimate with its limits in the columns the estimators take. Made-up responses:

set.seed(7)
n <- 600
resp <- data.frame(
  group = sample(c("First Nations", "Inuit", "Non-Indigenous"), n, TRUE),
  wt = runif(n, 50, 150)
)
resp$goodmh <- rbinom(n, 1, ifelse(resp$group == "Non-Indigenous", 0.6, 0.5))
des <- svydes(resp, w = "wt")
tab <- as.data.frame(
  svyest(des, resp, "goodmh", by = "group", stat = "proportion", lvls = "1")
)
indgap(tab, "group", "estimate", "lower", "upper")
#> ╭── Gap against a reference group ──────────────────────────────────────╮
#> │   Reference  Non-Indigenous                                           │
#> │   Groups     2                                                        │
#> │   Intervals  95%                                                      │
#> │                                                                       │
#> │   group          estimate  refest   ratio  ratio CI        difference │
#> │   -------------  --------  ------  ------  --------------  ---------- │
#> │   First Nations    0.5152  0.5710  0.9023  [0.747, 1.090]     -0.0558 │
#> │   Inuit            0.4227  0.5710  0.7402  [0.600, 0.913]     -0.1484 │
#> │                                                                       │
#> │   group          difference CI                                        │
#> │   -------------  ----------------                                     │
#> │   First Nations  [-0.158, 0.046]                                      │
#> │   Inuit          [-0.250, -0.047]                                     │
#> │                                                                       │
#> │   Below the reference (interval below 1): Inuit.                      │
#> │   Not distinguishable from 1 (interval includes it): First Nations.   │
#> │                                                                       │
#> │   ratio: the group's estimate over the reference's; difference: the   │
#> │   group's less the reference's, in the estimate's own units.          │
#> │   Report a gap beside each group's own level, not alone.              │
#> ╰───────────────────────────────────────────────────────────────────────╯

3. The other estimators, same layout. A table with an indicator column goes to welprof(), mhcont(), healchg() (add a period column), holbal() (a domain column) or commidx() (with the community’s own weights), named the same way:

prof <- data.frame(
  community = rep(c("Lakeside", "Riverbend"), each = 3),
  item = rep(c("Hopeful", "Satisfied", "Rarely hopeful"), 2),
  pct = c(64, 71, 12, 58, 66, 16),
  lo = c(59, 66, 9, 52, 60, 12),
  hi = c(69, 76, 15, 64, 72, 20),
  stringsAsFactors = FALSE
)
welprof(prof, "community", "item", "pct", "lo", "hi",
        well = c("Hopeful", "Satisfied"), distress = "Rarely hopeful")
#> ╭── Wellness profile ─────────────────────────────────────────────────────────╮
#> │   Groups               Lakeside; Riverbend                                  │
#> │   Wellness indicators  2                                                    │
#> │   Distress indicators  1                                                    │
#> │                                                                             │
#> │   Hopeful (wellness)                                                        │
#> │   group      estimate  interval                                             │
#> │   ---------  --------  --------                                             │
#> │   Lakeside         64  [59, 69]                                             │
#> │   Riverbend        58  [52, 64]                                             │
#> │                                                                             │
#> │   Satisfied (wellness)                                                      │
#> │   group      estimate  interval                                             │
#> │   ---------  --------  --------                                             │
#> │   Lakeside         71  [66, 76]                                             │
#> │   Riverbend        66  [60, 72]                                             │
#> │                                                                             │
#> │   Rarely hopeful (distress)                                                 │
#> │   group      estimate  interval                                             │
#> │   ---------  --------  --------                                             │
#> │   Lakeside         12  [9, 15]                                              │
#> │   Riverbend        16  [12, 20]                                             │
#> │                                                                             │
#> │   On none of the 3 indicators do the groups' intervals part: the data       │
#> │   cannot rank the groups on any of them.                                    │
#> │                                                                             │
#> │   Wellness indicators come first, so each group's strengths are read        │
#> │   before its distress (Stewart 2008).                                       │
#> ╰─────────────────────────────────────────────────────────────────────────────╯

4. Plain numbers. Counts by age group go to agestd() and age-specific percentages to stdprev(), as vectors, with your own standard population:

agestd(events = c(4, 9, 11, 6, 3), pop = c(9000, 8000, 6000, 4000, 2000),
       std = c(30000, 25000, 20000, 15000, 10000))
#> ╭── Directly age-standardised rate ─────────────────────────────────╮
#> │   Groups    1                                                     │
#> │   Events    33                                                    │
#> │   Rate per  100,000                                               │
#> │   Interval  gamma (Fay and Feuer 1997)                            │
#> │                                                                   │
#> │   group  events   rate  95% CI                                    │
#> │   -----  ------  -----  ---------------                           │
#> │   all        33  115.6  [79.24, 165.00]                           │
#> │                                                                   │
#> │   all: rate 115.6 [79.24, 165.00].                                │
#> │                                                                   │
#> │   Each rate is what the group's age-specific rates give for the   │
#> │   standard population's ages, so the groups compare without them. │
#> │   Compare only rates standardised to the same population.         │
#> ╰───────────────────────────────────────────────────────────────────╯

What the package does not compute

Some instruments are individual questionnaires, and aggregate tables cannot reproduce them:

Two-Eyed Seeing (Bartlett, Marshall and Marshall 2012) is a way of bringing knowledges together, not a formula. The package follows it by reporting results beside one another without merging them into one index (dualcont(), welprof()).

Acknowledgements

The estimators for holistic balance, wellness beside distress, the four directions, healing as a process and connection across places put to quantitative use the frameworks of Professor Suzanne Stewart (Yellowknives Dene First Nation), Director of the Waakebiness-Bryce Institute for Indigenous Health at the University of Toronto’s Dalla Lana School of Public Health, and Professor Angela Mashford-Pringle (Timiskaming First Nation), its Associate Director. This work stands on theirs. Any error in turning their qualitative work into estimators is the package’s, not theirs.

References

Search and rescue: detection, coverage and the best use of effort

Search theory began with Koopman’s work for the U.S. Navy (1946) and is how search and rescue services plan a search: how likely a sensor is to see the object at a given distance, how much ground a search covers, how likely the object is to be in each area, and where a limited effort does the most good. Frost’s The Theory of Search: A Simplified Explanation (1996) explains it for search planners in four chapters of figures. This vignette walks those chapters and draws every one of the book’s 36 figures with srchfig(), which builds each from the package’s functions with the book’s parameters. srchfig() with no argument lists them and says how each is drawn: from the model, as a schematic, or as an approximation where the book’s plate cannot be recovered exactly.

table(srchfig()$how)
#> 
#> approximation         model     schematic 
#>             3            29             4

The functions behind the figures work on any inputs: ltrlrng() (lateral range curves), swpwdth() (sweep width), srchpod() and pswppod() (probability of detection), srchscs() (probability of success and the update after a search), datumpmap() (where the object may be), opteffort(), optcrcsrch() and optsqr() (the best use of effort). The book’s figures come first; Your own search then runs each function on numbers of your own.

Datums and probability maps

Figure 3-1. A point datum is a circular normal distribution of possible positions. The probability that the object is within R standard deviations of the datum is POC = 1 - exp(-R^2 / 2).

srchfig("3-1")

Figures 3-2, 3-3 and 3-5. The square circumscribed about a circle of radius R holds (2 Phi(R) - 1)^2 of the distribution: 47% for one standard deviation, 58% for the total probable error of position E = sqrt(2 ln 2) = 1.1774 standard deviations (the radius holding 50%), and 65% for the simplified method’s first search radius of 1.1 E.

srchfig("3-2")

srchfig("3-3")

srchfig("3-5")

Figure 3-4. A probability map lays a grid over the distribution. With cells 2E wide the centre holds 57.91%, a side cell 9.08% and a corner 1.42%.

srchfig("3-4")

Figures 3-6 and 3-7. The simplified method extended a point datum to a line by drawing circles at the two end points and boxing them in. The density this implies is drawn as a circular normal spread along the segment (datum = "line"), an approximation of the book’s plate.

srchfig("3-6")

srchfig("3-7")

Figures 3-8 to 3-10. A simpler line datum is normal across the line and uniform along it, with nothing beyond its ends (datum = "corridor"). A strip one standard deviation either side holds 68%, a strip of 2E holds 76.1%, and dividing strips into cells along the line gives the completed map of figure 3-10.

srchfig("3-8")

srchfig("3-9")

srchfig("3-10")

Figures 3-11 and 3-12. An area datum spreads the probability evenly: each of 20 equal cells holds 5%.

srchfig("3-11")

srchfig("3-12")

Figures 3-13 and 3-14. A generalized datum takes the planner’s own cell probabilities (datum = "general", wts =). The book’s surface is illustrative; here it is drawn from the cells of figure 3-14.

srchfig("3-13")

srchfig("3-14")

Lateral range curves and sweep width

Figure 4-1. A lateral range curve p(x) is the probability of detecting an object that passes at distance x from the track. The sweep width is its area, W = integral of p(x) dx; the dashed rectangles have the same area.

srchfig("4-1")

Figures 4-2 and 4-3. A definite range sensor sees everything within W / 2 of its track and nothing beyond: it sweeps a strip exactly W wide.

srchfig("4-2")

srchfig("4-3")

Figures 4-4 and 4-5. The M-Beta curve has height M over a width of W / M; with M = 0.5 and M = 0.25 the sweep width is unchanged.

srchfig("4-4")

srchfig("4-5")

Figures 4-6 and 4-7. Koopman’s inverse cube law takes the chance of a glimpse to be proportional to the solid angle of the wake, k h / r^3, which gives p(x) = 1 - exp(-W^2 / (4 pi x^2)). Its sweep width is W: the objects missed inside the sweep width (A) balance those seen outside it (B).

srchfig("4-6")

srchfig("4-7")

swpwdth(ltrlrng(seq(-500, 500, by = 0.01), W = 1))
#> [1] 0.9996817

Effort, coverage and probability of detection

Figure 5-1. Search effort is sweep width times track length, Z = W L, and coverage is effort over area, C = Z / A. A search whose path is random gives POD = 1 - exp(-C).

srchfig("5-1")

Figures 5-2 to 5-5. A parallel sweep at track spacing S has coverage C = W / S. With a definite range sensor the POD is min(C, 1): half the gap is seen at S = 2W, all of it at S = W.

srchfig("5-2")

srchfig("5-3")

srchfig("5-4")

srchfig("5-5")

Figures 5-6 and 5-7. M-Beta curves overlap once the tracks are closer than their width: POD = C below C = M and 1 - (1 - M)^(C / M) above it, 0.75 and 0.6836 at coverage 1 for M = 0.5 and M = 0.25, nearer the random search curve as M falls.

srchfig("5-6")

srchfig("5-7")

srchpod(c(0.5, 1, 2), "mbeta", m = 0.25)
#> ╭── Probability of detection ─────────────────────────────────────────╮
#> │   Model      M-Beta                                                 │
#> │   M          0.25                                                   │
#> │   Coverages  3                                                      │
#> │                                                                     │
#> │   coverage       pod                                                │
#> │   --------  --------                                                │
#> │    0.50000  0.437500                                                │
#> │    1.00000  0.683594                                                │
#> │    2.00000  0.899887                                                │
#> │                                                                     │
#> │   M-Beta model, M = 0.25                                            │
#> │   Coverage is effort over area (sweep width x track length / area); │
#> │   pod is the chance a search at that coverage finds the object.     │
#> ╰─────────────────────────────────────────────────────────────────────╯

Figure 5-8. The inverse cube parallel sweep gives POD = erf(sqrt(pi) C / 2), 0.7899 at coverage 1, between the definite range and random search curves. pswppod() computes the same sweep directly, for any curve.

srchfig("5-8")

pswppod(W = 1, S = 1)
#> [1] 0.7899086

Figure 5-9. Poorer conditions shrink the sweep width, and a cross-track navigation error of standard deviation sigma flattens the curve: the expected curve is the convolution phi(x) = integral of p(x - t) Z(t, sigma) dt, whose area is still W (ltrlrng(..., navsd = )). The plate is blank in the published file; it is drawn from the text.

srchfig("5-9")

Your own search, step by step

The same functions plan a real search from your own numbers, one call per step. Use one distance unit throughout (kilometres here; nautical miles at sea): the grid, the positions, the sweep width and the track length are all in it, and areas in its square.

1. The search area and its grid. A hiker is overdue; the last known position (LKP) is 2 km east and 3 km north of the base. Lay a grid of cells over the area you could search: here 12 km by 12 km in half-kilometre cells. xedg and yedg are the cell edges, west to east and south to north; cells need not be equal.

xe <- seq(-4, 8, by = 0.5)
ye <- seq(-3, 9, by = 0.5)

2. The probability map. datumpmap() puts the chance the hiker is in each cell on the grid. For a point datum give the LKP as center and the position error as sigma, one standard deviation (about 1.5 km, from the time and terrain since the LKP). A route between two points is datum = "line" or "corridor" with end, a region known only by its bounds datum = "area", and a map of your own (cell shares from the planners’ scenarios, say) is datum = "general" with wts, one row per y cell and one column per x cell.

m <- datumpmap(xe, ye, center = c(2, 3), sigma = 1.5)
m
#> ╭── Datum probability map (point datum) ──────────────────────────────────╮
#> │   Datum   point                                                         │
#> │   Centre  2, 3                                                          │
#> │   Sigma   1.5                                                           │
#> │   Grid    24 x 24 cells (x by y)                                        │
#> │                                                                         │
#> │   The most likely cells                                                 │
#> │       x      y      poc                                                 │
#> │   -----  -----  -------                                                 │
#> │   1.750  2.750  0.01705                                                 │
#> │   1.750  3.250  0.01705                                                 │
#> │   2.250  2.750  0.01705                                                 │
#> │   2.250  3.250  0.01705                                                 │
#> │   1.250  2.750  0.01527                                                 │
#> │                                                                         │
#> │   24 x 24 cells (x by y); 0.0% of the probability lies outside the grid │
#> ╰─────────────────────────────────────────────────────────────────────────╯
plot(m)

Two scenarios weighted by the planners (70% on the LKP, 30% that the hiker took the trail to the north-east) are mixed cell by cell, and the mixture is a general datum (its weights are scaled to sum to 1 over the grid):

trail <- datumpmap(xe, ye, "corridor", c(2, 3), sigma = 0.5, end = c(2, 8))
mix <- datumpmap(xe, ye, "general", wts = 0.7 * m$poc + 0.3 * trail$poc)
plot(mix)

3. The sweep width and the detection model. The sweep width W comes from a published table for the sensor, the object and the conditions, or from a field trial of your own with swpwdth() (above). The detection model says how detection falls off with distance: "invcube" for good conditions and evenly spaced tracks, "random" for poor conditions or uneven tracks (the conservative choice), "definite" for a sensor that sees everything within W / 2. Here the ground teams’ sweep width is 0.6 km.

4. Effort, and the probability of detection. The effort is the track length walked in each cell. Say today’s teams walked 2 km of track in each of the four cells round the LKP. srchpod() with W, effort and area turns that into the coverage and the POD of each cell; area = m takes the map’s cell areas, and effort is laid out like the map (one row per y cell):

walked <- matrix(0, nrow(m$poc), ncol(m$poc))
walked[12:13, 12:13] <- 2
pd <- srchpod(W = 0.6, effort = walked, area = m, model = "random")
range(pd)
#> [1] 0.0000000 0.9917703

5. The probability of success, and the map after a search that found nothing. srchscs() takes the map and the PODs: pos is the chance today’s search found the hiker, and map the probability map updated by what was searched (cells searched hard become less likely, the rest more likely):

day1 <- srchscs(m, pd)
day1
#> ╭── Probability of success ────────────────────────────────────────────╮
#> │   Areas        576                                                   │
#> │   Searches     1                                                     │
#> │   Overall POS  0.06762                                               │
#> │   Map          24 x 24 cells (x by y), updated in $map               │
#> │                                                                      │
#> │   area      poc  cumpod  posarea   pocupd                            │
#> │   ----  -------  ------  -------  -------                            │
#> │    276  0.01705  0.9918  0.01691  0.00015                            │
#> │    277  0.01705  0.9918  0.01691  0.00015                            │
#> │    300  0.01705  0.9918  0.01691  0.00015                            │
#> │    301  0.01705  0.9918  0.01691  0.00015                            │
#> │    252  0.01527  0.0000  0.00000  0.01638                            │
#> │    253  0.01527  0.0000  0.00000  0.01638                            │
#> │    275  0.01527  0.0000  0.00000  0.01638                            │
#> │    278  0.01527  0.0000  0.00000  0.01638                            │
#> │    299  0.01527  0.0000  0.00000  0.01638                            │
#> │    302  0.01527  0.0000  0.00000  0.01638                            │
#> │    324  0.01527  0.0000  0.00000  0.01638                            │
#> │    325  0.01527  0.0000  0.00000  0.01638                            │
#> │    251  0.01368  0.0000  0.00000  0.01467                            │
#> │    254  0.01368  0.0000  0.00000  0.01467                            │
#> │    323  0.01368  0.0000  0.00000  0.01467                            │
#> │    326  0.01368  0.0000  0.00000  0.01467                            │
#> │    228  0.01225  0.0000  0.00000  0.01314                            │
#> │    229  0.01225  0.0000  0.00000  0.01314                            │
#> │    274  0.01225  0.0000  0.00000  0.01314                            │
#> │    298  0.01225  0.0000  0.00000  0.01314                            │
#> │                                                                      │
#> │   POS 0.06762 after 1 search                                         │
#> │   poc: chance the object is in the area; cumpod: chance the searches │
#> │   detect it there; pocupd: the poc after the searches found nothing. │
#> │                                                                      │
#> │   The 20 most likely of 576 cells.                                   │
#> ╰──────────────────────────────────────────────────────────────────────╯
plot(day1$map)

6. Plan the next search. opteffort() spreads tomorrow’s 40 km of track over the updated map to make the chance of finding the hiker as large as possible (Stone 1975), under random search detection. Its result goes straight back into srchscs() for the day after:

day2 <- opteffort(day1$map, W = 0.6, effort = 40)
day2
#> ╭── Optimal allocation of effort ─────────────────────────────╮
#> │   Total effort    40                                        │
#> │   Sweep width W   0.6                                       │
#> │   Cells searched  108 of 576                                │
#> │   POS             0.5467                                    │
#> │                                                             │
#> │   cell      poc  effort  coverage     pod                   │
#> │   ----  -------  ------  --------  ------                   │
#> │    252  0.01638  0.7356     1.766  0.8289                   │
#> │    253  0.01638  0.7356     1.766  0.8289                   │
#> │    275  0.01638  0.7356     1.766  0.8289                   │
#> │    278  0.01638  0.7356     1.766  0.8289                   │
#> │    299  0.01638  0.7356     1.766  0.8289                   │
#> │    302  0.01638  0.7356     1.766  0.8289                   │
#> │    324  0.01638  0.7356     1.766  0.8289                   │
#> │    325  0.01638  0.7356     1.766  0.8289                   │
#> │    251  0.01467  0.6898     1.655  0.8090                   │
#> │    254  0.01467  0.6898     1.655  0.8090                   │
#> │    323  0.01467  0.6898     1.655  0.8090                   │
#> │    326  0.01467  0.6898     1.655  0.8090                   │
#> │    228  0.01314  0.6439     1.545  0.7868                   │
#> │    229  0.01314  0.6439     1.545  0.7868                   │
#> │    274  0.01314  0.6439     1.545  0.7868                   │
#> │    298  0.01314  0.6439     1.545  0.7868                   │
#> │    279  0.01314  0.6439     1.545  0.7868                   │
#> │    303  0.01314  0.6439     1.545  0.7868                   │
#> │    348  0.01314  0.6439     1.545  0.7868                   │
#> │    349  0.01314  0.6439     1.545  0.7868                   │
#> │                                                             │
#> │   POS 0.5467 from effort 40 over 108 of 576 cells (W = 0.6) │
#> │                                                             │
#> │   ... and 88 more cells                                     │
#> ╰─────────────────────────────────────────────────────────────╯
day3 <- srchscs(day1$map, day2)
day3$pos
#> [1] 0.546653

7. Plot. Every result has a plot() method: the allocation as a map of the effort in each cell, the map as an image, a surface or its cells with their percentages.

plot(day2)

plot(day1$map, view = "both")

For a single point datum with no grid, optcrcsrch() gives Koopman’s optimal circle directly, and optsqr() the optimal square, search after search:

optcrcsrch(sigma = 1.5, W = 0.6, effort = 40)
#> ╭── Optimal search of a circular normal datum ────────────────────────╮
#> │   sigma    1.5                                                      │
#> │   W        0.6                                                      │
#> │   effort   40                                                       │
#> │   radius   2.88                                                     │
#> │   pos      0.5497                                                   │
#> │   cntrdns  3.071                                                    │
#> │                                                                     │
#> │   Spread the effort evenly over a circle of radius 2.88: the chance │
#> │   of finding the object is then 55.0%.                              │
#> │                                                                     │
#> │   cntrdns: the coverage at the centre of the circle.                │
#> ╰─────────────────────────────────────────────────────────────────────╯
optsqr(c(40, 40) * 0.6, E = 1.1774 * 1.5)
#> ╭── Optimal square searches ────────────────────────────────────────────────╮
#> │   Searches                 2                                              │
#> │   Total effort             48                                             │
#> │   Chance found by the end  0.8047                                         │
#> │                                                                           │
#> │   search  effort  releff  cumeff  factor   side     poc  coverage     pod │
#> │   ------  ------  ------  ------  ------  -----  ------  --------  ------ │
#> │        1      24   7.695    7.69   1.270  4.485  0.7484    1.1929  0.8651 │
#> │        2      24   7.695   15.39   1.696  5.991  0.7460    0.6686  0.5980 │
#> │                                                                           │
#> │   search     pos  cumpos                                                  │
#> │   ------  ------  ------                                                  │
#> │        1  0.6475  0.6475                                                  │
#> │        2  0.4461  0.8047                                                  │
#> │                                                                           │
#> │   Each search covers the square (of side 'side') that gives it the best   │
#> │   chance of success after the searches before it. poc is the chance the   │
#> │   object is in that square, pod the chance of detecting it there, pos     │
#> │   their product and cumpos the running chance it has been found.          │
#> ╰───────────────────────────────────────────────────────────────────────────╯

References

Frost, J. R. (1996). The Theory of Search: A Simplified Explanation. Soza and Company and U.S. Coast Guard Office of Search and Rescue. Chapters 3-6 and their figures.

Koopman, B. O. (1946). Search and Screening. OEG Report 56, Office of the Chief of Naval Operations.

Koopman, B. O. (1956-1957). The theory of search, parts I-III. Operations Research 4(3), 324-346; 4(5), 503-531; 5(5), 613-626.

Stone, L. D. (1975). Theory of Optimal Search. Academic Press.