The answer seems to be positive and actually at least in the context of regular (locally) noetherian Deligne--Mumford stacks. (Actually, Artin stack should also be enough as we can compute the Brauer group also as the fppf-cohomology of $\mathbb{G}_m$.)
Let $p\colon X \to \mathcal{X}$ be an \'etale cover. We obtain a descent spectral sequence
$$H^i(X^{\times_{\mathcal{X}}j}; \mathbb{G}_m) \Rightarrow H^{i+j}(\mathcal{X}; \mathbb{G}_m).$$
Likewise we obtain a spectral sequence
$$H^i((X-p^{-1}(\mathcal{Z}))^{\times_{\mathcal{X}}j}; \mathbb{G}_m) \Rightarrow H^{i+j}(\mathcal{X}-\mathcal{Z}; \mathbb{G}_m).$$
As \'etale maps are smooth of relative dimension zero, $X - p^{-1}(\mathcal{Z})$ has still codimension at least $2$ in $X$ and likewise for all higher fiber products. Moreover, under your regularity hypotheses purity also holds for the group of units and for the Picard group. Thus, the map $\mathcal{X} - \mathcal{Z} \hookrightarrow \mathcal{X}$ induces an isomorphism of $E_2$-terms of spectral sequence for $i\leq 2$ and all $j$. Thus, it also induces an isomorphism on convergenda for $i+j\leq 2$, which is exactly the statement we want.