Recent gravitational wave observations show evidence for the presence of higher harmonics, thus possibly indicating that these waves were generated in the inspiral of compact objects with asymmetric mass ratios. Signals with higher harmonics contain a trove of information that can lead to a better estimation of system parameters and possibly to more stringent tests of general relativity. Gravitational wave model that include higher harmonics, however, have only been developed within general relativity, while models to test theory-agnostic deviations from general relativity have been purely based on the signal's dominant mode. We here extend the parameterized post-Einsteinian framework to include the ℓ = 2, 3 and 4 higher harmonics to first post-Newtonian order, therefore providing a ready-to-use Fourier-domain waveform model for tests of general relativity with higher harmonics. We find that the deformations to the higher harmonics of the Fourier phase can be easily mapped to the deformation of the dominant harmonic, while the deformations to the higherharmonics of the Fourier amplitude in general cannot in a theory-agnostic way. Nonetheless, we develop a simple ansatz for the deformations of the waveform amplitude (through a re-scaling deformation of the time-domain amplitude) that both minimizes the number of independent amplitude deformations parameters and captures the predictions of all known modified theories to date.