It is usually said that Dirac got his equation by looking for the square root of the 4-momentum norm (see Dirac’s coop here). The relativistic 4-momentum norm is
$$(E)^2-(\mathbf{p}c)^2=(mc^2)^2 \tag{1}$$
A quick look at this equation, one immediately spots that the left-hand side is similar to the square of a binomial, so, the simplest choice for the square root of this equation is
$$aE+\mathbf{b \cdot p}c=mc^2 \tag{3}$$
As long as the operators $a$ and $\mathbf{b}$ obeys the following rules: (a) $a^2=1, \, b_i^2=-1$ ($i=1,2,3$ is 3D-space index), (b) they all anti-commute and (c) are constant. Replacing $E$ and $\mathbf{p}$ in eq. (4) by their respective operators, one immediately gets the Dirac equation
$$i \hbar a ^ \mu \partial _\mu \Psi - mc^2 \Psi =0. \tag{4}$$
Where we changed the notation for $a$ and $\mathbf{b}$ by the new operators $a ^\mu$ with $\mu=0,1,2,3$ a spacetime index. Using the new notation, the rules above can be summarized by
$$\{a^\mu,a^\nu \}=2\eta_{\mu \nu} \mathbf{1} \tag{2}$$
To my knowledge it’s because of the Dirac's choice of these operators (he choose the $4 \times 4$ Dirac-matrices) that equation (4) is deemed to be valid only for spin-$\frac{1}{2}$ particles. However, from a deeper perspective $a^\mu$ (a.k.a $\gamma ^ \mu$) are just any set of objects as long as they obey conditions (2) which do not include any restrictions on form, nature or dimensions on them, but just demand that cross-terms of eq. (3) to vanish and squared term of $\mathbf{p}$ to be negative.
So, why do we say that Dirac equation is only valid for spin-$\frac{1}{2}$ particles if one could plug a different set of operators (maybe matrix representations of Clifford algebras) representing another spin?. Are the Dirac’s $\gamma$-matrices the only set of operators that when inserted into equation (3) leads to equation (1)?