--- jupytext: text_representation: extension: .md format_name: myst kernelspec: display_name: Python 3 language: python name: python3 --- ```{contents} :local: :depth: 2 ``` # Problems: Model Validation :::{admonition} Get this problem set :class: seealso {download}`Download everything (Topic2.2-Model_Validation_Problems.zip) ` — the notebook and `henry_law.csv`, in a folder that is ready to run as-is. The Download badge at the top of the page will also give you the notebook on its own. On Vocareum everything is already set up for you. ::: :::{admonition} Before you start :class: tip This problem set accompanies {doc}`Model Validation `. It is worth **100 points**: - **Part A — Skill Checks (30 pts)** — short answers, auto-graded, resubmit as often as you like until they pass. - **Part B — Visualization (35 pts)** — plots plus written interpretation, peer graded. - **Part C — Open Ended (35 pts)** — one synthesis problem, peer graded. The parts build on each other: Part A works out the syntax you need for Part B, and Part B produces the evidence you argue from in Part C. Do them in order. ::: ## Setup Henry's law constant $H$ describes how much of a gas dissolves in water at equilibrium. Its temperature dependence is set by the enthalpy of solvation $\Delta H_\text{sol}$ through the van 't Hoff relation. A **linear free-energy relationship** predicts that species which dissolve more exothermically should also be more soluble — that is, that $\ln H$ should fall roughly linearly with $\Delta H_\text{sol}$. You will test that claim on Sander's compilation of Henry's law constants (version 5.0.0), which gathers values from the published literature for thousands of species. Each row is **one species as reported by one literature reference**: | column | meaning | |---|---| | `species`, `formula`, `casrn` | chemical identity | | `ref` | which literature reference reported this value | | `htype` | how the value was obtained — `M` measured, `L` literature review, `V` from VLE data, `T` from thermodynamic data, `Q` estimated by a QSAR | | `H_mol_m3_Pa`, `lnH` | Henry's law solubility constant at 298.15 K, and its log | | `mindHR_K`, `dHsol_kJ_mol` | van 't Hoff temperature coefficient, and $\Delta H_\text{sol}$ | The important structural feature: **many species appear more than once**, because several groups measured them independently. That gives you something rare — a direct measurement of how much independent laboratories disagree. ```{code-cell} ipython3 %matplotlib inline import numpy as np import pandas as pd import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression from sklearn.model_selection import KFold from sklearn.metrics import r2_score, mean_absolute_error, root_mean_squared_error try: plt.style.use('../settings/plot_style.mplstyle') # available inside the book except OSError: pass # downloaded notebook: use defaults df = pd.read_csv('data/henry_law.csv') print(f"{len(df)} rows, {df.species.nunique()} species, {df.ref.nunique()} references") df.head() ``` Because that repetition matters later, here is what it looks like. To pull out every row for one species, build a **boolean array** by comparing the species column against a name, then use it to index — the same trick works on any array of the same length: ```{code-cell} ipython3 sel = (df['species'].values == 'carbon dioxide') # a True/False array, one per row print(f"{sel.sum()} entries for carbon dioxide, from {df['ref'].values[sel].size} references") print(df['lnH'].values[sel][:8].round(4)) ``` Those are the first eight of twenty-nine independent reports of the same constant, and they do not all agree. Keep that in mind — Part C is about exactly this. Run this once. It defines the `check` helper used after each Part A answer. It will not run inside the book — use the downloaded notebook. ```{code-cell} ipython3 :tags: [skip-execution] # Run this once. It defines `check`, used after each Part A answer below. import hashlib _EXPECTED = { 'q1': (4, ['57ac221ae7d3', '91faaf5e510b', '67558e22c994', '67bd4bdfe127', 'cbfa01c1d052', 'e0c027c6edd0', '92039333635d', '2d9023de4e90', '46458459d3f5', '27edf60435b1', 'ba87409b22ea', '973f75322414', 'cb3d629732e4', '65f45eafdb5d', 'ed40404ab65d']), 'q2': (4, ['0b08c753b20b', 'df11515bf4c3', '59e63a565041', '27ba3fc28a0e', '2a24e4750169', 'fce97f5ebcdb', '65f5bfbea45a', '0541eb18b7ca', '64bebad620f8', '43b0aca226b1', 'f7fd34bfaa05', '44d794ada822', 'ff769759650c']), 'q3': (3, ['b4be85884507', 'd3b0f9a5ae0b', 'd66dd59b43f5', 'fc8ca000190b', '8754e0fea022', 'f59185661213', '255bb6100d46', '2c7e26c6ba57', '69689d8f0e6f']), } def check(qid, value): """Compare an answer against a stored digest. Practice only -- nothing is recorded.""" places, digests = _EXPECTED[qid] if value is Ellipsis: print(f"{qid}: not answered yet") return got = hashlib.sha256(f"{float(value):.{places}f}".encode()).hexdigest()[:12] print(f"{qid}: {'correct' if got in digests else 'not the expected value'}") ``` --- ## Part A — Skill Checks (30 pts) Three questions, 10 points each. Throughout, the model is an ordinary least-squares fit of `lnH` (target) on the single feature `dHsol_kJ_mol`, using **all 4,344 rows**. ### A1. Fit the relationship (10 pts) :::{exercise} :label: pr-reg-henry-lfer Fit `sklearn.linear_model.LinearRegression` with `dHsol_kJ_mol` as the only feature and `lnH` as the target, using the whole dataset. Report the $r^2$ **on the same data you fit on**. Assign it to `r2_full`. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE r2_full = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q1", r2_full) ``` ### A2. Validate it honestly (10 pts) :::{exercise} :label: pr-reg-henry-cv The A1 number is optimistic — it scores the model on the data used to fit it. Redo the evaluation with **5-fold cross-validation**, built the way the chapter builds it: create ``` kf = KFold(n_splits=5, shuffle=True, random_state=0) ``` and loop over `kf.split(X)`, fitting on the training indices and scoring on the test indices each time. Report the **mean** of the five fold $r^2$ values as `r2_cv`. Do not change the seed. A shuffled split is a random quantity, so a reported cross-validation score means nothing unless the split is reproducible. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE r2_cv = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q2", r2_cv) ``` ### A3. Put the error in physical units (10 pts) :::{exercise} :label: pr-reg-henry-rmse An $r^2$ alone does not say how wrong a prediction really is. Collect the **out-of-fold predictions** — inside the same loop as A2, store each fold's predictions in the right positions of a full-length array, so that every row ends up predicted by a model that never saw it. Then report the root-mean-squared error of those predictions against `lnH`, using `root_mean_squared_error`. Assign it to `rmse_cv`. The units are ln units of $H$. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE rmse_cv = ... ``` ```{code-cell} ipython3 :tags: [skip-execution] check("q3", rmse_cv) ``` --- ## Part B — Visualization (35 pts) :::{exercise} :label: pr-reg-henry-diagnostics Three plots, using the three diagnostics from the chapter. 1. **Parity plot.** Out-of-fold prediction against actual `lnH`, with the 1:1 line drawn. Color the points by whether `htype` is `'Q'` (a QSAR estimate — a value produced by a correlation rather than measured) or anything else. You can select rows of a NumPy array with a boolean array, e.g. `is_q = (df['htype'].values == 'Q')` and then `y[is_q]`. The `'Q'` points sit closer to the 1:1 line than the rest. In **two sentences**, explain why that is expected given how a QSAR value is produced, and say what it implies about using this pooled dataset as evidence that the relationship is physically real. 2. **Error histogram.** The out-of-fold errors $y - \hat{y}$. Mark the mean and, using `np.percentile`, the 5th and 95th percentiles. In **one sentence**, say whether the histogram is centered and symmetric, and what its shape and width add that $r^2$ alone does not. 3. **How stable is the estimate?** Recompute the mean cross-validated $r^2$ for `n_splits` in $\{2, 5, 10, 20\}$, each for at least 10 different `random_state` values, and show the spread as a boxplot per `n_splits`. In **one sentence**, say whether the choice of `n_splits` matters much here and why, referring to the model's size relative to the dataset. Label all axes with units. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE fig, ax = plt.subplots() ``` --- ## Part C — Open Ended (35 pts) :::{exercise} :label: pr-reg-henry-headroom The cross-validated RMSE you computed in A3 sounds terrible. But no model can do better than the data it is checked against, and this dataset has something unusual that lets you find out how good "perfect" would even be: **many species were measured independently by several different groups**, so the disagreement among those repeated values is measurement scatter that no model could ever predict. Decide whether this model is worth improving, and defend it with evidence. Your answer should include: 1. For **at least four** species that appear many times, the standard deviation of `lnH` across their repeated entries. Select a species' rows with a boolean array, e.g. `sel = (df['species'].values == 'carbon dioxide')`, then use `np.std(..., ddof=1)`. Some species with many entries: carbon dioxide, methane, benzene, ammonia, oxygen, trichloromethane, methylbenzene. 2. A comparison of your `rmse_cv` from A3 against those spreads. How many times larger is the model's error than the disagreement between laboratories? 3. A statement, in **two sentences**, of how much of the model's error could *in principle* be removed, with your assumptions stated. Your four species will not agree with each other — say what that does to the argument. 4. At least one concrete attempt to improve things, evaluated with the **same** cross-validation procedure as A2. A second feature, a different model from {doc}`Non-parametric Models `, or a restriction of the dataset are all fair game. 5. A verdict, in **two sentences**, on whether the effort is warranted, citing your own numbers. Be careful with step 4: at least one obvious-looking move makes the reported score *worse*, for a reason that is a property of the data rather than a failure of the model. If you hit it, explain it in **one sentence**. There is more than one defensible conclusion. You are graded on the reasoning and the evidence, not on reaching a particular verdict. ::: ```{code-cell} ipython3 :tags: [skip-execution] # YOUR CODE HERE ``` --- ## Additional Reading 1. R. Sander, [Compilation of Henry's law constants (version 5.0.0) for water as solvent](https://doi.org/10.5194/acp-23-10901-2023), *Atmos. Chem. Phys.* **23**, 10901–12440 (2023). 2. [`sklearn.model_selection.KFold`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.KFold.html) and [`sklearn.metrics.root_mean_squared_error`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.root_mean_squared_error.html).