Doing Analysis With a Model

“Vibe science” is the practice --- increasingly common and rarely written down --- of letting a language model carry substantial parts of the scientific workflow: proposing hypotheses, writing the analysis code, summarizing a literature, interpreting a table of results, and drafting the paragraph that reports them. The name is deliberately unflattering, and it points at the failure mode rather than the tool: the output reads plausibly, the reader’s sense that it is right substitutes for a check that it is right, and no step of the process leaves the kind of record that would let anyone establish afterwards what was actually done. A statistician has an unusually good vocabulary for what goes wrong here, because every one of the problems has a classical name. Iterating on a prompt until the answer looks convincing is the garden of forking paths. Treating a model’s fluent summary of a literature as evidence is measurement without validation. Grading model output with another model is a rater with unknown and uncalibrated bias. Reporting one sampled completion is reporting one draw from a distribution as if it were a point estimate. The useful posture is neither refusal nor enthusiasm but instrumentation: decide in advance what the model is allowed to decide, keep the prompt and the seed with the result, and make every model-generated claim answer to a check that does not itself come from a model. This section works through where these assistants genuinely help, where they quietly manufacture confidence, and what a defensible protocol for using one looks like.

1Prompt iteration is the garden of forking paths¶

The characteristic move of vibe science is to run an analysis, dislike the answer, rephrase the request, and keep the version that reads well. Each rephrasing is a fork, and none of them is logged. The simulation below fixes a null effect and lets an assistant choose among three plausible outcomes and two plausible subgroups.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
R, n = 4000, 60
x = rng.integers(0, 2, (R, n))          # treatment; truly inert
Y = rng.normal(size=(R, n, 3))          # three plausible outcomes
sub = rng.integers(0, 2, (R, n))        # one plausible subgroup split

def pval(y, keep, x):                   # two-sample t on the kept rows
    w1, w0 = keep & (x == 1), keep & (x == 0)
    n1, n0 = w1.sum(1), w0.sum(1)
    m1, m0 = (y * w1).sum(1) / n1, (y * w0).sum(1) / n0
    ss = ((w1 * (y - m1[:, None]) ** 2).sum(1)
          + (w0 * (y - m0[:, None]) ** 2).sum(1))
    s = ss / (n1 + n0 - 2)
    t = np.abs(m1 - m0) / np.sqrt(s * (1 / n1 + 1 / n0))
    return 2 * stats.t.sf(t, n1 + n0 - 2)

masks = (np.ones((R, n), bool), sub == 0, sub == 1)
P = np.column_stack([pval(Y[:, :, j], k, x) for j in range(3) for k in masks])
print(f"the analysis you would prespecify   Pr(p<.05) = "
      f"{np.mean(P[:, 0] < 0.05):.3f}")
print(f"best of 9 outcome-by-subgroup forks Pr(p<.05) = "
      f"{np.mean(P.min(1) < 0.05):.3f}")
print(f"independent-forks bound 1-(1-.05)^9 = {1 - 0.95 ** 9:.3f}")

the analysis you would prespecify   Pr(p<.05) = 0.054
best of 9 outcome-by-subgroup forks Pr(p<.05) = 0.309
independent-forks bound 1-(1-.05)^9 = 0.370

The prespecified analysis holds its nominal size; the best of nine forks rejects almost a third of the time. The garden here is unusually treacherous because the forks are cheap, the discarded branches leave no trace, and no correction is even computable after the fact --- you cannot Bonferroni over a number of attempts nobody counted.

1.1One completion is one draw, not a point estimate¶

The same request issued twice can return different analyses, because the assistant is sampling over specifications as well as over words. Reporting the one you received is reporting a single Monte Carlo draw as though it were the estimate.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 120
x = rng.integers(0, 2, n)
y = 0.35 * x + rng.standard_t(3, n)     # heavy tails invite "cleaning"

