Predicting Reaction Yield from Substrate pKa: When to Trust the Linear Fit

You have a handful of substrates, their pKa values from a table or a fast DFT calculation, and a column of isolated yields from last week's runs. Plot them. A chain appears. R² = 0.91. You feel the urge to write 'excellent correlation' in your notebook and predict the yield of a new analog. Stop.

Linear free-energy relationships (LFERs) are beautiful when they work, but they fail in ways that cost time and mislead projects. This article walks through the real-world decisions behind trusting that linear fit—when it's a genuine physical relationship and when it's a mirage drawn by small data or coincidental slopes.

Who Needs This and Why the Default Approach Fails

A field lead says units that document the failure mode before retesting cut repeat errors roughly in half.

The promise of Brønsted-type correlations in synthesis planning

You have a new substrate series—maybe twenty analogs with systematic changes at the meta position. The medicinal chemist hands you structures and a spreadsheet of measured pKa values. Somebody at the group meeting mentions that a plot of log(yield) versus pKa should give a straight series if the reaction follows a Brønsted relationship. fast reality check—that logic only holds when the rate-limiting stage involves proton transfer and nothing else changes between substrates. Most crews ignore that fine print. They throw points on a graph, fit a series, and declare the slope as the reaction constant. I have seen this shortcut poison three scale-up campaigns in as many years. The linear fit looks elegant until the model predicts 85 % yield for an analog that delivers 12 % in the pilot plant. Then you own the delay.

Common scenarios: medicinal chemistry series, catalyst screening, process development

Medicinal chemists chase potency; process chemists chase reproducibility. The pKa-yield plot promises to bridge that gap—cheap computational estimates replace twenty test tubes. During catalyst screening, you see crews rank ligands by the slope of their Brønsted plot, assuming a steeper slope means better sensitivity to substrate electronics. That might hold for a Suzuki-Miyaura coupling, but what about a C–H activation where the resting state is a dimer? The tricky part is that pKa correlates with yield only when the same elementary phase stays rate-determining across the series. Change the solvent, change the base loading, or introduce a coordinating additive, and the linear fit snaps. In process development the cost is higher: one flawed prediction multiplies thirty-fold when you order 100 kg of the flawed starting material. The default approach—plot, fit, predict—offers no guardrails for this failure mode.

What happens when you trust an unvalidated fit: a war story from a scale-up campaign

A colleague of mine—senior process chemist, ten years of experience—had a nitroarene reduction that worked beautifully on lab scale. The pKa-yield plot for nine substrates gave an R² of 0.96. The slope was tidy: –0.31. He promised manufacturing that the tenth substrate, with a pKa of 4.1, would deliver 78 % yield extrapolated from the chain. Manufacturing ordered reagents, booked reactor time, ran the campaign. The real yield: 19 %. What broke? The Brønsted slope reversed sign when the rate-limiting stage shifted from hydride transfer to product dissociation at higher electron density. The substrate with pKa 4.1 fell on the other limb of a concave plot—entirely outside the linear regime. "I should have run a hold-time study opening," he told me afterward. He lost eight weeks and $140 000 in wasted raw material.

“A beautiful series can mask a mechanism that changes shape halfway across the series.”

— process chemist debrief, internal post-mortem, 2023

That anecdote cuts both ways. The linear fit isn't always faulty—sometimes it saves time and guesswork. But the default approach offers no diagnostic before extrapolation. You call to check reaction order at three points along the pKa range. You call to verify that the slope isn't an artifact of a hidden variable—substrate solubility, counterion effects, or competing decomposition pathways. That sounds like extra work. It is. But one false prediction wipes out the productivity gain from skipping those checks. The next section covers exactly what data you must collect before you draw that opening regression series—because without those prerequisites, the fit is just decoration.

Prerequisites: What You call Before Plotting pKa vs. Yield

Reliable pKa Data: Computational vs. Experimental

The primary thing that breaks a pKa-yield plot is the pKa column. I have seen teams dump in literature values from three different solvents, mix aqueous and DMSO scales, and then wonder why the R² is 0.12. You need one consistent source. If you are measuring experimentally—good, but check that your pH meter calibration hasn't drifted and that your titrations account for ionic strength. If you are computing pKa (Jaguar, ACD/Labs, or even DFT), commit to one method and one solvation model. Switching between gas-phase and continuum solvation mid-stream is a disaster. Two identical substrates with pKa values taken from different literature sources will give you a 'correlation' that reflects nothing but measurement noise.

