Let me first remind you of (or perhaps introduce you to) a couple of aspects of quantum mechanics in general as a model for physical systems. It seems to me that many of your questions can be answered with a better understanding of these general aspects followed by an appeal to how spin systems emerge as a special case.
General remarks about quantum states and measurement.
The state of a quantum system is modeled as a unit-length element $|\psi\rangle$ of a complex Hilbert space $\mathcal H$, a special kind of vector space with an inner product. Every observable quantity (like momentum or spin) associated with such a system whose value one might want to measure is represented by a self-adjoint operator $O$ on that space. If one builds a device to measure such an observable, and if one uses that device to make a measurement of that observable on the system, then the machine will output an eigenvalue $\lambda$ of that observable. Moreover, if the system is in a state $|\psi\rangle$, then the probability that the result of measuring that quantity will be the eigenvalue of the observable is
\begin{align}
p(\lambda) = |\langle \lambda|\psi\rangle|^2
\end{align}
where $|\lambda\rangle$ is the normalized eigenvector corresponding to the eigenvalue $\lambda$.
Specialization to spin systems.
Suppose, now, that the system we are considering consists of the spin of a particle. The Hilbert space that models the spin state of a system with spin $s$ is a $2s+1$ dimensional Hilbert space. Elements of this vector space are often called "spinors," but don't let this distract you, they are just like any other vector in a Hilbert space whose job it is to model the quantum state of the system.
The primary observables whose measurement one usually discusses for spin systems are the cartesian components of the spin of the system. In other words, there are three self-adjoint operators conventionally called $S_x, S_y, S_z$ whose eigenvalues are the possible values one might get if one measures one of these components of the system's spin. The spectrum (set of eigenvalues) of each of these operators is the same. For a system of spin $s$, each of their spectra consists of the following values:
\begin{align}
\sigma(S_i) = \{m_i\hbar\,|\, m_i=-s,-s+1,\dots, s-1,s\}
\end{align}
where in my notation $i=x,y,z$. So for example, if you build a machine to measure the $z$ component of the spin of a spin-$1$ system, then the machine will yield one of the values in the set $\{-\hbar, 0, \hbar\}$ every time. Corresponding to each of these eigenvalues, each spin component operator has a normalized eigenvector $|S_i, m_i\rangle$. As indicated by the general remarks above, if the state of the system is $|\psi\rangle$, and one wants to know the probability that the measurement of the spin component $S_i$ will yield a certain value $m_i\hbar$, then one simply computes
\begin{align}
|\langle S_i, m_i |\psi\rangle|^2.
\end{align}
For example, if the system has spin-$1$, and if one wants to know the probability that a measurement of $S_y$ will yield the eigenvalue $-\hbar$, then one computes
\begin{align}
|\langle S_y, -1|\psi\rangle|^2
\end{align}
Spinors.
In the above context, spinors are simply the matrix representations of states of a particular spin system in a certain ordered basis, and the Pauli spin matrices are, up to a normalization, the matrix representations of the spin component operators in that basis specifically for a system with spin-$1/2$. Matrix representations often facilitate computation and conceptual understanding which is why we use them.
More explicitly, suppose that one considers a spin-$1/2$ system, and one chooses to represent states and observables in the basis $B =(|S_z, -1/2\rangle, |S_z, 1/2\rangle)$ consisting of the normalized eigenvectors of the $z$ component of spin, then one would find the following matrix representations in that basis
\begin{align}
[S_x]_B &= \frac{\hbar}{2}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \frac{\hbar}{2}\sigma_x\\
[S_y]_B &= \frac{\hbar}{2}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} = \frac{\hbar}{2}\sigma_y\\
[S_z]_B &= \frac{\hbar}{2}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} =\frac{\hbar}{2}\sigma_z\\
\end{align}
Notice that these representations are precisely the Pauli matrices up to the extra $\hbar/2$ factor. Moreover, each state of the system would be represented by a $2\times 1$ matrix, or "spinor"
\begin{align}
[|\psi\rangle]_B = \begin{pmatrix} a \\ b\end{pmatrix}.
\end{align}
And one could use these representations to carry out the computations referred to above.
There is a general expression in my article A Compact Formula for Rotations
as Spin Matrix Polynomials, SIGMA 10 (2014), 084, to the effect that, e.g., for the doublet,
\begin{gather*}
e^{i(\theta/2)(\hat{\boldsymbol{n}}\cdot\boldsymbol{\sigma})}=I_{2}\cos{\theta/2}+i(\hat{\boldsymbol{n}}\cdot\boldsymbol{\sigma})\sin{\theta/2},
\end{gather*}
and the triplet, $j=1$, so $J_{3}=\mathrm{diag}(1,0,-1)$,
\begin{gather*}
e^{i\theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}=I_{3}+i(\boldsymbol{\hat
{n}}\cdot\boldsymbol{J})\sin{\theta}+(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})^{2}(\cos\theta-1)
\\
\phantom{e^{i\theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}}
=I_{3}+(2i\hat{\boldsymbol{n}}\cdot\boldsymbol{J}\sin(\theta/2))\cos
(\theta/2)+\tfrac{1}{2}(2i\hat{\boldsymbol{n}}\cdot\boldsymbol{J}\sin (\theta/2))^{2}.
\end{gather*}
For the quartet, $j=3/2$,
\begin{gather}
e^{i \theta (\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} = I_4 \cos (\theta/2)\left(1+\tfrac{1}{2}\sin^2
(\theta/2)\right)+(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin (\theta/2))\left(1+\tfrac{1}{6} \sin^2 (\theta/2)
\right)
\nonumber
\\
\phantom{e^{i \theta (\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}=}{}
+\frac{1}{2!} \bigl (2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J}\sin (\theta/2) \bigr)^2 \cos (\theta/2)+\frac {1}{3!}
\bigl (2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J}\sin (\theta/2) \bigr)^3.
\label{quartet}
\end{gather}
For the quintet, $j=2$,
\begin{gather*}
e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} = I_5+(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta/2))
\cos(\theta/2)\left(1+\tfrac{2}{3}\sin^2(\theta/2)\right)
\\
\phantom{e^{i \theta (\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}=}{}
+\frac{1}{2!} {(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta/2))^2}\left(1+\tfrac{1}{3} \sin^2 (\theta/2)\right)
\\
\phantom{e^{i \theta (\hat{\boldsymbol{n}}\cdot\boldsymbol{J})}=}{}
+\frac{1}{3!} {(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta/2))^3} \cos(\theta /2) +\frac{1}{4!} (2i
\hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin (\theta/2))^4.
\end{gather*}
For the sextet, $j=5/2$,
\begin{gather*}
e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} = I_6 \cos(\theta/2)\left(1+ \tfrac{1}{2} \sin^2
(\theta/2+\tfrac{3}{8} \sin^4 (\theta/2)\right)
\\
\phantom{e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} =}{}
+(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta /2))\left(1+\tfrac{1}{6}\sin^2(\theta/2)
+\tfrac{3}{40}\sin^4(\theta /2)\right)
\\
\phantom{e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} =}{}
+\frac{1}{2!} {(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta/2))^2} \cos(\theta /2)
\left(1+\tfrac{5}{6}\sin^2(\theta/2)\right)
\\
\phantom{e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} =}{}
+\frac{1}{3!} {(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin (\theta/2))^3}\left(1+\tfrac{1}{2}\sin^2(\theta/2)
\right)
\\
\phantom{e^{i \theta(\hat{\boldsymbol{n}}\cdot\boldsymbol{J})} =}{}
+\frac{1}{4!} {(2i \hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin(\theta/2))^4}\cos(\theta /2) +\frac{1}{5!} {(2i
\hat{\boldsymbol{n}}\cdot\boldsymbol{J} \sin (\theta/2))^5}.
\end{gather*}
etc...
There is a simple pattern and compact formula for arbitrary spin detailed in that paper.
Best Answer
As a matter of fact, there is a much easier way to derive those matrices. What you are really after is the matrix representation of the $\mathfrak{su}(2)$ Lie algebra generators $J_x$, $J_y$ and $J_z$ in the irreducible $\mathrm{SU}(2)$ representation with spin $j$, given by
$$ J_z|j,m\rangle = m |j,m\rangle,\quad J_\pm|j,m\rangle = C_\pm(j,m) |j,m\pm 1\rangle, $$ with $$ J_\pm=J_x \pm i J_y $$ and $$ C_\pm(j,m)=\sqrt{(j\mp m)(j\pm m+1)}. $$ The matrix representations of $J_z$ and $J_\pm$ in this basis are $$ J_z=\begin{pmatrix} j\\ & j-1\\ && \ddots\\ &&& -j+1\\ &&&& -j \end{pmatrix} $$ $$ J_+= \begin{pmatrix} 0 & C_+(j,j-1)\\ & 0 & C_+(j,j-2)\\ && \ddots\\ &&& 0 & C_+(j,-j)\\ &&&& 0 \end{pmatrix} $$ $$ J_-= \begin{pmatrix} 0 \\ C_-(j,j) & 0\\ && \ddots\\ && C_-(j,-j+2) &0\\ &&& C_-(j,-j+1) &0 \end{pmatrix}, $$ and you can recover $J_x$ and $J_y$ as $$ J_x=\frac{J_+ + J_-}{2},\quad J_y=\frac{J_+ - J_-}{2i}. $$ For example, with $j=1$ you would get $$ J_z=\begin{pmatrix}1 & 0 & 0\\ 0 & 0 & 0\\ 0& 0 & -1 \end{pmatrix},\quad J_x=\frac{\sqrt{2}}{2}\begin{pmatrix}0 & 1 & 0\\ 1 & 0 & 1\\ 0& 1 & 0 \end{pmatrix},\quad J_y=\frac{\sqrt{2}}{2}\begin{pmatrix}0 & -i & 0\\ i & 0 & -i\\ 0& i & 0 \end{pmatrix}. $$