---
title: "School attendance, Philadelphia and Brooklyn — 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.
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", error = FALSE)
options(stringsAsFactors = FALSE)
```

Philadelphia and New York both publish school attendance, and they publish
mirrored halves of the same fact: Philadelphia the share of students attending
**at least 90%** of days, New York the share **chronically absent** — missing
10% or more. This document normalizes everything to the second, which is the
federal standard metric, and then reaches for the one source that measures all
three cities on one definition in one year.

It is one of three implementations of this calculation. The others are the
original TypeScript in `src/` and `calc/attendance.py`; all three write the
contract described in `calc/CONTRACT.md` and the three files are compared byte
for byte.

```{r constants}
ROOT <- Sys.getenv("SDP_ROOT", ".")
A <- file.path(ROOT, "attendance")

nyc_years <- c("2016-17", "2017-18", "2018-19", "2019-20", "2020-21")
crdc_leas <- c("4218990", "0103360", "3600091", "3600092")
# 2017-18 is the last clean pre-pandemic collection; 2022-23 the most recent.
# The 2020-21 collection covers the remote year and is left out.
crdc_years <- c("2017", "2022")
tuscaloosa_seda_lea <- "103360"

k8 <- data.frame(
  org = c("1300", "1280", "1340", "1260", "1390", "1470", "2380", "2590", "2470", "2480", "2210"),
  label = c("Harrington", "Penn Alexander", "Lea", "Comegys", "Powel", "Alain Locke",
            "Meredith", "Nebinger", "Greenfield", "Marian Anderson", "Bache-Martin"),
  nces = c("421899006528", "421899000619", "421899006530", "421899003782", "421899003783",
           "421899004732", "421899003649", "421899003650", "421899005186", "421899003805",
           "421899003618")
)

nyc_schools <- data.frame(
  dbn = c("13K011", "13K009", "15K321", "15K107", "15K039", "15K051",
          "15K154", "15K010", "15K118", "13K133", "13K282"),
  nces = c("360009102074", "360009102069", "360009202818", "360009202403", "360009202175",
           "360009202226", "360009202557", "360009202072", "360009206406", "360009102490",
           "360009102793"),
  label = c("PS 11 Purvis J. Behan", "PS 9 Teunis G. Bergen", "PS 321 William Penn",
            "PS 107 John W. Kimball", "PS 39 Henry Bristow", "MS 51 William Alexander",
            "PS 154 Windsor Terrace", "PS 10 Magnet School", "PS 118 Maurice Sendak",
            "PS 133 William A. Butler", "PS 282 Park Slope")
)

# Each source, the short name it goes by in the contract, and the columns every
# implementation reads out of it. The digest covers those columns only: hashing
# a column nothing uses would fail the check for noise rather than for signal.
csv_sources <- list(
  ada_archive = list(file = "ex/Student_ADA_Yearly.csv",
                     cols = c("School Year", "ULCS Code", "Average Daily Attendance")),
  ada_current = list(file = "Student_ADA_School.csv",
                     cols = c("School Year", "ULCS Code", "Average Daily Attendance (YTD)")),
  attending90 = list(file = "Student_Attendance_90_School.csv",
                     cols = c("School Year", "ULCS Code", "Category", "% with 90%+ Attendance (Yearly)")),
  attending95_archive = list(file = "ex/SDP_95_Attendance_School_S.csv",
                             cols = c("School Year", "ULCS Code", "Category", "% with 95%+ Attendance (Yearly)")),
  attending95_current = list(file = "Student_Attendance_95_School.csv",
                             cols = c("School Year", "ULCS Code", "Category", "% with 95%+ Attendance (Yearly)")),
  teacher_ada = list(file = "ex/Archive_Teacher_95_ADA_2012-13_2022-23/OD_Teacher_Attendance_Schools_Suppressed_20240328.csv",
                     cols = c("school_year", "ulcscode", "pct_att_ada")),
  seda_covariates = list(file = "../seda/seda_cov_school_pool_4.1.csv",
                         cols = c("sedasch", "sedalea", "schnam"))
)
```

## Reading the files

The Philadelphia files are latin-1, and one of them carries a byte-order mark.
Left in place, that mark becomes part of the first column's name, so every
lookup against that column silently returns nothing — a failure that shows up as
an empty chart rather than as an error. It is stripped from the raw bytes,
before any decoding, which is the only place it is unambiguous.

```{r read}
read_extract <- function(rel) {
  path <- file.path(A, rel)
  if (!file.exists(path)) stop("missing ", path, " - run ./fetch-attendance.sh")
  raw <- readBin(path, "raw", file.size(path))
  # UTF-8 byte-order mark: EF BB BF.
  if (length(raw) >= 3 && identical(as.integer(raw[1:3]), c(239L, 187L, 191L))) raw <- raw[-(1:3)]
  txt <- iconv(rawToChar(raw), from = "latin1", to = "UTF-8")
  read.csv(text = txt, colClasses = "character", check.names = FALSE, na.strings = NULL)
}