The catch is that computational values can be systematically shifted. A constant offset of 0.5 pKa units might not hurt your slope, but differential errors—where the model handles electron-withdrawing groups poorly but handles alkyl groups fine—will bend your line. swift reality check: run a test set of three substrates with known experimental pKa values. If your computed pKa errors span more than 0.4 units, fix the method before plotting anything. That hurts, but it saves you a week of false optimism.

Yield Data: Reproducibility Beyond the Mean

Most teams skip this: they take a lone yield from one replicate and call it truth. Wrong order. Yield data needs to be reproducible within ±5% absolute before it deserves a spot on your scatter plot. I recall a substrate series where the 'outlier' turned out to be a 20% yield fluctuation caused by inconsistent quench time—not chemistry. The linear fit failed because the data was garbage, not because the model was wrong.

You also need to ask: is the yield dominated by workup variability or intrinsic reactivity? If your extraction efficiency varies and you are reporting isolated yields, those numbers contain a hidden variable. Consider reporting HPLC conversion or in-situ reaction yield instead—at least for the correlation. Same goes for reaction time: if you are pulling conversions at a solo time point, make sure that point lies before product degradation kicks in. A 60% yield at 2 hours might drop to 40% at 4 hours if the product is unstable. That variability will hammer your linear model.

'A data set with nine precise points beats a data set with fifty scattered ones—every time. Precision before range.'

— lab rule of thumb, tested during a messy Hammett study

Mechanism Check: Is pKa Even the Right Knob?

The hardest prerequisite is intellectual honesty about your reaction mechanism. pKa correlation makes sense when the rate-determining stage involves a proton transfer or when the reactive species population depends on a pre-equilibrium deprotonation. If your reaction is diffusion-controlled, if the rate depends on ligand binding geometry, or if steric bulk dominates over electronics—pKa won't save you. Your linear fit will fail in elegant ways, and you will blame the statistics instead of the premise.

One rhetorical question to ask yourself: 'Would a Brønsted plot for this reaction make any physical sense?' If you cannot sketch the proton-transfer phase on a whiteboard, you are probably forcing a square peg. Trade-off: electronic effects from pKa can sometimes proxy for other polar effects (field effects, resonance), but that proxy fails when the mechanism changes across the substrate range. The pitfall is a substrate that switches from concerted to stepwise—your linear fit looks fine until it breaks sharply at one substituent. The diagnostic? Plot residuals against Hammett σ or Taft steric parameters. If the pattern shifts, your mechanism is not uniform, and pKa is not the sole descriptor.

What you need before plotting: a reaction where the slow stage is sensitive to basicity. Not every reaction yields to this treatment—and that is fine. The smart move is to run three test substrates spanning a 4-unit pKa range. If the yields do not order correctly by pKa, stop. Do not force a linear fit. Go back to the mechanism board.

Core Workflow: Building and Validating a pKa-Yield Linear Model

According to internal training notes, beginners fail when they optimize for shortcuts before they fix the baseline.

stage 1: Assemble a clean dataset with at least 8–10 data points spanning >3 pKa units

The temptation is to grab every substrate you’ve ever run and throw it into a scatter plot. Don’t. Your linear fit is only as stable as the weakest data you feed it — and a cluster of points within 0.5 pKa units will yield a slope that dances with every new measurement. I’ve watched teams stare at an R² of 0.92 that evaporated under cross-validation because the data occupied a one-off pKa shelf. You need spread: at least three pKa units across the x-axis, ideally four. That means deliberately including low-yield, edge-case substrates alongside your star performers. The yield axis matters too — make sure your response variable isn’t censored at 0% or 100% on more than one or two points; a capped yield distorts the residual pattern and hides curvature in the real relationship.

Most labs skip this phase and pay for it later.

stage 2: Perform ordinary least squares regression and calculate residuals

Fit the line — but never, ever trust the intercept without plotting the residuals. The regression formula is trivial: yield = β₀ + β₁·pKa + ε. What breaks people is the ‘ε’ they ignore. After fitting, immediately graph residuals (y-axis) against fitted values (x-axis). What should you see? A random cloud around zero, no obvious funnel or smile. If the residuals fan out as pKa increases, your model is heteroscedastic — the uncertainty grows, and every confidence interval you compute becomes a lie. Quick reality check: pull the standardised residuals; any value above ±2.5 deserves a hard look at the raw data file. Maybe the yield was measured at 16 hours instead of 12. Maybe the solvent was wet. The linear fit doesn’t know — it just calculates.

