The basic purpose of site response analysis is to evaluate possible peak and spectral accelerations on the ground surface to estimate probable earthquake damage for the existing building stock and for the design of new structures. The basic issues in the site response analysis are the uncertainties in source characteristics, soil profile, soil properties, and site response analysis procedure. In addition, characteristics of the building inventories would introduce critical uncertainties associated with these analyses. Recent advances, with growing computational capacity, emphasize probabilistic frameworks to capture these uncertainties. The probability distribution of the related peak and spectral accelerations on the ground surface may be determined considering all possible input acceleration time histories, site profiles, and dynamic soil properties. One option to account for the variability in earthquake source and path effects may be to consider using large number of acceleration records compatible with the site-dependent earthquake hazard partially based on hazard deaggregation for the investigated site. Likewise, stochastic soil profiles generated via Monte Carlo simulations can be used to account for the site condition variability. A seismic microzonation methodology is proposed based on the probabilistic assessment of these factors involved in site response analyses. The second important issue is the selection of microzonation parameters. The selection of microzonation parameters such as Cumulative Absolute Velocity (CAV) and Housner Intensity (HI) is emphasized for their stronger empirical correlation with structural damage. The main approach is to develop a microzonation procedure for ground shaking intensity accounting for variability in ground motion and soil response. The third issue is the reliability and correctness of the site response analysis procedure. The adopted methodology advances traditional site response analysis by integrating frequency- and stress-dependent soil behavior models to achieve a more accurate numerical model and proposing the use of 3D site response analysis to reflect the multi-directional nature of seismic excitations. It also highlights the need for representative time histories with known exceedance probabilities. Even though the selected representative acceleration time histories may be scaled with respect to probabilistic acceleration spectrum or peak ground acceleration obtained based on probabilistic site response analysis; the probability of the selected acceleration time histories are not known. The only possible option is to estimate probabilistic acceleration time histories with predetermined exceedance probabilities to enhance fully probabilistic site response analysis. The proposed methodology is demonstrated through case studies to underline the importance of fully probabilistic site response analysis in seismic microzonation, aiming to improve the reliability of ground motion predictions and support informed decision-making in earthquake engineering.