The wurtzite structure In1−x−yGaxAlyN quaternary system is studied with respect to the unstable region in mixing. The composition in the unstable region is calculated from the free energy of mixing by using the strictly regular solution model. The interaction parameter used in this calculation is obtained by using the delta-lattice-parameter method. Here, the proportionality constant connecting the lattice constants and the band-gap energy is determined by fitting the calculation to the composition data obtained experimentally from InGaN grown by metallorganic vapor phase epitaxy. From this calculation, the ternary alloys of InAlN, InGaN and GaAlN are predicted to always, sometimes, and hardly ever, respectively, have an unstable mixing region. The essential mismatch in thermal equilibrium between the strictly regular solution approximation and the growth conditions in MOVPE is removed by using a fitting calculation and experimental data. Also, the mismatch between the zinc-blende structure and the wurtzite structure is corrected. As a result, this prediction of the phase separation in In1−x−yGaxAlyN becomes more reliable.