As bhobba already said, in nonrelativistic QT (NOT in relativistic!) the spin commutes with the position and momentum variables. A complete set of free single-particle eigenstates are, e.g., the position-spin eigenstates ##|\vec{x},s,\sigma_z \rangle##, where the spin quantum-number ##s \in \{0,1/2,1,\ldots \}## is fixed. Both, position and spin, are represented by linear self-adjoint operators on an abstract separable Hilbert space. There are no matrices here!
You come to a matrix-differential-operator algebra on the Hilbert space ##L^2(\mathbb{R}^3,\mathbb{C}^{2s+1}##, i.e., the Hilbert space of square Lebesgue-integrable ##\mathbb{C}^{2s+1}##-valued functions when choosing the position representation. So let ##|\Psi \rangle## be a normalizable (true) Hilbert-space vector, then the mapping from the abstract Hilbert space ##\mathcal{H}## to this function-Hilbert space is given by
$$|\Psi \rangle \mapsto \Psi_{\sigma}(\vec{x})=\langle \vec{x},\sigma|\Psi \rangle, \quad \vec{x} \in \mathbb{R}^3,\sigma \in \{-s,-s+1,\ldots,s-1,2 \}.$$
In this representation operators are represented by "matrices". The matrix elements for position and spin-##z## components are very simple to calculate, because we've chosen the common (generalized) eigenvectors as a (generalized) basis for our position-spin wave-mechanics representation:
$$\hat{\vec{x}} \mapsto \langle \vec{x}_1,\sigma_1|\hat{\vec{x}} \vec{x}_2,\sigma_2 \rangle=\vec{x}_2 \langle \vec{x}_1,\sigma_1|\vec{x}_2,\sigma_2 \rangle =\vec{x}_2 \delta^{(3)} (\vec{x}_1-\vec{x}_2) \delta_{\sigma_1 \sigma_2},$$
$$\hat{\sigma}_z \mapsto \langle \vec{x}_1,\sigma_1|\hat{\sigma}_z \vec{x}_2,\sigma_2 \rangle=\sigma_2 \langle \vec{x}_1,\sigma_1|\vec{x}_2,\sigma_2 \rangle=\sigma_2 \delta^{(3)} (\vec{x}_1-\vec{x}_2) \delta_{\sigma_1 \sigma_2}.$$
Thus in the position-spin representation one has
$$\hat{\vec{x}},\sigma |\Psi \rangle \mapsto \sum_{\sigma'=-s}^s \int_{\mathbb{R}^3} \langle \vec{x},\sigma|\hat{\vec{x}} \vec{x}',\sigma' \rangle \Psi_{\sigma'}(\vec{x}') = \vec{x} \Psi_{\sigma}(\vec{x}),$$
$$\hat{\sigma}_z,\sigma |\Psi \rangle \mapsto \sum_{\sigma'=-s}^s \int_{\mathbb{R}^3} \langle \vec{x},\sigma|\hat{\sigma_z} \vec{x}',\sigma' \rangle \Psi_{\sigma'}(\vec{x}') = \sigma \Psi_{\sigma}(\vec{x}).$$
For other operators one has to use the commutator algebra (which follows from the corresponding representation of the Galilei group, which is characterized by the mass as a central charge of this group and the spin-quantum number ##s## of the particle, which is a Casimir operator of the group). E.g., from these considerations you come to
$$\hat{\vec{p}} |\Psi \rangle \mapsto -\mathrm{i} \hbar \vec{\nabla} \Psi_{\sigma}(\vec{x}).$$