Determination of Personalized IOL-Constants for the Haigis Formula under Consideration of Measurement Precision

The capabilities of a weighted least squares approach for the optimization of the intraocular lens (IOL) constants for the Haigis formula are studied in comparison to an ordinary least squares approach. The weights are set to the inverse variances of the effective optical anterior chamber depth. The effect of random measurement noise is simulated 100000 times using data from N = 69 cataract patients and the measurement uncertainty of two different biometers. A second, independent data set (N = 33) is used to show the differences that can be expected between both methods. The weighted least squares formalism reduces the effect of measurement error on the final constants. In more than 64% it will result in a better approximation, if the measurement errors are estimated correctly. The IOL constants can be calculated with higher precision using the weighted least squares method.


Introduction
Cataract is the major cause of blindness world wide [1]. It is characterized by a clouding of the natural human lens. Its treatment requires surgical removal of the lens. The lens can then be replaced by an artificial intraocular lens (IOL). Proper preoperative calculation of the refractive power of the IOL is important in order to achieve the desired visual outcome. State of the art IOL calculation requires at least three parameters: the axial length of the eye, the refractive power of the cornea and an estimation of the IOL position after surgery [2][3][4]. If a distinct refractive outcome is desired, the target refraction will also be required as a parameter. The Haigis formula requires the distance between the corneal epithelium and the natural lens (the phakic anterior chamber depth) as an additional parameter [5,6].
The estimation of the postoperative IOL position based on preoperative biometry (refractive power of the cornea, axial length, anterior chamber depth) is a major source of error in the calculation of IOL power [7]. Several formulae can be used to estimate the appropriate IOL power [2-5, 8, 9]. The design parameters provided by the manufacturers (lens constants) allow for an initial estimation of the best suited IOL power, but have to be refined in order to reduce the deviations between planned and achieved post-surgical refraction [10,11]. These lens constants are optimized according to the postoperative outcomes observed with a specific IOL, in order to minimize the numerical errors for the sample.
Some measurements may include larger measurement errors than others. Measurements that have larger errors should contribute less to the refinement of the parameters than those with smaller errors. The statistical measurement uncertainty results from natural variations in the shape of the eye, uncontrolled lens accommodation and pupil size, variations in the axis of fixation during data acquisition, and measurement resolution [12]. The uncertainty of the axial length L measurement contributes approximately 2.8m −1 mm −1 (additional refractive power in m −1 at spectacle plane per measurement error in mm), and the uncertainty of the anterior chamber depth ACD measurements −1.4m −1 mm −1 to the error of the IOL power calculation [7]. Measurement uncertainty is significantly reduced through the proper use of low coherence interferometry for the measurement of the axial length [7,13,14]. Nevertheless, it can be beneficial to prepare the data by adequate screening [9] for outliers originating from erroneous data transfer or inaccurate measurements. An approximation method that takes the probability-distribution of measurement errors into account such as a maximum likelihood approach might also increase the predictive potential of personalized IOL-constants.
The Haigis formula [5,6] is based on a simplified thin lens model of the cornea using only the keratometry values of the anterior cornea to calculate the effective corneal refractive power K m defined as the average over both keratometry measurements using the keratometer index of n c = 1.332. The postoperative optical anterior chamber depth d does not necessarily correspond to the physiological post-surgical anterior chamber depth. It is defined as the parameter d in the the prediction of the lens power D [5,6] that on average results in the smallest difference between the planned refraction R x and the achieved post surgical refractive outcome F. The index of refraction of the aqueous humour is accepted as n = 1.336 [15]. The phakic anterior chamber depth ACD and the axial length L are used to predict the postoperative anterior chamber depth d of the thin lens calculation (Eq 1). The optical anterior chamber depth d can be estimated with the Haigis formula [5,6] given by The constants a 0 , a 1 , and a 2 in the Haigis formula characterize each IOL-type. The postoperative optical anterior chamber depth can be derived by solving Eq 1 for d given the post surgical refractive outcome F. The planned refraction R x in Eq 1 is replaced by F. Substitution of the effective refractive power allows to write the solution of Eq 1 in the simple form In clinical practice the effective optical anterior chamber depth d is calculated from Eq 4 for a large number of patients and the corrected IOL-constants for the Haigis formula are calculated by ordinary linear least squares fitting to Eq 2. The ordinary least squares fitting procedure is sensitive to outliers and does not take into account the measurement uncertainties for L, ACD, K m , D, or F. In order to improve the robustness of the optimization of lens constants against measurement uncertainty, we propose to use a weighted least squares method. In this paper, the weighted least squares method is tested and compared to the ordinary least squares approximation with the help of simulated measurement noise for two different biometry devices. The method is demonstrated for the optimization of the IOL constants for the Haigis formula, but can be applied to other IOL power formulae as well.

