Extreme mass-ratio inspirals~(EMRIs) detectable by the Laser Inteferometric Space Antenna~(LISA) are unique probes of astrophysics and fundamental physics. Parameter estimation for these sources is challenging, especially because the waveforms are long, complicated, known only numerically, and slow to compute in the most relevant regime, where the dynamics is relativistic. We perform a time-consuming Fisher-matrix error analysis of the EMRI parameters using fully-relativistic numerical waveforms to leading order in an adiabatic expansion on a Kerr background, taking into account the motion of the LISA constellation, higher harmonics, and also including the leading correction from the spin of the secondary in the post-adiabatic approximation. We pay particular attention to the convergence of the numerical derivatives in the Fisher matrix and to the numerical stability of the covariance matrix, which for some systems requires computing the numerical waveforms with approximately $90$-digit precision. Our analysis confirms previous results (obtained with approximated but much more computationally efficient waveforms) for the measurement errors on the binarys parameters. We also show that the inclusion of higher harmonics improves the errors on the luminosity distance and on the orbital angular momentum angles by one order and two orders of magnitude, respectively, which might be useful to identify the environments where EMRIs live. We particularly focus on the measurability of the spin of the secondary, confirming that it cannot be measured with sufficient accuracy. However, due to correlations, its inclusion in the waveform model can deteriorate the accuracy on the measurements of other parameters by orders of magnitude, unless a physically-motivated prior on the secondary spin is imposed.