Statistical thermodynamics of residue fluctuations of native proteins in a temperature, pressure, and force reservoir is formulated. The general theory is discussed in terms of harmonic and anharmonic fluctuations of residues. The two elastic network models based on the harmonic approximation, the anisotropic network and the Gaussian network models are discussed as the limiting cases of the general theory. The heat capacity and the correlations between the energy fluctuations and residue fluctuations are obtained for the harmonic approximation. The formulation is extended to large fluctuations of residues in order to account for effects of anharmonicity. The fluctuation probability function is constructed for this purpose as a tensorial Hermite series expansion with higher order moments of fluctuations as coefficients. Evaluation of the higher order moments using the proposed statistical thermodynamics model is explained. The formulation is applied to a hexapeptide and the fluctuations of residues obtained by molecular dynamics simulations are characterized in the framework of the model developed. In particular, coupling of two different modes in the nonlinear model is discussed in detail.