SPECS = {"as collected": np.ones(n, bool),
         "drop |y|>3": np.abs(y) <= 3,
         "drop |y|>2": np.abs(y) <= 2,
         "every 10th row dropped": np.arange(n) % 10 != 0}
for name, m in SPECS.items():
    d = y[m & (x == 1)].mean() - y[m & (x == 0)].mean()
    p = stats.ttest_ind(y[m & (x == 1)], y[m & (x == 0)]).pvalue
    print(f"{name:<23} n={m.sum():3d}  diff {d:+.3f}  p = {p:.3f}")

draws = rng.choice(list(SPECS), size=2000)
sig = np.array([stats.ttest_ind(y[SPECS[d] & (x == 1)],
                                y[SPECS[d] & (x == 0)]).pvalue < 0.05
                for d in draws])
print(f"\nPr(assistant reports p<.05) over its choice of spec: "
      f"{sig.mean():.3f}")

as collected            n=120  diff +0.523  p = 0.056
drop |y|>3              n=111  diff +0.468  p = 0.017
drop |y|>2              n=105  diff +0.301  p = 0.074
every 10th row dropped  n=108  diff +0.478  p = 0.087

Pr(assistant reports p<.05) over its choice of spec: 0.243

The four specifications are all defensible and they straddle $p = 0.05$. The last line is the honest summary of the procedure: the probability that this workflow hands you significance is a property of the assistant’s sampling distribution, not of the data. Re-running with a different seed and reporting the spread is the cheapest sensitivity analysis available, and it is almost never done.

1.2Grading model output with a model¶

Scoring generated answers with a second model is convenient and it is also a measurement decision: you have introduced an unblinded rater whose bias is unknown. The problem is that the judge and the candidate typically share training data, so they share blind spots.

import numpy as np

rng = np.random.default_rng(0)
n = 3000
truth = rng.integers(0, 2, n)
hard = rng.random(n) < 0.30                    # items both models find hard
gen = np.where(hard & (rng.random(n) < 0.7), 1 - truth, truth)
shared = hard & (rng.random(n) < 0.6)          # judge inherits the blind spot
judge = np.where(shared, gen,
                 np.where(rng.random(n) < 0.9, truth, 1 - truth))
human = np.where(rng.random(n) < 0.95, truth, 1 - truth)

acc = (gen == truth).mean()
print(f"true accuracy of the candidate      {acc:.3f}")
print(f"accuracy as scored by the LLM judge {(judge == gen).mean():.3f}  "
      f"(bias {(judge == gen).mean() - acc:+.3f})")
print(f"accuracy as scored by the human     {(human == gen).mean():.3f}")
po = (judge == human).mean()
pe = sum((judge == v).mean() * (human == v).mean() for v in (0, 1))
print(f"judge vs human: agreement {po:.3f}  chance {pe:.3f}  "
      f"kappa {(po - pe) / (1 - pe):.3f}")

true accuracy of the candidate      0.785
accuracy as scored by the LLM judge 0.847  (bias +0.062)
accuracy as scored by the human     0.758
judge vs human: agreement 0.758  chance 0.500  kappa 0.515

The judge reports an accuracy six points above the truth, and it does so in the direction that flatters the system being evaluated, because the items it gets wrong are the items the candidate also gets wrong. Note that the judge--human kappa of 0.52 would pass a casual reliability check while the bias is still there: agreement with a human on easy items does not certify grading on hard ones.

1.3Voting reduces variance, not a shared bias¶

Self-consistency --- sample $m$ answers and take the majority --- is bagging, and it inherits bagging’s limitation exactly. Averaging kills the independent component of the error and leaves everything the samples have in common.

import numpy as np

rng = np.random.default_rng(0)
R, m = 40000, 11
p_item = np.clip(0.60 + 0.7 * (rng.random(R) - 0.5), 0, 1)   # shared bias

