Derive sex-specific reference limits from the study's percentiles

1library(yamaa)23yamaa_spec("spec.yaml")$run()$output

1import yamaa23yamaa.spec("spec.yaml").run().output

1# Reference solution for the yamaa benchmark adam-adlb-nonparametric-limits2# (R track).3#4# Written from the benchmark's full prompt and its inputs alone, as an agent5# would write it. Harbor's oracle agent runs it to check that the prompt can6# be solved and that the grader scores a correct answer 1.78adsl <- read.csv("/app/input/adsl.csv", na.strings = "",9                 stringsAsFactors = FALSE)10lb <- read.csv(11  "/app/input/lb.csv",12  colClasses = c(LBSTRESN = "numeric"),13  na.strings = "",14  stringsAsFactors = FALSE15)1617adlb <- merge(lb, adsl[, c("STUDYID", "USUBJID", "SEX")],18              by = c("STUDYID", "USUBJID"), all.x = TRUE, sort = FALSE)19adlb <- adlb[order(adlb$STUDYID, adlb$USUBJID, adlb$LBTESTCD), ]2021# Each parameter and sex's 2.5th and 97.5th percentiles (R's default, type 7)22# and sample variance, a test not done left out.23group <- interaction(adlb$LBTESTCD, adlb$SEX, drop = TRUE)24percentile <- function(p) {25  function(x) {26    if (all(is.na(x))) return(NA_real_)27    unname(quantile(x, p, na.rm = TRUE, type = 7))28  }29}3031out <- data.frame(32  STUDYID = adlb$STUDYID,33  USUBJID = adlb$USUBJID,34  SEX = adlb$SEX,35  PARAMCD = adlb$LBTESTCD,36  AVAL = adlb$LBSTRESN,37  ANRLO = ave(adlb$LBSTRESN, group, FUN = percentile(0.025)),38  ANRHI = ave(adlb$LBSTRESN, group, FUN = percentile(0.975)),39  AVALVAR = ave(adlb$LBSTRESN, group, FUN = function(x) var(x, na.rm = TRUE))40)4142dir.create("/app/output", showWarnings = FALSE)43write.csv(out, "/app/output/adlb.csv", row.names = FALSE, na = "")

1# Reference solution for the yamaa benchmark adam-adlb-nonparametric-limits2# (Python track).3#4# Written from the benchmark's full prompt and its inputs alone, as an agent5# would write it. Harbor's oracle agent runs it to check that the prompt can6# be solved and that the grader scores a correct answer 1.78import csv9import statistics1011ADSL_CSV = "/app/input/adsl.csv"12LB_CSV = "/app/input/lb.csv"13OUTPUT_CSV = "/app/output/adlb.csv"141516def number(text):17    text = text.strip()18    return float(text) if text else None192021def cell(value):22    return "" if value is None else repr(value)232425def type7(values, p):26    """Hyndman and Fan's type 7 quantile of the values at probability p."""27    x = sorted(values)28    h = (len(x) - 1) * p + 129    j = int(h)30    if j >= len(x):31        return x[-1]32    return x[j - 1] + (h - j) * (x[j] - x[j - 1])333435def main():36    with open(ADSL_CSV, newline="") as f:37        sex = {38            (r["STUDYID"], r["USUBJID"]): r["SEX"] for r in csv.DictReader(f)39        }40    with open(LB_CSV, newline="") as f:41        lb = list(csv.DictReader(f))4243    # Each parameter and sex's results, a test not done left out.44    results = {}45    for row in lb:46        value = number(row["LBSTRESN"])47        subject_sex = sex.get((row["STUDYID"], row["USUBJID"]), "")48        group = (row["LBTESTCD"], subject_sex)49        results.setdefault(group, [])50        if value is not None:51            results[group].append(value)5253    summary = {}54    for group, values in results.items():55        low = type7(values, 0.025) if values else None56        high = type7(values, 0.975) if values else None57        variance = statistics.variance(values) if len(values) > 1 else None58        summary[group] = (low, high, variance)5960    fieldnames = [61        "STUDYID", "USUBJID", "SEX", "PARAMCD", "AVAL", "ANRLO", "ANRHI",62        "AVALVAR",63    ]64    with open(OUTPUT_CSV, "w", newline="") as f:65        writer = csv.DictWriter(f, fieldnames=fieldnames)66        writer.writeheader()67        for row in lb:68            subject_sex = sex.get((row["STUDYID"], row["USUBJID"]), "")69            low, high, variance = summary[(row["LBTESTCD"], subject_sex)]70            writer.writerow(71                {72                    "STUDYID": row["STUDYID"],73                    "USUBJID": row["USUBJID"],74                    "SEX": subject_sex,75                    "PARAMCD": row["LBTESTCD"],76                    "AVAL": cell(number(row["LBSTRESN"])),77                    "ANRLO": cell(low),78                    "ANRHI": cell(high),79                    "AVALVAR": cell(variance),80                }81            )828384if __name__ == "__main__":85    main()

