Parallel form filter design
A few points first


Since you are using digital filters, you should be using residuez() not residue()
Matlab's residuez() is not particularly well implemented: in case of higher order or steep filters, you can get a lot of numerically noise or results that are just totally wrong. Always check the result against your original target.


Back to the question: each residue comes as a conjugate complex pair (at least if you are lucky and don't have to deal with real or multiple poles). Just add these to form a single real second order section.

$$H(z) = \frac{r}{1-p \cdot z^{-1}} + \frac{r^{*}}{1-p^* \cdot z^{-1}} $$
$$= \frac{r \cdot (1-p^* \cdot z^{-1})+r^* \cdot (1-p^* \cdot z^{-1})}{(1-p \cdot z^{-1}) \cdot (1-p^* \cdot z^{-1})}$$

$$= \frac{(r + r^*) - (r \cdot p^* + r^* \cdot p  ) \cdot z^{-1} }{1 -(p+p^*)\cdot z^{-1}+ p \cdot p^* \cdot z^{-2}}$$

The last equation  gives you directly the real valued biquad coefficients (if I did the math correctly,  that is).