The compressive kinking behavior of non-slender unidirectional glass fiber-reinforced polymer (GFRP) specimens has been analyzed by finite element (FE) models. The experimentally observed imperfections, including the initial fiber waviness throughout the entire specimen volume and the scattered resin/interface defects, were taken into account in the FE models. The birth-and-death method was employed to simulate the progressive damage to the material. The consideration of the coexistence of initial fiber waviness and initial resin/interface defects was found to be essential for accurate modeling of the kinking failure process. Kinking was initiated due to the disproportional increase of the fiber microbuckling at the locations of initial defects. The numerically obtained peak load, fiber microbuckling amplitudes, kink band angle and width, and compressive strain concentrations at the kink band edges were well predicted compared to the experimental results. The number of defects was less significant than the fact that defects existed that served as initiation points of the kink band formation.