---
title: "SDP teacher retention — the calculation"
output: md_document
---

```{r setup, include = FALSE}
# error = FALSE so that knitting *fails* on an error rather than embedding it in
# the document. A calculation that quietly renders its own stack trace and
# carries on is worse than one that stops.
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", error = FALSE)
options(stringsAsFactors = FALSE)
```

This document is both the description of the retention calculation and the
program that performs it. Every number the site publishes about teacher
retention is produced by running this file; the prose between the chunks is
where the definitions are argued for, which is the part that usually goes
missing when a calculation is written as a script.

It is one of three implementations. The other two — the original TypeScript in
`src/`, and `calc/retention.py` — compute the same facts from the same files and
write them out in the same format, described in `calc/CONTRACT.md`. The three
output files are compared byte for byte. Where a definition below is stated
oddly, it is usually because it is stating exactly what the other two do.

```{r constants}
ROOT <- Sys.getenv("SDP_ROOT", ".")
DATA <- file.path(ROOT, "data", "extracted")

snapshots <- data.frame(
  year = 2015:2025,
  file = c(
    "employee_information_Oct2015.csv", "employee_information_Oct2016.csv",
    "employee_information_Oct2017.csv", "employee_information_Oct2018.csv",
    "employee_information_Oct2019.csv", "employee_information_Oct2020.csv",
    "employee_information_Oct2021.csv", "employee_information_Oct2022.csv",
    "Employee_Data_2023-2024/employee_information_October2023.csv",
    "employee_information_Oct2024.csv",
    "Oct 2025/employee_information_Oct2025.csv"
  )
)

years <- snapshots$year
transitions <- data.frame(from = head(years, -1), to = tail(years, -1))

# Peer levels, and the token each becomes in a contract key. Keys carry no
# spaces, so that every implementation sorts them the same way.
levels_tbl <- data.frame(
  level = c("ELEMENTARY", "MIDDLE SCHOOL", "HIGH SCHOOL"),
  token = c("ELEMENTARY", "MIDDLE_SCHOOL", "HIGH_SCHOOL")
)
scopes <- c("teacher", "teacherWithAlt", "all")
MIN_TEACHERS <- 15

k8_orgs <- c("1300", "1280", "1340", "1260", "1390", "1470",
             "2380", "2590", "2470", "2480", "2210")
schools <- rbind(
  data.frame(org = k8_orgs, level = "ELEMENTARY"),
  data.frame(org = c("2050", "2650", "2680", "6040"),
             level = c("MIDDLE SCHOOL", "HIGH SCHOOL", "HIGH SCHOOL", "HIGH SCHOOL"))
)
```

## Reading the extracts

The district publishes an "Employee Information" extract every quarter. Only the
October ones are used: October is the district's own census month, so each pair
of consecutive Octobers is one school year.

The files are **latin-1**, not UTF-8 — a handful of accented names contain bytes
that are not valid UTF-8 and would otherwise decode to replacement characters.
They are read as raw bytes and converted explicitly, rather than left to
whatever the session's locale happens to be.

```{r read}
read_extract <- function(path) {
  raw <- readBin(path, "raw", file.size(path))
  txt <- iconv(rawToChar(raw), from = "latin1", to = "UTF-8")
  read.csv(text = txt, colClasses = "character", check.names = FALSE, na.strings = NULL)
}

raw_snap <- lapply(seq_len(nrow(snapshots)), function(i) {
  path <- file.path(DATA, snapshots$file[i])
  if (!file.exists(path)) stop("missing snapshot: ", path, " - run ./fetch-data.sh")
  read_extract(path)
})
names(raw_snap) <- as.character(years)

sapply(raw_snap, nrow)
```

## Who counts as a teacher

A classroom teacher needs **both** a PFT teaching local and a TEACHER title.
Requiring the union membership as well as the title is what excludes the
`NON REP` rows — student teachers on placement, and `CONTRACTED TEACHER`
positions bought in from outside — who are not district hires and turn over by
design. It also excludes the PFT-represented counsellors, nurses and therapists,
who are not classroom staff.

The title test is "TEACHER followed by a word boundary", written out here as an
explicit character test rather than as `\b`, because the three languages
disagree about whether `\b` is Unicode-aware and this is not a place to find
that out later.

