Nielsen and Chuang's textbook, Equation 8.101 (section 8.3.4 'Depolarizing Channel') shows that applying a random Pauli to a density matrix representing one qubit equals the identity (times one half):
$$\frac{I}{2}=\frac{\rho+X \rho X+Y \rho Y+Z \rho Z}{4}$$
Does this generalize to density matrices with more than one qubit?
You can generalize the question as follows: given a set of unitaries $X\subset\mathbf U(d)$ and weights $p_U$ normalized as $\sum_{U\in X}p_U=1$, when do we have $$\sum_{U\in X} p_U U\rho U^\dagger=\operatorname{tr}(\rho)\frac{I}{d}\tag1$$ for all states (and more generally Hermitian operators) $\rho$?
There's a few different ways to look at this problem:
As discussed in this answer, if you look at the relation componentwise, you see it amounts to $$\sum_{U\in X} p_U U_{ij} \bar U_{k\ell} = \frac1 d\delta_{j\ell}\delta_{ik},\,\, \forall i,j,k,\ell.$$ In the case where you have $|X|=d$ and $p_U=1/d$, you can see this relation as the completeness relation for the set of unitaries. That is, just like the completeness relation for a set of vectors $\{ v_k\}$ is $\sum_k v_k v_k^\dagger=I$, if you think of the unitaries as vectors in the vector space of matrices with respect to the standard $L_2$ inner product, you get the equation above.
If you pass through the Choi of the channel $\rho\mapsto \sum_U p_U U\rho U^\dagger$, the connection with the completeness relation of the unitaries becomes even clearer, because the Choi of a map $\rho\mapsto A\rho B^\dagger$ is $\operatorname{vec}(A)\operatorname{vec}(B)^\dagger$, and then (1) translates into asking the vectors $\operatorname{vec}(U)$ to satisfy the completeness relation. The relation with the Choi is discussed e.g. in this answer.
Note that these arguments do not characterise the set of solutions to (1), but tell you that a sufficient condition is that the unitaries form (modulo scaling factor) an orthonormal basis of operators. This is precisely the case for Pauli matrices for a single qubit, or products of Pauli matrices for $n$-qubit states.
A direct generalisation of the above argument is to observe that we can go beyond complete orthonormal bases via the formalism of frame theory, which essentially amounts to studying how to do decompositions with possibly overcomplete sets of vectors. This is relevant here because (1) corresponds to asking that the frame of vectors $(\sqrt{p_U}\operatorname{vec}(U))_{U\in X}$ forms a tight frame, which is equivalent to saying that any other vector admits a specific type of simple decomposition in terms of this "overcomplete basis".
Another perspective discussed eg in this answer is to observe that (1) corresponds to saying that the ensemble $\{(p_U,U)\}_{U\in X}$ forms a weighted unitary 1-design. And then a natural generalisation of the notion is to consider unitary $t$-designs. There's plenty to say about these things, starting from the examples cited in the above answer. This perspective is very similar to the one in terms of frames, as unitary 1-designs are essentially the same as tight frames.
Yet another perspective is to look at the relation with the tools of group theory. For this, see e.g. Show that, averaging over uniformly random unitaries, $\mathbb{E}[UXU^\dagger]=\operatorname{Tr}(X)\frac{I}{d}$ and Why does the twirl of a quantum channel give a depolarizing channel?.
Yes, this generalizes to density matrices with more than one qubit as well. But you would need to apply all the elements from the Pauli group corresponding to the number of qubits. For example, if you consider density matrices on two qubits then you need to apply all the elements from the set $\{II,XI,YI,ZI,IX,IY,IZ,XX,XY,XZ,YX,YY,YZ,ZX,ZY,ZZ\}$.