cases = [("independent errors  ", (rng.random((R, m)) < 0.60).sum(1)),
         ("shared per-item bias",
          (rng.random((R, m)) < p_item[:, None]).sum(1))]
for label, v in cases:
    right = v > m / 2
    consensus = np.maximum(v, m - v) / m       # how unanimous the vote looked
    print(f"{label}  single sample 0.60 -> majority of 11: {right.mean():.3f}")
    print(f"{' ' * 20}  mean consensus when the vote is WRONG "
          f"{consensus[~right].mean():.3f}")
    print(f"{' ' * 20}  Pr(>= 10 of 11 agree | wrong) "
          f"{np.mean(consensus[~right] >= 10 / 11):.3f}")

independent errors    single sample 0.60 -> majority of 11: 0.756
                      mean consensus when the vote is WRONG 0.596
                      Pr(>= 10 of 11 agree | wrong) 0.003
shared per-item bias  single sample 0.60 -> majority of 11: 0.643
                      mean consensus when the vote is WRONG 0.669
                      Pr(>= 10 of 11 agree | wrong) 0.064

Independent errors give the expected gain, from 0.60 to 0.756 at $m = 11$; a shared per-item bias yields 0.643 and stalls there. Worse, the vote’s apparent unanimity moves in the wrong direction: when the correlated ensemble is wrong, it agrees with itself more strongly than when it is right, so consensus is being manufactured rather than earned. A confident majority is evidence about correlation, not about truth.

1.4A retrieved citation is provenance by similarity¶

Attaching references to a generated claim looks like sourcing it. Retrieval ranks passages by surface similarity to the claim, and similarity is not support.

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics.pairwise import cosine_similarity

corpus = [
    "statin therapy reduced LDL cholesterol in a randomized trial "
    "of 400 adults",
    "statin therapy showed no effect on all-cause mortality in this cohort",
    "LDL cholesterol was measured by an enzymatic assay in 400 adults",
    "aspirin reduced all-cause mortality in a randomized trial of 400 adults",
]
supports = {                        # does passage i really support the claim?
    "statins reduce LDL cholesterol": [True, False, False, False],
    "statins reduce all-cause mortality": [False, False, False, False],
}
V = TfidfVectorizer().fit(corpus + list(supports))
C = V.transform(corpus)
for claim, sup in supports.items():
    s = cosine_similarity(V.transform([claim]), C)[0]
    top = int(s.argmax())
    print(f"claim: {claim}")
    print(f"  retrieved passage {top} at cosine {s[top]:.3f}")
    print(f"  that passage supports the claim: {sup[top]}")
    print(f"  any passage in the corpus supports it: {any(sup)}")

claim: statins reduce LDL cholesterol
  retrieved passage 0 at cosine 0.252
  that passage supports the claim: True
  any passage in the corpus supports it: True
claim: statins reduce all-cause mortality
  retrieved passage 3 at cosine 0.340
  that passage supports the claim: False
  any passage in the corpus supports it: False

The second claim is false and the corpus contains nothing supporting it, yet retrieval returns a real, on-topic, high-similarity passage --- about a different drug --- with a higher score than the correct match got for the true claim. Every element of the provenance chain is genuine except the one that matters, and the failure is invisible unless someone reads the cited passage. This is precisely the failure that deep-research assistants produce at scale, and Section that section develops the sampling design that measures its rate rather than asserting it.

1.5The benchmark may already be in the training corpus¶

Public evaluations are the usual evidence that an assistant is fit for a task. They are external validation studies whose test set may have leaked into the instrument being validated.

import numpy as np
from statsmodels.stats.proportion import proportion_confint

rng = np.random.default_rng(0)
n = 200
seen = rng.random(n) < 0.30      # items also present in the training corpus
correct = np.where(seen, 1, rng.random(n) < 0.55).astype(int)
clean = correct[~seen]
lo, hi = proportion_confint(clean.sum(), clean.size, method="wilson")

