--- 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 to load the autograder. It will not run inside the book — use the downloaded notebook. ```{code-cell} ipython3 :tags: [skip-execution] import otter grader = otter.Notebook() ``` --- ## 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] grader.check("q1") ``` ### 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] grader.check("q2") ``` ### A3. Put the error in physical units (10 pts) :::{exercise} :label: pr-reg-henry-rmse An $r^2$ of 0.59 does not say how wrong a prediction typically 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] grader.check("q3") ``` --- ## Part B — Visualization (35 pts) :::{exercise} :label: pr-reg-henry-diagnostics Three panels, 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]`. 2. **Error histogram.** The out-of-fold errors $y - \hat{y}$. Mark the mean and, using `np.percentile`, the 5th and 95th percentiles. 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`. Then, in **4–6 sentences**: - The `'Q'` points sit closer to the 1:1 line than the rest. 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. - Is the error histogram centered and symmetric? Say what its shape and width add that $r^2$ alone does not. - Does the choice of `n_splits` matter much here? Say 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, axes = plt.subplots(1, 3, figsize=(14, 4)) ``` --- ## Part C — Open Ended (35 pts) :::{exercise} :label: pr-reg-henry-headroom An RMSE of 3.68 ln units 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 of how much of the model's error could *in principle* be removed, with your assumptions stated. Note that 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 short written argument (**one paragraph**) about 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. 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 ``` --- ## Summary - Part A produced three numbers that answer three different questions: how well the model fits data it has already seen (0.602), how well it predicts data it has not (0.594), and how large a typical prediction error actually is (3.68 ln units, a factor of about 40 in $H$ itself). - Part B looked at the same model three ways the chapter recommends — a parity plot, an error histogram, and the spread of the cross-validated score — and found that part of the apparent fit quality comes from rows that are model output rather than measurement. - Part C set the model's error against the disagreement between independent laboratories, and found both that the headroom is large and that the floor itself is not a single number. ## 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).