I'm trying to find a parameterised affine transform that will transform the unit circle into the family of ellipses shown.
The unit circle passes through (0,1) and (1,0) and is tangential to the lines (1,x) and (x,1). In this question, the circle is transformed into an ellipse that retains those properties. Because the properties only fix 4 parameters and a unique ellipse requires 5, there will be a family of ellipses, defined by a single parameter "p" ("p" being arbitrary within some range), which all share these same 4 properties with the unit circle.
Since a transform from a unit circle to an ellipse (however rotated or translated) is affine, I'd expect to be able to represent the transform from the unit circle to a member of this family of ellipses in this kind of form:
x -> A.x + B.y + C
y -> D.x + E.y + F
where A,B,C,D,E,F are functions of "p".
(It isn't important what "p" actually represents, so long as it parameterises the family of ellipses.)
I'm having difficulty getting a suitable transform for this.
The approach I tried was to look at what happens to 3 easy points under the required transform: (0,1) and (1,0) stay fixed at (0,1) and (1,0), and (0,0) must be mapped due to symmetry, to some other point on the diagonal (Z, Z) where Z < 1/sqrt(2) and Z will often be negative. This gives 6 trivial equations in 6 variables, and as there should be just one solution for any given value of the parameter this should be it.
But when I do this, I get an ellipse that isn't tangential to the 2 lines, I think that's because the lines that the ellipse I get is tangent to at the fixed points, have also been transformed and are no longer horizontal/vertical.
Has my approach inadvertently loosened up the parameters so that I'm getting multiple families? If not, then what's going on, and what do I have to do to fix it?

Your assumption that the two tangent points of the circle will remain fixed is leading you astray, as you’ve discovered. As MvG explains in his answer, the only affine transformation that leaves these points fixed and maps the lines $x=1$ and $y=1$ to themselves is the identity. Instead, we’ll build up the transformation in stages. It will consist of a dilation, a 45° rotation and a translation.
To make the calculations simpler, consider an ellipse in standard position. We want its normal at the point in the first quadrant with $y=1/\sqrt2$ to be parallel to the line $x=y$. From this condition, we get $${x\over a^2}=\frac1{\sqrt2b^2}\quad\text{and}\quad{x^2\over a^2}+{1\over2b^2}=1.$$ Combining these two equations gives a relationship between the major and minor axes of the ellipse: $$a^2=b^2(2b^2-1).$$ Solving for $b$ looks a bit messy, so we’ll take it (the semi-minor axis length) as our parameter. (Note that we must have $b\gt1/\sqrt2$ for an ellipse.) We thus get the dilation matrix $$S=\begin{bmatrix}b\sqrt{2b^2-1}&0&0\\0&b&0\\0&0&1\end{bmatrix}.$$ Next, we compute the required translation. For our point of interest to end up at $(0,1)$ after rotation, we must translate it to $(1/\sqrt2,1/\sqrt2)$, which can be accomplished by translating by $$\frac1{\sqrt2}-\frac1{\sqrt2}{a^2\over b^2} = \sqrt2(1-b^2)$$ along the $x$-axis. Alternatively, we can translate last, in which case the translation will be along the line $x=y$, giving the translation matrix $$T=\begin{bmatrix}1&0&1-b^2\\0&1&1-b^2\\0&0&1\end{bmatrix}.$$
The rotation matrix is, of course, $$R=\begin{bmatrix}\frac1{\sqrt2}&-\frac1{\sqrt2}&0\\\frac1{\sqrt2}&\frac1{\sqrt2}&0\\0&0&1\end{bmatrix}.$$ Composing these three transformations produces the required family of affine transformations: $$M=TRS=\begin{bmatrix}b\sqrt{b^2-\frac12}&-\frac1{\sqrt2}b&1-b^2\\b\sqrt{b^2-\frac12}&\frac1{\sqrt2}b&1-b^2\\0&0&1\end{bmatrix}.$$ If we take $p=2b^2-1$, this becomes $$\begin{bmatrix}-\frac12\sqrt{p(1+p)}&-\frac12\sqrt{1+p}&\frac12(1-p) \\ -\frac12\sqrt{p(1+p)}&\frac12\sqrt{1+p}&\frac12(1-p) \\ 0&0&1\end{bmatrix}$$ for $p>0$, which is a perhaps more aesthetically pleasing range.
The equation for this family of ellipses is $$\begin{bmatrix}x&y&1\end{bmatrix}(M^{-1})^T\begin{bmatrix}1&0&0\\0&1&0\\0&0&-1\end{bmatrix}M^{-1}\begin{bmatrix}x\\y\\1\end{bmatrix}=0.$$ $M^{-1}$ is easily computed by inverting its components individually, resulting in $$\begin{bmatrix}x&y&1\end{bmatrix}\begin{bmatrix}\frac1p&-{p-1\over p(p+1)}&{p-1\over p(p+1)} \\ -{p-1\over p(p+1)}&\frac1p&{p-1\over p(p+1)} \\ {p-1\over p(p+1)}&{p-1\over p(p+1)}&-{3p-1\over p(p+1)}\end{bmatrix}\begin{bmatrix}x\\y\\1\end{bmatrix}=0$$ or, after multiplying through by $p(p+1)$, $$(p+1)x^2-2(p-1)xy+(p+1)y^2+2(p-1)(x+y)-3p+1=0.$$ Alternatively, we can multiply the $M^{-1}\begin{bmatrix}x&y&1\end{bmatrix}^T$ terms out first, yielding $${(y-x)^2\over p+1}+{(x+y+p-1)^2\over p(p+1)}=1.$$ The parameter $p$ turns out to have an interesting meaning: the resulting ellipse will be tangent to the lines $x=-p$ and $y=-p$ at $(-p,1-p)$ and $(1-p,-p)$, respectively, i.e., it gives the other two sides of the bounding box of the ellipse.