When assessing liver iron content using relaxometry, an average relaxation rate (R1, R2 or R2*) is usually determined from a region of interest or the entire liver. This is commonly performed by fitting the signal decay in individual voxels to an appropriate relaxation function. The voxel-level parameters resulting from the fits are combined to determine the average relaxation rate, and an empirically derived calibration curve is used to convert this single value to iron content. The goal of this study was to compare the precision and accuracy of this voxel-wise fitting to an alternative method that relies on first averaging the signals from all voxels within the region of interest and then determining the relaxation rate from a single fit. Systematic differences were observed when both methods were applied to clinical images. Mathematical simulations were employed to determine which method provided more robust estimates of the true relaxation rate. The mathematical simulations were then expanded to include a range of conditions expected in typical relaxometry images. The results show that voxel-wise fitting skews the relaxation rate estimates and increases variance, particularly when the true relaxation rate is moderate to fast, as it would be in liver with high iron content. The potential impact of these results on clinical decisions is discussed.