# ED Atlas model files This dataset holds the Bayesian logistic regressions of the ED Atlas (30-day mortality and ICD-10 diagnoses, given sex, age and presenting complaint). The files are NetCDF (`.nc`), so any language with a NetCDF/HDF5 reader works. See the GitHub repository (https://github.com/npielawski/ed-atlas) for the full details. This readme explains the logic, and should support any other implementation. ## Files | File | Content | |---|---| | `trace.nc` | mortality model | | `trace_base.nc` | mortality base model (sex and age only; used for relative risk) | | `traces_diag_.zip` | the diagnosis models whose code starts with the letter `L`, as `trace_diag_.nc` files (one per ICD-10 chapter, block or three-character code). `traces_diag_num.zip` holds the numeric codes. Unzip before use | | `columns_mortality.csv` | column order of the mortality models (no header, one name per line) | | `columns_diag.csv` | column order shared by all diagnosis models (same format) | Each `.nc` file has a group `posterior` with one variable `betas` of shape `(chain, draw, dims)`: the posterior samples of the regression coefficients. `dims` has one entry per coefficient. A prediction is `sigmoid(row · betas)` for one design row, evaluated for every sample. Every model is `logit(p) = row · betas`, with `row` built from sex, age and complaint as below. ## 1. Read the file Python (ArviZ), the same idea works with any NetCDF4/HDF5 reader (R `ncdf4`/`tidync`, Julia `NCDatasets`, ...): ```python import arviz as az betas = az.from_netcdf("trace.nc").posterior.betas.to_numpy() # (chain, draw, n_coef) betas = betas.reshape(-1, betas.shape[-1]) # (n_samples, n_coef) ``` Reading raw NetCDF: open group `posterior`, read variable `betas`. ## 2. Inputs | Input | Coding | |---|---| | sex | `1` male, `0` female | | age | years, number. Trained on adults (18 and up) | | complaint | one SNOMED CT concept id, e.g. `29857009` (chest pain), or `nosct:other` | The complaint must be one of the complaints the model knows. The list is in the `columns_*.csv` file of the model (`columns_mortality.csv` or `columns_diag.csv`): it is every line that starts with `symptoms_`. Complaint number `j` is the `j`-th such line, counting from 0, in file order. Example: the line `symptoms_29857009` is complaint `29857009`. The line `symptoms_nosct:other` is the catch-all. Notes: - The first two lines of `columns_*.csv` are `ed_age` and `is_male`. - All diagnosis models share `columns_diag.csv`. The mortality models use `columns_mortality.csv`, whose complaint list differs (two extra complaints). Always look the complaint up in the file that belongs to the model you use. - The models only take SNOMED CT ids, so any free text must be mapped to ids. - The model expects exactly one complaint selected per patient, not several. ## 3. Build the design row The design matrix was built using the python `patsy` package, very similar to R formulas. It was built using the formula `sex * cr(age) * symptoms`. The coordinate labels of `dims` in the `.nc` file spell it out: `sex`, `cr(...)[m]`, `symptoms[j]`, `sex:cr(...)[m]:symptoms[j]`, ... If `len(dims)` is not `10 + 10*S`, stop and check the columns file. For the `_base` model there are no complaint terms: the row is just `row[0..9]`. `cr(age)` means that the age is transformed to a natural cubic spline. See below for details. Let `S` be the number of complaints in the model (the number of `symptoms_` lines in its `columns_*.csv`), `j` the index of the chosen complaint, `b = B(age)` (4 numbers), `x = sex` (0 or 1). The row has `10 + 10*S` entries, all zero except: ``` row[0] = 1 # intercept row[1] = x # sex row[2..5] = b # age row[6..9] = x * b # sex * age row[10 + j] = 1 # complaint row[10 + S + j] = x # sex * complaint row[10 + 2S + 4j .. +3] = b # age * complaint row[10 + 6S + 4j .. +3] = x * b # sex * age * complaint ``` ### The age spline `B(age)` Age enters through a natural cubic regression spline with 4 knots, fixed for all models: ``` knots = [22, 43, 65, 88] # lower bound, two inner knots, upper bound ``` This is the `cr` spline of patsy/mgcv with those knots. It maps an age to 4 numbers `B[0..3]`. (They always sum to 1; below 22 and above 88 the spline continues as a straight line.) Pseudocode: ``` t = [22, 43, 65, 88]; K = 4; h[i] = t[i+1] - t[i] # i = 0..2 # One-time setup: K x K matrix F, with zero first and last rows Bm[i][i] = (h[i] + h[i+1]) / 3 # i = 0..1 Bm[i][i+1] = Bm[i+1][i] = h[i+1] / 6 # i = 0 Dm[i][i] = 1 / h[i] Dm[i][i+1] = -(1/h[i] + 1/h[i+1]) Dm[i][i+2] = 1 / h[i+1] F[1..2][*] = solve(Bm, Dm) # rows 1,2 of F; rows 0,3 are zero # Value of the 4 basis functions at x inside interval j (t[j] <= x <= t[j+1]): function segment(x, j): m = t[j+1] - x; p = x - t[j] a_m = m / h[j]; a_p = p / h[j] c_m = (m^3 / h[j] - h[j] * m) / 6 c_p = (p^3 / h[j] - h[j] * p) / 6 r = c_m * F[j] + c_p * F[j+1] # vectors of length 4 r[j] += a_m r[j+1] += a_p return r # Slope of the basis at the two ends (vectors of length 4): slope_lo = (-h[0]/3) * F[0] + (-h[0]/6) * F[1]; slope_lo[0] += -1/h[0]; slope_lo[1] += 1/h[0] slope_hi = ( h[2]/6) * F[2] + ( h[2]/3) * F[3]; slope_hi[2] += -1/h[2]; slope_hi[3] += 1/h[2] function B(x): if x < 22: return segment(22, 0) + (x - 22) * slope_lo if x > 88: return segment(88, 2) + (x - 88) * slope_hi j = 0 if x < 43, 1 if x < 65, else 2 return segment(x, j) ``` Test values, to check your implementation (`B(age)`): | age | B[0] | B[1] | B[2] | B[3] | |---|---|---|---|---| | 18 | 1.240090 | -0.299026 | 0.070009 | -0.011073 | | 22 | 1 | 0 | 0 | 0 | | 30 | 0.534220 | 0.566545 | -0.119698 | 0.018933 | | 50 | -0.084319 | 0.796394 | 0.335970 | -0.048045 | | 65 | 0 | 0 | 1 | 0 | | 88 | 0 | 0 | 0 | 1 | | 100 | -0.039849 | 0.226578 | -0.850696 | 1.663966 | ## 4. Predict ``` logit = betas @ row # one value per posterior sample p = 1 / (1 + exp(-logit)) # one probability per sample estimate = mean(p) # point estimate interval = [quantile(p, 0.025), quantile(p, 0.975)] # 95% credible interval ``` The mean is taken over the probabilities, not over the logits. ## Worked check Verify your model predictions against https://ed-atlas.org/ to make sure you have implemented things correctly.