Concentrated hard-sphere suspensions in the liquid and crystal coexistence region show a unique nonlinear behavior under a large amplitude oscillatory shear flow, the so-called strain stiffening, in which the viscosity or modulus suddenly starts to increase near a critical strain amplitude. Even though this phenomenon has been widely reported in experiments, its key mechanism has never been investigated in a systematic way. To have a good understanding of this behavior, a numerical simulation was performed using the lattice Boltzmann method (LBM). Strain stiffening was clearly observed at large strain amplitudes, and the critical strain amplitude showed an angular frequency dependency. The distortion of the shear stress appeared near the critical strain amplitude, and the nonlinear behavior was quantified by the Fourier transformation (FT) and the stress decomposition methods. Above the critical strain amplitude, an increase in the global bond order parameter Ψ(6) was observed at the flow reversal. The maximum of Ψ(6) and the maximum shear stress occurred at the same strain. These results show how strongly the ordered structure of the particles is related to the stress distortion. The ordered particles maintained a bond number of "two" with alignment with the compressive axis, and they were distributed over a narrow range of angular distribution (110°-130°). In addition, the ordered structure was formed near the lowest shear rate region (the flow reversal). The characteristics of the ordered structure were remarkably different from those of the hydroclusters which are regarded as the origin of shear thickening. It is clear that strain stiffening and shear thickening originate from different mechanisms. Our results clearly demonstrate how the ordering of the particles induces strain stiffening in the liquid and crystal coexistence region.