The quaternions can be defined as $$\mathbb{R}\langle X,Y\rangle/(X^2+1,Y^2+1,XY+YX)$$ From these relations, it is relatively easy to prove that $1,X,Y,XY$ span the quaternions over $\mathbb{R}$. But I cannot find any way to prove that this is a basis.
The quaternions could alternatively be constructed as the set $\mathbb{R}^4$ together with the product $$(a_1,b_1,c_1,d_1)\cdot(a_2,b_2,c_2,d_2)=(a_1a_2-b_1b_2-c_1c_2-d_1d_2,a_1b_2+b_1a_2+c_1d_2-d_1c_2,a_1c_3-b_1d_2+c_1a_1+d_1b_2,a_1d_2+b_1c_2-c_1b_2+d_1a_2)$$ From this point you can prove the product is distributive and associative to show that it forms a ring. You can also prove the identities used in constructing $\mathbb{H}$ as a quotient of the free algebra on 2 generators. Hence, using the universal property, you could show that it is a quotient of $\mathbb{R}\langle X,Y\rangle/(X^2+1,Y^2+1,XY+YX)$, so $\mathbb{R}\langle X,Y\rangle/(X^2+1,Y^2+1,XY+YX)$ must have dimension at least $4$. So this is technically an answer.
However, this feels very much like going the long way around, so I want to know if there is a neater way to show that $\mathbb{R}\langle X,Y\rangle/(X^2+1,Y^2+1,XY+YX)$ is $4$ dimensional.
The quaternions $\mathbb{H}$ have both a $4$-dimensional real representation embedding them as a subalgebra of $M_4(\mathbb{R})$ as well as a $2$-dimensional complex representation embedding them as a subalgebra of $M_2(\mathbb{C})$. Both of these are just the regular representation of $\mathbb{H}$ acting on the left on itself (where for the complex embedding we need to pick a copy of $\mathbb{C}$ in $\mathbb{H}$ acting on the right), but the point is that you can write these embeddings down explicitly, and then verify that the matrices of $1, X, Y, XY$ are linearly independent. This saves you from having to check associativity.
I'll rewrite the quaternions as $\mathbb{R} \langle i, j \rangle / (i^2 + 1, j^2 + 1, ij + ji)$ and pick the copy of $\mathbb{C}$ corresponding to $\mathbb{R}[i]$ to get the $2 \times 2$ complex representation. Pretend that we already know that the quaternions are $4$-dimensional, or equivalently $2$-dimensional over $\mathbb{C}$, with basis $\{ 1, j \}$. Then multiplication on the left by $i$ sends this basis to $\{ i, ij = -ji \}$ so the corresponding $2 \times 2$ complex matrix is the diagonal matrix
$$\rho(i) = \left[ \begin{array}{cc} i & 0 \\ 0 & -i \end{array} \right].$$
(The $-i$ comes from the fact that we need to think of the copy of $\mathbb{C}$ as acting on the right so that it commutes with left multiplication by elements of $\mathbb{H}$, which is why we needed to rewrite $ij$ as $-ji$.) On the other hand, multiplication on the left by $j$ sends this basis to $\{ j, -1 \}$ so the corresponding $2 \times 2$ complex matrix is just the real rotation matrix
$$\rho(j) = \left[ \begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array} \right].$$
Now we stop pretending that we know the quaternions are $4$-dimensional, and instead we check that $\rho(i)^2 = \rho(j)^2 = -1$, which should be pretty clear, then that
$$\rho(ij) = \left[ \begin{array}{cc} 0 & -i \\ -i & 0 \end{array} \right]$$ $$\rho(ji) = \left[ \begin{array}{cc} 0 & i \\ i & 0 \end{array} \right]$$
so $\rho(ij) + \rho(ji) = 0$, meaning this is really a representation of $\mathbb{H}$. Finally we check that $1, \rho(i), \rho(j), \rho(ij)$ are linearly independent in $M_2(\mathbb{C})$ (over $\mathbb{R}$!) which is pretty straightforward.
We can characterize their real span as the subalgebra of $M_2(\mathbb{C})$ of matrices of the form $\left[ \begin{array}{cc} z & w \\ - \overline{w} & \overline{z} \end{array} \right]$, which should remind you of how $\mathbb{C}$ sits in $M_2(\mathbb{R})$ and can be used as an alternative definition of $\mathbb{H}$ (here the work goes into showing closure under multiplication). The determinant of this matrix is $|z|^2 + |w|^2$ and the inverse is another matrix of the same form so it's even not hard to see that this must be a division algebra.