TL;DR: Two steps. First, we find the entanglement fidelity $F_e(\mathcal{P}_{p_x,p_y,p_z})=1-(p_x+p_y+p_z)$. Next, we use Horodecki's formula that connects the average and entanglement fidelities
$$
\overline{F}(\mathcal{E})=\frac{NF_e(\mathcal{E})+1}{N+1}\tag1
$$
to get
$$
\overline{F}(\mathcal{E})=1-\frac23(p_x+p_y+p_z).\tag2
$$
There is also an alternative approach that exploits the fact that Pauli eigenstates form a spherical 3-design.
Entanglement fidelity
There are many ways to quickly find the entanglement fidelity $F_e(\mathcal{E})$ of a channel $\mathcal{E}$. For example, it can be shown directly from the definition that $F_e(\mathcal{E})$ can be read off of the channel's Pauli transfer matrix (aka process matrix) $\chi(\mathcal{E})$. Namely, $F_e(\mathcal{E})$ is the diagonal element corresponding to the identity. In the present case, we have $\chi(\mathcal{P}_{p_x,p_y,p_z})=\mathrm{diag}(1-(p_x+p_y+p_z), p_x, p_y, p_z)$, so $F_e(\mathcal{P}_{p_x,p_y,p_z})=1-(p_x+p_y+p_z)$.
Alternatively, we can use the formula for $F_e(\mathcal{E})$ in terms of the channel's Kraus operators
$$
F_e(\mathcal{E})=\frac{1}{N^2}\sum_i|\mathrm{tr}K_i|^2\tag4
$$
c.f. equation $(9.135)$ on page $421$ in Nielsen & Chuang. For a Pauli channel all but one operators in $(4)$ have zero trace and once again we obtain $F_e(\mathcal{P}_{p_x,p_y,p_z})=1-(p_x+p_y+p_z)$.
Horodecki's formula
You can find a beautiful proof of formula $(1)$ in this paper. The key fact established in the course of the proof is that twirling$^1$ any channel produces a depolarizing channel with the same entanglement fidelity as the original channel.
You can also find a more technical but shorter proof of $(1)$ in this paper.
Spherical 2-designs
Alternatively, we can use a spherical 2-design to simplify the evaluation of the integral
$$
\overline{F}(\mathcal{P}_{p_x,p_y,p_z})=\int\langle\psi|\mathcal{P}_{p_x,p_y,p_z}(|\psi\rangle\langle\psi|)|\psi\rangle d\psi.\tag5
$$
Tetrahedron is the smallest spherical 2-design, but it is also not the most convenient to work with in this context. Instead, we can use an octahedron which is the smallest spherical 3-design. See this paper. This is convenient because the eigenstates of the Pauli operators form an octahedron. We get
$$
\overline{F}(\mathcal{P}_{p_x,p_y,p_z})=\frac{1}{|S|}\sum_{|\psi\rangle\in S}\langle\psi|\mathcal{P}_{p_x,p_y,p_z}(|\psi\rangle\langle\psi|)|\psi\rangle\tag6
$$
where $S=\{|0\rangle, |1\rangle, |+\rangle, |-\rangle, |{+i}\rangle, |{-i}\rangle\}$. We then calculate
$$
\begin{align}
\overline{F}(\mathcal{P}_{p_x,p_y,p_z})&=\frac{1}{6}\left(6\cdot(1-(p_x+p_y+p_z)+2p_x+2p_y+2p_z\right)\\
&=\frac{1}{6}\left(6-4(p_x+p_y+p_z)\right)\\
&=1-\frac23(p_x+p_y+p_z)\tag7
\end{align}
$$
in agreement with $(2)$.
$^1$ Twirl takes a channel and produces another channel by averaging the effect of the conjugation of the input channel by a Haar-random unitary.