  The laser-interferometer space antenna (LISA) will be launched in the mid
2030s. It promises to observe the coalescence of massive black-hole (BH)
binaries with signal-to-noise ratios (SNRs) reaching thousands. Crucially, it
will detect some of these binaries with high SNR both in the inspiral and the
merger-ringdown stages. Such signals are ideal for tests of General Relativity
(GR) using information from the whole waveform. Here, we consider
astrophysically motivated binary systems at the high-mass end of the population
observable by LISA, and simulate their signals using the newly developed
multipolar effective-one-body model: \texttt{pSEOBNRv5HM}. The merger-ringdown
signal in this model depends on the binary properties, and also on parameters
that describe fractional deviations from the GR quasi-normal-mode (complex)
frequencies of the remnant BH. Performing full Bayesian analyses, we assess to
which accuracy LISA will be able to constrain deviations from GR in the
ringdown signal when using information from the whole signal. We find that
these deviations can typically be constrained to within $10\%$ and in the best
cases to within $1\%$. We also show that with this model we can measure the
binary masses and spins with great accuracy even for very massive BH systems
with low SNR in the inspiral. We also probe the accuracy of the
\texttt{SEOBNRv5HM} waveform family by performing synthetic injections of GR
numerical-relativity waveforms. Using a novel method that we develop here to
quantify the impact of systematic errors, we show that already for sources with
SNR $\mathcal{O}(100)$, we would measure erroneous deviations from GR due to
waveform model inaccuracies. One of the main sources of error is the
mismodelling of the relative alignment between harmonics. Our results confirm
the need for improving waveform models to perform tests of GR with high SNR
binary BHs observed by LISA.