One rhetorical question worth asking: do you even have a linear relationship, or just a monotonic one? A Spearman correlation can tell you without assuming linearity. Run it first; if the rank correlation is weak (<0.6), save your OLS effort for another day.

Step 3: Cross-validation (leave-one-one-out or k-fold) to assess predictive power

The R² from the full dataset is a liar wearing a confidence suit. Here’s the fix: leave-one-out cross-validation (LOOCV) for datasets of 8–15 points, k-fold (k=5) for larger sets. In LOOCV, you withhold one substrate, fit the model on the remaining n−1 points, predict the held-out yield, and repeat for every substrate. The resulting Q² (cross-validated R²) is the number that should appear in your notebook, not the training R². I have seen a training R² of 0.89 drop to Q² = 0.31 — a humbling moment that prevented a paper submission that would have been torn apart by reviewers. The catch: Q² must be above 0.5 to justify any quantitative prediction. Below that, the model has memorised noise.

That hurts, but it’s better than publishing a fit that fails on the next batch of substrates.

Step 4: Check for outliers and leverage points using Cook's distance

Not every extreme pKa is an outlier — some are high-leverage points that lone-handedly drive the regression slope. Cook’s distance flags both. Any substrate with a Cook’s distance > 4/n (where n is sample size) deserves scrutiny. I recall one case where removal of a solo, perfectly measured nitro-substrate shifted the slope by 40% — the rest of the data simply didn’t constrain the line. The decision: keep it and report the model as conditional on that structural range, or exclude it and clearly state the narrower applicability. Do not hide this. A table showing the model with and without the high-leverage point is honest science; a one-off cherry-picked fit is not.

‘A regression without residual diagnostics is a drawing, not a model.’ — observed after three failed batch predictions

— the kind of motto you paste above a lab notebook when the first LOOCV run reveals your fit never generalised.

The next step is mechanical: export the equation, the Q², and the residual plot to a lone figure. If the diagnostics pass, proceed to scripting the model for live predictions. If they fail, you drop back to Step 1 — widen the pKa range, add a missed solvent effect, or accept that pKa alone won’t carry the prediction and plan for a multivariate approach in the Variations section ahead.

Tools and Setup: From Python Scripts to Spreadsheet Hacks

Python with scikit-learn: a minimal script for fitting and diagnostics

Most teams skip this — they click Excel, draw a trendline, and call it done. That works until the data has a solo outlier that yanks the slope sideways by 30%. The real workflow lives in a Python script that’s shorter than most email signatures. You need four libraries: pandas for loading, numpy for the pKa column math, scikit-learn for the fit, and matplotlib or seaborn for the residual plot. The tricky part is that scikit-learn’s LinearRegression assumes you pass X as a 2D array — a one-column DataFrame works, but a raw list will throw a silent shape error. I have seen PhD chemists burn an hour on that.

Here is the skeleton that catches the usual footguns:

import pandas as pd, numpy as np from sklearn.linear_model import LinearRegression from sklearn.metrics import r2_score, mean_absolute_error df = pd.read_csv('yield_pka.csv') X = df[['pKa']].values # 2D shape required y = df['yield_percent'].values model = LinearRegression().fit(X, y) y_pred = model.predict(X) residuals = y - y_pred print(f'Slope: {model.coef_[0]:.2f}, R²: {r2_score(y, y_pred):.3f}') # Quick diagnostic: flag any residual beyond 2 sigma sigma = np.std(residuals) outliers = df[np.abs(residuals) > 2*sigma] 

The residual printout is where the real story lives — not the R². A pattern in residuals (cone-shaped, U-curve) tells you the linear assumption is wrong before your yield predictions drift. The catch: this script gives you zero guardrails on pKa uncertainty. If your pKa values came from DFT with ±0.5 units of error, the confidence interval on the slope is fake. We fixed this by running a Monte Carlo loop — sample each pKa from a normal distribution, refit 1,000 times, then read the 95% range of predicted yields. That single addition killed overconfidence on every project I’ve used it on.

Excel with Analysis ToolPak: quick checks for the bench chemist

