Problems: Model Validation

Get this problem set

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.

Before you start

This problem set accompanies 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.

%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()

4344 rows, 1519 species, 519 references

	species	formula	casrn	ref	htype	H_mol_m3_Pa	mindHR_K	lnH	dHsol_kJ_mol
0	(1,1-dimethylethyl)-benzene	C6H5C4H9	98-06-6	2891	M	0.00160	4700	-6.43775	-39.0780
1	(1,1-dimethylethyl)-benzene	C6H5C4H9	98-06-6	2904	M	0.00094	2400	-6.96963	-19.9547
2	(1,1-dimethylethyl)-methanoate	HCOOC4H9	762-75-4	2646	M	0.01400	3600	-4.26870	-29.9321
3	(1,1-dimethylethyl)-methanoate	HCOOC4H9	762-75-4	3518	L	0.01400	3600	-4.26870	-29.9321
4	(1-methylpropyl)-benzene	C6H5C4H9	135-98-8	2891	M	0.00130	4600	-6.64539	-38.2465

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:

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))

29 entries for carbon dioxide, from 29 references
[-8.0164 -7.9866 -7.9866 -7.9866 -8.0164 -8.0164 -8.0164 -8.0164]

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.

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 147

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.

# YOUR CODE HERE
r2_full = ...

grader.check("q1")

A2. Validate it honestly (10 pts)

Exercise 148

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.

# YOUR CODE HERE
r2_cv = ...

grader.check("q2")

A3. Put the error in physical units (10 pts)

Exercise 149

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\).

# YOUR CODE HERE
rmse_cv = ...

grader.check("q3")

---

Part B — Visualization (35 pts)

Exercise 150

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.

# YOUR CODE HERE
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

---

Part C — Open Ended (35 pts)

Exercise 151

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 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.

# 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, Atmos. Chem. Phys. 23, 10901–12440 (2023).

2. sklearn.model_selection.KFold and sklearn.metrics.root_mean_squared_error.