Abstract
Background
LaserDoppler imaging (LDI) of cutaneous blood flow is beginning to be used by burn surgeons to predict the healing time of burn wounds; predicted healing time is used to determine wound treatment as either dressings or surgery. In this paper, we do a statistical analysis of the performance of the technique.
Methods
We used data from a study carried out by five burn centers: LDI was done once between days 2 to 5 post burn, and healing was assessed at both 14 days and 21 days post burn. Randomeffects ordinal logistic regression and other models such as the continuation ratio model were used to model healingtime as a function of the LDI data, and of demographic and wound history variables. Statistical methods were also used to study the falsecolor palette, which enables the laserDoppler imager to be used by clinicians as a decisionsupport tool.
Results
Overall performance is that diagnoses are over 90% correct. Related questions addressed were what was the best blood flow summary statistic and whether, given the blood flow measurements, demographic and observational variables had any additional predictive power (age, sex, race, % total body surface area burned (%TBSA), site and cause of burn, day of LDI scan, burn center). It was found that mean laserDoppler flux over a wound area was the best statistic, and that, given the same mean flux, women recover slightly more slowly than men. Further, the likely degradation in predictive performance on moving to a patient group with larger %TBSA than those in the data sample was studied, and shown to be small.
Conclusion
Modeling healing time is a complex statistical problem, with random effects due to multiple burn areas per individual, and censoring caused by patients missing hospital visits and undergoing surgery. This analysis applies stateofthe art statistical methods such as the bootstrap and permutation tests to a medical problem of topical interest. New medical findings are that age and %TBSA are not important predictors of healing time when the LDI results are known, whereas gender does influence recovery time, even when blood flow is controlled for.
The conclusion regarding the palette is that an optimum threecolor palette can be chosen 'automatically', but the optimum choice of a 5color palette cannot be made solely by optimizing the percentage of correct diagnoses.
Background
The use of laserDoppler imaging (LDI) for predicting burn healing time is increasing, and there are several recent reviews that compare it favorably with other techniques ([14]). The basis of the methodology is that the proximate cause of healing is skin blood flow, which is measured by the laserDoppler imager in perfusion units (PU).
The medical decision to be made is whether to allow healing to occur naturally, or, if that would be too slow, to operate to excise the burn and apply a skin graft. As a decision aid, clinicians consider it desirable to predict healing time as being either less than 14 days, 14–21 days, or over 21 days. Hence model predictions are ordinal. There are many papers describing the benefits of LDI as a decision aid for clinicians, e.g. [58]. In essence, LDI can enable predictions of healing time that are correct over 90% of the time, whereas unaided clinician judgement is correct only about 70% of the time [5].
This paper uses data collected by five burn centers, two in the UK, one in Belgium and two in the USA, to develop and validate a probabilistic model of burn healing time, as a function of the laserDoppler flux measurements, and of demographic and observational variables. Besides confirming the performance of laserDoppler methodology, such modeling can address the question of its statistical sufficiency. This is, whether, given laserDoppler blood flow measurements, demographic and observational variables have any further predictive ability, or whether the blood flow measurement recorded as mean flux in PU is a sufficient statistic (see [9] for a definition of this term). The model also allows 'what if' questions to be asked, such as how far predictive ability might degrade on moving to a different patient mix than that seen in the data, for example for patients with a much higher percentage of total burned surface area (%TBSA). Current practice is for the clinician to view a falsecolor image of the region of interest, in which high blood flow areas are colored red and low blood flow regions blue. Figures 1 and 2 give an example. A suitable color palette must be chosen for this purpose, and this problem is also discussed here.
We may expect laserDoppler devices in the future to do more than display blood flow pictorially to aid decisions. In this paper, a probabilistic model is developed that would enable the predicted probabilities of time to heal within 14 days, 14–21 days or over 21 days to be displayed. The clinician would then have a prediction of healing time together with an indication of its probability of being correct. EngKean Yeong et al (2005) [10] have predicted burn healing times using artificial neural networks and a reflectance spectrometer. The work here is in that spirit, but using a 'traditional' probabilistic model, the proportional odds (PO) model.
To enter a caveat: it is not yet practicable to use a mean fluxbased methodology for clinical wound class prediction. This is because real burns have a distribution of flux values and the spatial distribution of these is also important. For these reasons two extra, overlap colors are used for the LDI image colorcoding, instead of discrete class boundaries, to aid clinical interpretation with flexibility. At present, it is necessary to eyeball the flux levels and their spatial distributions, in conjunction with recognising confounding factors such as undebrided dead tissue (the dead epidermis not removed), to make an accurate assessment of wound healing potential.
The next section briefly describes the medical methodology, then the dataset used is described, then the statistical methodology is described and the results of fitting models to the data are given. Next, a statistical study of the 5color palette currently in use is given, followed by conclusions.
Methods
Medical Methodology
The full methodology will be described in a separate clinical publication; briefly, the laserDoppler imager records flux measurements of skin blood flow, measured in PU, for each of typically thousands of pixels corresponding to a wound area. The mean flux measurement with which we shall be largely concerned is the average flux across these pixels.
Burn wounds were assessed only once by laserDoppler imager (moorLDI, Moor Instruments, UK) between 2 and 5 days post burn, and the burns were photographed at this time and also at 14 and 21 days post burn to record the burn at time of LDI and its subsequent healing.
Assessments were made at 14 and 21 days post burn because these times are important for clinical decisions on the burn treatment that is likely to result in least scarring: surgical (skin graft) or conservative (dressings only). A burn wound to most adults that heals within 21 days will do so with minimal scarring: 'Second degree burns that heal within 3 weeks are unlikely to leave scars' [11] ch5, p 70. Exceptions to this include infants, where the risk of infection is higher, and patients of ethnic origins predisposed to hypertrophic scarring. For these groups, surgical management is frequently performed on wounds that are expected to take longer than 14 days to heal [8].
Healing was assessed by clinicians and was defined as a wound with a continuous covering of epithelial cells. The boundary of healed wound at 14 days and 21 days post burn was assessed from the photographs taken at these times. Areas healed and not healed at 14 and 21 days were mapped back on to the original laser Doppler images for later analysis. This was necessary because a wound is usually heterogeneous, with parts healing at different times. These areas were then redrawn on a greyscale photo image that was pixelposition identical with the flux image (obtained simultaneously with the flux image); flux values within the corresponding regions of interest were then exported for analysis.
Exclusions were made where burn wound boundary selections would have had a significant effect on the result, for example a boundary at a steep flux gradient. Prior to analysis, the data used here were screened for clinical factors such as wound infection, tattoos, drugs and concomitant patient sickness; and technical factors such as patient movement, edge effects and reflections.
The day of the LDI scan was determined by clinical factors, staff availability and the study protocol: the protocol restricted the day of scan to within 2 to 5 days post burn based on previous observations by others on the reliability of the LD technique.
Besides time to healing, the usual demographic variables, age and gender, were recorded, also skin pigmentation, which could affect recovery time. The medical history variables relevant to burns studies were total burned surface area, burn site on the body, and burn cause. Treatment and procedural variables were burn center and day of LDI scan.
The clinical decision to perform surgery inevitably censored the healing time in some cases. Some patients went to surgery for some of their burns before day 14 or between days 14 and 21. If the surgery was before day 14, the wound was excluded; if after day 14 we have recorded healing time > 14 days (similar to a patient who did not attend at day 21) unless there was further biopsy evidence. Where possible, biopsies were taken from wounds at surgery and some of these results for healing time have been included in the current analysis. Where the histological analysis found a full thickness wound, these are known to take longer than 21 days to heal; where histology found the wound to be superficial dermal (at the wound edge), these were classified as healing time < 14. Wounds found to be deep dermal were not included because these could heal before or after 21 days and even before 14 days if adjacent tissue had very high flux.
Exploratory Data Analysis
In total, data on 768 burn areas were available, but of these only 310 had mean flux measured, of which 299 were complete, the other 11 having a missing outcome assessment because the patient did not attend one followup visit or because of surgery. The analysis here is focused on the 310 burn areas for which flux measurements were available. These areas came from 100 patients.
Table 1 summarises the univariate statistics, for the 310 cases for which flux data were available. Note that for one of the five burn centers, there were no cases with flux measured. Median age was 32 years, minimum 1 and maximum 88 years, and the %TBSA (total burned surface area) had median 6%, minimum 1% and maximum 68%. A larger sample of 581 burn areas, where flux measurements were not always available, was used to boost the statistics for an exploratory analysis of the relationships among the observational variables.
Table 1. Data summary, showing breakdowns of burn areas from 100 patients.
There are naturally many correlations among these variables. Demographic mix and %TBSA varied across the burn centers. Similarly, mean age varies with burn site.
There are some interesting gender differences. The type of burn site varies a little by gender, and this is statistically significant as shown by chisquared test (χ^{2}[3] = 11.7, p = 0.008), with men relatively more likely to be burnt on the face, and women on the torso. Burn cause varies considerably between the genders (table 2). It is mainly men who experience flash, chemical, electrical and contact burns, presumably because of gendered employment. This difference is very significant (χ^{2}[5] = 75.6, p < 0.001.)
Table 2. Gender differences in burn statistics.
The age and %TBSA distributions are similar. It is also clear that healing time is longer for females than for males. There is however the problem of confounding, for example the slower female healing time could be at least partly the result of genderspecific burn sites and burn causes. This issue will be addressed in the survivaltime modeling.
The slower female healing time is an important fact, and we briefly summarise previous findings in this area. There are a few studies on gender and morbidity, with mortality being the more frequent focus. The effect of gender, in these studies, is considered with other patient and wound variables: e.g. age, race, %TBSA and wound type, inhalation injuries and biochemical markers. Large %TBSA burns are included in the patient groups and these studies indicate that gender does influence outcome. Length of hospital stay and duration in intensive care have been used to assess morbidity. For adults, length of stay was greater for women [12] but in children, duration in the intensive care unit was found to be greater for boys [13]. The child gender difference was also found for mortality, higher for boys than for girls [14]. Adult mortality has been found to be greater in women by twofold [15,16], and there is debate over the age group at highest risk [17,18].
The findings of the current investigation of the effect of gender on healing time are therefore consistent with previous work. However, it must be stressed that there are no studies looking at the residual effect of gender once the laserDoppler measurement is known, which is the main focus of this study
Statistical Methodology
There were 433 burn areas for which clinician's predictions of healing time using the LDI were available. Of these a subset of 310 wounds was appropriate for computing average laserDoppler flux. The analysis here focuses on the use of flux as a predictor rather than using the clinician predictions, mainly because this shows the performance of the technology when quite divorced from clinical judgement. Also, the mean flux measurement data are probably better suited to answering questions about the role of covariates. Proportionalodds (PO) ordinal logistic regression is the most popular method of analyzing data such as these, where a dependent variable Y that takes ordered integer values is modeled as a function of a vector x of covariates. An introductory account is given in [19] and general descriptions in [2022]. We seek to predict a dependent variable (healing time) that is ordered, and model the logits of the cumulative probabilities of healing in under 14 days or in 14–21 days as a linear function of the covariates. In the 'parallel lines' model usually fitted to the kth (of 2) cumulative probabilities P_{1 }and P_{2}, the logit is a function
where only the intercept a is a function of k, and the vector of coefficients β of the covariates x is not. The probabilities of the three healing times, p_{1}, p_{2 }and p_{3 }are given by p_{1 }= P_{1}, p_{2 }= P_{2 } P_{1 }and p_{3 }= 1  p_{1 } p_{2 }= 1  P_{2}.
Surgeons occasionally desire to predict probabilities of healing over different time intervals than the ones used here. We suggest that this could be attempted by regarding α as a function of time t, for example α = γlog(t/t_{0}). Then we can determine γ and t_{0 }from α_{1 }= γlog(14/t_{0}), α_{2 }= γlog(21/t_{0}). With this parameterization, probabilities of healing over two or three different intervals could be found. This model corresponds to a loglogistic distribution of healing time for fixed covariates.
In a later section of this paper, on palette derivation (a discriminant problem) we follow the methodology of [23,24], which dispenses entirely with the need to choose an ordinal model. However, for assessing the significance of covariates using likelihoodbased inference, a model is required. Other models besides PO are also available ([2527]), and there is no overriding reason to choose the PO model. Hence two other models were also fitted, to see whether the PO model could be improved on, the probit and continuation ratio (CR) models [26,28]. The latter model is applicable where, as here, the three data groups are periods. As will be seen, neither of these models fitted better than the PO model. Other models, such as the multinomial model, are useful only when the categories are unordered, and are less efficient than ordinal models when used with ordinal data [29].
The covariates used in this study, besides the LD flux measurement, are those demographic variables that are usually important in medicine, age and gender, plus those known to be important in predicting healing time, such as %TBSA and burn site, and some that might possibly be thought to be relevant, such as skin pigmentation and burn cause. It is thought that these last two are probably not important. For example, natural skin pigmentation does not affect LDI flux from debrided wounds because the pigment is removed with the epidermis.
The proper selection of predictor variables for use in a model raises some statistical problems (e.g., [28], ch. 4. In particular, modeling choices made by the statistician after viewing the data are often not reflected in the final pvalues and confidence intervals produced. Results then appear unduly significant or accurate, and traditional approaches, such as stepwise regression, can give misleading results. This problem looms large when sample sizes are small and there are many variables. Here, fortunately, we have the opposite situation.
The data posed some problems that necessitated a purposewritten computer program. An example is the existence of a few censored cases, where patients did not attend a hospital visit. In these cases, it was known only that healing took place in under 21 days or after more than 14 days. Where surgery was performed, burn severity information was obtained from biopsy results: less than 14 days was assumed for superficial dermal wounds; more than 21 days was assumed for full thickness wounds. Fitting the model by maximumlikelihood estimation meant that these cases could be included. A more serious problem was the existence of multiple burn areas on the same individual. The 310 burn areas occurred on only 100 patients. Although demographic variables that could influence healing time such as age and gender were known, there might be a further 'frailty' that varied between individuals, so that healing times would not even be conditionally independent [30]. Obesity and smoking, for example, affect health, but were not tested here. To model frailty, we add a hypothetical frailty variable, and each individual has a random value of it. A normal distribution for the variable is the simplest choice, and as frailty turned out to be a small effect, no more elaborate modeling was done.
Taking the mean contribution from this variable as zero and the variance as σ^{2 }gave a likelihood function
where the jth burn area on the ith of N individuals has healing time Y_{j }and covariate vector x_{j}. The assumption of a zero mean is nugatory, as a nonzero mean would be absorbed into the constants α_{1 }and α_{2}. For censored data, where for example, healing is only known to occur in under 21 days, the probability in (2) is the probability of this observed event, p_{1 }+ p_{2 }or P_{2}.
The integration could be done by the GaussHermite method as recommended by [31]; in fact adaptive integration using the Numerical Algorithms Group (NAG) integration subroutine D01AMF was used here. The website [32] gives full details of the NAG routines mentioned here.
The loglikelihood function ℓ = ln ℒ was maximised using the NAG function minimisers E04UCF and E04CCF [32]. The former uses a sequential quadratic programming method and is relatively fast, and the latter uses the NelderMeade method and is slow but robust. The best of a large number of random restarts (20) from the current function maximum was used, to ensure that the global maximum of the loglikelihood function had been found.
The asymptotic covariance matrix for fitted model parameters is taken as the inverse of the Hessian matrix where the β_{i }are model parameters, and are the maximumlikelihood estimates. To obtain a more accurate estimate of parameter errors for small samples, the likelihood was maximised for 1000 bootstrapped samples, in which individuals were sampled from the dataset with replacement; see e.g. [33]. Bootstrap resamples that did not allow all parameters to be estimated were rejected. This can happen, if for example, the resample does not contain any examples of a particular burn cause.
Estimates of the model's predictive ability were found by taking the predicted healing time as that with the highest probability. The resulting percentage of correct assignments is however liable to be overoptimistic, because the model is evaluated insample. A recommended method of correcting for this is to calculate the 'optimism' of an insample estimate as the mean of a large number of differences D_{i }between estimates from a bootstrapped dataset evaluated on the bootstrapped dataset and on the original dataset. Finally, the mean optimism is subtracted from the sample estimate (e.g. [28], section 5.2.5, [33]). This procedure is better than simply splitting the sample into two, one for model fitting and one for validation [34]. The predictive ability of a linear model may be measured by the coefficient of determination, R^{2}. For instance, [35] and others proposed a pseudo R^{2 }for general models given by
where β = 0 denotes the 'null' model. Nagelkerke (1991) [36] proposed a correction, since the maximum value of R^{2 }attainable is less than 1. The correction required normalising R^{2 }to its maximum value of 1  exp{2n^{1}ℓ(0)}, and we use this corrected value of R^{2}.
Significance tests that a parameter β_{i }is zero would classically be done using the Wald statistic, i.e. using the estimated standard error of the parameter estimate (see e.g. [28], section 9.2.2). Alternatively, the increase in loglikelihood Δℓ on 'floating' the parameter can be used as a test statistic, when under H_{0 }we have that 2Δℓ ~ χ^{2}[1]. An alternative is to carry out a Wald test using the bootstrapped error estimate. However, an exact test can be obtained by retaining Δℓ as a test statistic, but obtaining its reference distribution under H_{0 }by permuting the relevant variable x_{i }among the cases. The pvalue of the test is the proportion of permutations for which the computed Δℓ exceeds the value for the original sample (see e.g. [37] and [33], chapter 15.). Variables such as gender must of course be permuted among individuals rather than among burn areas. The logic is that under H_{0 }gender is irrelevant, and so the permuted datasets generate the reference distribution for the test statistic. Note that this is not the more elaborate procedure used by [38].
Since the number of distinguishable permutations (combinations) is very large, a random sample of 1000 was generated to compute the pvalue. A random permutation of n labels held in an array was achieved by swapping each array element in turn with a random array element.
The main advantage of a permutation test over a bootstrap hypothesis test is that the permutation test correctly generates the reference distribution of the test statistic under H_{0}, whereas the bootstrap hypothesis test can only generate the distribution of the test statistic under H_{1}. Also, the pvalue obtained is exact, and not an asymptotic approximation. Permutation tests have a simple rationale; under H_{0}, gender labels are irrelevant to healing, and so we can obtain an equivalent dataset by shuffling them. The only drawback here is the extra computing time needed; a separate set of permutations must be carried out for each variable or related set of variables in turn.
In the analyses done here, variable selection was not a problem: given the laserDoppler flux, only gender was significant. This contrasts with many clinical studies, where many covariates are significant, and the choice of best model is not easy.
Results
Data Analysis: LDI performance
Predictors of healing time without LDI measurements
Before considering the laserDoppler measurements as predictors of healing time, it is worthwhile modeling healing time purely as a function of demographic and observational variables. This allows the partial confounding caused by correlation between covariates to be resolved, and in particular answers the question of whether the slower female healing time is solely due to gendered burn site type. Table 3 shows the results of a proportional odds model analysis, using the model from (1) based on 581 regions of interest from 130 individuals. The variables x were dummy variables corresponding to 3 of the 4 burn sites, gender, race, burn cause, and two variables for age effect. The reference categories were burn site 1 (limbs) and scald for burn cause. Analyses were also done including factors such as burn center and day of scan, but any effects due to these were not significant, and table 3 shows an analysis including only the more interesting variables. Note that the joint significance of related variables such as linear and quadratic age terms was also evaluated by likelihood ratio tests (not shown).
Table 3. Demographic and observational predictors of healing based on 581 burn areas, where se* is the bootstrapped standard error of the estimated parameter.
It can be seen that the bootstrapped standard errors of parameter estimates are always larger than the asymptotic theory values that are commonly used. This is particularly marked for variables such as race and some burn causes, where very few individuals may have a particular burn cause. This shows the wisdom of using bootstrapped estimates rather than relying on asymptotic theory. The pvalues in the tables are those of two simple bootstrap hypothesis tests. The first is the bootstrap percentile test, where the resampled distribution of fitted parameters is shifted left by so that it approximates to the distribution of the fitted parameter values under H_{0}. This gives a simple (nonparametric) twosided bootstrap hypothesis test. The second pvalue is derived using the normal approximation and the bootstrapped standard error se*. The analytic approximation sometimes gave much lower pvalues than either of these tests; for example, for gender, the analytic pvalue was 0.0004, as compared to the last two columns of table 3 (pvalues of 0.084 and 0.073).
The Cox and Snell R^{2 }= 0.138, the Nagelkerke R^{2 }= 0.158, and the proportion of correct assignments is only 61%, showing the poor predictive ability of this model.
Healing appears much faster at some burn sites such as extremities and face, and permutation tests showed that the variation in healing time between burn site types is very significant (p < 0.001). Gender was also significant (p = .026). Figure 3 shows the coefficient for gender and its reference distribution derived by permutation of gender labels, together with a fitted normal curve. Age was also significant (p = .024). The effect of age is in the table modeled using a linear and quadratic term, the quadratic term being (age  35)^{2 }rather than age^{2 }for numerical stability. These two terms give a healing time that improves slightly with age up to about age 30, and subsequently decreases. Modeling the effect of age using dummy variables corresponding to age ranges gave similar results. Burn cause was not significant (p = 0.194). Burn center and day of scan were also not significant factors.
Figure 3. The coefficient (vertical dotted line) and its permutation distribution under H_{0 }for gender, in the absence of a laserDoppler measurement, with fitted normal distribution.
Predictions of healing time with LDI measurements
Table 4 presents the results when the laserDoppler flux measurement is included. Permutation tests show that burn site is now not significant (p = 0.78), and neither is age (p = .24). However, gender is still a significant predictor (p = 0.002), and its statistical significance has increased. Aside from gender, the laserDoppler measurement is a sufficient statistic. The Cox and Snell R^{2 }= 0.931 and the Nagelkerke R^{2 }= 0.976, confirming the excellent predictive ability of this model.
Table 4. Predictors of healing, when LDI mean flux measurement is included, based on 310 burn areas for which flux measurements were available.
Including a random effect as in (2), when gender is omitted and only mean flux is used as a predictor, a drop in minus loglikelihood from 73.56 to 69.21 is obtained on including the random effect, which is thus clearly significant (pvalue of significance test is 0.0016). The standard error of the random effect was estimated as 1.91. However, on including gender as a covariate, minus the loglikelihood then fell from 63.44 to 62.10, giving a nonsignificant random effect (pvalue of 0.052). The standard error was 1.38, a little lower. It seems that once gender is included as a covariate, the random effect is no longer needed. This suggests that covariates that retard healing and were not available in the dataset, such as smoking and obesity, might also be irrelevant once blood flow measured as mean flux is known.
To test whether the gender effect varied with age, further terms were added to the model, such as an age by gender interaction term (equal to age for females, zero for males). A permutation test was carried out, by permuting this interaction variable among individuals and finding the proportion of such random permutations for which the loglikelihood was at least as large as the observed value, when this variable was included in the model. The pvalue was 0.105, showing that there is no evidence that the gender difference in healing times varies with age. Finally, we considered use of maximum and minimum logged flux, standard deviation of flux (coefficient of variation), and higher powers of the logged flux. None of these variables were statistically significant.
The model finally arrived at is very simple. Including gender besides mean flux as a predictor, table 5 shows the crosstabulation of predicted and observed healing times. As can be seen, the bootstrapcorrected percentage of correct classifications is 92%. This compares to the situation where gender is omitted from the model, when the correct classification rate is 91.3%, reducing to 90.9% when the bootstrap 'optimism' correction is applied. Figure 4 shows the predicted probabilities of the three healing periods as a function of mean flux.
Table 5. Classification performance of the PO model with gender and mean flux.
Figure 4. Model predictions for probabilities of slow, medium and fast healing for males and females under the PO model, as a function of mean laserDoppler flux.
Varying the model used from the proportionalodds model had very little effect. Using the probit link function instead of the logit decreased the loglikelihood function by 1.16, a slightly worse fit, and changing to the continuation ratio model decreased the loglikelihood by 0.207. The justification for taking the PO model with gender and mean flux as covariates as the definitive model was that it was the minimumAIC model (MAICE model).
The effect of total burned surface area (%TBSA)
Since the effects of total burned body area and skin pigmentation are likely to be important issues in the use of LDI, the statistical conclusions here will be spelt out. In the case of skin pigmentation, there are only 15 burn areas from blackskinned patients. It is not surprising that the coefficient for this factor is not statistically significant, and we really cannot draw any conclusion here from the data. There is however no reason to think that skin pigmentation would have an effect on healing for a given observed flux once skin is debrided.
The case of percentage total body surface area burned is slightly different. Here, apart from a few cases, the data do not extend much beyond 30% TBSA. In fact, the distribution of %TBSA fits well to a lognormal distribution, so that the natural logarithm of %TBSA has mean 1.956 and standard error 0.883.
However, this is a wide variation, and it allows the %TBSA coefficient in equation 1 to be determined with some accuracy; from table 4, the estimated coefficient is 0.0481, with a standard error of 0.051. Note that the sign of the effect is such that survival would improve with %TBSA, which would seem unlikely, so there is not a shred of evidence for any %TBSA effect.
A statistical task remains however, because the relevant question is slightly different. Given the possible (posterior) distribution of the %TBSA effect as determined from the data, would any clinically significant degradation in predictive performance result at high %TBSA from assuming a model with zero %TBSA effect?
To give some answer to this, a fortran95 simulation program was written. Fluxes were simulated by 'bootstrapping' the flux measurements (sampling with replacement), and a sample of 'patients' were generated with PU (perfusion unit) measurements from this distribution, and random %TBSA from the lognormal distribution of %TBSA. Healing times were simulated from the average flux model prediction for healing time, using flux and %TBSA as predictors. The 'data' healing times were generated randomly from the model with the coefficient of the %TBSA effect randomly chosen from its normal distribution. This can be done, as we know the estimated coefficient and its standard error. We thus imagine that the data will obey an unknown model that has a %TBSA coefficient consistent with our findings. The model used for prediction however has the %TBSA coefficient set to zero.
The misclassification rate obtained when the maximum likelihood estimate of healing time was used (i.e. we took the largest of the 3 probabilities) was 10.2% in this simulation. This agrees well with the 9% misclassification rate seen in the simple fit of the model to data. On moving to a scenario where the lognormal distribution of %TBSA is scaled up so that the mean %TBSA is 40% the misclassification rate increased to 13.6% and with 75% mean %TBSA, the misclassification rate increased to 17.7%.
This shows that, as expected, uncertainty in the amount of dependence of healing time (given a flux measurement) on %TBSA, does introduce extra error into healing time predictions as we extrapolate to high average %TBSA. Although we cannot reject the hypothesis that %TBSA has no effect on healing time once given a flux measurement, there could still be an effect large enough to degrade prediction ability when predictions are made ignoring this covariate. The misclassification rate could increase by a factor of 1.7. However, for healing time, this would only give a drop to maybe 85% accuracy. The conclusion is that, unless qualitatively different phenomena occur for patients with high %TBSA, the degree of accuracy of the model parameter estimates rules out any large increase in the rate of misclassification.
Goodness of Fit
This was explored in two ways: by fitting various models in the attempt to improve fit, and by diagnostic plots. The first approach led to the use of the logged mean flux as the best predictor in proportionalodds models.
Figures 5 and 6 show observed mean values for predictor variables for the three healing classes, and the predicted mean values from the PO model. Predicted mean values for a variable are calculated by including all cases, weighted by the predicted probability of group membership (e.g. [28], section 13.3.5). This relies on the use of Bayes' theorem, to turn the model into a predictive model for the covariates, given the healing time group. A bad fit is indicated by a strong divergence between the observed and predicted values, or an observed value that does not change monotonically with healing time. In figures 5 and 6, the agreement looks very good. The high value of R^{2 }and the high predictive ability confirm that the model fits the data well.
Figure 5. Observed and predicted proportions of females, under the proportional odds model.
Figure 6. Observed and predicted mean logged flux (PUs), under the proportional odds model.
Palette Construction
Recall that a false color image of the burn area is used by clinicians, as shown in figure 1. The 'palette' is the choice of PUcutpoints that defines whether the pixels of a burn region should be colored with color 1, 2, 3, 4 or 5, and does not specify which shades are actually used.
The 5color palette currently in use has PUcutpoints at 200, 260, 440 and 600 PU. The colors 1 to 5 are blue, green, yellow, pink and red, and in fact the blue region is divided into dark and light blue at 140 PU. The decision rule for healingtime prediction recommended to clinicians is that the 'buffer colors' green and pink should be lumped in with whichever of their neighbors has the larger area (rule 1).
Pixel data were available for two datasets, one used in palette construction, the other for palette validation. These were combined for the analysis here. The statistical analysis aimed to appraise the palette in current use, and to discover whether it could be optimised. Also, with an eye to the future, with individual pixel data available, we ask what is the best decision rule that can be constructed for 'automatic' use. Here, the device envisaged would show the probabilities of the three healing times for a region of interest. It is interesting to ask whether the mean flux that was used earlier as a healingtime predictor is a sufficient statistic, or whether, given the full distribution of flux density, a better predictive rule could be found. Using two flux cutpoints, the optimum was found by minimizing the mean percentage of wrong healingtime class assignments per healingtime group. The NAG function NelderMeade function minimizer E04CCF was used [32], with 50 sets of random starting values. It is also possible to convert this to a discrete optimization problem, by varying the two cutpoints in units of 1 PU. This gave the same results.
Appraising the palette
Table 6 shows how the current palette would perform, if the cutpoints were chosen in the center of each of the two 'buffer' colors. The decision rule based on the proportions of pixels below cut 1 (q_{1}), between the cuts (q_{2}) and above cut 2 (q_{3}) is
Table 6. Classification performance of: (A), the current palette, taking color boundaries in the center of 'buffer' colors at 230 and 520 PU.
When the five colors are used as designed, the classification of a case depends on the rule used for interpreting buffer colors. We can formulate at least four possible decision rules. Rule 1, 'lump buffer colors in with their neighboring color of larger area' is what is recommended to clinicians. Rule 2 is to ignore the buffer colors, and to predict healingtime from the largest area of the three colors blue, yellow and red. Rule 3 is to use the buffer colors to make the best case for each healing time. The visual decision made by a clinician would then be for example 'is the total area of color consistent with slow healing (colors 1 and 2) greater than the total area consistent with mediumtime healing (colors 2, 3 and 4)?'. This is illustrated in the pseudocode:
Rule 4 is 'regard the buffer colors as identical to the neighboring color corresponding to faster healing time'. Dark and light blue then mean slow healing, green or yellow mean medium healing time, and pink or red mean fast healing. This means that the 6color palette is effiectively used as a 3color palette.
Table 7 shows the performance of the palette under these four rules. It is interesting to see that the accuracy of diagnosis is very similar under rules 1, 2, and 3. This is reassuring, as no doubt, regardless of advice, different clinicians will interpret the falsecolor image in their individual ways. However, adoption of rule 4 would give worse results. From table 6 the optimized palette with 3 colors gives only a tiny improvement in the average percentage of wronglycategorized cases. The cutoff between fast (< 14 days) healing and medium (14–21 days) healing has moved down slightly, so that fewer fasthealing cases are misdiagnosed as mediumhealing time cases, but on the other hand, more medium healingtime cases are misdiagnosed as fast healing. In general, it was found by varying the cutoffs that many 5color palettes could be devised with similar performance. The buffer zones can be wider or narrower. The optimized palette in table 7 performs only slightly better than the palette in current use, and some of this improvement may be due to the double use of the data for fitting and validation. In fact, using the leaveoneout method of crossvalidation, the 11.6% misclassification rate increased to 13.7%, eroding the improvement completely.
Table 7. Classification performance of the current palette (A) under allocation rules 1–2, (B), under rule 3, and (C) under rule 4, and (D) of the optimized palette under rule 1.
Clearly, other criteria than the percentage of wrong diagnoses would be needed to optimize the 5color palette. A possible approach would be to penalize certainty in wrong diagnosis, and uncertainty in right diagnosis. It is arguably best to be right, and to be sure that you are right, less good to be right but unsure, worse to be wrong and unsure, and worst of all to be wrong and sure that you are right. The probability q_{1 }+ q_{3 }+ q_{5 }is a measure of certainty, and q_{2 }+ q_{4 }a measure of uncertainty. A loss function such as
penalizes complete certainty in wrong diagnosis as twice as bad as a wrong diagnosis with complete uncertainty, and a right diagnosis with complete uncertainty is 1/10 as a bad as the best wrong diagnosis. Minimizing this loss function, the PU cutoffs were 191, 243, 326 and 699 PU. The penalizing of uncertainty in diagnosis would allow an 'automatic' method of palette construction. Unfortunately, this resembles the problem of utility maximisation in healthcare and elsewhere; we could devise optimum policies if we could specify our utility function. Here, one could perhaps survey burns surgeons to elicit their preferred loss function, and only then could a palette be computergenerated. The value of the work here is that it has demonstrated that a 5 or 6color palette cannot simply be derived automatically, in the same way as can a 3color palette, but that further judgements are required.
Automatic decisions
Disregarding the clinical complications, for example ensuring that a region of interest is sufficiently homogeneous in blood flow for the entire region to heal at the same rate, we seek the best automatic decision rule. Table 6 shows that a rule based on the mean flux slightly outperforms rule 1. This is an example of the ODAO (Optimal Discrimant Analysis for Ordinal responses) method described by [23] and [24]. They found that this method outperforms other discriminant methods, because there is no need to model the functional form of the relation between (here) mean flux and healing time; one simply chops the mean into three ranges. Using the median flux or the mean logged flux gave slightly worse performance, and several other methods that were tried did not perform as well. It does indeed seem that the mean flux in a homogeneous burn region is the best predictor of healing time, and that the shape of the flux distribution is irrelevant. Clinically, this can make sense because mean blood flow in an area is proportional to total blood flow.
Discussion
The high accuracy of LDI, over 90%, relative to the 70% accuracy achievable by unaided clinical judgement, and the fact that the only demographic variable influencing healing time once the laserDoppler measurement is known is gender, are the main results of the analysis; it is significant that age and %TBSA did not aid prediction of healing time. Also, it has been established that the mean flux over an area is a 'sufficient' flux statistic; it is impossible to obtain better predictive power by substituting the median flux, by adding standard deviation, or by using quantiles of the flux distribution.
The accuracy of LDI seen in this study is lower than the levels reported by clinicians ([5,6,8] etc) and could be because the current assessment is based on only LDI whereas the other studies would have included clinical judgement. Also, instead of the 3way classification of healing time used here, frequently only a 2part classification is use: treat surgically or not. Only 6.5% of cases are misclassified by the methodology used here under that scheme.
Women tend to heal more slowly than men even when the blood flow as measured by laserDoppler flux is corrected for. It would therefore be possible to obtain a slight improvement (1.1%) in classification accuracy by using knowledge of gender in addition to the laserDoppler measurements.
Turning to the %TBSA results, this is an unusual statistical situation, where despite a nonsignificant result, we seek to explore the possible degradation of predictive ability at high %TBSA resulting from our uncertainty over the effect of %TBSA. Clearly, we have been able to transform a statement of 'no evidence for an effect' into a statement about the possible degradation of performance arising from the uncertainty over the size of the effect. However, any really believable evidence that LDI can be used up to high values of %TBSA must come from further clinical trials.
With an eye to the future, use of a probabilistic model and its incorporation into the laserDoppler imaging software would enable clinicians to see when the healing time prediction was very reliable and likely to be correct, and when the diagnosis was less certain. The best model obtained of the probabilities p_{1}, p_{2}, p_{3 }of fast, medium, and slow healing times is, for completeness:
where y = exp(6), is mean flux in PU, G is gender (0 for male, 1 for female), α_{1 }= .0222, α_{2 }= 5.563, β_{f }= 9.185, β_{g }= 2.164. The factor of exp(6) was added for numerical reasons.
Uncertainty over healing times is conveyed visually by using a 5color palette instead of the three colors needed for decision making. The statistical investigation into the derivation of an optimum palette has shown that the widths of the two 'buffer color zones' are not fixed by the requirement that the percentage of wrong healingtime diagnoses be minimized. They could be fixed by penalizing 'certain' wrong diagnoses more than uncertain ones. To construct an 'optimum' palette, it would be necessary to survey burn surgeons to elicit the best form of loss function to be used. Further work may be carried out along these lines.
Conclusion
This study has shown how contemporary statistical techniques can be used to assess performance of a medical imager and to reveal which clinical factors are important with and without use of the device. Ease of image interpretation is essential for such devices. In this respect, it has been shown that statistical techniques alone are incapable of providing a solution and that clinician's input is required.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
The first author carried out the statistical analysis, and the other authors set up the clinical investigation, obtained the clinical results and provided clinical information.
Acknowledgements
References

HemingtonGorse S: A comparison of laser Doppler imaging with other measurement techniques to assess burn depth.
J Wound Care 2005, 14:151153. PubMed Abstract  Publisher Full Text

Devgan L, Bhat S, Aylward S, RJ S: Modalities for the assessment of burn wound depth.
J Burns Wounds 2006, 5:e2. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chatterjee J: A critical evaluation of the clinimetrics of laser Doppler as a method of burn assessment in clinical practice.
J Burn Care Res 2006, 27:123130. PubMed Abstract  Publisher Full Text

Monstrey S, Hoeksema H, Verbelen J, Pirayesh A, Blondeel P: Assessment of burn depth and burn wound healing potential.
Burns 2008, 34:761769. PubMed Abstract  Publisher Full Text

Pape S, Skouras C, Byrne P: An audit of the use of laserDoppler Imaging (LDI) in the assessment of burns of intermediate depth.
Burns 2001, 27:233239. PubMed Abstract  Publisher Full Text

Jeng JC, Bridgeman A, Shivnan L, Thornton P, Alam H, Clarke T, Jablonski K, Jordan M: Laser Doppler imaging determines need for excision and grafting in advance of clinical judgment a prospective blinded trial.
Burns 2003, 29:665670. PubMed Abstract  Publisher Full Text

Holland A, Martin H, Cass D: Laser Doppler imaging prediction of burn wound outcome in children.
Burns 2002, 28:1117. PubMed Abstract  Publisher Full Text

La Hei E, Holland A, Martin H: Laser Doppler imaging of paediatric burns burn wound outcome can be predicted independent of clinical examination.
Burns 2006, 32:550553. PubMed Abstract  Publisher Full Text

Stuart A, Ord K: Kendall's Advanced Theory of Statistics. Volume 2. 5th edition. London: Edward Arnold; 1991.

Yeong EK, Mann R, Goldberg M, Engrav L, Heimbach D: Improved accuracy of burn wound assessment using laser Doppler.
Journal of traumainjury infection and critical care 1996, 40:956961. Publisher Full Text

JP B, DN H: Colour Atlas of Burn Care. Harcourt Publishers; 2001.

Attia A, Reda A, Mandil A, Arafa M, Massoud N: Predictive models for mortality and length of hospital stay in an Egyptian burns centre.
East Mediterr Health J 2000, 6:105561. PubMed Abstract

Jeschke M, Barrow R, Mlcak R, Herndon D: Endogenous anabolic hormones and hypermetabolism effect of trauma and gender differences.
Ann Surg 2005, 241:759767; discussion 7678. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Barrow R, Herndon D: Incidence of mortality in boys and girls after severe thermal burns.
Surg Gynecol Obstet 1990, 170:295298. PubMed Abstract

O'Keefe G, Hunt J, Purdue G: An evaluation of risk factors for mortality after burn trauma and the identification of genderdependent differences in outcomes.
Journal of the American college of surgeons 2001, 192:153160. PubMed Abstract  Publisher Full Text

McGwin GJ, George R, Cross J, Rei3 D, Chaudry I, Rue Lr: Gender differences in mortality following burn injury.
Shock 2002, 18:311315. PubMed Abstract  Publisher Full Text

George R, McGwin GJ, Schwacha M, Metzger J, Cross J, Chaudry I, Rue Lr: The association between sex and mortality among burn patients as modified by age.
J Burn Care Rehabil 2005, 26:416421. PubMed Abstract  Publisher Full Text

Kerby J, McGwin GJ, George R, Cross J, Chaudry I, Rue Lr: Sex differences in mortality after burn injury results of analysis of the National Burn Repository of the American Burn Association.
J Burn Care Res 2006, 27:452456. PubMed Abstract  Publisher Full Text

Bender R, Grouven U: Ordinal logistic regression in medical research.
Journal of the Royal College of Physicians of London 1997, 31:546551. PubMed Abstract

Agresti A: Categorical Data Analysis. 2nd edition. Hoboken, NJ: Wiley, Inc; 2002.

McCullagh P: Regression models for ordinal data (with discussion).
Journal of the Royal Statistical Society, Series B 1980, 42:109142.

McCullagh P, Nelder J: Generalized Linear Models. 2nd edition. London: Chapman & Hall; 1989.

Coste J, Walter E, Wasserman D, Venot A: Optimal discriminant analysis for ordinal responses.
Statistics in Medicine 1997, 16:561569. PubMed Abstract  Publisher Full Text

Le Teuff G, Quantin C, Venot A, Walter E, Coste J: Improving model robustness with bootstrapping Application to optimal discriminant analysis for ordinal responses (ODAO).
Methods of Information in Medicine 2005, 44:704711. PubMed Abstract  Publisher Full Text

Greenland S: Alternative models for ordinal logistic regression.
Statistics in Medicine 1994, 13:16651677. PubMed Abstract  Publisher Full Text

Ananth C, Kleinbaum D: Regression models for ordinal responses: a review of methods and applications.
International Journal of Epidemiology 1997, 26:13231333. PubMed Abstract  Publisher Full Text

Armstrong B, Sloan M: Ordinal regression models for epidemiologic data.
American Journal of Epidemiology 1989, 129:191204. PubMed Abstract  Publisher Full Text

Harrell FJ: Regression Modeling Strategies. New York: Springer; 2001.

Campbell M, Donner A: Classification efficiency of multinomial logistic regression relative to ordinal logistic regression.
Journal of the American Statistical Association 1989, 84:587591. Publisher Full Text

Vaupel J, Manton K, Stallard E: The Impact of Heterogeneity in Individual Frailty on the Dynamics of Mortality.
Demography 1979, 16:439454. PubMed Abstract  Publisher Full Text

Hedeker D, Gibbons R: A randomeffects ordinal regression model for multilevel analysis.
Biometrics 1994, 50:933944. PubMed Abstract  Publisher Full Text

NAG Plc: Numerical Algorithms Group Fortran Library – Supporting Documentation. [http://www.nag.co.uk/numeric/FL/FLdocumentation.asp] webcite

Efron B, Tibshirani R: An introduction to the Bootstrap. New York: Chapman & Hall; 1993.

Steyerberg E, Harrell FJ, Borsboon G, Eijkemans M, Vergouwe Y, Habbema J: Internal validation of predictive models: efficiency of some procedures for logistic regression analysis.
Journal of Clinical Epidemiology 2001, 54:774781. PubMed Abstract  Publisher Full Text

Cox D, Snell E: The analysis of Binary Data. London: Chapman and Hall; 1989.

Nagelkerke N: A note on a general definition of the coefficient of determination.
Biometrika 1991, 78:691692. Publisher Full Text

Edgington E, Onghena P: Randomization Tests. 4th edition. New York: Chapman & Hall/CRC; 2007.

Potter DM: A permutation test for inference in logistic regression with small and moderatesized datasets.
Statistics in Medicine 2005, 24:693708. PubMed Abstract  Publisher Full Text
Prepublication history
The prepublication history for this paper can be accessed here: