Performing a typed analytical plan

library(conductoR)

The question

Does nitrogen rate affect maize yield? agridat::hernandez.nitrogen holds 136 plots from five Minnesota sites, four replicates each, at nitrogen rates from 0 to 268 kg/ha. Two methods answer the question independently: the slope of a linear regression of yield on rate, and Kendall’s rank correlation. A plan puts a power-aware gate on the rank lens, reconciles the two lenses, and turns the consensus into a recommendation.

data("hernandez.nitrogen", package = "agridat")
d <- hernandez.nitrogen
table(d$site, d$rep)
#>     
#>      R1 R2 R3 R4
#>   S1  7  7  7  7
#>   S2  7  7  7  7
#>   S3  7  7  7  7
#>   S4  7  7  7  7
#>   S5  6  6  6  6

The tools

A tool is a function of (inputs, data, node) that returns an orchestraManifest::orchestra_manifest. lens_to_manifest() lifts a yes-or-no verdict into one; decision_to_manifest() does the same for a recommendation. The tools are registered by name so the plan can be written to a file.

ols_lens <- function(inputs, data, node) {
  fit <- lm(yield ~ nitro, data = data)
  co  <- summary(fit)$coefficients["nitro", ]
  p   <- unname(co["Pr(>|t|)"])
  lens_to_manifest(
    headline = if (p < 0.05) "effect" else "no_effect",
    effect_present = p < 0.05, emitter = "ols", method = "ols:slope_t",
    metrics = list(p_value = p, slope = unname(co["Estimate"]),
                   effect_mean = unname(co["Estimate"]) * diff(range(data$nitro)),
                   noise_sd = sigma(fit), n = nrow(data)))
}

rank_lens <- function(inputs, data, node) {
  ct <- suppressWarnings(cor.test(data$nitro, data$yield, method = "kendall"))
  lens_to_manifest(
    headline = if (ct$p.value < 0.05) "effect" else "no_effect",
    effect_present = ct$p.value < 0.05, emitter = "rank", method = "rank:kendall",
    metrics = list(p_value = ct$p.value, tau = unname(ct$estimate),
                   n = nrow(data)))
}

recommend <- function(inputs, data, node) {
  consensus <- inputs$consensus@summary$headline
  action <- switch(consensus,
    causal_effect_corroborated = "apply_nitrogen",
    no_effect_corroborated     = "withhold_nitrogen",
    "insufficient_evidence")
  decision_to_manifest(action, emitter = "agronomy", method = "agronomy:rule",
                       metrics = list(consensus = consensus),
                       consumed = list(inputs$consensus))
}

register_tool("ols_lens", ols_lens)
register_tool("rank_lens", rank_lens)
register_tool("recommend", recommend)

The plan

Five nodes. Each edge names the inferential target it expects from its source; the gate reads the rank lens and abstains when a null verdict comes from fewer than 25 plots or from a design that could not have resolved the implied effect.

score <- oap_plan(
  id = "nitrogen_response",
  meta = list(question = "Does nitrogen rate affect maize yield?"),
  nodes = list(
    node("lens_ols",  "ols_lens",  emits = "treatment_effects", role = "regression"),
    node("lens_rank", "rank_lens", emits = "treatment_effects", role = "rank"),
    node("gate", operator = "abstain_gate", emits = "decisions",
         consumes = list(edge("lens_rank", "treatment_effects")),
         params = list(n_min = 25L, z_floor = 2, mode = "flag")),
    node("consensus", operator = "triangulate", emits = "decisions",
         consumes = list(edge("lens_ols", "treatment_effects"),
                         edge("lens_rank", "treatment_effects"),
                         edge("gate", "decisions"))),
    node("decision", "recommend", emits = "decisions", role = "agronomy",
         consumes = list(edge("consensus", "decisions")))
  ))
plan_order(score)
#> [1] "lens_ols"  "lens_rank" "gate"      "consensus" "decision"

Scenario A: the full trial

run_a <- conduct(score, data = d)
#>   [compute     ] lens_ols                            -> treatment_effects  
#>   [compute     ] lens_rank                           -> treatment_effects  
#>   [abstain_gate] gate           lens_rank            -> decisions          
#>   [triangulate ] consensus      lens_ols+lens_rank+gate -> decisions          
#>   [compute     ] decision       consensus            -> decisions
run_a$terminal@summary$headline
#> [1] "apply_nitrogen"
run_a$bus$consensus@summary$metrics[c("agree", "n_lenses")]
#> $agree
#> [1] TRUE
#> 
#> $n_lenses
#> [1] 2

Both lenses find the effect, the gate passes, and the consensus is corroborated.

Scenario B: one site, one replicate, three rates

small <- d[d$site == "S1" & d$rep == "R1" & d$nitro <= 67.2, ]
nrow(small)
#> [1] 3
run_b <- conduct(score, data = small)
#>   [compute     ] lens_ols                            -> treatment_effects  
#>   [compute     ] lens_rank                           -> treatment_effects  
#>   [abstain_gate] gate           lens_rank            -> decisions          [ABSTAIN]
#>   [triangulate ] consensus      lens_ols+lens_rank+gate -> decisions          [ABSTAIN]
#>   [compute     ] decision       consensus            -> decisions          [ABSTAIN]
run_b$bus$gate@summary$abstain_reason
#> [1] "small-sample no_effect: n = 3 < 25, so the design cannot rule out the model-implied effect; abstaining"
run_b$terminal@summary$headline
#> [1] "insufficient_evidence"

With three plots the rank test cannot reject, and the gate refuses to let that null through. In flag mode the consensus abstains and the recommendation is insufficient_evidence. Rerouting instead drops the gated lens and reconciles what survives.

gate <- plan_node(score, "gate")
gate@params$mode <- "reroute"
rerouted <- score
rerouted@nodes[[3]] <- gate
run_b2 <- conduct(rerouted, data = small, verbose = FALSE)
run_b2$bus$consensus@summary$metrics[c("rerouted", "n_lenses")]
#> $rerouted
#> [1] TRUE
#> 
#> $n_lenses
#> [1] 1
run_b2$terminal@summary$headline
#> [1] "insufficient_evidence"

One lens on three plots is not a corroborated effect, so the recommendation is still to wait for evidence.

The plan as a file

Because every compute node names its tool, the plan writes to YAML and reads back to the same run.

f <- tempfile(fileext = ".oap.yaml")
plan_to_yaml(score, f)
cat(readLines(f, n = 12), sep = "\n")
#> id: "nitrogen_response"
#> meta:
#>   question: "Does nitrogen rate affect maize yield?"
#> nodes:
#> - id: "lens_ols"
#>   role: "regression"
#>   operator: "compute"
#>   emits: "treatment_effects"
#>   consumes: []
#>   params: []
#>   tool: "ols_lens"
#> - id: "lens_rank"
again <- plan_from_yaml(f)
run_c <- conduct(again, data = d, verbose = FALSE)
identical(run_c$terminal@run_id, run_a$terminal@run_id)
#> [1] TRUE

Writing a plan

A plan can also be written directly as a file, by a person or by an AI model. Version 0.2 of the plan language adds four optional fields that let a plan act on its own results:

Field	What it does
`when`	runs the node only if a condition on upstream results is true
`foreach`	runs the node once per item, from a list or from a column of an upstream result
`on_refusal`	says what an abstaining result does: `continue`, `halt-branch` or `stop`
`operator: select`	keeps the rows of an upstream table that meet a rule

The plan below asks which interactions in a two-way model of yield hold up. The first tool fits yield ~ nitro * site and returns one row per model term. A select node keeps the interaction terms with p < 0.05, and a check runs once for each kept term, but only if there is at least one.

term_table <- function(inputs, data, node) {
  a <- anova(lm(yield ~ nitro * site, data = data))
  terms <- rownames(a)[rownames(a) != "Residuals"]
  tab <- data.frame(term = terms,
                    order = lengths(strsplit(terms, ":", fixed = TRUE)),
                    p_value = a[terms, "Pr(>F)"])
  summ <- orchestraManifest::manifest_summary("terms", metrics = list(n = nrow(data)))
  orchestraManifest::orchestra_manifest(
    emitter_package = "example", emitter_version = "0.1.0",
    inferential_target = "parameters", run_id = "terms-1",
    method = "lm:anova", params = tab, summary = summ, timestamp = Sys.time(),
    data_hash = orchestraManifest::manifest_data_hash(tab, NULL, NULL, NULL,
                                                      NA_integer_, summ))
}

The check is a stand-in for an external validation: it refits the model on replicates R1 and R2, then on R3 and R4, and asks whether the term is significant in both halves. When a half has too few plots it abstains with insufficient_evidence, which on_refusal: continue records and passes on.

check_term <- function(inputs, data, node) {
  term <- node@params$item
  p <- vapply(list(c("R1", "R2"), c("R3", "R4")), function(half) {
    a <- anova(lm(yield ~ nitro * site, data = data[data$rep %in% half, ]))
    a[term, "Pr(>F)"]
  }, numeric(1))
  verdict <- if (anyNA(p)) "insufficient_evidence"
             else if (all(p < 0.05)) "replicated" else "not_replicated"
  decision_to_manifest(verdict, emitter = "split_half",
                       method = paste0("split_half:", term),
                       metrics = list(p_first = p[1], p_second = p[2]))
}

register_tool("term_table", term_table)
register_tool("check_term", check_term)

The plan as YAML. Every string is in quotes: YAML 1.1 reads an unquoted yes, no, on or off as a logical, so conductoR refuses those words unquoted rather than guess.

plan_text <- '
id: "interactions"
plan_version: "0.2"
meta:
  question: "Which nitrogen-by-site interactions hold up in both halves of the trial?"
nodes:
- id: "terms"
  tool: "term_table"
  emits: "parameters"
- id: "flagged"
  operator: "select"
  emits: "parameters"
  consumes:
  - from: "terms"
    contract: "parameters"
  params:
    rule: "order > 1 & p_value < 0.05"
- id: "check"
  tool: "check_term"
  emits: "decisions"
  consumes:
  - from: "flagged"
    contract: "parameters"
  when: "flagged.n_rows > 0"
  foreach:
    from: "flagged"
    column: "term"
  foreach_max: 10
  on_refusal: "continue"
'
interactions <- plan_from_yaml(plan_text)
run_i <- conduct(interactions, data = d)
#>   [compute     ] terms                               -> parameters         
#>   [select      ] flagged        terms                -> parameters         
#>   [compute     ] check          flagged              -> decisions
run_i$conditions
#>    node          condition value
#> 1 check flagged.n_rows > 0  TRUE
run_i$items[, c("item", "status")]
#>         item status
#> 1 nitro:site    ran
run_i$bus$check@outputs[["check[nitro:site]"]]@summary$headline
#> [1] "not_replicated"

One interaction, nitro:site, meets the rule on all 136 plots, and the check finds it significant in the first half but not in the second. The run record lists the condition, its value and every item with its status, so a blinded run and the unblinded run that follows it can be compared item by item.

Conditions are read by conductoR, never run as R code. They may refer to <node>.headline, <node>.abstained, <node>.n_rows and <node>.metrics.<name>, and use plain numbers, quoted strings, comparisons, %in% with c(...), &, |, ! and brackets. Anything else is refused when the plan is read:

res <- try(plan_from_yaml(sub("flagged.n_rows > 0", "system(cmd) == 0",
                              plan_text, fixed = TRUE)), silent = TRUE)
cat(conditionMessage(attr(res, "condition")))
#> The plan has 1 problem:
#> * node 'check', when: `system` is not allowed in a condition. Fix: a condition may use <node>.headline, <node>.abstained, <node>.n_rows, <node>.metrics.<name>, literals and comparisons (== != < <= > >=), %in% with c(...), & | ! and brackets.

The schema and the repair loop

plan_schema() returns a JSON Schema (draft 2020-12) of the language. An editor can check a plan against it as it is typed, and an AI model can use it to write a plan as structured output. The loop is: write the plan as JSON against the schema, check it with plan_validate(), and fix every problem the table lists. plan_validate() reports all the problems at once, each with the field, the node, the rule broken and the fix.

schema <- jsonlite::fromJSON(plan_schema(), simplifyVector = FALSE)
names(schema$`$defs`$node$properties)
#>  [1] "id"          "role"        "tool"        "operator"    "emits"      
#>  [6] "consumes"    "params"      "when"        "foreach"     "foreach_max"
#> [11] "on_refusal"

