Geometric least squares means ratios for the analysis of Plasmodium falciparum in vitro susceptibility to antimalarial drugs

Background The susceptibility of microbes such as Plasmodium falciparum to drugs is measured in vitro as the concentration of the drug achieving 50% of maximum effect (IC50); values from a population are summarized as geometric means. For antimalarial drugs, as well as for antibiotics, assessing changes in microbe susceptibility over time under drug pressure would help inform treatment policy decisions, but no standard statistical method exists as yet. Methods A mixed model was generated on loge-transformed IC50 values and calculated geometric least squares means (GLSM) with 90% confidence intervals (CIs). In order to compare IC50s between years, GLSM ratios (GLSMR) with 90%CIs were calculated and, when both limits of the 90%CIs were below or above 100%, the difference was considered statistically significant. Results were compared to those obtained from ANOVA and a generalized linear model (GLM). Results GLSMRs were more conservative than ANOVA and resulted in lower levels of statistical significance. The GLSMRs approach allowed for random effect and adjustment for multiple comparisons. GLM was limited in the number of year-to-year comparisons by the need for a single reference year. The three analyses yielded generally consistent results. Conclusion A robust analytical method can palliate inherent limitations of in vitro sensitivity testing. The random effects GLSMRs with adjustment for multiple comparisons and 90%CIs require only assumptions on the mixed model to be applied. Results are easy to display graphically and to interpret. The GLMSRs should be considered as an option for monitoring changes in drug susceptibility of P. falciparum malaria and other microbes.


Background
Plasmodium falciparum, the species causing most of the malaria burden in the world, can be grown in culture, which makes it possible to measure parasite drug susceptibility in vitro [1]. Results are generally expressed as IC 50 (the concentration achieving 50% of maximum effect). Turnidge et al [2] presented a new method to evaluate cut-off values to define antibacterial resistance. However, while methodologies are well established for antimicrobials [3], it is not possible, for the majority of drugs, to confidently classify a strain as resistant or sensitive because validated thresholds are not available, unlike the situation with antibiotics [4].
Of particular interest in P. falciparum drug susceptibility is not the sensitivity profile of an isolate from a given patient, but rather the monitoring of the sensitivity patterns to drugs when submitted to drug pressure in a given area [5]. This information may contribute to inform decision as to choice of drugs for the treatment of malaria. However, evaluating trends over time poses methodological problems, mainly because (i) the high variability of results makes conventional statistical methods, such as the analysis of variance (ANOVA), inadequate; (ii) the statistical unit is not the individual patient isolate but the composite population of parasites within an individual. Due to this distribution variability, these data are usually presented as geometric mean. A key issue in the analysis of these data is detecting trends over time. However, there is still no standard statistical approach for this [6], without loss of information.
Methods used in pharmacokinetic studies deal with similar problems in dealing with datasets with respect to drug disposition in an individual or a population of individuals [7][8][9]. Bioequivalence between drug formulations is typically assessed by using the analysis of variance of log etransformed data from a cross-over design with the null hypothesis expressed by: H 0 : μ T = μ R where μ T and μ R represent the log-transformed expected bioavailability parameters of the test and reference formulations respectively [10,11]. A 90% confidence interval (CI) for the ratio test/reference for bioavailability parameters is constructed by using the equation: where S is the square root of the mean square error from the analysis of variance, n is the number of subjects per period, t 0.05 (1) is the critical value of t at α = 0.05 with υ the degrees of freedom.
In the case of susceptibility testing conducted over time, when subjects are grouped by the year of sampling, the variance in the ANOVA of parasite susceptibility (expressed as IC 50 ) can be separated out into the contributions of the parasite, the subject and the time of measurement. The ANOVA will allow for the year of measurement and test whether the variability between years occurs at random or not. The ratio of the mean sum of squares of the parameter (e.g. IC 50 ) to the error mean sum of squares in the ANOVA will give an F-statistic to test the null hypothesis: H 0 : μ year = μ year0 . This will provide a test of whether the arithmetic mean of IC 50 s measured from a given year is identical to the mean of IC 50 s obtained in the reference year.
An underlying assumption in order to use the ANOVA is the normality of the residuals (the difference between an individual value and the mean of the sample it belongs to, x i -.) However, verifying the null hypothesis of identity between means may not be possible with distributions of residuals that are generally neither normal nor log-normal. Furthermore, the sample size is most often too small with respect to the variance for the comparison of means. The high variability of data leads to large error variability in terms of error sum of squares. Thus, the detection of a difference will be difficult to interpret since it will be a function of the variability of data and sample size for each year. Mixed models can be used when there are different levels of clustering in the observations. One can assume that there is a grouping per year and a random part of measurements within years due to the subject and the parasites strains the subject is infected with. It allows the user to analyse samples (here: years) with unequal sample sizes and to relax the assumption of independently and identically distributed residuals while accounting for the data structure in a more flexible way [12].
For malaria, with the introduction of new treatment regimens such as the Artemisinin-containing Combination Therapies (ACTs) [13], it is important to evaluate whether, with parasites being exposed to drug pressure, the amounts of drug needed to inhibit parasite growth departs from that of a reference year, prior to and during deployment to monitor the evolution of drug susceptibility. In addition, appropriate statistical methods are needed to account for the variability in the determination of IC 50 s in order to properly inform treatment policy decisions.
Several approaches were explored to describe the trends over time of parasite in vitro drug susceptibility of the parasite using a dataset collected in Casamance (south-western Senegal) during 2000-2004 slightly adapted for the purpose from Agnamey et al [14].