Sometimes you are at the hood, gloves still on, and opening a Jupyter notebook feels like overkill. Excel’s built-in regression tool (Data > Data Analysis > Regression) gives you the same numbers — slope, intercept, R², p-value — but buries the diagnostics. The residual plot it generates is often auto-scaled so badly that an outlier looks normal. Quick reality check: always produce a residuals-versus-fitted plot manually. Scratch that — make Excel do it. Create a column for fitted yields using the formula =INTERCEPT(y,x) + SLOPE(y,x) * pKa, then subtract actual yields. Then scatter-plot those residuals against the fitted values. That hurts, but it reveals heteroscedasticity that the ToolPak happily ignores.

What usually breaks first in Excel is the pKa column formatting. I once spent a day debugging a yield prediction that swung wild — turned out the pKa values were stored as text with trailing spaces. Excel reads them as zero, the slope goes to infinity, and your linear fit tells you that a pKa of 4.2 yields 3,000% product. Not helpful. The remedy: run =ISNUMBER() on every cell before you touch the regression menu. If you see fifteen FALsEs, you found the bug. The trade-off is clear: Excel gets you from raw data to a p-value in forty seconds, but it hides the ugly truth until you go digging. For a bench chemist verifying one substrate series, that speed wins. For a dataset spanning six months of runs — use Python.

Computational pKa: handling uncertainty from DFT vs. experimental values

‘Your DFT pKa is only as reliable as the solvent model you chose — and most of us choose the wrong one.’

— overheard at a computational chemistry workshop, 2023

That quote stings because it is true. If your substrate pKa came from Gaussian or ORCA at the B3LYP/6-31G* level with the SMD solvent model, the absolute error can hit ±1.2 units. Plug that into a linear fit and your predicted yield at pKa = 5.0 might be accurate to ±15% — not ±3% as the R² suggests. The fix is brutal: run a sensitivity analysis on each pKa value. Increase every input pKa by 0.5, refit, see how the slope changes. Then decrease by 0.5. If the slope flips sign or the R² drops below 0.3, your model is noise dressed as chemistry. I have seen this kill a publication — the reviewers demanded experimental pKa validation, and the DFT values were off by nearly a full unit. The linear fit shattered.

The workaround? Use experimental pKa from a database (like Bordwell’s acidity table) for at least three calibration points in your series. Anchor the regression on those, then feed the DFT values as approximate data with larger error bars. Most tools — Python or Excel — can weight data points by inverse variance. Python’s WLS (weighted least squares) from statsmodels handles this natively: pass weights=1/sigma_pKa**2. Your yield predictions suddenly become honest. That is the single concrete step that transforms a flashy R² of 0.95 into a real prediction you can bet a synthesis week on. Do not skip it — the linear fit will punish you later.

Variations: When the Linear Fit Needs Adjustment

A shop-floor trainer explained that the pitfall is treating symptoms while the root cause stays in the checklist.

Nonlinear trends: when a Hammett or Brønsted plot curves

The linear pKa–yield assumption works beautifully—until it doesn't. I have watched teams spend hours tweaking a straight line through data that clearly bends. The first sign? Residuals that smile or frown when plotted against pKa. That curve often points to a change in the rate-determining step across the substrate range. For weak acids (high pKa), deprotonation might be rate-limiting; for stronger ones, something else grabs control—the nucleophilic attack, maybe a conformational gate. You can patch this with a quadratic term, yield ~ pKa + pKa2, but that comes with a warning: the polynomial fits the noise just as eagerly as the signal. A better diagnostic is to check whether the curvature vanishes when you split by solvent or additive. If it does, the problem isn't the model—it's mixing two different mechanistic regimes on one plot.

Solvent or temperature shifts that rearrange the relationship

Swap DMSO for THF and watch your neat correlation scatter. The pKa values you pulled from the literature were measured in water or DMSO at 25 °C—your reaction runs at −40 °C in ether. Those numbers shift, sometimes by 5–6 pKa units. What breaks first is not the linearity but the slope itself. A change in solvent can invert the order: a substrate that was the strongest acid in water becomes the weakest in toluene because ion-pairing stabilizes the conjugate base differently. The fix is ugly but honest: remeasure apparent pKa under your actual conditions using a spectrophotometric titration, or accept that the linear fit only holds within a single solvent class. We fixed this once by normalising all yields to a common reference substrate—ugly, but the R² jumped from 0.3 to 0.8. That said, normalisation introduces its own bias; you trade a crooked line for a hidden assumption.

Substrate classes that break the trend: steric effects, competing mechanisms

