We present a statistical study of infrared variability using the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) database for a sample consisting of 1085 high-mass young stellar objects (YSOs) related to 6.7GHz methanol masers. A total of 383 maser sources were identified as NEOWISE variables and classified in two variability behavior classifications: 204 secular (linear, curved, and periodic) and 179 stochastic (burst, drop, and irregular) variables. Statistical analysis of the properties of these variables (e.g., the dust temperature, bolometric luminosity, hydrogen column density, W4 luminosity, and W1-W2 color) has revealed a potential evolutionary sequence among different light-curve types of variables. There is a possible general evolutionary (from less to more evolved) trend between the three variable types from secular to stochastic to nonvariable. For the specific classifications, the evolutionary trend for secular variables is linear to sin to sin+linear, and for stochastic variables it is burst to irregular to drop. These sequences may reflect the evolution of the envelope or accretion disk of high-mass YSOs, from large to small radii due to gravitational collapse. Although no significant variability correlation was established between the 6.7GHz methanol maser and the W2-band emission based on the data collected so far, a number of candidates were found for further investigating the accretion burst events via future variability monitoring programs of both mid-infrared emission and masers.