Why do the degrees of freedom of a density operator not match up with the degrees of freedom of a state vector? The density operator $\rho$ of a mixed 2-qubit system has $4^2-1=15$ degrees of freedom. We can require Tr[$\rho^2$] $ =1$ so that the system is in a pure state. Now we have 14 degrees of freedom.
If we describe the system with a complex 4-dimensional vector we have with the norm restriction $2 \times4-1=7$ degrees of freedom. The state vector even shows a redundant degree of freedom because $|\psi\rangle$ is the same state as $e^{i\phi}|\psi\rangle$. 
So what part of my counting is wrong, or what does it mean, when the density operator and the state vector are supposed to describe the same thing but the density operator has twice the degrees of freedom, although it is the vector state which has an obvious redundant degree of freedom?
 A: For one thing, you are neglecting positivity.
The other - and in some sense, more central - problem is that you pretend that non-linear constraints such as positivity and the non-linear $\mathrm{tr}\,\rho^2=1$ are independent constraints which remove one variable each.
However, neither non-linear constraints nor positivity are so simple. Quadratic equations (in a single variable) can e.g. have different numbers of solutions.
Positivity, on the other hand, is a condition on all eigenvalues, so it adds more than one constraint, similar to a matrix-valued equation $A\vec x=\vec b$ (but not the same!)
Essentially, what you have is conditions on the eigenvalues of the density matrix. If you don't insist on positivity, the eigenvalues $p_n$ have to satisfy $\sum p_n=1$, and $\sum p_n^2=1$.  These equations have many solutions if the $p_n$s can be negative as well.  So you see that positivity is essential to get the unique solution $p_1=1$ and the other $p_n=0$.  At the same time, you can also see that positivity does not just reduce the number of parameters by one, since if all but one $p_n$ are zero, we don't need to specify 3 out of 4 eigenvectors, and thus, many parameters suddenly disappear.
Essence: You can't say that every equation removes one parameter if the equations are not linear, or if they are matrix-valued.
A: Essentially the reason why you see the discrepancy in your calculation is that purity imposes not one constraint but of order of $N^2$.
Let me list some relevant spaces together with their respective dimensions. 


*

*Set of complex $N\times N$ matrices. $N^2$ complex dimension, $2 N^2$ real dimension.

*Subset of Hermitian matrices. These matrices are made of $N$ real numbers on the diagonal and $(N^2-N)/2$ complex number off the diagonal. The total dimension is $N^2$ real dimension.

*Subset of positive matrices. This is also a positive cone. This condition does not change the dimension of the set. Hence this set has still dimension $N^2$ (real). 

*Subset of density matrices. This is the intersection of the positive cone with the hyperplane $\mathrm{tr}(\rho) =1$. Note that this is only one real equation. Hence the dimension of this set is $N^2-1$. A common representation is the generalized Bloch form which makes this explicit. 

*Subset of pure states. Here is where the complication happens. This is due to the fact that purity imposes several constraint, since it is the requirement that  $\rho^2=\rho$. It can be shown that this is equivalent to the fact that $\rho = |\psi\rangle \langle \psi |$.  In this latter form it is much easier to compute its dimension and one sees that this set has dimension $2 N -2$. Note that this is  only a subset of the $N^2 -2$ dimensional boundary of the set of general density matrices. 
A: For a pure state $|\psi\rangle$, the corresponding density operator is the projector
$$\hat\rho=|\psi\rangle\langle\psi|.\tag1$$
The trace property is a necessary, but not a sufficient condition for a density operator to represent a pure state. That's why imposing this condition doesn't reduce enough the number of degrees of freedom for an arbitrary operator, even in a finite-dimensional Hilbert space.
The easiest way to see that, for a finite-dimensional Hilbert space, $\hat\rho$ has exactly the same number of degrees of freedom (up to the phase factor you've already mentioned) as $|\psi\rangle$, is the definition $(1)$. From it you can infer that, if you take a matrix representation of $\hat\rho$, all the rows in this matrix will be linearly dependent (its rank is 1).