Mixed model
A fixed effect model can be expressed as y ij = μ + t j + e ij where j is the year of measurement, y ij the observation on year j for patient I, μ the overall mean (also referred to as intercept in statistical softwares), t j the relative effect of year j, μ + t j is the mean effect for the year j and e ij the residual variance for year j on the i th patient [12]. The model allows for the patient effect, in which case the formula becomes y ij = μ + p i + t j + e ij where p i represents the i th patient effect. Instead of defining some effects as constants in the model, one could consider them as arising from independent samples with a normal distribution [12], i.e. x as random effects. The model containing both fixed and random effects can then be referred to as a mixed model [12,15]. The underlying assumptions to using mixed models are: normally distributed residuals, normally distributed random effects and residuals independent of the random effects.

GLSMRs calculations
A mixed linear model of log e -transformed values was estimated whereby the year was considered as fixed and the intercept as a random effect. The isolates were from different subjects each year. GLSM ratios (GLSMRs) were calculated for each betweenyear comparison. An adjustment for multiple comparisons was done in order to control for the overall type 1 error rate using the Tukey-Kramer method (chosen because it allows for unequal sample size between years). GLSMRs were considered statistically different if both bounds of the CIs fell on either side of the value of 1 (or 100% in percentage values). Previously [14], we had used GLSMRs calculated without this adjustment and evaluated the 95%CIs.

Standard statistical methods
Standard methods such as the one-way ANOVA were also used with the year as fixed factor to analyse the variations of IC 50 s over time. For non-normally distributed data (significant Kolmogorov-Smirnov test), a log e -transformation was applied. The Levene test for homogeneity of the variance was used and, in case of non-homogeneity, a Welch adjustment for the ANOVA. A non-parametric Kruskall-Wallis sign rank test was used when parametric analyses were not suitable. Pair-wise mean comparisons between years for each treatment were carried out following ANOVA with a Tukey adjustment. Normality of residuals was checked with a non significant Shapiro-Wilk test and normal probability plots. Concurrently a generalized linear model (GLM) was also estimated with the year as fixed factor using a normal probability function and an identity link function parameterization [16].
A p-value of < 0.05 was considered statistically significant. All tests were two-tailed. Statistical analyses were carried out with the statistical package SAS ® System version 9.1.3 (SAS Institute, Cary, NC, USA).

