As mentioned by Deane Yang in the comments and his (deleted) answer, one can estimate the components of the metric in normal coordinates using a transport ODE (I know it from Dolgov-Khriplovich (1983) Eq.(34), and also in more implicit form from Florides-Synge (1971), but perhaps there is a more canonical reference). The way I see the question, it does not ask to estimate the derivatives of the metric, only its pointwise deviation from the Euclidean metric $\delta_{ij}$. So I see this strategy as completely viable.
The radial ODE satisfied by the metric $g_{ij}(y)$ in normal coordinates is
$$\begin{aligned}
r \partial_r (r\partial_r +1) (g(ry) - \delta)_{ij}
&= \frac{1}{2} g^{kl}(ry) [r\partial_r(g(ry) - \delta)_{ik}] [r\partial_r (g(ry)-\delta)_{jl}]
- 2r^2 y^k y^l R_{ikjl}(ry) \\
&= \frac{1}{2} g^{kl}(ry) [r\partial_r g_{ik}(ry)] [r\partial_r g_{jl}(ry)]
- 2r^2 y^k y^l R_{ikjl}(ry)
=: r^2 F_{ij}(r, g(ry)-\delta, \partial_r g(ry)) ,
\end{aligned}$$
where $y=(y^i)$ is any fixed coordinate vector. The linear differential operator has the following forward Green function/fundamental solution:
$$
G(s,r) = \Theta(s-r) \frac{1}{r} \left(1-\frac{r}{s}\right) .
$$
To solve the equation with initial condition $g_{ij}(0) = \delta_{ij}$, we can adapt the usual argument using Picard iteration and the Banach fixed point theorem to the integral equation
$$\begin{aligned}
(g(sy)-\delta)
&= \int_0^s dr\, r\left(1-\frac{r}{s}\right) \, F(r, g(ry)-\delta, \partial_r g(ry))
=: \Gamma[g(ry)-\delta] , \\
\partial_s g(sy)
&= \int_0^s dr \, \frac{r^2}{s^2} F(r,g(ry)-\delta,\partial_r g(ry)) ,
\end{aligned}$$
where the second equation is a consequence of the first, which we'll need shortly. Some uniform estimates on $|F|$, coupled with the integrals $\int_0^s dr\, r(1-r/s) = s^2/6$ and $\int_0^s dr\, r^2/s^2 = s/3$, will allow us to find lower bounds on the interval of existence of the solution to the transport ODE.
The original question didn't specify what norm $|g(y) - \delta|$ refers to. I'll just suppose that it is the Frobenius norm with respect to $\delta_{ij}$, $|h| = (h_{ik} h_{jl} \delta^{ij}\delta^{kl})^{1/2}$, which is quite natural for dealing with tensors in normal coordinates. To get anywhere with estimating $|F|$, it helps to assume some a priori uniform bounds on $|g(ry) - \delta| \le D < 1$, $|\partial_r g(ry)| \le C$. The constants $C$, $D$ are for now arbitrary, but they can be tuned to specific values later. The first useful uniform estimate is
$$
|g^{kl}(ry) [r\partial_r g_{ik}(ry)] [r\partial_r g_{jl}(ry)]| \le \frac{C^2}{1-D} .
$$
Next, the curvature term satisfies
$$
|y^k y^l R_{.k.l}(ry)| \le |R(ry)| \|y\|^2 ,
$$
where $\|y\|^2 = |yy| = y^k y^l \delta_{kl}$ and
$$
|R(ry)| = (R_{ijkl}(ry) R_{i'k'j'l'}(ry) \delta^{ii'} \delta^{kk'} \delta^{jj'} \delta^{ll'})^{1/2}
$$ is the Frobenius norm of the endomorphism $h^{kl} \mapsto h^{kl} R_{ikjl}(ry)$. It would be nice to express the bound on the Riemann tensor in terms of some geometric quantity, but $|R(ry)|$ only has meaning inside the normal coordinate chart. However, replacing each $\delta^{ij}$ by $g^{ij}(ry)$ we get an invariant curvature scalar and we can essentially bound one with the other. Writing $\delta^{ij} = g^{ij}(ry) - (g^{ij}(ry) - \delta^{ij})$, noting that $|g^{..}(ry)-\delta^{..}| \le \frac{D}{1-D}$, and repeatedly using the Cauchy-Schwarz inequality, we get
$$\begin{aligned}
|R(ry)| &\le \frac{\|R\|(ry) \|y\|^2}{(1-D)^2} \le \frac{\mathcal{R} \|y\|^2}{(1-D)^2} , \\
\|R\| &= \max\{R^{(1)}, R^{(2)}, R^{(3)}, R^{(4)}\} , \\
R^{(1)} &= (R_{ijkl} R^{ijkl})^{1/2} , \\
R^{(2)} &= (R_{ijkl} R^{i'jkl} R_{i'j'k'l'} R^{ij'k'l'})^{1/4} , \\
R^{(3)} &= (R_{ijkl} R^{kli'j'} R_{i'j'k'l'} R^{k'l'ij})^{1/4} , \\
R^{(4)} &= (R_{ikjl} R^{ki'lj'} R_{i'k'j'l'} R^{k'il'j})^{1/4} , \\
\mathcal{R} &= \sup_{x\in M} \|R\|(x) .
\end{aligned}$$
The global geometric invariant $\mathcal{R}$ of $(M,g)$ is what we will eventually use to determine the constants sought in the original question.
The main upper bound on the size of the interval $s\in [0,t]$ is that the iteration must map $(g(ry)-\delta) \mapsto \Gamma(g(ry)-\delta)$ to the space satisfying the same a priori uniform bounds. Using the above estimates on $|F|$ these inequalities are expressed as
$$\begin{aligned}
\frac{t^2}{6} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R \|y\|^2}}{(1-D)^2}\right) &\le D , \\
\frac{t}{3} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R} \|y\|^2}{(1-D)^2}\right) &\le C ,
\end{aligned}$$
which has a consistent solution $D=1/4$, $C=(2\mathcal{R}\|y\|^2/3)^{1/2}$, $t=(3\mathcal{R}\|y\|^2/8)^{-1/2}$ (no claim about the optimality of $t$). So, the integral equation $(g(ry)-\delta) = \Gamma[g(ry)-\delta]$ and the above parameters for the a priori estimates gives use the uniform bound
$$
|g(ry)-\delta| \le \frac{r^2}{6} \left(\frac{1}{2} \frac{C^2}{1-D} + 2\frac{\mathcal{R} \|y\|^2}{(1-D)^2}\right) \le \frac{2}{3} \mathcal{R} \|ry\|^2
$$
for $\|ry\| \le (3\mathcal{R}/8)^{-1/2}$. This section on Wikipedia explains how to go from here prove the contraction property of a sufficiently high power of $\Gamma$ to prove existence.
Since now we can just replace the normal coordinate vector $ry^i$ by $y^i$ we get what was hopefully the desired estimate valid for any normal coordinate chart on $M$:
$$
\|y\| \le \min\left\{\sqrt{\frac{8}{3\mathcal{R}}}, \rho_I\right\}
\implies
|g(y) - \delta| \le \frac{2}{3} \mathcal{R} \|y\|^2 ,
$$
where $\rho_I$ is the injectivity radius of $(M,g)$ (positive by assumption), which has to be there because some normal coordinate charts cannot exceed the radius $\rho_I$ anyway.
As expected, we get a larger coefficient $2/3 > 1/3$ than in the Taylor formula. It might be possible to optimize this coefficient and the bound on $\|y\|$ a bit more, if needed. Also changing the Frobenius norm to a different one could also affect the constant and the definition of $\mathcal{R}$ in terms of global geometric invariants of $(M,g)$.