[proofplan]
The proof rewrites Wilks' lambda $\det W/\det T$ by factoring $T = W+B$ through the invertible matrix $W$. This reduces the determinant ratio to the reciprocal of $\det(I_p+W^{-1}B)$. The nonzero eigenvalues of $W^{-1}B$ are the canonical roots $\theta_1,\dots,\theta_s$, while all remaining eigenvalues are $0$, so the determinant is $\prod_{k=1}^s(1+\theta_k)$. Finally, the defining relation between canonical roots and squared sample canonical correlations converts each factor $(1+\theta_k)^{-1}$ into $1-\hat\rho_k^2$.
[/proofplan]
[step:Factor the total sums-of-products determinant through $W$]
Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $W$ is symmetric positive definite, $W$ is invertible and $\det W > 0$. Using $T = W+B$, we factor
\begin{align*}
T = W + B = W(I_p + W^{-1}B).
\end{align*}
Taking determinants and using multiplicativity of the determinant gives
\begin{align*}
\det T
= \det W \cdot \det(I_p + W^{-1}B).
\end{align*}
Because $\det W \neq 0$, division by $\det T$ yields
\begin{align*}
\frac{\det W}{\det T}
= \frac{1}{\det(I_p + W^{-1}B)}.
\end{align*}
[/step]
[step:Express $\det(I_p+W^{-1}B)$ through the canonical roots]
Let $A := W^{-1}B \in \mathbb{R}^{p \times p}$. By hypothesis, the nonzero eigenvalues of $A$ are $\theta_1,\dots,\theta_s$, counted with algebraic multiplicity. Since $B$ has rank $s$ and $W^{-1}$ is invertible, $A$ has rank $s$, so its remaining $p-s$ eigenvalues are equal to $0$.
Therefore the eigenvalues of $I_p + A$ are
\begin{align*}
1+\theta_1,\dots,1+\theta_s,
\underbrace{1,\dots,1}_{p-s\ \text{times}}.
\end{align*}
The determinant of a square matrix is the product of its eigenvalues counted with algebraic multiplicity, hence
\begin{align*}
\det(I_p + W^{-1}B)
= \prod_{k=1}^s (1+\theta_k).
\end{align*}
[guided]
Define $A := W^{-1}B \in \mathbb{R}^{p \times p}$. The role of this matrix is to encode the generalized eigenvalue problem for the MANOVA pair $(B,W)$ in ordinary eigenvalue form: if $Bv = \theta Wv$ with $v \neq 0$, then multiplying by $W^{-1}$ gives $Av = \theta v$.
The hypothesis says that the nonzero eigenvalues of $A$ are exactly $\theta_1,\dots,\theta_s$, counted with algebraic multiplicity. Also, since $W^{-1}$ is invertible, multiplication by $W^{-1}$ does not change rank, so
\begin{align*}
\operatorname{rank}(A)
= \operatorname{rank}(W^{-1}B)
= \operatorname{rank}(B)
= s.
\end{align*}
Thus $A$ has exactly $s$ nonzero eigenvalues and its remaining $p-s$ eigenvalues are $0$.
Now pass from $A$ to $I_p+A$. If $Av = \lambda v$ for some nonzero vector $v \in \mathbb{R}^p$, then
\begin{align*}
(I_p+A)v = v + \lambda v = (1+\lambda)v.
\end{align*}
Thus each eigenvalue $\lambda$ of $A$ contributes the eigenvalue $1+\lambda$ of $I_p+A$, with the same algebraic multiplicity. Applying this to the list of eigenvalues of $A$, the eigenvalues of $I_p+A$ are
\begin{align*}
1+\theta_1,\dots,1+\theta_s,
\underbrace{1,\dots,1}_{p-s\ \text{times}}.
\end{align*}
Since the determinant of a square matrix is the product of its eigenvalues counted with algebraic multiplicity, we obtain
\begin{align*}
\det(I_p + W^{-1}B)
= \prod_{k=1}^s (1+\theta_k)\cdot 1^{p-s}
= \prod_{k=1}^s (1+\theta_k).
\end{align*}
[/guided]
[/step]
[step:Convert canonical roots into squared sample canonical correlations]
For each $k \in \{1,\dots,s\}$, the relation between the canonical root $\theta_k$ and the squared sample canonical correlation $\hat\rho_k^2$ is
\begin{align*}
\hat\rho_k^2 = \frac{\theta_k}{1+\theta_k}.
\end{align*}
Since $B$ is positive semidefinite and $W$ is positive definite, each canonical root satisfies $\theta_k \geq 0$, so $1+\theta_k > 0$. Rearranging gives
\begin{align*}
1-\hat\rho_k^2
= 1 - \frac{\theta_k}{1+\theta_k}
= \frac{1}{1+\theta_k}.
\end{align*}
Multiplying over $k = 1,\dots,s$ gives
\begin{align*}
\prod_{k=1}^s \left(1-\hat\rho_k^2\right)
= \prod_{k=1}^s \frac{1}{1+\theta_k}
= \frac{1}{\prod_{k=1}^s(1+\theta_k)}.
\end{align*}
[/step]
[step:Combine the determinant factorization with the canonical-correlation factors]
From the determinant factorization,
\begin{align*}
\frac{\det W}{\det T}
= \frac{1}{\det(I_p+W^{-1}B)}.
\end{align*}
From the eigenvalue product formula,
\begin{align*}
\det(I_p+W^{-1}B)
= \prod_{k=1}^s(1+\theta_k).
\end{align*}
Therefore
\begin{align*}
\frac{\det W}{\det T}
= \frac{1}{\prod_{k=1}^s(1+\theta_k)}
= \prod_{k=1}^s \left(1-\hat\rho_k^2\right).
\end{align*}
This is the stated Wilks' lambda product formula.
[/step]