sources <- lapply(csv_sources, function(s) read_extract(s$file))
sapply(sources, nrow)
```

## What a published cell means

Suppressed cells are published as `s`, not as a blank, and a school with too few
students to report is not a school with zero attendance. Anything that is not a
number is dropped rather than coerced.

```{r parse}
parse_value <- function(v) {
  t <- trimws(v)
  out <- rep(NA_real_, length(t))
  ok <- !(t == "" | tolower(t) == "s" | tolower(t) == "n/a")
  out[ok] <- suppressWarnings(as.numeric(t[ok]))
  out[!is.finite(out)] <- NA_real_
  out
}

# SDP school codes appear both zero-padded and unpadded across files.
normalize_ulcs <- function(code) {
  t <- sub("^0+", "", trimws(code))
  ifelse(t == "", "0", t)
}
```

## One series per school

Every Philadelphia file is in long format — one row per school per year, and in
some files per demographic group as well. They collapse to the same shape: a
school, a year, a number. Where a school-year appears twice, the last row in the
file wins, which is what an assignment does in the other two implementations.

```{r series}
series_from_rows <- function(df, year_col, code_col, value_col, keep = NULL) {
  if (!is.null(keep)) df <- df[keep(df), , drop = FALSE]
  v <- parse_value(df[[value_col]])
  code <- normalize_ulcs(df[[code_col]])
  year <- trimws(df[[year_col]])
  ok <- !is.na(v) & code != "" & year != ""
  v <- v[ok]; code <- code[ok]; year <- year[ok]

  out <- list()
  for (cd in split(seq_along(code), code)) {
    yrs <- year[cd]
    vals <- v[cd]
    last <- !duplicated(yrs, fromLast = TRUE)
    out[[code[cd[1]]]] <- setNames(vals[last], yrs[last])
  }
  out
}

all_students <- function(df) df$Category == "All Students"

# Two releases of the same series; the newer file wins year by year.
merge_series <- function(old, new) {
  out <- old
  for (cd in names(new)) {
    have <- out[[cd]]
    add <- new[[cd]]
    if (is.null(have)) out[[cd]] <- add
    else out[[cd]] <- c(add, have[setdiff(names(have), names(add))])
  }
  out
}
```

```{r philadelphia}
phl_ada <- merge_series(
  series_from_rows(sources$ada_archive, "School Year", "ULCS Code", "Average Daily Attendance"),
  series_from_rows(sources$ada_current, "School Year", "ULCS Code", "Average Daily Attendance (YTD)")
)

attending90 <- series_from_rows(sources$attending90, "School Year", "ULCS Code",
                                "% with 90%+ Attendance (Yearly)", all_students)
# Philadelphia publishes the attending share; the standard metric is its complement.
phl_chronic <- lapply(attending90, function(y) 100 - y)

