Since we are working within the injectivity radius, we can use the exponential map to define the geodesic normal coordinate system. More precisely, let $S^{n-1}$ denote the unit sphere in $T_pM$, a neighborhood of $p$ (minus the point $p$ itself) can be identified with $(0,r_*)\times S^{n-1}\ni (r,\omega)$ via $(r,\omega) \mapsto \exp_p(r\omega)$. The Gauss Lemma states that the metric $g$ takes the form
$$ g = dr^2 + r^2 \gamma $$
where $\gamma = \gamma(r)$ is a one parameter family of Riemannian metrics on $S^{n-1}$, such that it suitably limits to the standard metric as $r\to 0$.
$\newcommand\dvol{\mathrm{dvol}}$ I will write $\dvol_\gamma$ for the surface measure on $S^{n-1}$ corresponding to the metric $\gamma$, and $\dvol$ for the volume measure of the metric $g$, which we note is equal to $r^{n-1} \mathrm{d}r ~\dvol_\gamma$.
Finally, we observe that $\partial_r$ (defined in the obvious way) is the unit outward normal for the ball $B_\rho$, which is just the set $B_{\rho} = \{(r,\omega): r< \rho\}$.
By the divergence theorem, we have
$$ \int_{B_\rho} \triangle \phi ~\dvol = r^{n-1}\int_{S^{n-1}} \partial_r\phi ~\dvol_\gamma \Big|_{r = \rho} $$
"Differentiation by parts" gives us
$$ = r^{n-1} \partial_r \Big( \int_{S^{n-1}} \phi ~\dvol_\gamma \Big) - r^{n-1} \int_{S^{n-1}} \phi ~\partial_r (\dvol_{\gamma} ) \Big|_{r = \rho}$$
Assume we have the claim that $\partial_r(\dvol_{\gamma}) = h\cdot \dvol_{\gamma}$ for some $h < 0$. Then under the assumption of your question, we have that $\triangle \phi \leq 0$ and $\phi \geq 0$ yielding
$$ 0 \geq r^{n-1} \partial_r \Big( \int_{S^{n-1}} \phi ~\dvol_\gamma \Big) $$
Using now that
$$ \lim_{\rho\to 0^+} \Big( \int_{S^{n-1}} \phi ~\dvol_{\gamma} \Big) \Big|_{r = \rho} = |S^{n-1}| \phi(p) $$
we can conclude that
$$ \int_{S^{n-1}} \phi ~\dvol_{\gamma} \Big|_{r = \rho} \leq |S^{n-1}| \phi(p) $$
Now you just need to multiply both sides by $\rho^{n-1}$ and integrate from $0$ to $r$ to get the desired result.
The key step is to prove that $h \leq 0$.
We will use the following standard results. (See Peter Petersen, Riemannian Geometry, second edition, Chapter 9, Proposition 39.)
Proposition Let $r$ (as above) be the Riemannian distance function from a point $p$, we have that
- The Lie derivative $L_{\partial_r} \dvol = \triangle r\cdot \dvol$
- $\partial_r (\triangle r) + \frac{1}{n-1} (\triangle r)^2 \leq - \mathrm{Ric}(\partial_r, \partial_r)$.
(As an aside, the second statement is a version of what is called the Raychaudhuri equation in the relativity literature.)
Since $\dvol_{\gamma} = r^{1-n} \iota_{\partial_r} \dvol$, we have that the quantity $h$ can be expressed as
$$ h = \triangle r - \frac{n-1}{r} $$
by the first point of the proposition.
Next, observe that
$$ \triangle (r^2) = 2 r\triangle r + 2 \nabla r \cdot \nabla r $$
As we approach $r = 0$, we have that $\triangle (r^2)$ converges to $2n$, using that the Hessian of $r^2$ evaluated at $p$ is exactly the metric.
This means that
$$ r \triangle r \to n-1$$
as $r\to 0$.
Now, apply point 2 of the proposition. In view of the Ricci sign condition assumed, we can rewrite it as
$$ - \partial_r \frac{1}{\triangle r} + \frac{1}{n-1} \leq 0 $$
Integrating from $r = 0$ where $\frac{1}{\triangle r}$ tends to $0$ we find
$$ \triangle r \leq \frac{n-1}{r}$$
and hence $h \leq 0$ as claimed.
The result above is essentially the "differential" version of the Bishop comparison theorem.