If you're OK with "almost every", we can use a $L^2$ method and multiplicative character sum bounds to do better.

Let's estimate $$ \sum_{c \in \mathbb F_p^\times} \left( \left[\sum_{a,b \in I | ab=c } 1\right] - p^{1-2\epsilon}\right)^2 $$

by Plancherel $$= \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} \left | \sum_{a,b \in I} \chi(ab) \right|^2 $$

and multiplicativity$$= \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} \left | \sum_{a\in I} \chi(a) \right|^4 $$

Now assume we have a nontrivial character sum bound $\left| \sum_{a \in I} \chi(a) \right| \leq N$ for all nontrivial characters $\chi$. Then we continue:

$$ \leq \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} N^2 \left | \sum_{a\in I} \chi(a) \right|^2 \leq \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee} } N^2 \left | \sum_{a\in I} \chi(a) \right|^2 = N^2 p^{1-\epsilon} $$

Then we see that

$$\left[\sum_{a,b \in I | ab=c } 1\right] \leq 2 p^{1-2\epsilon} $$

for all but

$$\frac{ N^2 p^{1-\epsilon}}{ p^{2-4 \epsilon} } = N^2 p^{3\epsilon-1}$$

This is $\ll p$ as long as $N \ll p^{1 - 3\epsilon/2}$. Using the Pólya-Vinagradov bound where $N \approx p^{1/2}$, this works for $\epsilon <1/3$. Using the Burgess bound we can do slightly better, but not get all the way to $\epsilon = 1/2$, which would require $N= p^{1/4}$ - full square root cancellation, outside of reach.