A generalization of the Majorana representation to $N+1$ dimensional density matrices can be performed as follows:
(By a generalization I mean a representation of the density matrix by means of a certain number of points on the Bloch sphere (not necessarily independent) + a polarization vector belonging to the $N+1$ dimensional probability simplex)
I'll consider first the generic case (multiplicity free eigenvalues) described in the question for completeness.
Every $N+1$ dimensional matrix can be written as:
$\rho = \sum_{i=0}^{N}p_i\theta_i$
where $p$ is the eigenvalue vector $p\in \Delta^n$ (the probability simplex in $\mathbb{R}^{n+1}$)
and $\theta_i$ are the one dimensional projectors on the the eigenvectors:
$\theta_i^2=\theta_i$
satisfying the orthonormality constraints.
$tr(\theta_i\theta_j)=\delta_{ij}$
The Majorana representation
$\mathrm{Sym}^N(\mathbb{C}P^1) \approx \mathbb{C}P^N$
allows to express the first projector in terms of $N$ points on the Bloch sphere, and the second one in terms of $N-1$ points because the orthonormality constraints and the third in terms of $N-2$ points etc.
Dimension count
$2 \times (0+1+ . . .+N) + N =(N+1)^2-1$
The case of an eigenvalue of multiplicity $1<M<N$
In this case the projector on the degeneracy eigenspace is of dimension greater than 1. The orbit of these projectors is the Grassmannian $Gr(M,N)$.
In order to perform a Majorana representation of this projector, we first embed the Grassmannian into a complex projective space by means of a Plucker embedding:
$Gr(M,N) \rightarrow \mathbb{C}P^{{N \choose M}-1}$
and then perform the Majorana map.
Now, both the Majorana map and the Plucker embedding are known explicitely, which makes the the whole construction possible. In the degenerate cases, the dimension count is less than that of the space of all density matrices.
Update:
This update is intended to provide a partial answer for Piotr's comment
Remark: I should have remarked that material on the classification of the density matrix orbits according to their eigenvalue degeneracy in this answer is based on
Geometry of quantum states by Bengtsson and Życzkowski (Mainly on chapter 8)
The case of the Grassmannian $Gr(M,N)$ (which is the unitary orbit of $N$ dimensional density matrices with two distinct eigenvalues, one of the of multiplicity $M
The rigid $SU(2)$ which acts on the qubits of the Majorana representation of $\mathbb{C}P^{{N \choose M}-1}$ is a subgroup of the isometry group $SU({N \choose M})$ of the complex projective space. One should not expect apriori that it is a subgroup of the isometry group $SU(N)$ of the Grassmannian.
However, I checked the the simplest case $Gr(2,4) \rightarrow \mathbb{C}P^5$. This Grassmannian is given by the Plucker relation given by the quadric $w \wedge w = 0$ in the homogeneous coordinates $w = \sum_{i=j=1, i< j}^4w_{ij} e_i\wedge e_j$ of $\mathbb{C}P^5$, ($e_i$ form a basis of $\mathbb{C}^4$).
After the identification of the homogeneous coordinates $w$ with the homogeneous coordinates of the Majorana representation such that on both the action of the $SU(2)$ generator $\sigma_3$ is diagonal, it turned out that the Plucker relation is invariant under the whole $SU(2)$ action. In other words there exists an $SU(4)$ transformation $U$ such that this action can be implemented as: $U^{-1} \theta U$, where $\theta$ is the projector onto the two dimensional degeneracy space of the density matrix. This result is new (and very interesting) to me and I don't have a deep understanding of its reason. I hope it generalizes to all cases.
The generic case of distinct eigenvalues:
In this case, the orbit of the quantum states is the flag manifold $Fl(N)$. In this case there are $N-1$ distinct Bloch spheres corresponding to the hierarchy of eigenvectors as given in the main answer.
A rigid $SU(2)$ action on the qubits of the Majorana representation of each of the $n$-th eigenvector is realized as a spin $\frac{N-n+1}{2}$ representation on the subspace orthogonal to the higher in hierarchy eigenvalues.
In addition there will be a corresponding action on all lower in hierarchy eigenvectors. I think that in order to feel comfortable with the construction, one can work out explicitely the three dimensional case with the flag manifold $Fl(3)$ expressed as a $\mathbb{C}P^1$ bundle over $\mathbb{C}P^2$
This post has been migrated from (A51.SE)