I am extremely late to the party, but here is another more direct approach to this problem.
Consider the map $$N:\mathbb{H}\rightarrow \mathbb{H}:q=x+yi+zj+tij\mapsto \overline{q}=x-yi-zj-tij.$$ It's easy to see that $N$ is an anti-algebra automorphism of $\mathbb{H}$ of order two, i.e. $N$ is an involution on $\mathbb{H}$. Recall that if $A$ is an algebra, $A^{\text{op}}$ is an algebra by reversing the multiplication. Using the map $N$ we get an algebra isomorphism $\mathbb{H}\cong \mathbb{H}^{\text{op}}$.
It's well known that $A\otimes_k A^{\text{op}}\cong \text{End}_k(A)\cong M_n(k)$ for finite-dimensional simple algebras $A$, where $n=\dim_k(A)$. Thus $\mathbb{H}\otimes_{\mathbb R} \mathbb{H}\cong \mathbb{H}\otimes_{\mathbb R} \mathbb{H}^{\text{op}}\cong \text{End}_{\mathbb{R}}(\mathbb{H})=M_4(\mathbb{R})$. This shows that the quaternions have order two in the Brauer group of $\mathbb{R}$.
The above approach is very natural and in fact this technique shows that the inverse of an element in a Brauer group is given by the opposite algebra. For the quaternions the opposite happens to be isomorphic to the original algebra. This approach is also very natural since the map $N$ is simply the norm on $\mathbb{H}$. It can be used to calculate inverses, indeed, the inverse of $q$ is $\frac{\overline{q}}{N(q)}$.