This paper describes statistical procedures for developing earthquake damage fragility functions. Although fragility curves abound in earthquake engineering and risk assessment literature, the focus has generally been on the methods for obtaining the damage data (i.e., the analysis of structures), and little emphasis is placed on the process for fitting fragility curves to this data. This paper provides a synthesis of the most commonly used methods for fitting fragility curves and highlights some of their significant limitations. More novel methods are described for parametric fragility curve development (generalized linear models and cumulative link models) and non-parametric curves (generalized additive model and Gaussian kernel smoothing). An extensive discussion of the advantages and disadvantages of each method is provided, as well as examples using both empirical and analytical data. The paper further proposes methods for treating the uncertainty in intensity measure, an issue common with empirical data. Finally, the paper describes approaches for choosing among various fragility models, based on an evaluation of prediction error for a user-defined loss function.