If $x_1,\dotsc,x_{n-1} \in \mathbb{R}^n$, one defines $x_1 \times \cdots \times x_{n-1} \in \mathbb{R}^n$ to be the unique vector such that
$$
\forall y \in \mathbb{R}^n, \quad \langle x_1 \times \cdots \times x_{n-1},y \rangle = \operatorname{det}(x_1,\dotsc,x_{n-1},y),
$$
where the determinant is being viewed as a function of the rows or columns of the usual matrix argument, i.e., as the unique antisymmetric $n$-form $\operatorname{det} : \mathbb{R}^n \times \cdots \times \mathbb{R}^n \to \mathbb{R}$ such that $\det(e_1,\dotsc,e_n) = 1$ for $\{e_k\}$ the standard ordered basis of $\mathbb{R}^n$.
Now, suppose that $x_1,\dotsc,x_{n-1} \in \mathbb{R}^n$ are linearly independent, and hence span a hyperplane $H$ ($n-1$-dimensional subspace) in $\mathbb{R}^n$. Then, in particular, $x_1 \times \cdots \times x_{n-1} \neq 0$ is orthogonal to each $x_k$, and hence defines a non-zero normal vector to $H$; write $$x_1 \times \cdots \times x_{n-1} = \|x_1 \times \cdots \times x_{n-1}\|\hat{n}$$ for $\hat{n}$ the corresponding unit normal. Let $y \notin H$. Then $x_1,\dotsc,x_{n-1},y$ are linearly independent and span an $n$-dimensional parallelopiped $P$ with $n$-dimensional volume
$$
|\operatorname{det}(x_1,\dotsc,x_{n-1},y)| = |\langle x_1 \times \cdots x_{n-1},y\rangle| = \|x_1 \times \cdots \times x_{n-1}\||\langle \hat{n},y\rangle|.
$$
Now, with respect to the decomposition $\mathbb{R}^n = H^\perp \oplus H$, let
$$
T = \begin{pmatrix} I_{H^\perp} & 0 \\ M & I_{H} \end{pmatrix}
$$
for $M : H^\perp \to H$ given by $$M(c \hat{n}) = -c \langle \hat{n},y \rangle^{-1} P_H y = -c\langle \hat{n},y\rangle^{-1}(y-\langle\hat{n},y\rangle\hat{n}),$$ where $P_H(v)$ denotes the orthogonal projection of $v$ onto $H$. Then $T(P)$ is a $n$-dimensional parallelepiped with with vertices $Tx_1 = x_1,\dotsc,Tx_{n-1}=x_{n-1}$, and
$$
Ty = \langle \hat{n},y \rangle \hat{n} = P_{H^\perp} y = y - P_H y,
$$
with the same volume as $P$. On the one hand, since $Ty = y - P_H y$ for $P_H y \in H = \{x_1 \times \cdots \times x_{n-1}\}^\perp$,
$$
\operatorname{Vol}_n(T(P)) = |\operatorname{det}(Tx_1,\dotsc,Tx_{n-1},Ty)|\\ = |\operatorname{det}(x_1,\dotsc,x_{n-1},y-P_H y)|\\ = |\operatorname{det}(x_1,\dotsc,x_{n-1},y)|\\ = \|x_1 \times \cdots \times x_{n-1}\||\langle \hat{n},y\rangle|.
$$
On the other hand, since $Ty \in H^\perp$, $T(P)$ is an honest cylinder with height $\|Ty\| = |\langle \hat{n},y\rangle|$ and base the $(n-1)$-dimensional parallelopiped $R$ spanned by $x_1,\dotsc,x_{n-1}$, so that
$$
\operatorname{Vol}_n(T(P)) = \operatorname{Vol}_{n-1}(R)|\langle \hat{n},y\rangle|.
$$
Thus,
$$
\operatorname{Vol}_{n-1}(R)|\langle \hat{n},y\rangle| = \operatorname{Vol}_n(T(P)) = \|x_1 \times \cdots \times x_{n-1}\||\langle \hat{n},y\rangle|,
$$
so that
$$
\operatorname{Vol}_{n-1}(R)| = \|x_1 \times \cdots \times x_{n-1}\|,
$$
as required.
EDIT: Theoretical Addendum
Let's see what $\phi x_1 \times \cdots \times \phi x_n$ is in terms of $x_1 \times \cdots \times x_{n-1}$ for $\phi$ a linear transformation on $\mathbb{R}^n$.
Define a linear map $T : (\mathbb{R}^n)^{\otimes(n-1)} \to (\mathbb{R}^n)^\ast$ by
$$
T : x_1 \otimes \cdots \otimes x_{n-1} \mapsto \operatorname{det}(x_1,\cdots,x_{n-1},\bullet),
$$
so that if $S : \mathbb{R}^n \to (\mathbb{R}^n)^\ast$ is the isomorphism $v \mapsto \langle v,\bullet \rangle$, then
$$
x_1 \times \cdots \times x_n = (S^{-1}T)(x_1 \otimes \cdots \otimes x_n).
$$
Now, since the determinant is antisymmetric, so too is $T$, and hence $T$ descends to a linear map $T : \bigwedge^{n-1} \mathbb{R}^n \to (\mathbb{R}^n)^\ast$,
$$
x_1 \wedge \cdots \wedge x_{n-1} \mapsto \operatorname{det}(x_1,\cdots,x_{n-1},\bullet);
$$
indeed, if $\operatorname{Vol} = e_1 \wedge \cdots \wedge e_n$ for $\{e_k\}$ the standard ordered basis for $\mathbb{R}^n$, then for any $y \in \mathbb{R}^n$,
$$
\langle x_1 \otimes \cdots \otimes x_{n-1},y \rangle \operatorname{Vol} = \operatorname{det}(x_1,\cdots,x_{n-1},y)\operatorname{Vol} = x_1 \wedge \cdots \wedge x_{n-1} \wedge y,
$$
which, in fact, shows that
$$
x_1 \times \cdots \times x_{n-1} = \ast (x_1 \wedge \cdots \wedge x_{n-1}),
$$
where $\ast : \wedge^{n-1} \mathbb{R}^n \to \mathbb{R}^n$ is the relevant Hodge $\ast$-operator. Thus, a cross product is really an $(n-1)$-form in the orientation-dependent disguise given by the Hodge $\ast$-operator; in particular, it will really transform as an $(n-1)$-form, as we'll see now.
Now, let $\phi : \mathbb{R}^n \to \mathbb{R}^n$ be linear. Observe that the adjugate matrix $\operatorname{Adj}(\phi)$ of $\phi$ can be invariantly defined as the unique linear transformation $\operatorname{Adj}(\phi) : \mathbb{R}^n \to \mathbb{R}^n$ such that for any $\omega \in \bigwedge^{n-1} \mathbb{R}^n$ and $y \in \mathbb{R}^n$,
$$
(\wedge^{n-1})\omega \wedge y = \omega \wedge \operatorname{Adj}(\phi) y,
$$
e.g., in our case,
$$
x_1 \wedge \cdots \wedge x_{n-1} \wedge \operatorname{Adj}(\phi) y = (\wedge^{n-1}\phi)(x_1 \wedge \cdots \wedge x_{n-1}) \wedge y = \phi x_1 \wedge \cdots \wedge \phi x_{n-1} \wedge y,
$$
and that, as a matrix, $\operatorname{Adj}(\phi) = \operatorname{Cof}(\phi)^T$, where $\operatorname{Cof}(\phi)$ denotes the cofactor matrix of $\phi$. Then for any $y$,
$$
\langle \phi x_1 \times \cdots \times \phi x_{n-1},y \rangle \operatorname{Vol} = \operatorname{det}(\phi x_1,\cdots,\phi x_{n-1},y)\operatorname{Vol}\\ = \phi x_1 \wedge \cdots \wedge \phi x_{n-1} \wedge y\\ = (\wedge^{n-1}\phi)(x_1 \wedge \cdots \wedge x_{n-1}) \wedge y\\ = (x_1 \wedge \cdots \wedge x_{n-1}) \wedge \operatorname{Adj}(\phi)y\\ = \langle x_1 \times \cdots \times x_{n-1},\operatorname{Adj}(\phi)y \rangle \operatorname{Vol}\\ = \langle \operatorname{Cof}(\phi)(x_1 \times \cdots \times x_{n-1}),y \rangle \operatorname{Vol},
$$
and hence, since $y$ was arbitrary,
$$
\phi x_1 \times \cdots \times \phi x_{n-1} = \operatorname{Cof}(\phi)(x_1 \times \cdots \times x_{n-1}) = (\ast \circ \wedge^{n-1}\phi \circ \ast^{-1})(x_1 \times \cdots \times x_{n-1}),
$$
in terms of the Hodge $\ast$-operation and the invariantly defined $\wedge^{n-1}\phi$.
Best Answer
My question has gathered several useful comments, which together have given me enough pieces of information that I think I now have the answer. I think this is correct, but please don't hesitate to correct me if I'm making any mistakes:
Throughout this answer, the letter $R$ will denote an element of (the fundamental representation of) $SO(7)$, and $\times_1$ and $\times_2$ will denote 7D cross products. We'll use square brackets to denote the rotated vector $R[a]$ and subscripts to denote any additional quantities that the rotation $R$ might depend on.
Since the definition of the cross product fixes its magnitude, we obviously have that $|a \times_1 b| \equiv |a \times_2 b|$ and so $$ a \times_2 b \equiv R_{\times_1, \times_2, a, b} [a \times_1 b]. $$ (I.e. the map is norm-preserving but not a priori linear in $a$ or $b$.)
However, it turns out that any two 7D cross products are related by a single rigid rotation of the entire vector space: that is, for any two cross products $\times_1$ and $\times_2$, we have $$ a \times_2 b \equiv \pm R_{\times_1, \times_2} \left[ R_{\times_1, \times_2}^{-1}[a] \times_1 R_{\times_1, \times_2}^{-1}[b] \right]. $$ (We can drop the "$\pm$" if we allow $R$ to be a general element of $O(7)$ instead of $SO(7)$.) This nontrivial result is what Massey means when he says that the 7D cross product is "unique up to isomorphism"; by "isomorphism" he means a single rigid rotation of the entire vector space.
Now let's think about this from a different perspective. Above, we started with two arbitrary fixed cross products $\times_1$ and $\times_2$ and considered a rotation $R_{\times_1, \times_2}$ that related them. Now let's instead start with a particular cross product $\times_1$ and imagine a generic "rotated version" $$ f_{\times_1,R}(a, b) := R\left[ R^{-1}[a] \times_1 R^{-1}[b] \right], \qquad R \in SO(7). $$ We can show that $f_{\times_1, R}$ is also a cross product $\times_2$ (with the same orientation as $\times_1$ if $R \in SO(7)$). This is in some sense the converse of the previous result; the previous result states that any two cross products are rotations of each other, and this result states that the rotation of any cross product is another cross product.
In 3D, it turns out that there is a single cross product $\times$ (up to orientation reversal) and $f_{R,\times} \equiv \times$ for all $R \in SO(3)$. But this is not true in 7D. In 7D it turns out that for any fixed cross product $\times_1$, the claim $f_{\times_1, R} \equiv \times_1$ is only true for a set of rotations $R$ that forms a proper subgroup of $SO(7)$ that is isomorphic to the exceptional Lie group $G_2$. (As Jason DeVito explains in the comments, the exact subgroup of $SO(7)$ depends on $\times_1$, but it is always isomorphic to $G_2$.)
The set of distinct 7D cross products is therefore $X := SO(7)/G_2$. (More precisely, this is the set of transformations that maps an arbitrary initial cross product into any other cross product.) Since $SO(n)$ is always a simple Lie group, $G_2$ is not a normal subgroup of $SO(7)$. So $X$ does not have the full structure of a Lie group, but only of a homogeneous space. In fact, it turns out that this homogeneous space $SO(7)/G_2$ is isometric to the real projective space $\mathbb{R}P^7$.
Now let's fix an orthonormal basis of $\mathbb{R}^7$ and consider the set of cross products that map basis vectors to basis vectors. One way to do this is to choose an arbitrary such cross product $\times_1$ and then consider the cosets $R\, G_2 \in X$ that transform $\times_1$ to another cross product with the same property of mapping basis vectors to basis vectors. We find that this is a finite subset of $X \cong \mathbb{R}P^7$ with 480 elements. (But the other elements of $X$ are still perfectly valid cross products.)
TLDR: There are an uncountably infinite number of distinct 7D cross products, and they form a Reimannian manifold isometric to the real projective space $\mathbb{R}P^7$.