1/* opensas reference solution for the yamaa benchmark2   adam-adlb-nonparametric-limits.3   Reads the benchmark inputs and derives the requested output.4   Execute with the existing opensas runtime: opensas /app/output/result.sas.5   The limits follow R's default definition (type 7), which none of the6   percentile definitions of the procedures gives, so a DATA step computes7   them from the sorted values. */89data lb;10  length STUDYID DOMAIN USUBJID LBTESTCD LBTEST LBSTRESU LBSTAT $40;11  infile "/app/input/lb.csv" dsd dlm="," firstobs=2 truncover lrecl=32767;12  input STUDYID : $char40. DOMAIN : $char40. USUBJID : $char40. LBSEQ13        LBTESTCD : $char40. LBTEST : $char40. LBSTRESN LBSTRESU : $char40.14        LBSTAT : $char40.;15run;1617data adsl;18  length STUDYID USUBJID SEX $40;19  infile "/app/input/adsl.csv" dsd dlm="," firstobs=2 truncover lrecl=32767;20  input STUDYID : $char40. USUBJID : $char40. SEX : $char40.;21run;2223/* Each record with its subject's sex. */24proc sql;25  create table withsex as26    select l.STUDYID, l.USUBJID, a.SEX, l.LBTESTCD, l.LBSTRESN27    from lb as l left join adsl as a28      on a.STUDYID = l.STUDYID and a.USUBJID = l.USUBJID;29quit;3031/* The results of each parameter and sex in ascending order, numbered from32   1; a test not done has none and takes no part. */33proc sort data=withsex(where=(not missing(LBSTRESN))) out=sorted;34  by LBTESTCD SEX LBSTRESN;35run;3637data ranked;38  set sorted;39  by LBTESTCD SEX;40  if first.SEX then RANK = 0;41  RANK + 1;42  keep LBTESTCD SEX RANK LBSTRESN;43run;4445proc sql;46  create table counts as47    select LBTESTCD, SEX, count(*) as N, mean(LBSTRESN) as M48    from ranked group by LBTESTCD, SEX;49quit;5051/* Type 7: the quantile at proportion P sits at position H = (N - 1) * P + 1,52   between the values ranked J = floor(H) and K = J + 1, at fraction G. */53data positions;54  set counts;55  STAT = 1; P = 0.025; H = (N - 1) * P + 1; J = floor(H); G = H - J;56  K = min(J + 1, N); output;57  STAT = 2; P = 0.975; H = (N - 1) * P + 1; J = floor(H); G = H - J;58  K = min(J + 1, N); output;59run;6061/* The limits, and the sample variance with divisor N - 1. */62proc sql;63  create table quantiles as64    select p.LBTESTCD, p.SEX, p.STAT,65           lo.LBSTRESN + p.G * (hi.LBSTRESN - lo.LBSTRESN) as Q66    from positions as p67      inner join ranked as lo68        on lo.LBTESTCD = p.LBTESTCD and lo.SEX = p.SEX and lo.RANK = p.J69      inner join ranked as hi70        on hi.LBTESTCD = p.LBTESTCD and hi.SEX = p.SEX and hi.RANK = p.K;71  create table squares as72    select r.LBTESTCD, r.SEX, c.N,73           (r.LBSTRESN - c.M) * (r.LBSTRESN - c.M) as SQ74    from ranked as r inner join counts as c75      on c.LBTESTCD = r.LBTESTCD and c.SEX = r.SEX;76  create table variances as77    select LBTESTCD, SEX,78           case when N > 1 then sum(SQ) / (N - 1) else . end as AVALVAR79    from squares group by LBTESTCD, SEX, N;80  /* A query of one table names its columns bare: opensas reads a table81     alias's qualified names as missing in GROUP BY and inside a summary82     function, which would put every row in one group with no limits. */83  create table limits as84    select LBTESTCD, SEX,85           max(case when STAT = 1 then Q else . end) as ANRLO,86           max(case when STAT = 2 then Q else . end) as ANRHI87    from quantiles group by LBTESTCD, SEX;88  create table final as89    select w.STUDYID, w.USUBJID, w.SEX, w.LBTESTCD as PARAMCD length=40,90           w.LBSTRESN as AVAL, s.ANRLO, s.ANRHI, v.AVALVAR91    from withsex as w92      left join limits as s on s.LBTESTCD = w.LBTESTCD and s.SEX = w.SEX93      left join variances as v on v.LBTESTCD = w.LBTESTCD and v.SEX = w.SEX94    order by 1, 2, 4;95quit;9697proc export data=final outfile="/app/output/adlb.csv" dbms=csv replace; run;