phl_95 <- merge_series(
  series_from_rows(sources$attending95_archive, "School Year", "ULCS Code",
                   "% with 95%+ Attendance (Yearly)", all_students),
  series_from_rows(sources$attending95_current, "School Year", "ULCS Code",
                   "% with 95%+ Attendance (Yearly)", all_students)
)

teacher_ada <- series_from_rows(sources$teacher_ada, "school_year", "ulcscode", "pct_att_ada")

length(phl_ada)
```

## Reference lines

A school's absence rate means little without knowing what the district looked
like that year, so each measure also gets a districtwide median across every
school SDP reports — not across the eleven on the page.

```{r medians}
med <- function(v) {
  v <- v[!is.na(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
}

district_medians <- function(series) {
  years <- unlist(lapply(series, names), use.names = FALSE)
  vals <- unlist(series, use.names = FALSE)
  split(vals, years)
}
```

## New York

New York publishes one row per school per year, already as a chronic-absence
percentage, so nothing has to be converted — only read.

```{r newyork}
nyc_raw <- lapply(nyc_years, function(y) {
  path <- file.path(A, sprintf("nyc_attendance_%s.json", y))
  if (!file.exists(path)) stop("missing ", path, " - run ./fetch-attendance.sh")
  jsonlite::fromJSON(path, simplifyDataFrame = TRUE)
})
names(nyc_raw) <- nyc_years

nyc_by_year <- lapply(nyc_raw, function(df) {
  ada <- parse_value(df$attendance)
  chronic <- parse_value(df$chronically_absent_1)
  ok <- !is.na(df$dbn) & !is.na(ada) & !is.na(chronic)
  keep <- ok & !duplicated(df$dbn, fromLast = TRUE)
  list(dbn = df$dbn[keep], ada = ada[keep], chronic = chronic[keep])
})

sapply(nyc_by_year, function(x) length(x$dbn))
```

## One year, three cities, one definition

Philadelphia, New York and Tuscaloosa each publish attendance their own way,
over different windows, so their own figures cannot be set side by side. The
federal Civil Rights Data Collection asks every district for the same count in
the same year, which makes it the one apples-to-apples measure available.

It publishes counts rather than rates, so enrolment comes from the Common Core
of Data as the denominator — and it returns **two rows per school** even with
every other subgroup pinned to "all": one for all students and one for the
homeless subgroup. Keeping the wrong one produces single-digit rates that look
entirely plausible.

```{r crdc}
read_results <- function(rel) {
  path <- file.path(A, rel)
  if (!file.exists(path)) stop("missing ", path, " - run ./fetch-crdc.sh")
  jsonlite::fromJSON(path, simplifyDataFrame = TRUE)$results
}

read_crdc_year <- function(year) {
  absent <- list(); enrolled <- list()
  for (lea in crdc_leas) {
    a <- read_results(sprintf("crdc_absent_%s_%s.json", year, lea))
    if (!is.null(a) && nrow(a)) {
      keep <- (is.na(a$homeless) | a$homeless == 99) & !is.na(a$ncessch) &
        !is.na(a$students_chronically_absent)
      a <- a[keep, , drop = FALSE]
      last <- !duplicated(a$ncessch, fromLast = TRUE)
      absent[a$ncessch[last]] <- a$students_chronically_absent[last]
    }
    e <- read_results(sprintf("crdc_enrolled_%s_%s.json", year, lea))
    if (!is.null(e) && nrow(e)) {
      e <- e[!is.na(e$ncessch) & !is.na(e$enrollment), , drop = FALSE]
      last <- !duplicated(e$ncessch, fromLast = TRUE)
      enrolled[e$ncessch[last]] <- e$enrollment[last]
    }
  }
  list(absent = absent, enrolled = enrolled)
}

crdc <- lapply(crdc_years, read_crdc_year)
names(crdc) <- crdc_years

# CRDC encodes suppressed or not-applicable cells as negative numbers (-3, -5,
# -9). Those must not be divided; nor may a rate above 100 be reported, which
# would mean the count and the denominator came from mismatched records.
federal_rate <- function(absent, enrolled) {
  if (is.null(absent) || is.na(absent) || absent < 0) return(NULL)
  if (is.null(enrolled) || is.na(enrolled) || enrolled <= 0) return(NULL)
  rate <- (absent / enrolled) * 100
  if (rate > 100) NULL else rate
}

sapply(crdc, function(x) length(x$absent))
```

## Writing the contract

```{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) == 0 || is.na(x)) return(NA_VALUE)
  if (!is.finite(x)) stop("contract: ", x, " is not finite")
  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

```{r facts-input}
for (name in names(csv_sources)) {
  df <- sources[[name]]
  cols <- csv_sources[[name]]$cols
  dump <- do.call(paste, c(lapply(cols, function(c) df[[c]]), sep = "\x1f"))
  add("input.rows", name, fmt_int(nrow(df)))
  add("input.digest", name, fmt_text(sha256(paste(dump, collapse = "\n"))))
}

# The JSON sources are hashed as raw bytes rather than as parsed columns: there
# is no cross-language question about how to read a JSON string, and what each
# implementation makes of the numbers is checked by the facts below.
json_files <- c(
  setNames(sprintf("nyc_attendance_%s.json", nyc_years), sprintf("nyc_%s", nyc_years)),
  setNames(as.vector(outer(crdc_years, crdc_leas, function(y, l) sprintf("crdc_absent_%s_%s.json", y, l))),
           as.vector(outer(crdc_years, crdc_leas, function(y, l) sprintf("crdc_absent_%s_%s", y, l)))),
  setNames(as.vector(outer(crdc_years, crdc_leas, function(y, l) sprintf("crdc_enrolled_%s_%s.json", y, l))),
           as.vector(outer(crdc_years, crdc_leas, function(y, l) sprintf("crdc_enrolled_%s_%s", y, l))))
)
for (name in names(json_files)) {
  add("input.filedigest", name,
      fmt_text(digest::digest(file = file.path(A, json_files[[name]]), algo = "sha256")))
}
```

```{r facts-philadelphia}
series <- list(ada = phl_ada, chronic = phl_chronic, a95 = phl_95, teacherAda = teacher_ada)

for (i in seq_len(nrow(k8))) {
  org <- k8$org[i]
  code <- normalize_ulcs(org)
  for (name in names(series)) {
    v <- series[[name]][[code]]
    if (is.null(v)) next
    for (y in names(v)) add(sprintf("phl.%s", name), paste0(org, "|", y), fmt_real(v[[y]]))
  }
}

for (name in names(series)) {
  if (name == "teacherAda") next
  by_year <- district_medians(series[[name]])
  for (y in names(by_year)) {
    add(sprintf("phlMedian.%s", name), y, fmt_real(med(by_year[[y]])))
    add(sprintf("phlCount.%s", name), y, fmt_int(length(by_year[[y]])))
  }
}
```

```{r facts-newyork}
for (i in seq_len(nrow(nyc_schools))) {
  dbn <- nyc_schools$dbn[i]
  for (y in nyc_years) {
    hit <- match(dbn, nyc_by_year[[y]]$dbn)
    if (is.na(hit)) next
    add("nyc.ada", paste0(dbn, "|", y), fmt_real(nyc_by_year[[y]]$ada[hit]))
    add("nyc.chronic", paste0(dbn, "|", y), fmt_real(nyc_by_year[[y]]$chronic[hit]))
  }
}

for (y in nyc_years) {
  add("nycMedian.ada", y, fmt_real(med(nyc_by_year[[y]]$ada)))
  add("nycMedian.chronic", y, fmt_real(med(nyc_by_year[[y]]$chronic)))
  add("nycMedian.n", y, fmt_int(length(nyc_by_year[[y]]$dbn)))
}
```

```{r facts-federal}
target_id <- c(k8$nces, nyc_schools$nces)
target_city <- c(rep("Philadelphia", nrow(k8)), rep("Brooklyn", nrow(nyc_schools)))
target_label <- c(k8$label, nyc_schools$label)

# SEDA drops Alabama's leading zero from the NCES id; the federal collections
# keep all twelve digits.
cov <- sources$seda_covariates
tusc <- cov[cov$sedalea == tuscaloosa_seda_lea & cov$sedasch != "", , drop = FALSE]
tid <- ifelse(nchar(tusc$sedasch) == 11, paste0("0", tusc$sedasch), tusc$sedasch)
target_id <- c(target_id, tid)
target_city <- c(target_city, rep("Tuscaloosa City", length(tid)))
target_label <- c(target_label, trimws(tusc$schnam))

# A repeated id is one entry, resolved to the last one seen.
last <- !duplicated(target_id, fromLast = TRUE)
target_id <- target_id[last]; target_city <- target_city[last]; target_label <- target_label[last]

for (i in seq_along(target_id)) {
  add("federal.city", target_id[i], fmt_text(gsub(" ", "_", target_city[i])))
  add("federal.label", target_id[i], fmt_text(target_label[i]))
}

city_rates <- list()
for (year in crdc_years) {
  rates_by_city <- list()
  for (i in seq_along(target_id)) {
    rate <- federal_rate(crdc[[year]]$absent[[target_id[i]]], crdc[[year]]$enrolled[[target_id[i]]])
    if (is.null(rate)) next
    add("federal.rate", paste0(year, "|", target_id[i]), fmt_real(rate))
    rates_by_city[[target_city[i]]] <- c(rates_by_city[[target_city[i]]], rate)
  }
  for (city in names(rates_by_city)) {
    key <- paste0(year, "|", gsub(" ", "_", city))
    add("federalMedian.median", key, fmt_real(med(rates_by_city[[city]])))
    add("federalMedian.min", key, fmt_real(min(rates_by_city[[city]])))
    add("federalMedian.max", key, fmt_real(max(rates_by_city[[city]])))
    add("federalMedian.n", key, fmt_int(length(rates_by_city[[city]])))
  }
  city_rates[[year]] <- rates_by_city
}

do.call(rbind, lapply(crdc_years, function(y) data.frame(
  year = y,
  city = names(city_rates[[y]]),
  n = sapply(city_rates[[y]], length),
  median = round(sapply(city_rates[[y]], med), 1),
  row.names = NULL)))
```

## Does a school's absence rate track its teachers' attendance?

The two are measured independently — students by the district's own attendance
system, teachers by payroll — so a correlation between them is not an artefact
of one measure feeding the other.

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.

```{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)
}

chronics <- numeric(0); teachers <- numeric(0)
for (i in seq_len(nrow(k8))) {
  org <- k8$org[i]
  code <- normalize_ulcs(org)
  ch <- phl_chronic[[code]]["2024-2025"]
  ta <- teacher_ada[[code]]
  # Year keys are compared as strings, exactly as the original does; they are
  # all ASCII digits, so every language orders them the same way.
  tv <- if (is.null(ta)) NULL else med(ta[names(ta) >= "2016"])
  if (is.null(ch) || is.na(ch) || is.null(tv)) next
  chronics <- c(chronics, ch); teachers <- c(teachers, tv)
  add("correlation.chronic", org, fmt_real(unname(ch)))
  add("correlation.teacher", org, fmt_real(tv))
}
r <- pearson(chronics, teachers)
add("correlation.n", "-", fmt_int(length(chronics)))
add("correlation.chronicVsTeacherAttendance", "-", fmt_real(r))

data.frame(measure = "chronic absence vs teacher attendance",
           r = round(r, 3), n = length(chronics))
```

## Output

```{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, "attendance-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.