print(f"headline accuracy over all {n} items      {correct.mean():.3f}")
print(f"accuracy on the {seen.sum():3d} contaminated items  "
      f"{correct[seen].mean():.3f}")
print(f"accuracy on the {clean.size:3d} clean items         "
      f"{clean.mean():.3f}  95% CI ({lo:.3f}, {hi:.3f})")
print(f"contamination inflates the headline by "
      f"{100 * (correct.mean() - clean.mean()):.1f} points")

headline accuracy over all 200 items      0.640
accuracy on the  53 contaminated items  1.000
accuracy on the 147 clean items         0.510  95% CI (0.430, 0.590)
contamination inflates the headline by 13.0 points

Thirty per cent contamination moves the headline from a true 0.51 to a reported 0.64, and the contaminated subset scores a perfect 1.000 --- the tell, when you can see it. Usually you cannot, which is why an assistant’s benchmark numbers are not a substitute for a small validation on your own held-out data.

1.6Reproducibility with a moving dependency¶

Everything above assumes you can say afterwards what produced a number. That requires recording the model identifier alongside the prompt and the seed, because the model is a dependency the vendor updates on their schedule. (An experiment tracker such as MLflow or Weights & Biases already stores one row per run with arbitrary parameters attached, so the practical obstacle is remembering to log the model version rather than the absence of a place to put it.)

import hashlib
import json

def run_id(record, keys):
    payload = json.dumps({k: record[k] for k in keys}, sort_keys=True)
    return hashlib.sha256(payload.encode()).hexdigest()[:12]

march = {"prompt": "summarise the effect", "seed": 0, "temperature": 0.0,
         "model": "assistant-v1.2", "answer": 0.31}
july = {**march, "model": "assistant-v1.4", "answer": 0.44}   # vendor updated

for keys, label in [(["prompt", "seed", "temperature"], "prompt+seed only"),
                    (["prompt", "seed", "temperature", "model"],
                     "prompt+seed+model")]:
    a, b = run_id(march, keys), run_id(july, keys)
    print(f"{label:<18} march={a} july={b}  flags it: {a != b}")
print(f"the two answers differ by {abs(july['answer'] - march['answer']):.2f}")

prompt+seed only   march=29946d68cd17 july=29946d68cd17  flags it: False
prompt+seed+model  march=a6ace2a00035 july=5ecf8e5bb41d  flags it: True
the two answers differ by 0.13

Logging the prompt and the seed alone gives two runs the same identifier while their answers differ by 0.13 --- a provenance record that actively certifies a result it cannot reproduce. The model version belongs in the hash, and the resulting posture is the ordinary one for any analysis with a moving dependency: record the version, expect deprecation, and treat exact reproducibility as having a shelf life.

1.7Tools in practice¶

Nothing in this section argues against using these systems; it argues that using one without a record is an experiment you cannot report. The tools below are grouped by which part of the record they supply, because that is the property that determines whether the diagnostics in this section can be run at all. A system that logs the prompt, the model identifier, the decoding settings and the returned completion turns vibe science into ordinary science with an unusual instrument in the loop; a chat window that keeps none of these does not. As elsewhere in this chapter, capabilities move fast and the current documentation is the only reliable source for them --- what is stable is the distinction between an interface that produces an artefact and one that produces a conversation.

• [Deep-research assistants] Multi-step search and synthesis. Run an extended retrieve-read-write loop and return a long report with inline citations. Fits: orienting yourself in an unfamiliar literature before committing reading time. Watch: the output has the shape of a review without its provenance discipline; see Section that section for the two distinct hallucination modes and how to measure each.

• [Elicit / Consensus / Semantic Scholar] Literature screening and extraction. Search a structured corpus of papers and pull specified fields into a table with a link to each source. Fits: the extraction step of a review, where the unit of output is a cell traceable to a document rather than a paragraph. Watch: extraction accuracy is a measurable quantity and these tools do not know theirs on your corpus. Double-extract a subsample by hand and report agreement before using the table.