Materials and Methods
The ordinary least squares method searches for the minimum of The values of ACD i , L i , and d i are subject to measurement noise. In a weighted least squares method, the values with higher uncertainties are given smaller weights by multiplying Eq 5 with the reciprocal of the variances s 2 The variances can be calculated by error-propagation The ordinary least square solution that minimizes w 2 olsq (Eq 5) is used as an estimate for the IOL-constants in the calculation of the variances.
The performance of the approximation model can be tested with the help of simulated measurement noise. The program work flow is illustrated in Fig 1: An ordinary least-squares approximation (Eq 5) is performed on a real data-set. An artificial system of realistic measurements is obtained by replacing the ACD i values with using the expression for d i given in Eq 4. The original ordinary least-squares approximation is thus the ideal solution for this artificial system and has no approximation error. Then, noise with a Gaussian distribution whose width is given by the corresponding statistical uncertainty of each measurement is added onto each of the terms ACD i , L i , D i , K i m , and F i . The effective optical anterior chamber depths d i are recalculated with these perturbed values according to Eq 4. The approximation of the IOL constants for the Haigis formula by the ordinary (Eq 5) and the weighted least squares (Eq 6) method is conducted 100000 times. Each time the result is affected by the random noise and one method may achieve better results than the other. The comparison of both methods therefore requires many repetitions.
One method can be considered more effective than the other when its results deviate less from the ideal values in more than half of all 100000 cases. This means that the variances where the a k 0 , a k 1 , a k 2 refer to the solution of the perturbed system, are smaller. Another important measure is the square-root of the mean of Eq 9 over all iterations k. Robust fitting reveals itself in a narrow distribution of constants a k 0 , a k 1 , a k 2 (that are not constant in this case, but can vary each iteration k).
These simulations are performed using the repeatability values for two different biometers; the IOLMaster 700 (Carl Zeiss Meditec AG, Jena, Germany) and the Aladdin (Topcon Corp., Tokio, Japan). The repeatability for the Zeiss instrument was studied by Kunert et al. [16] and Srivannaboon et al. [17]. The weighted averages of both values according to the number of studyparticipants are used in this study. They are σ ACD = 9.4 μm, σ L = 10 μm, and σ K m = 0.11 m −1 . The values for the Topcon instrument were derived by Huang et al. [18]. The average values of the measurements for the cataract group are used here: σ L = 20 μm, σ ACD = 45 μm, and σ K m = 0.11 m −1 . The statistical uncertainty for the assessment of the final refraction F is assumed to be σ F = 0.2 m −1 . In accordance with the principles of modern quality assurance, the tolerance interval for the IOL power should be 6 standard deviations [19]. The width of the Gaussian noise σ D for IOL-power D is thus set to 1/6 th of the total tolerance interval defined by the International Organization for Standardization (ISO) [20] which depends on the absolute value of D (Table 1). This is a collaborative work between Saarland University and the University Hospital in Vienna (AKH Wien). The data sets for this study were recorded routinely from patients before and after cataract surgery at AKH Wien by Leydolt and Menapace using the Zeiss IOLMaster 700 device. All data were anonymized and de-identified before they were sent to Saarland University for analysis. Since this is a retrospective study of routinely recorded data, no ethics committee approval was required.
The simulations were performed using N = 69 measurements of the EyeCeeOne NS-60YG lens (Croma Pharma GmbH, Leobendorf, Austria). Finally, the method was applied to an independent data set consisting of N = 33 measurements on patients that had the Acrysof SN60WF (Alcon Inc., Ft. Worth, TX, USA) implanted. This is intended to show the differences that can be expected between both methods. The refraction was measured six months postoperatively. Mean values and standard deviations of the relevant measurements in both data sets are shown in Tables 2 and 3.

Results
The first data set contains data of patients which having an EyeCeeOne NS-60Y implant. It comprises N = 69 valid equations. The original ordinary least squares approximation results for the IOL-constants is the ideal solution given in Table 4.
Firstly, the uncertainty values of the IOLMaster 700 were used. In 64.9% of all 100000 cases the weighted least squares procedure proved to be superior to the ordinary least squares fitting procedure according to Eq 9. The development of this ratio with iteration number is shown in Fig 2a. The root mean square error (RMSE) given by the square-root of the mean of the 100000 s 2 d k given by Eq 9 is RMSE wlsq = 40.97 μm ± 0.05μm, which is better than the RMSE value for an ordinary least squares approach RMSE olsq = 44.34 μm ± 0.06 μm. The RMSE fitting error of the weighted fit approach is smaller, and its distribution (Fig 3a) has a smaller width compared with the ordinary least squares fitting approach. The median value of the square-root of the s 2 d k is 36.0 μm with the weighted least squares, and 38.8μm with the ordinary least square method. The distributions of the IOL-constants for the Haigis formula with the weighted approach have smaller widths than those calculated with the ordinary least squares method (Fig 4). The mean and standard deviation of the calculated constants are shown in Table 4. The calculations are repeated using the uncertainty values of the Aladdin. In 64.3% of all 100000 cases the weighted least squares procedure proved to be superior to ordinary least squares fitting. The development of this ratio with iteration number is shown in Fig 2b. The RMSE errors are RMSE wlsq = 41.72 μm ± 0.056 μm, and RMSE olsq = 45.03 μ ± 0.061 μm for the weighted and ordinary least squares method respectively. The RMSE fitting error of the weighted fit approach is smaller, and its distribution (Fig 3b) has a smaller width compared with the ordinary least squares fitting approach. The corresponding median values of the square-root of the s 2 c k are 36.7 μm for the weighted fits and 39.4 μm for the ordinary least squares fits. The distributions of the IOL-constants for the Haigis formula with the weighted approach have smaller widths than those calculated with the ordinary least squares method (Fig 4). The mean and standard deviation of the calculated constants are shown in Table 4.  The second data set comprises 33 valid equations. The resulting IOL-constants are shown in Table 5 together with their corresponding weighted RMSE (wRMSE) values given by the root mean square of (d i − a 0 − a 1 Á ACD i − a 2 Á L i )/σ i . The uncertainty σ z of the effective refractive power z makes the greatest contribution to the weights s À2 i (Eq 7). The mean value and standard deviation of each contribution are shown in Table 6.

