The survey of Lyman $alpha$ emitters (LAEs) with Subaru Hyper Suprime-Cam, called SILVERRUSH (Ouchi et al.), is producing massive data of LAEs at $zgtrsim6$. Here we present LAE simulations to compare the SILVERRUSH data. In 162$^3$ comoving Mpc$^3$ boxes, where numerical radiative transfer calculations of reionization were performed, LAEs have been modeled with physically motivated analytic recipes as a function of halo mass. We have examined $2^3$ models depending on the presence or absence of dispersion of halo Ly$alpha$ emissivity, dispersion of the halo Ly$alpha$ optical depth, $tau_alpha$, and halo mass dependence of $tau_alpha$. The unique free parameter in our model, a pivot value of $tau_alpha$, is calibrated so as to reproduce the $z=5.7$ Ly$alpha$ luminosity function (LF). We compare our model predictions with Ly$alpha$ LFs at $z=6.6$ and $7.3$, LAE angular auto-correlation functions (ACFs) at $z=5.7$ and $6.6$, and LAE fractions in Lyman break galaxies at $5<z<7$. The Ly$alpha$ LFs and ACFs are reproduced by multiple models, but the LAE fraction turns out to be the most critical test. The dispersion of $tau_alpha$ and the halo mass dependence of $tau_alpha$ are essential to explain all observations reasonably. Therefore, a simple model of one-to-one correspondence between halo mass and Ly$alpha$ luminosity with a constant Ly$alpha$ escape fraction has been ruled out. Based on our best model, we present a formula to estimate the intergalactic neutral hydrogen fraction, $x_{rm HI}$, from the observed Ly$alpha$ luminosity density at $zgtrsim6$. We finally obtain $x_{rm HI}=0.5_{-0.3}^{+0.1}$ as a volume-average at $z=7.3$.