  General relativity (GR) has been well tested up to solar system scales, but
it is much less certain that standard gravity remains an accurate description
on the largest, that is, cosmological, scales. Many extensions to GR have been
studied that are not yet ruled out by the data, including by that of the recent
direct gravitational wave detections. Degeneracies among the standard model
($\Lambda$CDM) and modified gravity (MG) models, as well as among different MG
parameters, must be addressed in order to best exploit information from current
and future surveys and to unveil the nature of dark energy. We propose various
higher-order statistics in the weak-lensing signal as a new set of observables
able to break degeneracies between massive neutrinos and MG parameters. We have
tested our methodology on so-called $f(R)$ models, which constitute a class of
viable models that can explain the accelerated universal expansion by a
modification of the fundamental gravitational interaction. We have explored a
range of these models that still fit current observations at the background and
linear level, and we show using numerical simulations that certain models which
include massive neutrinos are able to mimic $\Lambda$CDM in terms of the 3D
power spectrum of matter density fluctuations. We find that depending on the
redshift and angular scale of observation, non-Gaussian information accessed by
higher-order weak-lensing statistics can be used to break the degeneracy
between $f(R)$ models and $\Lambda$CDM. In particular, peak counts computed in
aperture mass maps outperform third- and fourth-order moments.