• [ASReview] Active-learning title/abstract screening. An open-source screener that reorders records by predicted relevance so eligible studies surface early Schoot et al., 2021. Fits: systematic reviews, where the design is explicitly two-stage and the stopping rule is the object of study. Watch: the ranking is a stage-one instrument with unknown sensitivity on your topic; estimate it from a double-screened subsample rather than assuming it.

• [Weights & Biases / MLflow] Experiment tracking. Log parameters, artefacts and outputs of every run under an immutable run identifier. Fits: making the prompt, the seed, the model identifier and the completion into a row you can count, which is the precondition for every diagnostic in this section. Watch: logging only what the tool captures by default omits the model version, which is the field that mattered most in the provenance example above.

• [promptfoo / Inspect] Prompt evaluation harnesses. Run a fixed prompt set against one or more models repeatedly and score the outputs against declared checks. Fits: turning “it seemed to work” into a rate with an interval, and re-running it when the vendor updates the model underneath you. Watch: a harness whose scorer is itself a model inherits the judge-bias problem measured earlier; keep a human-scored subsample.

• [Jupyter / Quarto] Literate computing. Bind narrative, code and output in one re-executable document. Fits: the honest place to put an assistant-generated analysis, since the generated code, not the conversation, is the artefact. Watch: a rendered document proves the code ran, not that anyone read it.

Once runs are logged, the cheapest available sensitivity analysis becomes possible: re-run the same request and count how often the conclusion changes. The example below treats that as an estimation problem. A flip rate $\pi$ governs how often a rerun reverses the reported conclusion, and the question is how many reruns are needed before the answer is worth reporting.

import numpy as np
from statsmodels.stats.proportion import proportion_confint

rng = np.random.default_rng(0)
PI = 0.30                       # unknown rerun-to-rerun flip rate

print(" reruns  Pr(see >= 1 flip)  Wilson half-width for pi-hat")
for n in (1, 3, 5, 10, 30, 100):
    lo, hi = proportion_confint(round(PI * n), n, method="wilson")
    print(f"{n:7d}  {1 - (1 - PI) ** n:16.3f}  {(hi - lo) / 2:26.3f}")

for n in (5, 20, 50):
    flips = rng.random((4000, n)) < PI
    k = flips.sum(1)
    lo, hi = proportion_confint(k, n, method="wilson")
    cover = np.mean((lo <= PI) & (PI <= hi))
    print(f"\n{n:3d} reruns: median flips {int(np.median(k))}, "
          f"Pr(report a clean 0) {np.mean(k == 0):.3f}, "
          f"CI coverage {cover:.3f}")

 reruns  Pr(see >= 1 flip)  Wilson half-width for pi-hat
      1             0.300                       0.397
      3             0.657                       0.365
      5             0.832                       0.326
     10             0.972                       0.248
     30             1.000                       0.156
    100             1.000                       0.088

  5 reruns: median flips 1, Pr(report a clean 0) 0.174, CI coverage 0.979

 20 reruns: median flips 6, Pr(report a clean 0) 0.001, CI coverage 0.982

 50 reruns: median flips 15, Pr(report a clean 0) 0.000, CI coverage 0.954

Two thresholds fall out, and they are far apart. Detecting that the workflow is unstable at all is cheap: five reruns see at least one flip $83\%$ of the time, and $17\%$ of the time they see none and certify a procedure that reverses itself almost a third of the time. Estimating how unstable it is costs far more --- at five reruns the interval on $\pi$ has a half-width of 0.33, which is the whole scale, and getting under ten points requires around a hundred. The practical reading is that a handful of reruns is a screening test and not a measurement, and that the distinction should appear in what you write: reporting “we re-ran it a few times and it was stable” is reporting a negative result from a design with almost no power.