```{r definitions}
is_teacher <- function(rep, title) {
  t <- toupper(title)
  startsWith(toupper(rep), "PFT") &
    startsWith(t, "TEACHER") &
    !grepl("^[A-Za-z0-9_]", substr(t, 8, 8))
}

# Instructional staff who are not district hires.
is_alt <- function(title) {
  t <- toupper(title)
  grepl("STUDENT TEACHER", t, fixed = TRUE) | grepl("CONTRACTED TEACHER", t, fixed = TRUE)
}

# The building principal, not assistant principals or retired fill-ins.
is_principal <- function(title) trimws(toupper(title)) == "PRINCIPAL"
```

## Identity across years

The extracts carry no employee number, so a person is matched between two
Octobers by name alone. Case, punctuation and internal spacing are normalized
away, which stops a cosmetic spelling change from reading as a departure.

Two limits survive this and are stated on the site: somebody who changes their
surname reads as one departure and one arrival, and two people who share a name
read as one person.

```{r identity}
norm_name <- function(s) trimws(gsub("\\s+", " ", gsub("[^A-Z ]", "", toupper(s))))
emp_key <- function(last, first) paste0(norm_name(last), "|", norm_name(first))

staff <- lapply(names(raw_snap), function(y) {
  df <- raw_snap[[y]]
  data.frame(
    key       = emp_key(df$LAST_NAME, df$FIRST_NAME),
    org       = df$HOME_ORGANIZATION,
    level     = toupper(df$ORGANIZATION_LEVEL),
    teacher   = is_teacher(df$TYPE_OF_REPRESENTATION, df$TITLE_DESCRIPTION),
    alt       = is_alt(df$TITLE_DESCRIPTION),
    principal = is_principal(df$TITLE_DESCRIPTION),
    name      = paste(df$FIRST_NAME, df$LAST_NAME)
  )
})
names(staff) <- names(raw_snap)

data.frame(year = years,
           rows = sapply(staff, nrow),
           teachers = sapply(staff, function(e) sum(e$teacher)),
           principals = sapply(staff, function(e) sum(e$principal)))
```

## Indexing

Every school is compared against roughly 150 peers across ten transitions and
three staff definitions. Grouping the rows by school once, up front, is what
keeps that from being quadratic.

`keys_by[[year]][[scope]]` holds the identity keys of each school's staff under
one of the three definitions; the row counts are kept separately, because a
headcount counts rows and a retention rate counts distinct people.

```{r index}
scope_mask <- function(e, scope) {
  if (scope == "teacher") e$teacher
  else if (scope == "teacherWithAlt") e$teacher | e$alt
  else rep(TRUE, nrow(e))
}

keys_by <- lapply(staff, function(e) {
  out <- lapply(scopes, function(s) {
    m <- scope_mask(e, s)
    split(e$key[m], e$org[m])
  })
  names(out) <- scopes
  out
})

# Rows, not distinct people: this is the headcount used for the alternative-staff
# share and for the peer-set threshold.
rows_by <- lapply(keys_by, function(y) lapply(y, lengths))

at <- function(bucket, org) {
  v <- bucket[[org]]
  if (is.null(v)) character(0) else v
}
count_at <- function(bucket, org) {
  # `rows_by` holds named integer vectors, where indexing a name that is not
  # there is an error rather than a NULL; a school with no staff of some kind in
  # some year is an absent name, not a zero.
  v <- bucket[org]
  if (length(v) == 0 || is.na(v)) 0L else as.integer(v)
}
```

## Retention between two Octobers

Retention is the share of one October's people who are still there the next
October, counting distinct people rather than rows. It is deliberately *not*
`1 - leavers/headcount`: those differ whenever the headcount changes, and the
question is about the starting cohort.

```{r retention}
retention <- function(before, after) {
  start <- length(unique(before))
  stayed <- if (start == 0) 0L else length(intersect(before, after))
  list(stayed = as.integer(stayed), start = as.integer(start),
       rate = if (start == 0) NULL else stayed / start)
}
```

## Peer sets

A retention figure means little on its own — 2021 was a bad year everywhere. Each
school is therefore also ranked against schools of its own level with at least
`r MIN_TEACHERS` teachers that October.