Dataset
The dataset is an update of the one described in Agnamey et al [14]. Briefly, in vitro susceptibility of local isolates to chloroquine (CQ), quinine (QN), artemisinin (ART) and the amodiaquine metabolite monodesethylamodiaquine (MdAQ) were monitored the using the DELI test [17] before (1997) and during the deployment (2000)(2001)(2002)(2003)(2004) of artesunate+amodiaquine combination in Mlomp, a village in the district of Oussouye in Casamance, Southern Senegal, where malaria is mesoendemic (25 infective bites/person-year) and transmission occurs year-round with a peak during the rainy season (July-December). Samples for the in vitro assay were from consecutive patients recruited as part of an observational study [18] with a P. falciparum mono-infection and parasitaemia ≥ 0.2% [19].

Results
IC 50 s from 242 subjects for CQ, 250 subjects for QN, 236 subjects for MdAQ and 183 subjects for ART were used in these statistical analyses [14]. The number of subjects with in vitro results was different among years and products tested (Table 1). Means with two standard deviations along with data distributions are plotted in Figure 1. Mean, Geometric means and GLSMs are presented together in Table 1. For all products except log-transformed ART, values were not normally distributed (Kolmogorov-Smirnov test p < 0.05).

Results of the analyses using the different methods
The ANOVA pair wise means comparisons (  Using the GLM with 1997 as the reference year (Table 3), no relationship was found for CQ and ART, a significant positive estimate for MdAQ (IC 50 s decreased) and a significant positive estimate for QN in 2001, followed by a negative estimate (IC 50 s increased) in 2002.
GLSMR (Table 4) of CQ IC 50 s showed no significant differences, consistent with a stable response to CQ over the study period (

Discussion
In this study, three different statistical methods to assess changes of IC 50 over time (ANOVA, GLM, GLSMRs) were compared. The use of data from a single site of moderate to high transmission (25 infecting bites per person-year), and with consistent treatment policies and practices, meant that all patients were expected to be infected with parasites having been under the same degree of drug pressure.

GLSMR
The GLSMR had broader applicability than the other methods because, even if a mixed model is used to obtain the LSMs, it does not require either normally distributed IC 50 s or variance homogeneity. However, one needs to verify the assumptions needed for a mixed model such as the normality of the residuals, the normality of the random effects, and the independence of the residuals and the random effects. As log e -transformation serves the purpose of deriving GLSMs of the results, the linear mixed   model allowed to relax the assumption of independence of the model residuals and to account for the inherent variability of the data structure in a more flexible way. The in vitro test reflects the susceptibility of the whole parasite population in one isolate, and thus cannot separate the effects of the various parasite clones in a given infected subject. Therefore one cannot parameterize the withinsubject effect or within-parasite population effect in the linear mixed model. However, it can be argued that these effects are taken into account in the overall mean, which was defined as random in the model. GLSMRs calculated at a 5% level without an adjustment for multiple comparisons are more likely to wrongly detect a significant difference because of the multiplicity of the statistical tests performed than the GLSMRs at a 1% level and with an adjustment for multiplicity.

Linear models
ANOVA and pair wise means comparisons between years were good indicators of a difference between years. Noteworthy, in bioequivalence studies GLSMR are generally computed using an ANOVA model, but effects specified as random in a linear model are treated as a fixed factor as they serve the sole purpose of producing the corresponding expected mean squares [20]. These expected mean squares lead to the traditional ANOVA components without accounting for the random effect in the variance.
In the GLSMRs calculations, the mixed linear model was computed using restricted maximum likelihood to evaluate variance parameters, which are in general preferred to ANOVA estimates [20]. Furthermore, mixed models are commonly used when there are different levels of clustering in the observations. The sole level of grouping was the year (treated as fixed in the model), so no other particular variable (or grouping level) was defined as random. The reasons for treating the time variable as discrete (i.e. by calendar year and not as a continuous variable) are: (i) there was a 3-year gap between 1997 (the baseline) and 2000 (the first of a series of five consecutive years.); (ii) only qualitative or discrete variables allowed for the model to simply extract estimates for the different categories of the year effect; (iii) the majority of malaria cases and treatments cluster between July-November during the wet season. Hence, subjects were grouped by year and it was assumed that there was a random part of measurements within years due to the contribution of the subject and the parasite strains subjects were infected with. Specifically, the number of isolates for each year was not the same and IC 50 s varied considerably from year to year ( Table 1).
The GLM appears not to be suitable to compare IC 50 s because it treated values as if there were repeat measures from the same subjects, while isolates came from different individuals. In addition, IC 50 s were non-normally distributed despite log e -transformation.

Estimates of the IC50
It is clear that customary approaches are not satisfactory as the difficulty in the analyses is that the data are not time series or longitudinal data. They are also not normally distributed -a necessary condition to use parametric statistical tests, and a critical point of this work. In a recently published paper, Kaddouri et al [21] developed a new inhibitory sigmoid Emax statistical model to estimate more precisely the IC 50 of a given subject. However, it requires cut-off values for resistance of the studied parasites strains to a range of treatment, an element which is not easy to derive for antimalarial drugs. A Bayesian approach was also proposed recently to provide a correction of the estimate of the true IC 50 [22]. This work was based on the assumption that resistance is systematically overestimated because: (i) the precision of the estimated IC 50 value of the most resistant isolate will usually be the poorest of all the isolates assayed, (ii) sigmoid curve fitting or probit analysis of a unique isolate takes no account of other isolates in the series tested. This approach requires the distribution of measured IC 50 to be estab-  lished before using it to produce a large number of points and estimate the median as the true IC 50 . In doing so, one is inevitably faced with the problem of normality and the choice to apply data transformations such as log or Box-Cox [23]. This work offers an alternative way to deal with data transformation and normality condition in the context of the evolution of microbial susceptibility to drugs.

Expression of results
GLSMRs were found to be more intuitive, as results are expressed as percentage difference (increase or decrease) between two years, while with LSMs comparisons increases and decreases from the reference appear as inverted (they are marked with a negative and a positive sign, respectively). As geometric means are generally used to express IC 50 s of a pool of isolates, GLSMR are naturally easier to understand and interpret than the other statistical methods tested here.
Two different plots were also produced in order to illustrate the difficulties in interpreting trends over time in drug susceptibility. Figure 1 is a traditional way of plotting the distribution of IC 50 s with cut-off values above which a strain is resistant to a given drug. This makes the reader falsely interpret means against the cut-off, while there is a great variability in the data and no validated thresholds for the majority of antimalarial drugs. In addition, this display does not provide any indication of the significance (or lack thereof) of change between years.

41
LGHQWLW\ parison of GLSM between years; it does not need the cutoff values; it depicts visually a statistical test as the confidence intervals of each GLSMR against the line of identity between sets of data (here: years).

Conclusion
Treatment policy decisions would benefit from reliable information on changes in susceptibility of parasite or bacterial isolates to drugs over time. This entails an adequate statistical method, which can also account for the inherent variability of in vitro drug susceptibility tests ( Figure 2). This is particularly important for antimalarial drugs and cases alike where validated thresholds for resistance are not available. The underlying linear mixed model of GLSMRs allowed accounting for this variability and for unequal number of isolates collected during field testing. Based on these data GLSMRs appear to be more accurate and to offer advantages over other tests for the "longitudinal" analysis of IC 50 s. We used a simple statistical model which produces easily interpretable results and can be found in any statistical software. Finally, the utility of GLSMRs in monitoring drug susceptibility of not only malaria parasites but also other microbes should be further tested.