--- title: "Philadelphia and Tuscaloosa on a common scale — the calculation" output: md_document --- Pennsylvania and Alabama use different tests with different cut scores, so the two states' published "percent proficient" figures cannot be set side by side. The Stanford Education Data Archive links every state's tests to one national scale, which is what makes a West Philadelphia K-8 and a Tuscaloosa City school comparable at all. Both measures here are deviations from the national average, **in grade levels**, and both can be negative: * **level** (`cs_mn_avg_eb`) — how far above or below the national average a school's students score. Zero is the national average. * **growth** (`cs_mn_grd_eb`) — how much they gain per grade compared with typical. Zero is a typical year's progress. The `_eb` columns are the empirical-Bayes estimates, shrunk toward the mean — the right choice at school level, where small samples make raw estimates noisy. Three places are on this scale: eleven West Philadelphia K-8 schools, the schools of Tuscaloosa City, and the Brooklyn schools around Park Slope and Prospect Heights. New York needs an id crosswalk first — SEDA rolls its 32 geographic districts into one synthetic citywide LEA, where the federal collections keep them apart. Brooklyn is computed twice. The eleven schools around Park Slope were named individually, and a median across schools somebody named is not a median for the borough they sit in, so the whole borough is carried alongside them. One of three implementations; the others are the original TypeScript in `src/` and `calc/seda.py`. All three write the contract in `calc/CONTRACT.md`. ```r ROOT <- Sys.getenv("SDP_ROOT", ".") SEDA <- file.path(ROOT, "seda") # Where the low-income measure runs out of room; see calc/CONTRACT.md. poverty_ceiling <- 0.95 phila_lea <- "4218990" tusc_city_lea <- "103360" tusc_county_lea <- "103390" district_names <- c("4218990" = "Philadelphia City SD", "103360" = "Tuscaloosa City", "103390" = "Tuscaloosa County") # SEDA files all of New York City under one synthetic citywide LEA rather than # the 32 geographic-district LEAs the federal collections use. The five-digit # school number is the same in both, so an id crosswalks by swapping the prefix. nyc_seda_lea <- "3620580" nyc_seda_id <- function(nces) { if (!grepl("^[0-9]{12}$", nces)) stop("nyc_seda_id: not a 12-digit NCES id: ", nces) paste0(nyc_seda_lea, substr(nces, 8, 12)) } nyc <- data.frame( 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") ) k8 <- data.frame( nces = c("421899006528", "421899000619", "421899006530", "421899003782", "421899003783", "421899004732", "421899003649", "421899003650", "421899005186", "421899003805", "421899003618"), label = c("Harrington", "Penn Alexander", "Lea", "Comegys", "Powel", "Alain Locke", "Meredith", "Nebinger", "Greenfield", "Marian Anderson", "Bache-Martin") ) # LEAs the current federal enrolment counts are read for, and CCD's codes for # "all students" and "Black students". demographics_dir <- file.path(ROOT, "demographics") ccd_leas <- c("4218990", "0103360", "3600091", "3600092") ccd_all <- 99 ccd_black <- 2 sources <- list( achievement = list(file = "seda_school_pool_cs_4.1.csv", cols = c("sedasch", "sedaschname", "cs_mn_avg_eb", "cs_mn_grd_eb")), covariates = list(file = "seda_cov_school_pool_4.1.csv", cols = c("sedasch", "sedalea", "schcity", "schnam", "level", "perblk", "totenrl", "perecd")), districts = list(file = "seda2022_dist.csv", cols = c("sedaadmin", "sedaadminname", "subject", "subgroup", "ys_mn19_eb", "ys_mn22_eb", "ys_chg_eb")) ) ``` ## Reading ```r read_seda <- function(rel) { path <- file.path(SEDA, rel) if (!file.exists(path)) stop("missing ", path, " - run ./fetch-seda.sh") 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) } data <- lapply(sources, function(s) read_seda(s$file)) sapply(data, nrow) #> achievement covariates districts #> 82488 89350 32804 ``` ## Reading a SEDA number SEDA writes some values with no leading zero — `.075` rather than `0.075` — and uses `NA` for a suppressed cell. The shape of the string is checked before it is converted, rather than after: R's `as.numeric` and Python's `float` each accept strings that JavaScript's `Number` does not, and the three have to agree about what is not a number as well as about what is. ```r NUMERIC <- "^[+-]?([0-9]+\\.?[0-9]*|\\.[0-9]+)([eE][+-]?[0-9]+)?$" num <- function(v) { t <- trimws(v) ok <- t != "" & toupper(t) != "NA" & grepl(NUMERIC, t) out <- rep(NA_real_, length(t)) out[ok] <- as.numeric(t[ok]) out[!is.finite(out)] <- NA_real_ out } ``` ## Percentiles, and one place three languages disagree A school's percentile is its position in the distribution of **all** schools in the release, not among the few dozen on the page — the point is to place these schools nationally. The rounding needs saying out loud. JavaScript's `Math.round` sends halves up; R's `round` and Python's `round` both send halves to the nearest even number. A percentile that lands exactly on 46.5 would therefore come out 47 in the original and 46 in a naive port. This is the single most likely place for three implementations to disagree by one, so all three do the same explicit thing. ```r js_round <- function(x) floor(x + 0.5) percentile_of <- function(sorted, v) { if (is.na(v) || !length(sorted)) return(NA_integer_) js_round((sum(sorted < v) / length(sorted)) * 100) } 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 } ``` ## Joining achievement to covariates The achievement file carries the two estimates; the covariate file carries the district, the grade span, enrolment and the demographics. Where a school id appears twice in the covariates the last row wins, matching what an assignment does in the other two implementations. ```r ach <- data$achievement cov <- data$covariates cov_last <- !duplicated(cov$sedasch, fromLast = TRUE) cov_id <- cov$sedasch[cov_last] cov_lea <- cov$sedalea[cov_last] cov_city <- cov$schcity[cov_last] cov_name <- cov$schnam[cov_last] cov_level <- cov$level[cov_last] cov_blk <- num(cov$perblk)[cov_last] cov_enr <- num(cov$totenrl)[cov_last] a_level <- num(ach$cs_mn_avg_eb) a_growth <- num(ach$cs_mn_grd_eb) levels_all <- sort(a_level[!is.na(a_level)]) growths_all <- sort(a_growth[!is.na(a_growth)]) join_seda <- function(keep_idx, labels) { hit <- match(ach$sedasch, cov_id) ok <- !is.na(a_level) & !is.na(hit) ok[ok] <- keep_idx(hit[ok]) idx <- which(ok) lab <- labels[ach$sedasch[idx]] # `??` in the original coalesces on absence, not on emptiness: a blank school # name stays blank rather than falling through to the next candidate. nm <- ifelse(is.na(lab), trimws(cov_name[hit[idx]]), lab) out <- data.frame( id = ach$sedasch[idx], name = nm, district = trimws(cov_lea[hit[idx]]), gradeLevel = trimws(cov_level[hit[idx]]), level = a_level[idx], growth = a_growth[idx], pctBlack = cov_blk[hit[idx]], enrollment = cov_enr[hit[idx]] ) out$levelPct <- sapply(out$level, function(v) percentile_of(levels_all, v)) out$growthPct <- sapply(out$growth, function(v) percentile_of(growths_all, v)) # Descending by level. A radix order is stable, as the original's sort is, so # schools with the same level keep their file order. out[order(-out$level, method = "radix"), , drop = FALSE] } labels <- setNames(k8$label, k8$nces) nyc_ids <- vapply(nyc$nces, nyc_seda_id, character(1), USE.NAMES = FALSE) nyc_labels <- setNames(nyc$label, nyc_ids) phila <- join_seda(function(i) cov_id[i] %in% k8$nces, labels) tuscaloosa <- join_seda(function(i) cov_lea[i] == tusc_city_lea, character(0)) brooklyn <- join_seda(function(i) cov_id[i] %in% nyc_ids, nyc_labels) all_brooklyn <- join_seda( function(i) cov_lea[i] == nyc_seda_lea & toupper(cov_city[i]) == "BROOKLYN", character(0)) data.frame(group = c("phila", "tuscaloosa", "brooklyn", "allBrooklyn"), schools = c(nrow(phila), nrow(tuscaloosa), nrow(brooklyn), nrow(all_brooklyn))) #> group schools #> 1 phila 11 #> 2 tuscaloosa 18 #> 3 brooklyn 10 #> 4 allBrooklyn 400 ``` ## What these schools look like now The scores are pooled over 2009-2018, and so is the racial composition SEDA ships beside them. Intakes move: a school that was three quarters Black in that window may be about half now. The current federal count is read alongside, so neither figure is ever published without the other. Only the `sex = 99` rows are kept — the per-sex rows are the same students counted again — and CCD writes suppressed cells as negative numbers, which must not become a negative share. ```r ccd_rows <- list() ccd_digest <- character(0) for (lea in ccd_leas) { path <- file.path(demographics_dir, sprintf("enrollment_%s.json", lea)) if (!file.exists(path)) stop("missing ", path, " - run ./fetch-demographics.sh") ccd_digest[[sprintf("enrollment_%s", lea)]] <- digest::digest(file = path, algo = "sha256") ccd_rows[[lea]] <- jsonlite::fromJSON(path, simplifyDataFrame = TRUE)$results } ccd <- do.call(rbind, ccd_rows) keep <- ccd$sex == 99 & !is.na(ccd$ncessch) & !is.na(ccd$enrollment) & ccd$enrollment >= 0 ccd <- ccd[keep, , drop = FALSE] tot <- ccd[ccd$race == ccd_all, , drop = FALSE] blk <- ccd[ccd$race == ccd_black, , drop = FALSE] now_total <- setNames(tot$enrollment[!duplicated(tot$ncessch, fromLast = TRUE)], tot$ncessch[!duplicated(tot$ncessch, fromLast = TRUE)]) now_black <- setNames(blk$enrollment[!duplicated(blk$ncessch, fromLast = TRUE)], blk$ncessch[!duplicated(blk$ncessch, fromLast = TRUE)]) # SEDA drops Alabama's leading zero and files New York under a synthetic # citywide LEA; both have to be undone to reach the federal counts. nyc_federal <- setNames(nyc$nces, nyc_ids) federal_id <- function(sedasch) { if (!is.na(nyc_federal[sedasch])) return(unname(nyc_federal[sedasch])) if (nchar(sedasch) == 11) paste0("0", sedasch) else sedasch } length(unique(c(names(now_total), names(now_black)))) #> [1] 330 ``` ## Writing the contract ```r 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 for (name in names(sources)) { df <- data[[name]] cols <- 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")))) } for (name in names(ccd_digest)) add("input.filedigest", name, fmt_text(ccd_digest[[name]])) add("input.rows", "ccd_enrollment", fmt_int(sum(sapply(ccd_rows, nrow)))) add("input.schools", "ccd_enrollment", fmt_int(length(unique(c(names(now_total), names(now_black)))))) ``` One line per school, carrying exactly the values the per-school facts do. The borough is four hundred schools, and writing it out would add nearly five thousand lines and bury everything else, so it is checked by hashing this dump instead. One differing field in one school changes the hash. ```r school_dump <- function(df) { lines <- character(nrow(df)) for (i in seq_len(nrow(df))) { fed <- federal_id(df$id[i]) lines[i] <- paste( df$id[i], df$name[i], df$district[i], df$gradeLevel[i], fmt_real(df$level[i]), fmt_real(df$growth[i]), if (is.na(df$levelPct[i])) NA_VALUE else fmt_int(df$levelPct[i]), if (is.na(df$growthPct[i])) NA_VALUE else fmt_int(df$growthPct[i]), fmt_real(df$pctBlack[i]), fmt_real(df$enrollment[i]), if (is.na(now_total[fed])) NA_VALUE else fmt_int(unname(now_total[fed])), if (is.na(now_black[fed])) NA_VALUE else fmt_int(unname(now_black[fed])), fmt_int(i - 1), sep = "\x1f") } paste(lines, collapse = "\n") } emit_group <- function(group, df, per_school) { add("seda.n", group, fmt_int(nrow(df))) # The page draws a median rule on every row and plots growth for whichever # schools have it, so both are contract facts rather than page arithmetic. add("seda.median", group, fmt_real(med(df$level))) add("seda.withGrowth", group, fmt_int(sum(!is.na(df$growth)))) add("seda.digest", group, fmt_text(sha256(school_dump(df)))) if (!per_school) return(invisible(NULL)) for (i in seq_len(nrow(df))) { key <- paste0(group, "|", df$id[i]) add("seda.level", key, fmt_real(df$level[i])) add("seda.growth", key, fmt_real(df$growth[i])) add("seda.levelPct", key, if (is.na(df$levelPct[i])) NA_VALUE else fmt_int(df$levelPct[i])) add("seda.growthPct", key, if (is.na(df$growthPct[i])) NA_VALUE else fmt_int(df$growthPct[i])) add("seda.pctBlack", key, fmt_real(df$pctBlack[i])) add("seda.enrollment", key, fmt_real(df$enrollment[i])) add("seda.name", key, fmt_text(df$name[i])) add("seda.district", key, fmt_text(df$district[i])) add("seda.gradeLevel", key, fmt_text(df$gradeLevel[i])) # The page lists schools by level, descending. Emitting the position checks # the ordering too, not just the values. add("seda.rank", key, fmt_int(i - 1)) fed <- federal_id(df$id[i]) add("seda.nowEnrollment", key, if (is.na(now_total[fed])) NA_VALUE else fmt_int(unname(now_total[fed]))) add("seda.nowBlack", key, if (is.na(now_black[fed])) NA_VALUE else fmt_int(unname(now_black[fed]))) } } emit_group("phila", phila, TRUE) emit_group("tuscaloosa", tuscaloosa, TRUE) emit_group("brooklyn", brooklyn, TRUE) emit_group("allBrooklyn", all_brooklyn, FALSE) add("national.levels", "-", fmt_int(length(levels_all))) add("national.growths", "-", fmt_int(length(growths_all))) # How much of each district SEDA can say anything about. New York is the reason # this is in the contract: it has a growth estimate for about a third of its # schools, against most of Philadelphia's and all of Tuscaloosa City's. for (pl in list(c("Philadelphia", phila_lea), c("Tuscaloosa_City", tusc_city_lea), c("New_York_City", nyc_seda_lea))) { ids <- cov$sedasch[cov$sedalea == pl[2]] m <- ach$sedasch %in% ids add("coverage.schools", pl[1], fmt_int(sum(m))) add("coverage.level", pl[1], fmt_int(sum(m & !is.na(a_level)))) add("coverage.growth", pl[1], fmt_int(sum(m & !is.na(a_growth)))) } # New York is absent from the 2022 district release entirely. add("district.newYorkRows", "-", fmt_int(sum(data$districts$stateabb == "NY"))) ``` ## Districts, and the pandemic drop SEDA's 2022 release is at district level and carries two snapshots — 2019 and 2022 — with the change between them. The national median district gives the drop a reference point, so that a fall reads as bigger or smaller than everyone else's rather than simply as a fall. ```r dist <- data$districts subject_of <- function(s) ifelse(s == "mth", "math", "reading") sel <- dist$subgroup == "all" & dist$sedaadmin %in% names(district_names) d19 <- num(dist$ys_mn19_eb); d22 <- num(dist$ys_mn22_eb); dch <- num(dist$ys_chg_eb) for (i in which(sel)) { key <- paste0(dist$sedaadmin[i], "|", subject_of(dist$subject[i])) add("district.level2019", key, fmt_real(d19[i])) add("district.level2022", key, fmt_real(d22[i])) add("district.change", key, fmt_real(dch[i])) } for (id in names(district_names)) add("district.name", id, fmt_text(district_names[[id]])) for (subj in c("mth", "rla")) { m <- dist$subgroup == "all" & dist$subject == subj add("usMedian.level2019", subject_of(subj), fmt_real(med(d19[m]))) add("usMedian.change", subject_of(subj), fmt_real(med(dch[m]))) # Every district row for the subject, including the ones with no estimate: # this is the denominator the page prints. add("usMedian.n", subject_of(subj), fmt_int(sum(m))) } data.frame( district = district_names[dist$sedaadmin[sel]], subject = subject_of(dist$subject[sel]), level2019 = round(d19[sel], 2), change = round(dch[sel], 2), row.names = NULL ) #> district subject level2019 change #> 1 Tuscaloosa City math -0.43 -0.10 #> 2 Tuscaloosa City reading -0.30 -0.04 #> 3 Tuscaloosa County math -0.38 -0.03 #> 4 Tuscaloosa County reading -0.27 0.03 #> 5 Philadelphia City SD math -0.60 -0.26 #> 6 Philadelphia City SD reading -0.52 -0.08 ``` ## How much of "level" is just intake The page says that level tracks family income closely while growth barely does. That is a claim about the whole release, not about the schools on the page, so it is computed against every school in the file rather than cited. `perecd` is the share of a school's students who are economically disadvantaged. Pearson's r is an explicit loop rather than `cor()`, so all three implementations add their terms in the same order. ```r 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) } perecd <- num(cov$perecd) # Last-wins over the rows that qualify, which is what a repeated assignment does. q <- which(!is.na(cov$sedasch) & !is.na(perecd) & perecd >= 0 & perecd <= 1) pov_id <- cov$sedasch[q] pov_val <- perecd[q] keep <- !duplicated(pov_id, fromLast = TRUE) pov_id <- pov_id[keep]; pov_val <- pov_val[keep] add("poverty.schools", "-", fmt_int(length(pov_id))) # Where `perecd` stops measuring anything. It is built from free-and-reduced- # lunch eligibility, so a district on universal free meals books nearly every # school at the maximum whatever its families earn. The share of each set at the # ceiling is the only way to see which rows have run out of room, so it is a # contract fact rather than something the page works out. add("poverty.ceiling", "-", fmt_int(sum(pov_val >= poverty_ceiling))) add("poverty.ceilingAt", "-", fmt_real(poverty_ceiling)) hit <- match(ach$sedasch, pov_id) correlations <- list() for (nm in c("level", "growth")) { v <- if (nm == "level") a_level else a_growth ok <- !is.na(v) & !is.na(hit) xs <- pov_val[hit[ok]] ys <- v[ok] correlations[[nm]] <- list(r = pearson(xs, ys), n = length(xs)) add(sprintf("poverty.%s.r", nm), "-", fmt_real(correlations[[nm]]$r)) add(sprintf("poverty.%s.n", nm), "-", fmt_int(correlations[[nm]]$n)) } data.frame( measure = c("level vs share economically disadvantaged", "growth vs share economically disadvantaged"), r = c(round(correlations$level$r, 3), round(correlations$growth$r, 3)), n = c(correlations$level$n, correlations$growth$n) ) #> measure r n #> 1 level vs share economically disadvantaged -0.814 72714 #> 2 growth vs share economically disadvantaged -0.118 51016 ``` ## The sets the page compares The intake table reads five sets off one national line. Its numbers were the last part of the page not covered here, which is exactly the part a reader is most likely to want checked. ```r level_of <- setNames(a_level, ach$sedasch) level_of <- level_of[!is.na(level_of)] pov_map <- setNames(pov_val, pov_id) set_stats <- function(name, ids) { levels_ <- unname(level_of[ids]); levels_ <- levels_[!is.na(levels_)] povs <- unname(pov_map[ids]); povs <- povs[!is.na(povs)] add("set.n", name, fmt_int(length(levels_))) add("set.level", name, fmt_real(med(levels_))) add("set.poverty", name, fmt_real(med(povs))) # The denominator is the schools with a share at all — the same set the median # is taken over, not `set.n`, which counts level estimates. add("set.ceiling", name, fmt_int(sum(povs >= poverty_ceiling))) add("set.ceilingOf", name, fmt_int(length(povs))) } cov_ids_where <- function(keep) cov$sedasch[keep] set_stats("philaEleven", k8$nces) set_stats("philaAll", cov_ids_where(cov$sedalea == phila_lea)) set_stats("brooklynTen", nyc_ids) set_stats("brooklynAll", cov_ids_where(cov$sedalea == nyc_seda_lea & toupper(cov$schcity) == "BROOKLYN")) set_stats("tuscaloosaCity", cov_ids_where(cov$sedalea == tusc_city_lea)) ``` Least squares of level on the low-income share, across every US school — the line the table's "predicted" column is read off. Written as an explicit loop for the same reason Pearson's r is. ```r ok_fit <- !is.na(a_level) & !is.na(hit) xs <- pov_val[hit[ok_fit]] ys <- a_level[ok_fit] n_fit <- length(xs) sx <- 0; sy <- 0 for (i in seq_len(n_fit)) { sx <- sx + xs[i]; sy <- sy + ys[i] } mx <- sx / n_fit; my <- sy / n_fit sxy <- 0; sxx <- 0 for (i in seq_len(n_fit)) { dx <- xs[i] - mx sxy <- sxy + dx * (ys[i] - my) sxx <- sxx + dx * dx } slope <- sxy / sxx add("povertyFit.slope", "-", fmt_real(slope)) add("povertyFit.intercept", "-", fmt_real(my - slope * mx)) add("povertyFit.n", "-", fmt_int(n_fit)) data.frame(intercept = round(my - slope * mx, 3), slope = round(slope, 3), n = n_fit) #> intercept slope n #> 1 0.699 -1.334 72714 ``` ## Output ```r 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, "seda-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")))) #> facts: 570 #> sha256: 3c66ba511b07dc26a5871f2195b9bdb835f7f47d9fc3415b3185ccdc7935e669 ``` The TypeScript and Python implementations print the same hash.