A school's level is read from the **first** row that school contributes to the
file, matching the original implementation: a few organizations carry rows
labelled with more than one level, and taking the first row is a rule that all
three implementations can follow without ambiguity.

```{r peers}
peer_cache <- new.env(parent = emptyenv())

peer_orgs <- function(year, level) {
  ck <- paste(year, level)
  if (!is.null(peer_cache[[ck]])) return(peer_cache[[ck]])
  e <- staff[[as.character(year)]]
  first <- !duplicated(e$org)
  org_level <- e$level[first]
  orgs <- e$org[first]
  tcount <- lengths(split(e$key[e$teacher], factor(e$org[e$teacher], levels = orgs)))
  out <- orgs[grepl(level, org_level, fixed = TRUE) & tcount >= MIN_TEACHERS]
  peer_cache[[ck]] <- out
  out
}

rate_cache <- new.env(parent = emptyenv())

peer_rates <- function(from, to, level, scope) {
  ck <- paste(from, to, level, scope)
  if (!is.null(rate_cache[[ck]])) return(rate_cache[[ck]])
  b <- keys_by[[as.character(from)]][[scope]]
  a <- keys_by[[as.character(to)]][[scope]]
  orgs <- peer_orgs(from, level)
  out <- numeric(length(orgs))
  n <- 0L
  for (o in orgs) {
    r <- retention(at(b, o), at(a, o))
    if (!is.null(r$rate)) { n <- n + 1L; out[n] <- r$rate }
  }
  out <- out[seq_len(n)]
  rate_cache[[ck]] <- out
  out
}
```

The two summary statistics over a peer set are a median and a percentile rank.
Both are defined the plain way — the middle value, and the share of peers at or
below this school — with no interpolation and no tie correction, because the
point is to be reproducible rather than to be clever. The quartiles index
straight into the sorted rates for the same reason.

```{r stats}
med <- function(v) {
  if (!length(v)) return(NULL)
  s <- sort(v)
  mid <- length(s) %/% 2
  if (length(s) %% 2 == 1) s[mid + 1] else (s[mid] + s[mid + 1]) / 2
}

pct_rank <- function(v, target) if (!length(v)) NULL else sum(v <= target) / length(v)

quantile_at <- function(sorted, q) {
  if (!length(sorted)) return(NULL)
  sorted[floor(q * (length(sorted) - 1)) + 1]
}
```

## Where the leavers went

A school losing staff to other district schools has a different problem from one
whose teachers leave the district altogether, so each year's cohort is followed
into the next October and sorted into four outcomes.

Where the same name appears twice as a teacher in one year, the last row in the
file wins — not because that is right, but because it is what a dictionary
assignment does in the other two implementations, and the three have to agree.

```{r person-index}
teacher_index <- lapply(staff, function(e) {
  m <- e$teacher
  keep <- !duplicated(e$key[m], fromLast = TRUE)
  list(keys = e$key[m][keep], orgs = e$org[m][keep], payroll = unique(e$key))
})
```

## Principal tenures

The extracts list a principal per school per October, so a tenure is a run of
consecutive Octobers with the same name. Where a school lists two people at once
the names are joined, sorted, so that the order they appear in the file does not
create a spurious change of leadership. An October with nobody listed is its own
kind of run, not a gap to be smoothed over.

Retention is then pooled within each run. A transition belongs to whoever held
the building in its **starting** October — that is the year whose staff either
stayed or left.

```{r timeline}
compress <- function(names_by_year) {
  runs <- list()
  for (i in seq_along(names_by_year)) {
    nm <- names_by_year[i]
    yr <- years[i]
    if (length(runs) && runs[[length(runs)]]$names == nm) {
      runs[[length(runs)]]$to <- yr
    } else {
      runs[[length(runs) + 1]] <- list(names = nm, from = yr, to = yr)
    }
  }
  runs
}

principal_names <- function(year, org) {
  e <- staff[[as.character(year)]]
  nm <- e$name[e$principal & e$org == org]
  if (!length(nm)) return("(none listed)")
  # C-locale ordering, so that R, Python and JavaScript sort the same names the
  # same way regardless of the locale each happens to be running under.
  paste(nm[order(nm, method = "radix")], collapse = " + ")
}
```