draft <- '{
  "id": "interactions",
  "plan_version": 0.2,
  "nodes": [
    {"id": "terms", "tool": "term_table", "emits": "parameters"},
    {"id": "flagged", "operator": "select", "emits": "parameters",
     "consumes": [{"from": "terms", "contract": "parameters"}],
     "params": {"rule": "order > 1 & p_value < 5e-2"}},
    {"id": "check", "tool": "check_term",
     "consumes": [{"from": "flagged", "contract": "parameters"}],
     "foreach": {"from": "flagged", "column": "term"},
     "on_refuse": "continue"}
  ]
}'
plan_validate(draft)[, c("node", "field", "fix")]
#>      node        field
#> 1         plan_version
#> 2 flagged  params.rule
#> 3   check    on_refuse
#> 4   check        emits
#>                                                                                                         fix
#> 1                                                                       write plan_version: "0.2" in quotes
#> 2                                                           write numbers in plain decimal form: 0.05, 1000
#> 3 node fields are id, role, tool, operator, emits, consumes, params, when, foreach, foreach_max, on_refusal
#> 4                                                      set emits to the inferential target its tool returns

Each row says what to change. With the version quoted, the number written plainly, emits declared and the field name corrected, the plan is valid, and plan_to_yaml() writes the YAML a person reads.

fixed <- draft
fixed <- sub('"plan_version": 0.2', '"plan_version": "0.2"', fixed, fixed = TRUE)
fixed <- sub("5e-2", "0.05", fixed, fixed = TRUE)
fixed <- sub('"tool": "check_term",', '"tool": "check_term", "emits": "decisions",',
             fixed, fixed = TRUE)
fixed <- sub('"on_refuse"', '"on_refusal"', fixed, fixed = TRUE)
nrow(plan_validate(fixed))
#> [1] 0
cat(plan_to_yaml(plan_from_json(fixed)))
#> id: "interactions"
#> plan_version: "0.2"
#> meta: []
#> nodes:
#> - id: "terms"
#>   role: "compute"
#>   operator: "compute"
#>   emits: "parameters"
#>   consumes: []
#>   params: []
#>   tool: "term_table"
#> - id: "flagged"
#>   role: "compute"
#>   operator: "select"
#>   emits: "parameters"
#>   consumes:
#>   - from: "terms"
#>     contract: "parameters"
#>   params:
#>     rule: "order > 1 & p_value < 0.05"
#> - id: "check"
#>   role: "compute"
#>   operator: "compute"
#>   emits: "decisions"
#>   consumes:
#>   - from: "flagged"
#>     contract: "parameters"
#>   params: []
#>   tool: "check_term"
#>   foreach:
#>     from: "flagged"
#>     column: "term"

A mistyped edge

An edge that declares the wrong contract is refused before the downstream node runs. Here the decision node is told to expect parameters from the consensus, which emits decisions.

bad <- oap_plan(id = "mistyped", nodes = list(
  node("lens_ols", "ols_lens", emits = "treatment_effects"),
  node("decision", "recommend", emits = "decisions",
       consumes = list(edge("lens_ols", "parameters")))
))
res <- try(conduct(bad, data = d, verbose = FALSE), silent = TRUE)
cat(conditionMessage(attr(res, "condition")))
#> TYPED-EDGE VIOLATION -- node 'decision' consumes 'lens_ols' as 'parameters', but 'lens_ols' emits inferential_target 'treatment_effects'. The step is not run.

Provenance and a targets pipeline

Giving conduct() an out_dir writes run.rds, run_card.md and provenance.md, the last carrying the chain of payload hashes in execution order. run_provenance() returns the same table in the session.

run_provenance(run_a)[, c("node", "operator", "emits", "abstained")]
#>        node     operator             emits abstained
#> 1  lens_ols      compute treatment_effects     FALSE
#> 2 lens_rank      compute treatment_effects     FALSE
#> 3      gate abstain_gate         decisions     FALSE
#> 4 consensus  triangulate         decisions     FALSE
#> 5  decision      compute         decisions     FALSE

The same plan compiles to a _targets.R in which each node is one target calling conduct_step(), so a long-running plan inherits caching from the ‘targets’ package.

compile_to_targets(score, "_targets.R", preamble = c(
  "library(conductoR)",
  "source('tools.R')",
  ".score <- plan_from_yaml('nitrogen_response.oap.yaml')",
  ".data <- agridat::hernandez.nitrogen"))
targets::tar_make()
targets::tar_read_raw(terminal_target(score))