Discussion
The weighted least squares approximation method provided superior results in most cases compared with the ordinary least squares method. More than 64% of its results produce smaller deviations from data that have been cleaned from measurement error than the ordinary least squares method, and the distributions of the IOL-constants for the Haigis formula show smaller widths. Thus, the weighted least squares approximation method delivers results with higher precision than the ordinary least squares method. This is expected as the weighted fit reduces the influence of data with higher uncertainty, and the uncertainties were set correspondingly. This holds true for the uncertainty values of both the Aladdin and the IOLMaster biometer. The results are not limited to these lens types or biometers. Similar improvements can be expected for other IOLs and biometers.
A data set consisting of 69 valid equations was used in order to simulated the influence of Gaussian measurement noise on the approximation. Usually, the IOL constants are optimized on a larger data basis. The data set provides realistic values which were modified before each of the 100000 iterations. The height number of iteration enables us to report very precise results (< 0.1%) concerning which method is superior for the eyes in this data set. Their statistical data is shown in Tables 2 and 3. The differences in the resulting constants between both methods tested on an independent second data set are small compared to the statistical error. The wRMSE value of the weighted least squares approach is smaller for the weighted least squares method. Bigger sample size is required in order to study possible differences in the resulting constants.
In applying this method to calculate the personalized IOL constants it is important that the statistical uncertainties of axial length, anterior chamber depth, corneal refraction,  postoperative spherical equivalent, and IOL-power measurements are estimated. The uncertainties of the biometer have only minor influence. This can be seen when the distributions of the Gaussian noise is set according to the uncertainties of one biometer and the weights according to those of another. For comparison, we used data of the EyeCeeOne NS-60YG lens, and set the distributions of the Gaussian noise according to the statistical measurement uncertainties of the Aladdin, but the weights according to the statistical measurement uncertainties of the IOLMaster. The results were identical with those obtained using the statistical uncertainties of the Aladdin for both the Gaussian noise and the calculation of the weights. The weighted least squares method is superior in 64.3% of all cases. Consequently, the weighted least squares approximation approach can be helpful even when the exact values of the statistical uncertainties are unknown. The uncertainty is clearly dominated by the contribution of the effective corneal power z (Eq 3). The statistical uncertainty of the measurement of the postoperative refraction σ F has the strongest influence.
The statistical uncertainties of the measurements are assumed to have Gaussian distributions. The weighted least squares approximation method performs well in the presence of Gaussian noise. However, in the presence of outliers the weighted least squares method might be improved by eliminating them prior to the fit or implementing a robust approximation method that might include a more sophisticated description of the measurement uncertainties, for example by applying lower weights to particularly long eyes. A more realistic description of the uncertainties of the refractive power of the lens D, or the uncertainties of the pseudophakic refraction F, might improve the performance of the weighted least squares approximation on real data.

Conclusion
An ordinary and a weighted least squares approach for the refinement of the IOL-constants for the Haigis formula were compared. The differences in the resulting constants are small. However, using a weighted fitting method offers superior robustness against random measurement error and should be used with the correct uncertainties. If the uncertainties are estimated correctly the weighted least squares approximation will have a probability higher than 64% of delivering superior results compared to ordinary least squares fitting. Further studies could test the weighted least square approach for larger data sets and/or different formulae. The procedure can be extended by applying more sophisticated models for the statistical measurement uncertainties. The dominant contribution to the statistical uncertainty in this approximation model originated from the statistical uncertainty of the postoperative refraction measurement. Supporting Information S1 File. EyeCeeOne data. The biometric information for patients who got the EyeCeeOne NS-60YG implanted. Right eyes are indicated with OD, left eyes with OS. The table shows the steep  and flat keratometry values (D1, D2)