## Writing the contract

The output format is a sorted TSV of `metric`, `key`, `value`. Rates are written
as their numerator and denominator, as integers, rather than as a quotient —
integers cannot drift — and the reals that remain get a fixed twelve decimal
places, which is finer than any real disagreement and coarse enough to absorb a
last-bit difference between three languages' arithmetic.

```{r contract}
NA_VALUE <- "NA"

fmt_int <- function(n) {
  n <- as.numeric(n)
  if (!is.finite(n) || n != round(n)) stop("contract: ", n, " is not an integer")
  sprintf("%.0f", n)
}

fmt_real <- function(x) {
  if (is.null(x) || (length(x) == 1 && is.na(x))) return(NA_VALUE)
  if (!is.finite(x)) stop("contract: ", x, " is not finite")
  # -0 and 0 are the same number; only their printed forms differ.
  sprintf("%.12f", if (x == 0) 0 else x)
}

fmt_text <- function(s) {
  if (grepl("[\t\r\n]", s)) stop("contract: value contains a tab or newline: ", s)
  s
}

rows <- new.env(parent = emptyenv())
rows$metric <- character(0); rows$key <- character(0); rows$value <- character(0)
add <- function(metric, key, value) {
  rows$metric <- c(rows$metric, metric)
  rows$key <- c(rows$key, key)
  rows$value <- c(rows$value, value)
}

sha256 <- function(s) digest::digest(s, algo = "sha256", serialize = FALSE)
```

## The facts

Inputs first, so that a disagreement about the numbers can be told apart from a
disagreement about what was read off the disk. The digests are the part that
does the work: two implementations that hash their parsed extracts to the same
value have demonstrably read the same file the same way, which is otherwise the
hardest thing to be sure of across three CSV readers.

```{r facts-input}
for (y in years) {
  e <- staff[[as.character(y)]]
  df <- raw_snap[[as.character(y)]]
  dump <- paste(df$LAST_NAME, df$FIRST_NAME, df$TITLE_DESCRIPTION, df$HOME_ORGANIZATION,
                df$ORGANIZATION_LEVEL, df$TYPE_OF_REPRESENTATION, sep = "\x1f")
  uk <- unique(e$key)
  add("input.rows", as.character(y), fmt_int(nrow(e)))
  add("input.digest", as.character(y), fmt_text(sha256(paste(dump, collapse = "\n"))))
  add("input.keydigest", as.character(y),
      fmt_text(sha256(paste(uk[order(uk, method = "radix")], collapse = "\n"))))
  add("input.keys", as.character(y), fmt_int(length(uk)))
  add("input.teachers", as.character(y), fmt_int(sum(e$teacher)))
  add("input.principals", as.character(y), fmt_int(sum(e$principal)))
}
```

```{r facts-peers}
for (li in seq_len(nrow(levels_tbl))) {
  level <- levels_tbl$level[li]; tok <- levels_tbl$token[li]
  for (y in years) {
    add("peer.orgs", paste0(y, "|", tok), fmt_int(length(peer_orgs(y, level))))
  }
  for (ti in seq_len(nrow(transitions))) {
    from <- transitions$from[ti]; to <- transitions$to[ti]
    for (scope in scopes) {
      r <- sort(peer_rates(from, to, level, scope))
      key <- paste0(from, "-", to, "|", tok, "|", scope)
      add("peer.n", key, fmt_int(length(r)))
      add("peer.median", key, fmt_real(med(r)))
      add("peer.q1", key, fmt_real(quantile_at(r, 0.25)))
      add("peer.q3", key, fmt_real(quantile_at(r, 0.75)))
    }
  }
}
```

