  Supermassive black hole binaries (SMBHBs) are expected to be detected by the
future space-based gravitational-wave detector LISA with a large
signal-to-noise ratio (SNR). This prospect enhances the possibility of
differentiating higher harmonics in the inspiral-merger-ringdown (IMR)
waveform. In this study, we test the ability of LISA to identify the presence
of different modes in the IMR waveform from a SMBHB. We analyze the
contribution of each mode to the total SNR for different sources. Higher modes,
in particular the mode $(3, 3)$ and $(4, 4)$, can dominate the signal observed
through the LISA detector for SMBHB of the order of $10^8 M_\odot$. With
Bayesian analysis, we can discriminate models with different harmonics. While
spherical harmonics are often considered orthogonal, we observe it is not the
case in the merger-ringdown phase observed by LISA. Omitting harmonics not only
diminishes the SNR but can also lead to biased parameter estimates. We analyze
the bias for each model in a source example and quantify the threshold SNR
where we can expect the parameter bias to be comparable to the statistical
error. By computing the waveform model error with the Fisher approximation and
comparing it with the posterior distribution from our sampler results, we can
evaluate the veracity of the analytical bias, which converges with the sampler
results as more harmonics are introduced. To conclude, SMBHB events with SNR of
a few hundred, as expected in LISA, are required to use templates with at least
modes $(2, 2)$, $(2, 1)$, $(3, 3)$, $(3, 2)$, $(4, 4)$, $(4, 3)$ to estimate
all intrinsic parameters correctly. Our work highlights the importance of
higher modes to describe the gravitational waveform of events detected by LISA.
