  We combine the rough path theory and stochastic backward error analysis to
develop a new framework for error analysis on numerical schemes. Based on our
approach, we prove that the almost sure convergence rate of the modified
Milstein scheme for stochastic differential equations driven by multiplicative
multidimensional fractional Brownian motion with Hurst parameter
$H\in(\frac14,\frac12)$ is $(2H-\frac12)^-$ for sufficiently smooth
coefficients, which is optimal in the sense that it is consistent with the
result of the corresponding implementable approximation of the L\'evy area of
fractional Brownian motion. Our result gives a positive answer to the
conjecture proposed in [11] for the case $H\in(\frac13,\frac12)$, and reveals
for the first time that numerical schemes constructed by a second-order Taylor
expansion do converge for the case $H\in(\frac14,\frac13]$.