```{r facts-schools}
pooled_teacher <- list(); changes <- list(); tenure <- list()

for (si in seq_len(nrow(schools))) {
  org <- schools$org[si]; level <- schools$level[si]
  tok <- levels_tbl$token[match(level, levels_tbl$level)]
  add("school.level", org, fmt_text(tok))

  series <- list()
  for (scope in scopes) {
    st <- integer(0); sc <- integer(0)
    for (ti in seq_len(nrow(transitions))) {
      from <- transitions$from[ti]; to <- transitions$to[ti]
      r <- retention(at(keys_by[[as.character(from)]][[scope]], org),
                     at(keys_by[[as.character(to)]][[scope]], org))
      p <- if (is.null(r$rate)) NULL else pct_rank(peer_rates(from, to, level, scope), r$rate)
      key <- paste0(org, "|", from, "-", to)
      add(paste0("series.", scope, ".stayed"), key, fmt_int(r$stayed))
      add(paste0("series.", scope, ".start"), key, fmt_int(r$start))
      add(paste0("series.", scope, ".percentile"), key, fmt_real(p))
      st <- c(st, r$stayed); sc <- c(sc, r$start)
    }
    series[[scope]] <- list(stayed = st, start = sc)
  }

  keep2020 <- transitions$from >= 2020
  pools <- list(
    teacher = series$teacher,
    teacherWithAlt = series$teacherWithAlt,
    all = series$all,
    teacherSince2020 = list(stayed = series$teacher$stayed[keep2020],
                            start = series$teacher$start[keep2020])
  )
  for (nm in names(pools)) {
    s <- sum(pools[[nm]]$stayed); n <- sum(pools[[nm]]$start)
    add(paste0("pooled.", nm, ".stayed"), org, fmt_int(s))
    add(paste0("pooled.", nm, ".start"), org, fmt_int(n))
    if (nm == "teacher") pooled_teacher[[org]] <- if (n == 0) NULL else s / n
  }

  for (y in years) {
    pft <- count_at(rows_by[[as.character(y)]]$teacher, org)
    alt <- count_at(rows_by[[as.character(y)]]$teacherWithAlt, org) - pft
    add("alt.pft", paste0(org, "|", y), fmt_int(pft))
    add("alt.alt", paste0(org, "|", y), fmt_int(alt))
  }

  runs <- compress(sapply(years, principal_names, org = org))
  add("principal.runs", org, fmt_int(length(runs)))
  # Runs minus one: three runs over the window means two changes.
  add("principal.changes", org, fmt_int(max(0, length(runs) - 1)))
  for (i in seq_along(runs)) {
    add("principal.run", sprintf("%s|%02d", org, i - 1),
        fmt_text(sprintf("%d-%d|%s", runs[[i]]$from, runs[[i]]$to, runs[[i]]$names)))
  }
  changes[[org]] <- max(0, length(runs) - 1)
  tenure[[org]] <- max(sapply(runs, function(r) r$to - r$from + 1))

  since <- years[years >= 2020]
  base <- unique(at(keys_by[[as.character(since[1])]]$teacher, org))
  add("cohort2020.base", org, fmt_int(length(base)))
  for (y in since) {
    now <- at(keys_by[[as.character(y)]]$teacher, org)
    add("cohort2020.remaining", paste0(org, "|", y), fmt_int(length(intersect(base, now))))
  }

  add("teacherCount", org, fmt_int(count_at(rows_by[[as.character(tail(years, 1))]]$teacher, org)))

  dest <- c(stayed = 0L, movedSchool = 0L, changedRole = 0L, leftDistrict = 0L)
  for (ti in seq_len(nrow(transitions))) {
    from <- transitions$from[ti]; to <- transitions$to[ti]
    ks <- unique(at(keys_by[[as.character(from)]]$teacher, org))
    idx <- teacher_index[[as.character(to)]]
    where <- idx$orgs[match(ks, idx$keys)]
    on_payroll <- ks %in% idx$payroll
    dest["stayed"] <- dest["stayed"] + sum(!is.na(where) & where == org)
    dest["movedSchool"] <- dest["movedSchool"] + sum(!is.na(where) & where != org)
    dest["changedRole"] <- dest["changedRole"] + sum(is.na(where) & on_payroll)
    dest["leftDistrict"] <- dest["leftDistrict"] + sum(is.na(where) & !on_payroll)
  }
  for (nm in names(dest)) add(paste0("destination.", nm), org, fmt_int(dest[[nm]]))

  eras <- list()
  for (r in runs) {
    m <- transitions$from >= r$from & transitions$from <= r$to
    s <- sum(series$teacher$stayed[m]); n <- sum(series$teacher$start[m])
    if (n > 0) eras[[length(eras) + 1]] <- list(run = r, stayed = s, start = n)
  }
  add("era.n", org, fmt_int(length(eras)))
  for (i in seq_along(eras)) {
    k <- sprintf("%s|%02d", org, i - 1)
    add("era.span", k, fmt_text(sprintf("%d-%d|%s", eras[[i]]$run$from, eras[[i]]$run$to, eras[[i]]$run$names)))
    add("era.stayed", k, fmt_int(eras[[i]]$stayed))
    add("era.start", k, fmt_int(eras[[i]]$start))
  }
}
```