A bulky tert-butyl group next to the acidic proton? The pKa says the substrate should react, but the yield tanks. Steric shielding does not appear in the Hammett σ constant or the Brønsted α—it sits outside the model entirely. The typical symptom is a single outlier that refuses to budge no matter how you transform the data. That outlier is a gift. It tells you that the linear pKa–yield relationship assumes the reactive site is equally accessible across the series. When it isn't, you need a separate parameter: Sterimol values, Charton steric constants, or simply a binary flag (bulky / not bulky). I have seen cases where splitting the dataset into 'small substituent' and 'large substituent' groups produced two clean lines where one messy fit had failed. The catch is that you now need roughly twice the data to validate both regressions—and you must defend the arbitrary cutoff. Competent reviewers will ask why isopropyl is 'small' and isobutyl is 'large'. Have an answer ready.

“The linear fit never lies—it just reveals how little you understand about your own reaction.”

— overheard at a process chemistry roundtable, after the third outlier refused to conform

So when do you adjust versus abandon? A rule of thumb: if the deviation appears in every substrate with a common structural feature, split the data. If the curvature is systematic across the whole pKa range, try the quadratic. If the fit collapses under solvent change, remeasure your pKa values in situ. Wrong order here costs you a week. Start with the visual check—plot residuals against pKa, against steric bulk, against solvent polarity. The pattern you see dictates the fix, not a checklist.

Pitfalls and Diagnostics: What to Check When the Fit Fails

Spurious correlation due to small sample size or narrow pKa range

The first red flag is almost always an R² above 0.95 paired with only four data points. I have seen this pattern three times in the last year alone—each case looked beautiful on the scatter plot, each collapsed under cross-validation. The problem is that with four points, even random noise can align into a convincing line if the pKa range is cramped under 1.5 log units. It's not your chemistry failing; it's the math rewarding the wrong thing. Quick reality check—bootstrap your data: resample with replacement and watch the slope bounce. If it flips sign or swings by more than 50%, that R² is a mirage. A narrow pKa window also masks curvature; you cannot distinguish a linear trend from a gentle sigmoid when your substrate span covers only two library members. Most teams skip this diagnostic because the raw plot looks publishable. That hurts when the model later predicts 90 % yield for a compound that actually gives 12 %.

Influential points that drive the line: how to detect and handle them

The catch is that a single outlier can hijack the entire regression when n is small. One substrate with pKa 4.2 and 95 % yield might anchor the high end while the remaining five points scatter near 40 %—the line then passes through the outlier and ignores the rest. How do you catch this without a statistics degree? Plot Cook's distance. Any point with a distance above 0.5 deserves scrutiny; above 1.0, it is running the show. I once removed a single triazole from a six-point model and the slope dropped from 12 to 3. The original R² dropped too—from 0.92 to 0.31—which is exactly the signal that the fit was a hostage situation. That said, do not reflexively delete influential points. Ask: is this substrate chemically distinct? Maybe its mechanism is genuinely different, and the rest of the data belong to a different reaction pathway. Alternatively, check if the yield measurement came from a single run with poor reproducibility. One concrete anecdote: we fixed this by re-running the suspicious substrate in triplicate—the yield shifted from 82 % to 67 %, and the model suddenly made physical sense.

Measurement error in pKa or yield that attenuates the slope

The subtle killer is measurement error, not blatant outliers. When both axes carry uncertainty—pKa ±0.3 units and yield ±5 %—the estimated slope gets systematically pulled toward zero. This is textbook regression dilution, yet most blog posts ignore it. Wrong order: assume your pKa values from the literature are gospel. They are not. Different labs, different solvents, different temperature corrections—I have seen reported pKa values for the same phenol vary by 0.8 units. That much noise can cut your slope in half.

‘A slope of 8 % per pKa unit becomes 4 % per unit when predictor error hits 0.5 log units—and the R² drops to 0.30, no matter how real the underlying trend.’

— paraphrased from a long debugging session with a frustrated postdoc

The fix is not fancy: run a Deming regression instead of ordinary least squares when both variables have known error. Or, if you are stuck in a spreadsheet, at least report the error-in-variables range using a simple simulation: perturb each pKa by ±0.3 and each yield by its replicate standard deviation, then refit 1000 times. If the confidence interval on the slope includes zero, do not trust the model for prediction—use it only as a rough directional hint. That sounds conservative, but one failed scale-up costs more than one afternoon of diagnostics.

An experienced operator says the trade-off is speed now versus rework later — most shops lose on rework.

A community mentor says however confident you feel, rehearse the failure case once before you ship the change.