## Does principal churn predict teacher retention?

This is the claim the whole project started from: that a principal needs about
five years to move a school, and that stability at the top is what drives
everything else. It is testable directly — correlate how often a school changed
principal against how well it held on to its teachers.

Pearson's r is written as an explicit loop rather than with `cor()`, so that all
three implementations add their terms up in the same order and land on the same
bits. Both input series are published alongside the coefficient, so that a
disagreement can be traced to the inputs or to the arithmetic.

```{r correlation}
pearson <- function(xs, ys) {
  n <- min(length(xs), length(ys))
  if (n < 2) return(NULL)
  sx <- 0; sy <- 0
  for (i in seq_len(n)) { sx <- sx + xs[i]; sy <- sy + ys[i] }
  mx <- sx / n; my <- sy / n
  num <- 0; dx <- 0; dy <- 0
  for (i in seq_len(n)) {
    a <- xs[i] - mx; b <- ys[i] - my
    num <- num + a * b; dx <- dx + a * a; dy <- dy + b * b
  }
  if (dx == 0 || dy == 0) return(NULL)
  num / sqrt(dx * dy)
}

rankable <- k8_orgs[!sapply(k8_orgs, function(o) is.null(pooled_teacher[[o]]))]
for (o in rankable) {
  add("correlation.churn", o, fmt_int(changes[[o]]))
  add("correlation.tenure", o, fmt_int(tenure[[o]]))
}
churn_r <- pearson(sapply(rankable, function(o) changes[[o]]),
                   sapply(rankable, function(o) pooled_teacher[[o]]))
tenure_r <- pearson(sapply(rankable, function(o) tenure[[o]]),
                    sapply(rankable, function(o) pooled_teacher[[o]]))
add("correlation.n", "-", fmt_int(length(rankable)))
add("correlation.churnVsRetention", "-", fmt_real(churn_r))
add("correlation.tenureVsRetention", "-", fmt_real(tenure_r))

data.frame(
  measure = c("principal changes vs pooled teacher retention",
              "longest tenure vs pooled teacher retention"),
  r = c(round(churn_r, 3), round(tenure_r, 3)),
  n = length(rankable)
)
```

The claim predicts a clearly negative first row.

## Pooled retention, by school

```{r summary}
summary_tbl <- data.frame(
  org = k8_orgs,
  pooled_15_25 = round(sapply(k8_orgs, function(o) pooled_teacher[[o]]), 4),
  principals = sapply(k8_orgs, function(o) changes[[o]] + 1),
  longest_run = sapply(k8_orgs, function(o) tenure[[o]])
)
summary_tbl[order(-summary_tbl$pooled_15_25), ]
```

## Output

Sorted by metric then key, in byte order. Uppercase sorts before lowercase, and
`radix` is what pins that down regardless of the session's locale.

```{r write}
stopifnot(!any(duplicated(paste0(rows$metric, "\t", rows$key))))
ord <- order(paste0(rows$metric, "\t", rows$key), method = "radix")
lines <- c("metric\tkey\tvalue",
           paste(rows$metric[ord], rows$key[ord], rows$value[ord], sep = "\t"))

out_dir <- file.path(ROOT, "calc", "out")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
con <- file(file.path(out_dir, "r.tsv"), open = "wb")
writeLines(lines, con, sep = "\n", useBytes = TRUE)
close(con)

cat(sprintf("facts: %d\nsha256: %s\n", length(lines) - 1,
            sha256(paste0(paste(lines, collapse = "\n"), "\n"))))
```

The TypeScript and Python implementations print the same hash.