A. DAVIES et al

Mapping between measurement scales in meta-analysis, with application to
measures of body mass index in children

Corresponding author is Annabel L Davies.

Canynge Hall, Clifton, Bristol, BS8 2PN

# Mapping between measurement scales in meta-analysis, with application to
measures of body mass index in children

Annabel L Davies  A E Ades  Julian PT Higgins  Bristol Medical School,
University of Bristol, Bristol, United Kingdom<EMAIL_ADDRESS>A. L. Davies  A. E. Ades  J. P. T. Higgins

(Date Month Year; Date Month Year; Date Month Year)

###### Abstract

[Abstract]Quantitative evidence synthesis methods aim to combine data from
multiple medical trials to infer relative effects of different interventions.
A challenge arises when trials report continuous outcomes on different
measurement scales. To include all evidence in one coherent analysis, we
require methods to ‘map’ the outcomes onto a single scale. This is
particularly challenging when trials report aggregate rather than individual
data. We are motivated by a meta-analysis of interventions to prevent obesity
in children. Trials report aggregate measurements of body mass index (BMI)
either expressed as raw values or standardised for age and sex. We develop
three methods for mapping between aggregate BMI data using known relationships
between individual measurements on different scales. The first is an
analytical method based on the mathematical definitions of z-scores and
percentiles. The other two approaches involve sampling individual participant
data on which to perform the conversions. One method is a straightforward
sampling routine, while the other involves optimization with respect to the
reported outcomes. In contrast to the analytical approach, these methods also
have wider applicability for mapping between any pair of measurement scales
with known or estimable individual-level relationships. We verify and contrast
our methods using trials from our data set which report outcomes on multiple
scales. We find that all methods recreate mean values with reasonable
accuracy, but for standard deviations, optimization outperforms the other
methods. However, the optimization method is more likely to underestimate
standard deviations and is vulnerable to non-convergence.

, , and (2024), Mapping between aggregate level BMI data on different
measurement scales in a meta-analysis of childhood obesity prevention
interventions, 2024;00:1–21.

###### keywords:

Mapping, body mass index, standardization, meta-analysis, sampling,
optimization

††articletype: Research Article††footnotetext: Abbreviations: BMI, body mass
index; zBMI, body mass index z-score; SD, standard deviation; RMSE, root mean
squared error; MAE, mean absolute error.

##  1 Introduction

Evidence synthesis techniques, such as pairwise and network meta-analysis,
combine data from multiple medical trials in order to compare different
treatment options for the same condition.1, 2 The results of these analyses
are extremely influential, and are commonly used to inform national clinical
guidance regarding treatment recommendations.3 In the simplest scenario, each
trial provides an estimate of the effect of a specific treatment on the same
clinical outcome and the meta-analysis averages these values across the
trials. A challenge facing meta-analysis occurs when trials report continuous
outcomes on different measurement scales.4 That is, all trials measure the
same underlying construct but use different tools or test instruments. For
example, if the outcome of interest is temperature, trialists may choose to
measure this in Celsius, Fahrenheit or even Kelvin. In order to combine
evidence from all eligible trials, we require methods to ‘map’ the outcomes
onto a common scale.5 Unlike measurements of temperature, the relationships
between different clinical scales may be complex or even unknown. Mapping
between these scales is particularly challenging in the common scenario where
trial data is only available at the aggregate level (e.g. as means and
standard deviations) rather than per participant.

This paper is motivated by a systematic review of interventions to prevent
childhood obesity.6, 7 Body mass index (BMI), which measures a person’s weight
relative to their height, is one of the most widely used tools for assessing
obesity.8 However, in children, BMI varies according to the child’s age and
sex. Therefore, obesity in children is more appropriately assessed using age-
and sex- standardized BMI.9 A common standardization approach is to transform
BMI to a z-score (zBMI) based on the distribution of BMI in an appropriate
reference population.10 Using normality assumptions, the z-score can also be
expressed as a percentile.

In our motivating data set, trials report aggregate-level outcomes in terms of
either unstandardized BMI, zBMI or percentile BMI. Transformations between
these scales are well-defined for individual participant data, but cannot be
used for mapping between aggregate measurements. For example, the conversion
between BMI and zBMI relies on information about an individual’s age and sex,
as well as about the wider population, while percentile and zBMI are related
via a cumulative distribution function. Due to the lack of appropriate
aggregate mapping methods, previous analyses of this, and similar data have
been restricted to separate syntheses on each scale.6, 7, 11, 12 In this
scenario only a subset of the total data informs each analysis, limiting the
statistical power and precision of the results. A more desirable approach is
to include all the trials in a single analysis by expressing the outcomes on a
common scale. We choose zBMI as our primary outcome of interest, and therefore
aim to map aggregate measurements of BMI and percentile onto this scale.

In meta-analysis, a common approach to synthesising outcomes on different
scales is to measure treatment effects as ‘standardized mean differences’,
where each difference in means is divided by some scale-specific standard
deviation (SD). However, this method relies on the assumption that the
different scales are linear transformations of each other,13 which is not true
for the BMI outcomes. Another recommended method is the ‘ratio of means’
approach.14 Here, effect estimates are calculated by taking the ratio between
the mean outcomes in each arm, thus transforming them onto a unitless scale.
This assumes a multiplicative rather than an additive model and requires that
measurements on all scales take values greater than zero.4 Since zBMI is
defined to take positive and negative values, this approach is also not
suitable for our data. More recently, Lu et al (2014)5 developed a more
sophisticated approach that simultaneously synthesizes treatment effects from
multiple outcomes while estimating mapping coefficients between the scales.
However, this method assumes individual-level transformations are unknown and
makes strong assumptions about the form of aggregate relationships between the
scales. Based on our knowledge of the different BMI scales, we seek a method
that, rather than enforcing simplifying assumptions, exploits fully the known
individual relationships.

In this paper we develop three methods for mapping between aggregate-level BMI
data. The first two approaches involve making distributional assumptions for
one of the outcome scales and sampling individual participant data on which to
perform the conversions. One method performs the mapping based on sampling
alone, while the other makes use of a simple optimization routine. Although we
develop these methods in the context of BMI outcome scales, they also have
wider applicability to mapping between any aggregate outcomes where
individual-level relationships are available. In addition to these methods, we
also develop an analytic approach to mapping between zBMI and percentile based
on their mathematical definitions.

We begin by describing our motivating data set in Section 2. In Section 3 we
define the three BMI measurement scales and describe how to convert between
them for individual-level data. Next, we introduce our three aggregate mapping
methods in Section 4. We begin with the analytical approach for mapping
between percentile and zBMI. Then we define the more general sampling and
optimization methods, which can be used to map from either BMI or percentile.
We provide details of the specific algorithm in each case. In Section 5 we
compare the accuracy of the different methods using observations from trials
in our data set that report measurements on more than one scale. Finally, we
summarize and discuss our results in Section 6.

##  2 Motivating data set

Our work is motivated by a systematic review of nearly 250 randomized
controlled trials of interventions to prevent obesity in children aged 5 to 18
years.6, 7 We will refer to this as the obesity data set. 204 trials provide
usable data on BMI, zBMI or percentile BMI. The outcome of interest on each
scale is the mean difference (MD) between the intervention arm and reference
arm in change from baseline. This is known as contrast-level data because it
represents the comparison between two arms of a trial. For brevity, we will
refer these data points as ‘contrasts’. Arm-level data, on the other hand,
refer to the means and standard deviations of measurements (BMI, zBMI or
percentile) in each arm at each time point. We will perform our mapping
procedure on the arm-level data. The results can then be used to reconstruct
the relevant contrasts for analysis.

Figure 1: CDC (Centers for Disease Control and Prevention) reference
parameters $\\{\lambda,\mu,\sigma\\}$ as a function of age. Parameter values
obtained from Reference15.

The obesity data set includes multi-arm trials and observations at multiple
follow-up times, such that each trial can be associated with more than one
data point. Additionally, some trials report the same contrast on more than
one outcome scale. Overall, there are 182 contrasts measured on the BMI scale,
171 on the zBMI scale and 48 on the percentile scale. Due to some trials
reporting multiple scales, a fully mapped data set would include 286 observed
contrasts. An appropriate mapping procedure would therefore yield a much
larger sample compared with analysing any one of the outcomes alone.

To test the methods we develop, we will use trials that report arm-level data
on more than one outcome scale. From these trials we have 171 arm-level means
measured on both the BMI and zBMI scale, 20 means measured as percentile and
zBMI, and 34 values measured on all three scales. We will use these data to
investigate how well the mapping procedures recreate the observed zBMI means
from the reported BMI and percentile values. In these tests we will treat each
data point independently, though we note that this is not strictly true since
an individual trial is associated with multiple arm-level means.

##  3 BMI outcomes: individual-level definitions

We define $B$, $Z$, and $P$ as BMI, zBMI and percentile respectively. For each
arm and time point of a particular trial, we write
$\mathcal{D}_{X}=(\bar{X},\text{SD}(X))$ for the aggregate-level data on scale
$X\in\\{B,Z,P\\}$, where $\bar{X}$ is the reported mean and $\text{SD}(X)$ is
the standard deviation on that scale. In the following, we define the
relationships between these scales for a particular individual.

###  3.1 BMI and zBMI

The body mass index (BMI) of an individual is defined as

$\displaystyle B=\frac{M}{H^{2}},$ (1)

where $M$ is the person’s mass in kilograms (kg) and $H$ is their height in
metres (m). In adults, obesity is defined when an individual’s BMI crosses a
certain threshold. For example, the World Health Organization (WHO) define
obesity to be a BMI greater than or equal to $30\text{ kg}/\text{m}^{2}$.16 In
children, however, BMI varies with age and sex making it more difficult to
define meaningful thresholds. Instead, weight classifications in children are
defined by comparing a child’s BMI with the distribution of BMI in other
children of the same age and sex. This comparison is done by expressing BMI as
a z-score.10

For a normally distributed random variable, the z-score describes how many
standard deviations an individual measurement is away from the mean of the
distribution. It is calculated by subtracting the mean from the individual
measurement and dividing by the standard deviation. This defines the z-score
as a standard normal variable. However, measurements of BMI in the general
population exhibit a skewed distribution meaning z-scores defined in this way
do not follow the standard normal distribution.17 To deal with this problem,
zBMI is usually calculated using the so-called ‘LMS method’, an approach
developed by Cole et al.18, 19 as a way of defining z-scores for skewed
variables using age-related reference standards.

The LMS method characterises the distribution of the variable of interest in
terms of three parameters, $\\{\lambda,\mu,\sigma\\}$, where $\mu$ is the
median of the variable’s distribution, $\sigma$ is the coefficient of
variation, and $\lambda$ is the Box-Cox power. The method assumes that raising
the variable to the power $\lambda$ (or taking the logarithm of the variable
if $\lambda=0$) yields a normal distribution. For a BMI, $B$, the LMS method
defines the Z-score, $Z$, as a standard normal variable calculated via

$\displaystyle
Z=\begin{cases}\frac{1}{\lambda\sigma}((\frac{B}{\mu})^{\lambda}-1)&\lambda\neq
0\\\ \frac{1}{\sigma}\log\frac{B}{\mu}&\lambda=0.\end{cases}$ (2)

The distribution of BMI in children depends on age, sex, and the population of
interest.9 Therefore, to calculate zBMI for an individual, we require the age-
and sex-specific values of $\\{\lambda,\mu,\sigma\\}$ for an appropriate
reference population. These values have been calculated from large-scale
surveys of children’s heights and weights in various populations20, 21, 22,
23, 24 and are available from national and international public health
organization websites. For example, Figure 1 shows the LMS parameters for boys
and girls as a function of age for the population of the United States,
available through the Centers for Disease Control and Prevention (CDC).15
Here, we see that $\lambda<0$ for all children between the ages of 2 and 18
years. Reference charts based on surveys of international populations are
available via the World Health Organization (WHO)25 and the International
Obesity Task Force (IOTF).26 As shown in Section A of the Supplementary
Material, both these international references also indicate $\lambda<0$ across
all ages. Therefore, Equation (2) reduces to the first line. The procedure for
calculating a child’s zBMI is then as follows: (i) record the child’s BMI,
age, and sex, (ii) choose an appropriate LMS reference chart for the relevant
population, (iii) read off the values of $\\{\lambda,\mu,\sigma\\}$ for the
age and sex of the child, (iv) substitute the parameters and the measured BMI
into the first line of Equation (2).

###  3.2 zBMI and percentile

A percentile describes the probability that a random variable is less than a
particular value. It is related to the z-score via the standard cumulative
distribution function (cdf), $\Phi(.)$,

$\displaystyle
P(z)=\text{Prob}(Z<z)=\Phi(z)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{z}\text{e}^{-x^{2}/2}\text{d}x.$
(3)

Therefore, while there is a one-to-one mapping between zBMI and percentile
BMI, this relationship is non-linear.

##  4 Mapping between aggregate-level BMI outcomes

If we have measurements on one of the BMI scales for individuals of known age
and sex, we can use Equations (2) and (3) to convert the values to either of
the other two scales. However, in the obesity data set we do not have access
to individual participant data, only the mean and standard deviation of
measurements within each arm. Due to the complex non-linear relationships
between the scales, we cannot simply apply Equations (2) and (3) to these
aggregate data. Instead, in this section we develop more sophisticated methods
that build on the individual-level transformations. These allow us to map
aggregate BMI and percentile data onto the zBMI scale. We label these mappings
$\mathcal{D}_{B}\xrightarrow{}{}\mathcal{D}_{Z}$ and
$\mathcal{D}_{P}\xrightarrow{}{}\mathcal{D}_{Z}$ respectively.

###  4.1 Method 1: Analytical

####  4.1.1 Percentile to zBMI

In this section we develop an analytical method to perform the aggregate-level
mapping from percentile to zBMI,
$\mathcal{D}_{P}\xrightarrow{}{}\mathcal{D}_{Z}$. The method uses expressions
for the expectation and variance of continuous random variables along with
standard integral results from Owen (1980),27

$\displaystyle\int_{-\infty}^{\infty}\Phi(a+bx)\phi(x)\text{d}x$
$\displaystyle=\Phi\left(\frac{a}{\sqrt{1+b^{2}}}\right),$ (4)
$\displaystyle\int_{-\infty}^{\infty}[\Phi(a+bx)]^{2}\phi(x)\text{d}x$
$\displaystyle=\Phi\left(\frac{a}{\sqrt{1+b^{2}}}\right)$ $\displaystyle\hskip
10.0pt-2T\left(\frac{a}{\sqrt{1+b^{2}}},\frac{1}{\sqrt{1+2b^{2}}}\right),$ (5)

where $\phi(.)$ is the probability density function (pdf) of the standard
normal distribution and $T(h,a)$ is Owen’s T function.28

In Section 3.1 we showed that zBMI is defined such that it follows a standard
normal distribution in the general population. We assume it also follows a
normal distribution for samples from the general population such as
participants in randomized trials, $Z\sim\mathcal{N}(m_{z},s_{z}^{2})$, with
some sample-specific mean $m_{z}$ and standard deviation $s_{z}$.†††Throughout
this article, we will use $m$ and $s$ to represent the mean and standard
deviation of distributions to avoid confusion with the $\mu$ and $\sigma$
parameters in the LMS method. Given the relationship between percentile BMI
and zBMI in Equation (3), the expected value of percentile in this sample is
then

$\displaystyle\mathbb{E}_{Z}(P)=\mathbb{E}_{Z}(\Phi(Z))=\int_{-\infty}^{\infty}f_{Z}(x)\Phi(x)\text{d}x,$
(6)

where $f_{Z}(x)$ is the pdf of $Z$. Using the expression for the pdf of a
normal distribution, we find

$\displaystyle\mathbb{E}_{Z}(P)=\int_{-\infty}^{\infty}\frac{1}{s_{z}}\phi\left(\frac{x-m_{z}}{s_{z}}\right)\Phi(x)\text{d}x.$
(7)

To evaluate this integral, we use the substitution
$x^{\prime}=\frac{x-m_{z}}{s_{z}}$ so that Equation (7) takes the form of the
standard integral in Equation (4) with $a=m_{z}$ and $b=s_{z}$. This gives the
solution

$\displaystyle\mathbb{E}_{Z}(P)=\Phi\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}}\right).$
(8)

Next, we derive an expression for the variance of $P$ as a function of $m_{z}$
and $s_{z}$. In terms of the expectation, the variance is defined as

$\displaystyle\text{Var}(P)=\mathbb{E}_{Z}(P^{2})-\mathbb{E}_{Z}(P)^{2}.$ (9)

Given that we already have an expression for $\mathbb{E}_{Z}(P)$, we must now
evaluate

$\displaystyle\mathbb{E}_{Z}(P^{2})=\mathbb{E}_{Z}([\Phi(Z)]^{2})=\int_{-\infty}^{\infty}\frac{1}{s_{z}}\phi\left(\frac{x-m_{z}}{s_{z}}\right)[\Phi(x)]^{2}\text{d}x.$
(10)

Using the same substitution as before, $x^{\prime}=\frac{x-m_{z}}{s_{z}}$,
Equation (10) takes the form of the second standard integral, Equation (5),
with $a=m_{z}$ and $b=s_{z}$. From this, we obtain the result

$\displaystyle\mathbb{E}_{Z}(P^{2})=\Phi\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}}\right)-2T\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}},\frac{1}{\sqrt{1+2s_{z}^{2}}}\right).$
(11)

Finally, substituting Equations (8) and (11) into Equation (9) gives the
following expression for the variance of $P$,

$\text{Var}(P)=\Phi\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}}\right)-2T\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}},\frac{1}{\sqrt{1+2s_{z}^{2}}}\right)\\\
-\left(\Phi\left(\frac{m_{z}}{\sqrt{1+s_{z}^{2}}}\right)\right)^{2}.$ (12)

To perform the mapping $\mathcal{D}_{P}\xrightarrow{}{}\mathcal{D}_{Z}$, we
set the expectation of percentile in Equation (8) equal to the reported sample
mean, $\mathbb{E}_{Z}(P)=\bar{P}$, and the variance in Equation (12) equal to
the squared standard deviation, $\text{Var}(P)=\text{SD}(P)^{2}$. Then, we
obtain the mapped dataset, $\mathcal{D}_{Z}$, by solving the resulting
simultaneous equations for $m_{z}$ and $s_{z}$. This must be done using
numerical methods, meaning the method is not strictly analytical. However, we
retain the name ‘analytical method’ to distinguish it from the subsequent
methods based on sampling and optimization.

###  4.2 Method 2: Sampling

The relationship between BMI and zBMI in Equation (2) depends implicitly on
the individual’s age and sex via dependence on the LMS parameters. This
prevents us from deriving an analytic method to perform this mapping on the
aggregate level. Instead, we propose a sampling method which can be used to
map from both percentile and BMI.

The general idea, when mapping from some scale $X$ to another scale $Y$, is as
follows: (i) assume a sampling distribution for scale $X$, (ii) specify the
parameters of the distribution using the observed aggregate data on that
scale, $\mathcal{D}_{X}$, (iii) sample individual participant data from this
distribution, (iv) use the appropriate individual-level mapping relationship
to convert the sampled data on scale $X$ to scale $Y$, (v) calculate the mean
and SD of the individual-level data on scale $Y$ to define $\mathcal{D}_{Y}$.

In Sections 4.2.1 and 4.2.2 we outline the specific sampling routines to
perform the mappings $\mathcal{D}_{P}\xrightarrow{}{}\mathcal{D}_{Z}$ and
$\mathcal{D}_{B}\xrightarrow{}{}\mathcal{D}_{Z}$ respectively.

####  4.2.1 Algorithm 1: Percentile to zBMI via sampling

To map from percentile to zBMI using the sampling method, we must first assume
a distribution for $P$. As shown in Equation (3), percentile BMI is defined as
the standard normal cdf of $Z$. Given zBMI is a standard normal variable in
the general population, the corresponding distribution of percentile is
uniform between 0 and 1. However, to describe percentile BMI in a sample
population (e.g. in a particular trial), we require a more flexible
distribution. We choose the beta distribution, $\text{Beta}(\alpha,\beta)$,
because it has the correct $[0,1]$ support and contains the uniform
distribution as a special case when $\alpha=\beta=1$. Furthermore, as
discussed in Section B of the Supplementary Material, a beta distribution for
$P$ produces a normal or approximately normal distribution for $Z$ under
certain parameterisations. This choice therefore aligns with our previous
assumptions in Section 4.1.1 about the distribution of zBMI in a sample
population.

The mean and variance of a beta random variable,
$X\sim\text{Beta}(\alpha,\beta)$, with parameters $\alpha$ and $\beta$ are
given by,

$\displaystyle\mathbb{E}(X)$ $\displaystyle=\frac{\alpha}{\alpha+\beta},$ (13)
$\displaystyle\text{Var}(X)$
$\displaystyle=\frac{\alpha\beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}.$ (14)

To implement the mapping $\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$ via
sampling, we obtain the $\alpha$ and $\beta$ parameters of the percentile
sampling distribution by equating the observed mean and variance of the
aggregate data with these expectations. The method then proceeds as follows:

  1. 1.

For $i=1,\ldots,N$:

    1. (i)

Sample percentile BMI from $P_{i}\sim\text{Beta}(\alpha,\beta)$ where

$\displaystyle\alpha$
$\displaystyle=\bar{P}^{2}\left(\frac{1-\bar{P}}{\text{SD}(P)^{2}}-\frac{1}{\bar{P}}\right),$
(15) $\displaystyle\beta$
$\displaystyle=\alpha\left(\frac{1}{\bar{P}}-1\right).$ (16)

    2. (ii)

Calculate the corresponding value of zBMI by rearranging Equation (3),

$\displaystyle Z_{i}=\Phi^{-1}(P_{i}).$ (17)

  2. 2.

Calculate the mean and standard deviation of the $\\{Z_{i}\\}$ to define the
mapped dataset $\mathcal{D}_{Z}=(\bar{Z},\text{SD}(Z))$.

####  4.2.2 Algorithm 2: BMI to zBMI via sampling

Silverman and Lipscombe (2022) show that BMI in the adult population follows a
lognormal distribution. Given that $\lambda\neq 0$ in any BMI reference chart,
this is not strictly true for children. However, the lognormal distribution is
a practical choice as it is easy to sample from, provides flexibility and has
the property of skew. Therefore, to set-up the sampling method for
$\mathcal{D}_{B}\xrightarrow{}\mathcal{D}_{Z}$, we set the observed sample
mean and SD on the BMI scale, $\mathcal{D}_{B}=(\bar{B},\text{SD}(B))$, equal
to the mean and (square root of) variance of a lognormal variable.

To obtain zBMI from BMI for an individual sample $i$, we require the
appropriate LMS parameters, $\\{\lambda,\mu,\sigma\\}$. For a given reference
chart, these parameters are functions of the age, $a_{i}$, and sex, $x_{i}$,
of the individual. We assume trials provide aggregate information about the
age and sex of its participants. Therefore, at each iteration we sample sex
(male = 1, female = 0) from a Bernoulli distribution with probability equal to
the proportion of males in the trial, $p_{m}$. For age, we investigate both
normal and uniform distributions, specified in terms of the reported mean and
SD of age, $\mathcal{D}_{a}=(\bar{a},\text{SD}(a))$. The sampling method for
BMI is then as follows:

  1. 1.

For $i=1,\ldots,N$:

    1. (i)

Sample BMI from a lognormal distribution,
$B_{i}\sim\text{LN}(m_{B},s_{B}^{2})$, where

$\displaystyle m_{B}$
$\displaystyle=\log\left(\frac{\bar{B}}{\sqrt{\bar{B}^{2}+\text{SD}(B)^{2}}}\right),$
(18) $\displaystyle s_{B}^{2}$
$\displaystyle=\log\left(1+\frac{\text{SD}(B)^{2}}{\bar{B}^{2}}\right).$ (19)

    2. (ii)

Sample sex from a Bernoulli distribution, $x_{i}\sim\text{Bern}(p_{m})$.

    3. (iii)

Sample age from one of the following distributions:

      1. (a)

Normal: $a_{i}\sim\mathcal{N}(\bar{a},\text{SD}(a))$,

      2. (b)

Uniform: $a_{i}\sim\text{Unif}(\min(a),\max(a)),$ where
$\min(a)=\bar{a}-2\text{SD}(a)$ and $\max(a)=\bar{a}+2\text{SD}(a)$.

    4. (iv)

Read off the values of the LMS parameters at age $a_{i}$ and sex $x_{i}$ from
an appropriate reference chart,
$\lambda(a_{i},x_{i}),\mu(a_{i},x_{i}),\sigma(a_{i},x_{i})$.

    5. (v)

Calculate the corresponding value of zBMI via Equation (2) for $\lambda\neq
0$,

$\displaystyle
Z_{i}=\frac{1}{\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})}\left(\left(\frac{B_{i}}{\mu(a_{i},x_{i})}\right)^{\lambda(a_{i},x_{i})}-1\right).$
(20)

  2. 2.

Calculate the mean and standard deviation of the $\\{Z_{i}\\}$ to define the
mapped dataset $\mathcal{D}_{Z}=(\bar{Z},\text{SD}(Z))$.

###  4.3 Method 3: Optimization

To map from scale $X$ to scale $Y$, the sampling method assumes, and samples
from, a distribution for $X$. However, in some scenarios, it may be easier to
pick a distribution and sample from the $Y$ scale. Therefore, we develop a
second mapping method that uses an optimization scheme to map from $X$ to $Y$
based on sampling $Y$. In particular, we use optimization to perform the
mappings $\mathcal{D}_{B}\xrightarrow{}\mathcal{D}_{Z}$ and
$\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$ by sampling individual
participant data on the zBMI scale.

The general idea of the optimization method, for mapping from scale $X$ to
scale $Y$, is as follows: (i) Assume a distribution for scale $Y$, (ii)
initialise the parameters of this distribution to some sensible values, (iii)
sample individual data from this distribution, (iv) use the appropriate
individual-level mapping relationship to convert sampled data on scale $Y$ to
scale $X$, (v) calculate the mean and SD of the sampled data on scale $X$,
(vi) compare these values to the observed data, $\mathcal{D}_{X}$, and adjust
the parameters of the $Y$ distribution accordingly, (vii) repeat steps
(iii)-(vi) until the sampled means and SDs on scale $X$ are within some
specified threshold of the observed values, (viii) calculate the mean and SD
of the sampled $Y$ values at the last iteration and use these to define the
mapped dataset $\mathcal{D}_{Y}$.

For the obesity data set we are mapping to the zBMI scale and therefore, we
need to assume a sampling distribution for $Z$. As in Section 4.1.1 we model
zBMI in a sample population as following a normal distribution. We initialise
the mean and standard deviation to 0 and 1, representing the distribution of
zBMI in the general population.

To implement the method, we use a naïve optimization scheme. At iteration $t$,
we increase (or decrease) a parameter of the $Y$ distribution by a fixed
amount if the difference between the estimated parameter on the $X$ scale and
its corresponding observed value in $\mathcal{D}_{X}$ exceeds some specified
threshold. In particular, we adjust the mean and SD of the zBMI distribution
based on the corresponding sampled mean and SD of BMI or percentile at that
iteration. However, rather than resampling the zBMI values for each updated
set of parameters (step (iii) above), we use the properties of the normal
distribution to update the samples “smoothly”. First, assume that at iteration
$t$ we sample $N$ values of zBMI, $\\{Z_{i}^{(t)};i=1,\ldots,N\\}$, from a
normal distribution with mean $m_{z}^{(t)}$ and SD $s_{z}^{(t)}$. Then, as
shown in Appendix A, if the mean and SD are increased by $\delta_{m}$ and
$\delta_{s}$ respectively, the $i^{th}$ sample of zBMI can be updated via,

$\displaystyle
Z_{i}^{(t+1)}=\left(1+\frac{\delta_{s}}{s_{z}^{(t)}}\right)(Z_{i}^{(t)}-m_{z}^{(t)})+(m_{z}^{(t)}+\delta_{m}).$
(21)

Using this updating scheme means we only need to generate one set of $N$ zBMI
samples from the initialized distribution, thus speeding up the optimization.

####  4.3.1 Algorithm 3: Percentile to zBMI via optimization

We now specify the optimization method for the mapping
$\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$. First, we choose an
incremental step, $\Delta_{\text{step}}>0$, and a tolerance threshold,
$\Delta_{\text{tol}}>0$. The former defines the amount by which we increase or
decrease the parameters at each iteration, while the latter dictates the
accuracy required before successful convergence is reached. We discuss choices
for these values in Section C of the Supplement. The method then proceeds as
follows:

  1. 1.

Initialize the parameters: $t=0$, $m_{z}^{(t)}=0$, $s_{z}^{(t)}=1$.

  2. 2.

Generate the initial sample: For $i=1,\ldots,N$

    1. (i)

Sample zBMI from $Z_{i}^{(t)}\sim\mathcal{N}(m_{z}^{(t)},s_{z}^{(t)2})$.

    2. (ii)

Calculate the equivalent percentile, $P_{i}^{(t)}=\Phi(Z_{i}^{(t)})$.

  3. 3.

Calculate the means and standard deviations of the $\\{Z_{i}^{(t)}\\}$ and
$\\{P_{i}^{(t)}\\}$ samples. For scale $X\in\\{Z,P\\}$ we label these values
as $\mathcal{D}_{X}^{(t)}=(\bar{X}^{(t)},\text{SD}(X)^{(t)})$.

  4. 4.

Define $\delta_{m}$ and $\delta_{s}$ by comparing $\mathcal{D}_{P}^{(t)}$ with
the observed data, $\mathcal{D}_{P}$:

$\displaystyle\delta_{m}$
$\displaystyle=\begin{cases}\phantom{-}0&\text{if}\hskip
10.0pt|\bar{P}^{(t)}-\bar{P}|\leq\Delta_{\text{tol}}\\\
\phantom{-}\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\bar{P}^{(t)}-\bar{P}\phantom{|}<\Delta_{\text{tol}}\\\
-\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\bar{P}^{(t)}-\bar{P}\phantom{|}>\Delta_{\text{tol}}\end{cases}$
$\displaystyle\delta_{s}$
$\displaystyle=\begin{cases}\phantom{-}0&\text{if}\hskip
10.0pt|\text{SD}(P)^{(t)}-\text{SD}(P)|\leq\Delta_{\text{tol}}\\\
\phantom{-}\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\text{SD}(P)^{(t)}-\text{SD}(P)\phantom{|}<\Delta_{\text{tol}}\\\
-\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\text{SD}(P)^{(t)}-\text{SD}(P)\phantom{|}>\Delta_{\text{tol}}\end{cases}$

  5. 5.

Update the samples: For $i=1,\ldots,N$

    1. (i)

$Z_{i}^{(t+1)}=\left(1+\delta_{s}/s_{z}^{(t)}\right)(Z_{i}^{(t)}-m_{z}^{(t)})+(m_{z}^{(t)}+\delta_{m})$

    2. (ii)

$P_{i}^{(t+1)}=\Phi(Z_{i}^{(t+1)})$

  6. 6.

Update the parameters of the zBMI distribution:
$m_{z}^{(t+1)}=m_{z}^{(t)}+\delta_{m}$ and
$s_{z}^{(t+1)}=s_{z}^{(t)}+\delta_{s}$.

  7. 7.

Set $t=t+1$. Go to Step 3.

The method repeats steps 3 to 7 until one of the following conditions is met:

  1. C1.

$|\bar{P}^{(t)}-\bar{P}|\leq\Delta_{\text{tol}}$ and
$|\text{SD}(P)^{(t)}-\text{SD}(P)|\leq\Delta_{\text{tol}}$

  2. C2.

The number of iterations exceeds some specified maximum, $n_{\text{max}}$.

We then define the mapped zBMI data set using the mean and SD of the
$\\{Z_{i}^{(t)}\\}$ samples at the final iteration,
$\mathcal{D}_{Z}=\mathcal{D}_{Z}^{(t)}$. Alternatively, we could estimate
$\mathcal{D}_{Z}$ using the final values of the parameters $m_{z}^{(t)}$ and
$s_{z}^{(t)}$. We refer to the former as the sample estimates and the latter
as the distribution estimates. In Section 5.2.1 we investigate the behaviour
of both. For the mapping $\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$, the
sample and distribution estimates should be almost equivalent. However, as we
will see below, this is not necessarily the case for the BMI mapping due to
the extra approximations involved.

####  4.3.2 Algorithm 4: BMI to zBMI via optimization

The optimization method for the mapping
$\mathcal{D}_{B}\xrightarrow{}\mathcal{D}_{Z}$ follows the same steps as for
percentile, but with the additional sampling of age and sex in order to
specify the LMS parameters:

  1. 1.

Initialize the parameters: $t=0$, $m_{z}^{(t)}=0$, $s_{z}^{(t)}=1$.

  2. 2.

Generate the initial sample: For $i=1,\ldots,N$

    1. (i)

Sample zBMI from $Z_{i}^{(t)}\sim\mathcal{N}(m_{z}^{(t)},s_{z}^{(t)2})$.

    2. (ii)

Sample sex from a Bernoulli distribution, $x_{i}\sim\text{Bern}(p_{m})$.

    3. (iii)

Sample age from one of the following distributions:

      1. (a)

Normal: $a_{i}\sim\mathcal{N}(\bar{a},\text{SD}(a))$

      2. (b)

Uniform: $a_{i}\sim\text{Unif}(\min(a),\max(a)),$ where
$\min(a)=\bar{a}-2\text{SD}(a)$ and $\max(a)=\bar{a}+2\text{SD}(a)$.

    4. (iv)

Read off the values of the LMS parameters at age $a_{i}$ and sex $x_{i}$ from
an appropriate reference chart,
$\lambda(a_{i},x_{i}),\mu(a_{i},x_{i}),\sigma(a_{i},x_{i})$.

    5. (v)

Calculate BMI by rearranging Equation (2) for $\lambda\neq 0$,

$\displaystyle
B_{i}=\mu(a_{i},x_{i})(1+\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})Z_{i})^{\frac{1}{\lambda(a_{i},x_{i})}}.$
(22)

  3. 3.

Calculate the means and standard deviations of the $\\{Z_{i}^{(t)}\\}$ and
$\\{B_{i}^{(t)}\\}$ samples, labelled
$\mathcal{D}_{X}^{(t)}=(\bar{X}^{(t)},\text{SD}(X)^{(t)})$ for
$X\in\\{Z,B\\}$.

  4. 4.

Define $\delta_{m}$ and $\delta_{s}$ by comparing $\mathcal{D}_{B}^{(t)}$ to
the observed data, $\mathcal{D}_{B}$:

$\displaystyle\delta_{m}$
$\displaystyle=\begin{cases}\phantom{-}0&\text{if}\hskip
10.0pt|\bar{B}^{(t)}-\bar{B}|\leq\Delta_{\text{tol}}\\\
\phantom{-}\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\bar{B}^{(t)}-\bar{B}\phantom{|}<\Delta_{\text{tol}}\\\
-\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\bar{B}^{(t)}-\bar{B}\phantom{|}>\Delta_{\text{tol}}\end{cases},$
$\displaystyle\delta_{s}$
$\displaystyle=\begin{cases}\phantom{-}0&\text{if}\hskip
10.0pt|\text{SD}(B)^{(t)}-\text{SD}(B)|\leq\Delta_{\text{tol}}\\\
\phantom{-}\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\text{SD}(B)^{(t)}-\text{SD}(B)\phantom{|}<\Delta_{\text{tol}}\\\
-\Delta_{\text{step}}&\text{if}\hskip
10.0pt\phantom{|}\text{SD}(B)^{(t)}-\text{SD}(B)\phantom{|}>\Delta_{\text{tol}}\end{cases}$

  5. 5.

Update the samples: For $i=1,\ldots,N$

    1. (i)

$Z_{i}^{(t+1)}=\left(1+\delta_{s}/s_{z}^{(t)}\right)(Z_{i}^{(t)}-m_{z}^{(t)})+(m_{z}^{(t)}+\delta_{m})$

    2. (ii)

$B_{i}=\mu(a_{i},x_{i})(1+\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})Z_{i})^{\frac{1}{\lambda(a_{i},x_{i})}}$

  6. 6.

Update the parameters of the zBMI distribution:
$m_{z}^{(t+1)}=m_{z}^{(t)}+\delta_{m}$ and
$s_{z}^{(t+1)}=s_{z}^{(t)}+\delta_{s}$.

  7. 7.

Set $t=t+1$. Go to Step 3.

As before, steps 3 to 7 are repeated until both the sample mean and SD of BMI
are within $\Delta_{\text{tol}}$ of the observed data, or the maximum number
of repetitions is reached.

##### Aside: Limitations of the LMS method when calculating BMI from zBMI

The calculation of BMI from zBMI in Equation (22) involves the exponent
$1/\lambda(a_{i},x_{i})$. In the three major reference charts (CDC, WHO and
IOTF), $\lambda<-1$ for nearly all values of $(a_{i},x_{i})$. This means that
the exponent almost always represents a root. Since we can only take the root
of a positive number, the base of this exponent must be non-negative,

$\displaystyle 1+\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})Z_{i}>0,$ (23)

where $\sigma(a_{i},x_{i})>0$ and $Z_{i}\in[-\infty,\infty]$. Sampled values
of zBMI are then restricted to

$\displaystyle Z_{i}<\frac{-1}{\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})}.$ (24)

The minimum value of $-1/\lambda\sigma$ across all the reference charts is
2.785, which occurs in the CDC chart for females at $a_{i}=240.5$ months.
Therefore, for this chart, a sampled value of $Z_{i}$ greater than 2.785 will
give an invalid value of BMI when $x_{i}=0$ (female) and $a_{i}=240.5$. The
other reference charts have less restrictive cut-offs, but encounter the same
issue for large values of $Z_{i}$ at particular ages. To deal with this
problem, we truncate any problematic $Z_{i}$ samples to a value close to the
boundary. In particular, if $Z_{i}>-1/\lambda(a_{i},x_{i})\sigma(a_{i},x_{i})$
we set that sample equal to
$0.99\times(-1/\lambda(a_{i},x_{i})\sigma(a_{i},x_{i}))$. This sets the base
of the exponent equal to 0.01, which produces a large value of $B_{i}$.

Because we have to artificially restrict the distribution of $Z$ for certain
samples, the mean and SD of the sampled values at any iteration,
$\mathcal{D}_{Z}^{(t)}=(\bar{Z}^{(t)},\text{SD}^{(t)}(Z))$, will not
necessarily be the same as the distribution parameters,
$(m_{z}^{(t)},s_{z}^{(t)})$. In Section 5.2.2 we investigate the behaviour of
both. However, the restriction is only necessary in a minority of cases so we
do not expect it to have a major impact.

##  5 Application

###  5.1 Methods

To investigate the behaviour of the mapping methods, we compiled two test data
sets from the obesity data. The first includes results from all trials that
report arm-level means and SDs on both the BMI and zBMI scale (205
observations from 41 trials). The second includes all trials with arm-level
data for both percentile and zBMI (54 observations from 10 trials). We refer
to these as the $BZ$ and $PZ$ data sets respectively. In the following, we
apply the mapping methods to either the BMI or percentile data in each test
data set and compare the mapped values to the corresponding values reported on
the zBMI scale.

Figure 2: A plot comparing the different mapping methods for estimating
aggregate zBMI data from aggregate percentile data.

####  5.1.1 Reconstructing arm-level data

Most trials in our test data sets report arm-level data at baseline and
follow-up for each arm of the trial. Therefore, a trial with $A_{i}$ arms and
$T_{i}$ follow-up times, contributes $A_{i}(T_{i}+1)$ data points. Some
trials, on the other hand, only provide arm-level data at one time point (e.g.
baseline or follow-up alone) and therefore contribute $A_{i}$ data points.
Other trials report baseline values alongside mean change from baseline
(change scores). In these cases, we reconstructed follow-up means by adding
the change score to the baseline score, and set the follow-up SD on each arm
equal to the corresponding SD at baseline. The latter assumption is supported
by the observation in other trials that SDs at different time points are
approximately equal.

The obesity data include cluster-randomized trials. To construct our test data
sets, we extracted raw values that had not been adjusted for clustering. When
only adjusted values were available, we ‘unadjusted’ the SDs by dividing by
the square root of the design effect. Where available, we defined the design
effect in terms of the intra-class correlation coefficient (ICC) reported in
the trial. Otherwise, we assumed a value of $\text{ICC}=0.02$ based on ICCs
reported in other trials (see our previous analyses6, 7).

####  5.1.2 Age and sex

To perform the mapping $\mathcal{D}_{B}\xrightarrow{}\mathcal{D}_{Z}$, we
require information about the distribution of age and sex in the sample.
Therefore, to construct the $BZ$ data set, we extracted the mean and SD of age
(at baseline) and the proportion of male participants in each trial. For
follow-up measurements, we increased the mean age by the length of follow-up.
In some trials, age was reported as a range. Here, we assumed a normal
distribution such that the mean is given by the mid point, and the SD is
approximately a quarter of the range. When the proportion of males and females
was not reported, we assumed an even split ($p_{m}=0.5$).

####  5.1.3 Reference charts

Mapping from BMI to zBMI also requires that we have age- and sex-specific LMS
values for an appropriate reference population. Where reported, we extracted
which reference chart was used by the trial authors. In all but one trial,29
which used LMS values specific to the Flemish population,‡‡‡The Flemish
reference chart cited in Haerens et al (2006)29 is no longer accessible.
Instead we used the updated chart published by the same authors in 200930. the
trialists reported CDC, WHO or IOTF charts. Where a reference was not given,
we used the CDC chart for trials set in the US and the IOTF chart for trials
in any other country.

The sex-specific reference charts were downloaded as Microsoft Excel files
with age tabulated against $\lambda,\mu,$ and $\sigma$. For each sample
$(a_{i},x_{i})$, we identified the tabulated age closest to $a_{i}$ on the
appropriate reference chart for sex $x_{i}$ and read off the corresponding LMS
parameters. Each reference chart reports parameters for a different range of
ages, at different intervals. For example, the CDC chart reports LMS
parameters for ages 2 to 20 years in steps of 1 month, the IOTF chart goes
from 2 to 18 years in steps of 6 months and WHO ranges from 5 to 19 years in
steps of 1 month. When sampling age, we re-sampled any value of $a_{i}$ that
was outside the range of the relevant reference chart.

###  5.2 Results

####  5.2.1 zBMI from percentile

Table 1: The root mean squared error (RMSE) and mean absolute error (MAE) on the estimated aggregate data (mean and SD) from the percentile mapping methods. | Estimated mean zBMI | Estimated SD of zBMI
---|---|---
Method | RMSE | MAE | RMSE | MAE
Analytical | 0.147 | 0.090 | 0.303 | 0.138
Sampling | 0.196 | 0.122 | 0.204 | 0.135
Optimization (sample) | 0.119 | 0.079 | 0.146 | 0.089
Optimization (distribution) | 0.119 | 0.079 | 0.148 | 0.090

To perform the mapping $\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$, we
implemented the sampling method (Algorithm 1) with $N=10^{4}$ samples. For the
optimization method (Algorithm 3), we used $N=10^{3}$ samples and a maximum of
$n_{\max}=5\times 10^{3}$ iterations. To pick $\Delta_{\text{step}}$ and
$\Delta_{\text{tol}}$, we explored how changes in the zBMI distribution lead
to changes in the percentile distribution. We used this to propose a range of
sensible values for $(\Delta_{\text{step}},\Delta_{\text{tol}})$ and
investigated the corresponding convergence of the algorithm (see Section C of
the Supplementary Material). Based on these preliminary tests, we chose
$\Delta_{\text{step}}=0.002$ and $\Delta_{\text{tol}}=0.005$ which led to 100%
convergence.

Figure 2 shows estimated zBMI against reported zBMI for the different
$\mathcal{D}_{P}\xrightarrow{}\mathcal{D}_{Z}$ methods: analytical, sampling
and optimization. For the latter, we show both the sample estimates,
$\mathcal{D}_{Z}^{(t)}=(\bar{Z}^{(t)},\text{SD}(Z)^{(t)})$, and the
distribution estimates, $(m_{z},s_{z})$. The left panel of Figure 2 depicts
the results for mean zBMI and the right panel shows the corresponding standard
deviation.

From the Figure, it appears the different methods all exhibit similar
accuracy, though sampling has slightly higher estimates on average (for both
mean and SD) and optimization slightly lower. As expected, the sample and
distribution estimates from the optimization method are very similar. All
methods do a reasonable job of recreating mean zBMI, with most data points
approximately following the ‘true’ line. Estimated SDs are more scattered, but
all take values within a realistic range.

The right-most data point on the SD plot is underestimated by all methods and
represents a notable outlier. This value corresponds to the follow-up
measurement of the control arm in a trial by Brown et al. (2013).31 On
inspection we find that, on the zBMI scale, this SD measurement (2.21) is
approximately twice as large as the SD reported for the other arms and time
points (which range from 1.14 to 1.16). However, on the percentile scale this
SD (27.3) is very similar to the other values (26.5-27.4). It is therefore not
surprising that our mapping methods fail to recreate this value. This
potentially points to a reporting error in the trial.

In Table 1 we compare the percentile mapping methods using the root mean
squared error (RMSE) and the mean absolute error (MAE). These metrics quantify
how accurately the mapped percentile data recreate the corresponding zBMI data
reported in the trials. While all methods exhibit similar errors, the
optimization method (sample estimates) produces the most accurate measurements
of both mean and SD. According to these metrics, the sampling method performs
worst for the mean and the analytical method performs worst for the SD.

####  5.2.2 zBMI from BMI

To map from BMI to zBMI using algorithms 2 and 4, we chose the same values of
$N$ and $n_{\max}$ as we did for the percentile methods. As shown in Section C
of the Supplementary Material, we explored various values of
$(\Delta_{\text{step}},\Delta_{\text{tol}})$ for Algorithm 4 and they all led
to comparable estimates. We show results for $\Delta_{\text{step}}=0.01$ and
$\Delta_{\text{tol}}=0.1$ which gave the highest convergence rate (90%).

Figure 3 shows estimated zBMI against reported zBMI for the
$\mathcal{D}_{B}\xrightarrow{}\mathcal{D}_{Z}$ sampling and optimization
methods. Mean values are plotted in the left panel and SDs on the right. For
both methods, we present the results when age is sampled from a normal
distribution. As shown in Section D of the Supplement, we observe negligible
differences between the results for different age distributions.

Figure 3: A plot comparing the different mapping methods for estimating aggregate zBMI data from aggregate BMI data. Table 2: The root mean squared error (RMSE) and mean absolute error (MAE) on the estimated aggregate data (mean and SD) from the BMI mapping methods. | Estimated mean zBMI | Estimated SD of zBMI
---|---|---
Method | RMSE | MAE | RMSE | MAE
Sampling | 0.355 | 0.219 | 0.505 | 0.338
Optimization (sample) | 0.339 | 0.201 | 0.284 | 0.196
Optimization (distribution) | 0.341 | 0.204 | 0.283 | 0.195

On inspection of Figure 3, we find that the sample and distribution estimates
from the optimization method give comparable results, indicating that the
restriction placed on extreme values of zBMI has a negligible effect. For mean
zBMI, the sampling and optimization methods also appear to be similarly
accurate. Except for a cluster of outliers (indicated in the figure), they
recreate the mean values reasonably well. The circled outliers in Figure 3 are
all from the same trial, Papadaki et al. (2010).32 On inspection of this
study, there is no obvious cause of the issue. The trialists state that weight
status was determined using cut-offs from the IOTF, but they do not report how
zBMI was calculated. It is possible that they used the simple definition of a
z-score for normal distributions, rather than the LMS method. If so, this may
explain the discrepancy in these results.

As shown in the right panel of Figure 3, the optimization method produces more
accurate estimates of SD than the sampling method, which tends to
overestimate. However, the optimization method is more likely to underestimate
SD. In the context of using the data in an analysis, underestimating SD is
less desirable as it leads to an overweighting of the mapped values. It is
possible that optimization leads to smaller SDs than the sampling method
because of the truncation of extreme zBMI samples. However, the similarity of
the sample and distribution estimates from the optimization method suggests
that this does not have an important effect. The left-most data points on the
right panel reflect reported values which seem unrealistically small for
standard deviations. It is therefore not surprising they are estimated poorly
by both methods. On inspection we found no issue with data extraction so this
could point to a reporting error in the trial. The right-most data point is
the same outlier from Brown et al. (2013) that we observed in the percentile
results.

In the Supplementary Material (Section E) we explore three potential causes of
inaccurate results: (i) studies that don’t report a reference chart, (ii) SDs
that we unadjusted for clustering, and (iii) follow-up values imputed from
change scores. Excluding measurements corresponding to (i)-(iii) removes the
most significant scatter on the plots for both mean and SD. To further explore
the effect of reference charts, in Section F of the Supplement we plot the
results separately for trials which reported each chart. For mean zBMI, no
chart appears to perform systematically better or worse than any other. For
standard deviation, the results associated with the CDC chart are associated
with the largest scatter.

In Section G of the Supplementary Material we show the results of the
optimization method excluding any estimates which did not converge. The plot
is largely similar to Figure 3 but, notably, does not include the right-most
outlier on the SD plot.

Finally, Table 2 compares the BMI mapping methods in terms of RMSE and MAE.
According to these metrics, all methods exhibit similar accuracy in the
estimation of mean zBMI, but the sampling method has the largest error on SD.
While these results are in agreement with the observations from Figure 3, they
do not reflect the tendency of the optimization method to underestimate SD.

##  6 Discussion

In this paper, we devised a range of methods for mapping between aggregate
arm-level observations of BMI outcomes by leveraging known relationships
between individual measurements on the different scales. Specifically, we
proposed several methods for mapping percentile and un-standardized BMI onto
zBMI. When applied to our motivating dataset, all methods yielded reasonable
results, with estimated mean zBMI values aligning well with those reported in
the trials. While the methods performed less effectively in reproducing
standard deviations on zBMI, all estimated SDs fell within a realistic range.

Assuming a normal distribution for zBMI, we developed an analytical method
based on standard integral results to map observations of mean and SD from the
percentile to the zBMI scale. This deterministic approach, unlike numerical
methods, remains unaffected by sample size or convergence issues, making it
computationally less demanding. However, its specificity to the definition of
z-scores and percentiles limits its broader applicability. Conversely, both
the sampling and optimization methods offer versatility beyond our initial
dataset, enabling mapping between any pair of outcomes with known individual-
level relationships. To map from scale $X$ to scale $Y$, the sampling method
assumes and samples from a distribution for scale $X$, while the optimization
method samples from a distribution for $Y$. Therefore, the ease of sampling on
each scale can guide the selection of which method to use.

All the methods we developed are sensitive to the choice of distributions for
the different scales. Guided by theoretical and practical considerations, we
assumed a beta distribution for percentile, a lognormal distribution for BMI,
and a normal distribution for zBMI. Future work could explore the validity of
these assumptions, for example using individual participant data. However, we
expect that differences in trial demographics and eligibility criteria will
cause variability in distributions across trials. Therefore, any choice of
distribution will necessarily be an approximation. We believe the
distributions we assumed are sensible choices given the information available.

The relationship between BMI and zBMI relies on information about an
individual’s age and sex, thus introducing additional assumptions when mapping
from BMI. These include assumptions about the age distribution in the sample
and the reference chart used to define the LMS parameters. Furthermore, when
extracting the age-specific LMS values we used the age tabulated in the
reference chart that was closest to the sampled age. A more sophisticated
approach could interpolate between these values. However, we do not expect
that this method had a large impact on our results. Indeed, trial authors
themselves are likely to have used tabulated rather than interpolated values,
meaning our approach more closely resembles the probable data generating
process. Another important assumption we make is that all trial authors used
the LMS method to evaluate zBMI. While this is by far the most common
definition, it is possible that some authors used the standard definition of a
z-score assuming a normal distribution and some measure of mean BMI. We did
not find any reference to the latter, but several trials did not explain their
method for calculating zBMI. Considering the number of assumptions we had to
make, we believe the BMI mapping methods performed remarkably well.

Another related concern centres around potential ‘inaccuracies’ in the LMS
values themselves. Due to demographic differences and changes over time, the
reported parameters may not be appropriate for a particular trial or set of
trials. Our findings indicate that the BMI sampling method tends to
overestimate SD, while optimization is more likely to underestimate. Since the
two methods operate the LMS relation in the opposite direction
($B\xrightarrow{}Z$ for the sampling method and $Z\xrightarrow{}B$ for
optimization), the observed discrepancies may stem from systematic ‘errors’ in
the LMS parameters. Indeed, we found that removing data from trials that did
not report a reference chart eliminated some of the most inaccurate results.
This suggests that the accuracy of the mapping methods relies on appropriate
LMS parameters. Furthermore, we observed that the most extreme scatter on the
SD results was associated with the CDC reference chart, pointing to a
potential issue with these values. In future work, we plan to investigate
these ideas further, for example examining whether variations in accuracy are
related to the age or average BMI of the population. We will also explore
possible methods for adjusting the LMS parameters to use in specific
populations.

In our application to the motivating dataset, the RMSE and MAE metrics
slightly favoured the optimization method over the sampling method. However,
optimization is sensitive to convergence and may underestimate standard
deviation when mapping from BMI. While the sampling method gives less accurate
estimates of SD from BMI, it leans towards conservatism, reducing the risk of
overweighting mapped data in subsequent analyses. To map from BMI to zBMI, the
sampling method samples values of BMI while the optimization method samples
zBMI. Both methods also sample age and sex, and do this independently from the
sampling of the outcome scale. Given the dependence of BMI on age and sex, the
sampling method might be improved by constructing an appropriate joint
distribution for these variables. On the other hand, independence from age and
sex is a more reasonable assumption for zBMI, which may explain the superior
performance of the optimization method.

To validate and compare our methods, we relied on aggregate data reported in
trials, which are susceptible to misreporting and misinterpretation (as
highlighted by some outliers in our results). Thus, small differences in RMSE
and MAE should be interpreted cautiously. In future work, individual
participant data might facilitate more accurate assessments of the methods.

Our optimization scheme, while largely effective, employed a naïve approach of
incrementing or decrementing parameters by a fixed amount at each iteration.
Additionally, we assumed that the mean of the percentile or BMI samples was
only influenced by the mean of the zBMI distribution and similarly, that the
SD of the samples depended only on the SD of the distribution. Our
investigation in Section C of the Supplement revealed that changes in both
zBMI parameters could affect both the resulting sample statistics, suggesting
the method may benefit from a more sophisticated multivariate optimization
scheme. However, since we managed to achieve a high level of convergence using
our naïve scheme (100% for percentile and 90% for BMI), we anticipate that any
improvements would be limited to the efficiency and reliability of the method
rather than accuracy.

We have discussed a number of ways our methods can be improved and refined. In
future work, we plan to develop this work further and explore the utility of
our methods in other applications such as mapping between different depression
scales. However, we recognise that even the best mapping method is not a
substitute for data reported directly in trials. Our motivating example
highlights that, in order to synthesize evidence about childhood weight
status, trials are most useful when they report arm-level data on all BMI
scales (un-standardized, zBMI and percentile). Trial reports should state the
precise methods used to perform standardization, along with the relevant
reference chart. Since childhood BMI is correlated with age and sex, it would
be desirable for trials to describe the distribution of ages and the
proportion of males and females in each intervention arm. It would also be
useful for trials to report on the observed distribution of each outcome scale
in their population. This would improve our understanding of the different
scales and help to map between outcomes in historical trials.

##  7 Summary

To the best of our knowledge, this study is the first attempt to develop
methods to map between arm-level aggregate data on different measurement
scales. Previous techniques have involved simple re-scaling of contrast-level
data via standardization and ratio techniques, or by fitting approximate
models. In contrast, our methods leverage known transformations between scales
at the individual level. The sampling and optimization methods possess broad
applicability for mapping between any pair of measurement scales with known or
estimable relationships between individual measurements. These mapping
techniques enhance the potential for meta-analysis by facilitating the
creation of larger datasets, enabling a more comprehensive summary of evidence
from trials and improving the precision of estimates.

## Author contributions

Annabel L Davies: writing - original draft (lead), formal analysis (lead),
methodology (equal), software (lead), visualization (lead), writing - review &
editing (equal). A E Ades: conceptualization (equal), methodology (equal),
writing - review & editing (equal). Julian PT Higgins: conceptualization
(equal), funding acquisition (lead), writing - review & editing (equal).

## Acknowledgments

ALD and JPTH acknowledge funding from the National Institute for Health and
Care Research (NIHR) Public Health Research programme (NIHR131572). JPTH is an
NIHR Senior Investigator.

## Data availability statement

The data and code that support the findings of this paper are provided in the
GitHub repository here: https://github.com/AnnieDavies/Mapping_BMI

## Financial disclosure

None reported.

## Conflict of interest

The authors declare no potential conflict of interests.

## Supporting information

Supplementary material is provided online.

## Highlights

What is already known:

  * •

Meta-analysis of continuous data is difficult when trials report aggregate
outcomes on different measurement scales.

  * •

Existing methods to synthesize evidence on different scales rely on
simplifying assumptions such as linearity and non-negative outcomes.

What is new:

  * •

We develop three methods (analytical, sampling and optimization) for mapping
between arm-level aggregate data on different measurement scales.

  * •

We validate and compare the methods using trials that report childhood obesity
outcomes on different body mass index (BMI) scales.

Potential impact for readers:

  * •

The sampling and optimization methods proposed can be used to map between
aggregate data on any pair of measurement scales, provided we know, or can
estimate, the relationship between individual measurements on the different
scales.

  * •

Therefore, the methods facilitate the creation of larger, more comprehensive
datasets for meta-analysis.

## References

  * 1 Higgins JPT, Thomas J, Chandler J, et al. , eds.Cochrane Handbook for Systematic Reviews of Interventions.  Chichester, UK: John Wiley & Sons.  2 ed. 2019.
  * 2 Dias S, Ades AE, Welton NJ, Jansen JP, Sutton AJ. Network Meta-Analysis for Decision Making.  Oxford, UK: Wiley .  2018\.
  * 3 The National Institute of Health and Care Excellence . Developing NICE guidelines: the manual. Online; 2014.  Available from: https://www.nice.org.uk/process/pmg20.
  * 4 Daly C, Welton NJ, Dias S, Anwer S, Ades A. Meta-Analysis of Continuous Outcomes. Guideline Methodology Document 2. NICE Guidelines Technical Support Unit; 2021.  Online. Available from https://www.bristol.ac.uk/population-health-sciences/centres/cresyda/mpes/nice/guideline-methodology-documents-gmds/.
  * 5 Lu G, Kounali D, Ades AE. Simultaneous multioutcome synthesis and mapping of treatment effects to a common scale. Value Health 2014; 17(2): 280–287.
  * 6 Spiga F, Davies AL, Tomlinson E, et al. Supporting Evidence: Interventions to prevent obesity in children aged 5 to 11 years old. Overweight and obesity management (Guideline in development) [GID-NG10182], National Institute for Health and Care Excellence; 2023.  Available from: https://www.nice.org.uk/guidance/gid-ng10182/documents/supporting-evidence-3.
  * 7 Spiga F, Tomlinson E, Davies AL, et al. Supporting Evidence: Interventions to prevent obesity in children aged 12 to 18 years old. Overweight and obesity management (Guideline in development) [GID-NG10182], National Institute for Health and Care Excellence; 2023.  Available from: https://www.nice.org.uk/guidance/gid-ng10182/documents/supporting-evidence-4.
  * 8 Adab P, Pallan M, Whincup PH. Is BMI the best measure of obesity?. BMJ 2018; 360.  doi: 10.1136/bmj.k1274
  * 9 de Onis M, Lobstein T. Defining obesity risk status in the general childhood population: Which cut-offs should we use?. International Journal of Pediatric Obesity 2010; 5(6): 458-460.
  * 10 Must A, Anderson S. Body mass index in children and adolescents: considerations for population-based applications. Int. J. Obes. 2006; 30: 590-594.
  * 11 Brown T, Moore THM, Hooper L, et al. Interventions for preventing obesity in children. Cochrane Database of Systematic Reviews 2019(7).
  * 12 Hodder RK, O’Brien KM, Lorien S, et al. Interventions to prevent obesity in school-aged children 6-18 years: An update of a Cochrane systematic review and meta-analysis including studies from 2015-2021. EClinicalMedicine 2022; 54: 101635.
  * 13 Cummings P. Arguments for and Against Standardized Mean Differences (Effect Sizes). Archives of Pediatrics and Adolescent Medicine 2011; 165(7): 592-596.  doi: 10.1001/archpediatrics.2011.97
  * 14 Friedrich JO, Adhikari NK, Beyene J. The ratio of means method as an alternative to mean differences for analyzing continuous outcome variables in meta-analysis: a simulation study. y. BMC Medical Research Methodology 2008: 8:32.
  * 15 Centers for Disease Control and Prevention . CDC Growth Charts: Data Tables: Percentile Data Files with LMS Values. Online; .  Available from: https://www.cdc.gov/growthcharts/percentile_data_files.htm. Accessed September 2023.
  * 16 World Health Organization . Obesity: preventing and managing the global epidemic. Report of a WHO consultation. World Health Organization technical report series, 894; 2000.
  * 17 Silverman M, Lipscombe T. Exact Statistical Distribution of the Body Mass Index (BMI): Analysis and Experimental Confirmation. Open Journal of Statistics 2022; 12: 324–356.
  * 18 Cole TJ. Fitting Smoothed Centile Curves to Reference Data. Journal of the Royal Statistical Society. Series A (Statistics in Society) 1988; 151(3): 385–418.
  * 19 Cole TJ, Green PJ. Smoothing reference centile curves: The lms method and penalized likelihood. Statistics in Medicine 1992; 11(10): 1305-1319.
  * 20 Cole TJ, Bellizzi MC, Flegal KM, Dietz WH. Establishing a standard definition for child overweight and obesity worldwide: international survey. BMJ 2000; 320(7244): 1240.
  * 21 Kuczmarski RJ, Ogden CL, Grummer-Strawn LM, et al. CDC growth charts: United States. Adv Data 2000; 8(314): 1-27.
  * 22 Kuczmarski RJ, Ogden CL, Guo SS, et al. 2000 CDC Growth Charts for the United States: methods and development. Vital Health Stat 11 2002(246): 1-190.
  * 23 WHO Multicentre Growth Reference Study Group , de Onis M. WHO Child Growth Standards based on length/height, weight and age. Acta Paediatrica 2006; 95(S450): 76-85.
  * 24 de Onis M, Onyango AW, Broeck dJV, Chumlea WC, Martorell R. Measurement and Standardization Protocols for Anthropometry Used in the Construction of a New International Growth Reference. Food and Nutrition Bulletin 2004; 25(1_suppl_1): S27-S36.
  * 25 World Health Organization . BMI-for-age (5-19 years): Expanded tables for constructing national health cards. Online; .  Available from: https://www.who.int/tools/growth-reference-data-for-5to19-years/indicators/bmi-for-age. Accessed September 2023.
  * 26 Cole TJ, Lobstein T. Extended international (IOTF) body mass index cut-offs for thinness, overweight and obesity. Pediatric Obesity 2012; 7(4): 284-294.
  * 27 Owen D. A table of normal integrals. Communications in Statistics: Simulation and Computation 1980; 4(B9): 389–419.
  * 28 Owen DB. Tables for computing bivariate normal probabilities. Annals of Mathematical Statistics 1956; 27: 1075–1090.
  * 29 Haerens L, Deforche B, Maes L, Stevens V, Cardon G, De Bourdeaudhuij I. Body Mass Effects of a Physical Activity and Healthy Food Intervention in Middle Schools. Obesity 2006; 14(5): 847-854.
  * 30 Roelants M, Hauspie R, Hoppenbrouwers K. References for growth and pubertal development from birth to 21 years in Flanders, Belgium. Annals of Human Biology 2009; 36(6): 680-694.
  * 31 Brown B, Noonan C, Harris KJ, et al. Developing and Piloting the Journey to Native Youth Health Program in Northern Plains Indian Communities. The Diabetes Educator 2013; 39(1): 109-118.
  * 32 Papadaki A, Linardakis M, Larsen TM, et al. The Effect of Protein and Glycemic Index on Children’s Body Composition: The DiOGenes Randomized Study. Pediatrics 2010; 126(5): e1143-e1152.

##  Appendix A Smooth updating scheme for the optimization method

In this section we derive the updating scheme implemented in the optimization
method (Equation (21)). Let us assume that, at iteration $t$, zBMI is sampled
from a normal distribution with mean $m_{z}^{(t)}$ and standard deviation
$s_{z}^{(t)}$. Then, at iteration $t+1$, the mean and SD are increased by
$\delta_{m}$ and $\delta_{s}$ respectively. The resulting sampling
distributions at each iteration are

$\displaystyle Z_{i}^{(t)}$
$\displaystyle\sim\mathcal{N}(m_{z}^{(t)},s_{z}^{(t)2}),$ (25) $\displaystyle
Z_{i}^{(t+1)}$
$\displaystyle\sim\mathcal{N}(m_{z}^{(t)}+\delta_{m},(s_{z}^{(t)}+\delta_{s})^{2}).$
(26)

To write $Z_{i}^{(t+1)}$ as a function of $Z_{i}^{(t)}$ we transform both
variables to a standard normal and equate. If a random variable $X$ is
normally distributed via $X\sim\mathcal{N}(m,s^{2})$, then the transformed
variable $\frac{X-m}{s}$ follows a standard normal. Therefore,

$\displaystyle\frac{Z_{i}^{(t)}-m_{z}^{(t)}}{s_{z}^{(t)}}$
$\displaystyle\sim\mathcal{N}(0,1),$ (27)
$\displaystyle\frac{Z_{i}^{(t+1)}-(m_{z}^{(t)}+\delta_{m})}{s_{z}^{(t)}+\delta_{s}}$
$\displaystyle\sim\mathcal{N}(0,1).$ (28)

Equating the left hand side of Equations (27) and (28) leads to

$\displaystyle\frac{Z_{i}^{(t)}-m_{z}^{(t)}}{s_{z}^{(t)}}=\frac{Z_{i}^{(t+1)}-(m_{z}^{(t)}+\delta_{m})}{s_{z}^{(t)}+\delta_{s}}$
(29)

which, rearranging for $Z_{i}^{(t+1)}$, gives the result quoted in Section 4.3
of the main paper,

$\displaystyle
Z_{i}^{(t+1)}=\left(1+\frac{\delta_{s}}{s_{z}^{(t)}}\right)(Z_{i}^{(t)}-m_{z}^{(t)})+(m_{z}^{(t)}+\delta_{m}).$
(30)


#  Lifshitz model in the presence of spin-orbit coupling

M. E. Raikh  Department of Physics and Astronomy, University of Utah, Salt
Lake City, UT 84112

###### Abstract

Wave function of a localized state created by a short-range impurity in two
dimensions falls off with distance, $r$, from the impurity as
$\frac{1}{r^{1/2}}\exp(-\frac{r}{a})$, where $a$ is the localization radius.
With randomly positioned identical impurities with low concentration, $n\ll
a^{-2}$, the level smears into a band due to the overlap of the impurity wave
functions. This is the essence of the Lifshitz model. We demonstrate that,
upon incorporation of the spin-orbit coupling, the impurity wave functions
acquire oscillating factors which, subsequently, modify their overlap. As a
result of such modification, the density of states develops singularities at
certain energies.

###### pacs:

73.50.-h, 75.47.-m

##  I Introduction

Conventional Lifshitz modelLifshitz ; Gredeskul  describes single-electron
states in a system of randomly-positioned identical short-range impurities. It
is assumed that each impurity creates a localized level with a radius of wave
function, $a$, much smaller than a typical distance between the neighboring
impurities. In two dimensions this distance is $n^{-1/2}$, where $n$ is the
impurity concentration. Due to overlap of the wave functions of neighboring
impurities, individual levels smear into an impurity band. Presence of a small
parameter $na^{2}\ll 1$ allows to establish the form of the density of states
as a function of $E-E_{0}$, which is the electron energy measured from the
level position, $E_{0}$. The key argument summarized in the book Ref. book,
why this density can be found analytically is that the eigenstates of the
disordered system are composed of pairs of impurities which hybridize the
respective wave functions. As a result of this hybridization one of the levels
shifts up from $E_{0}$, while the other level shifts down. Thus, the density
of states in the Lifshitz model has the form of two peaks.

Assume now that the bare Hamiltonian of the 2D system contains a spin-orbit
term e.g. of the Rashba type.Rashba  This certainly does not have an effect on
the density of states if the spin splitting of the bare band is much smaller
than $E_{0}$. In the opposite limit, the localized states are classified
according to helicity. Importantly, the wave functions of the short-range
impurities get modified dramatically by spin-orbit coupling compared to a
simple exponential decay.Galstyan ; Chaplik ; Mkhitaryan ; Voloshin ;
Hutchinson  The main message of the present paper is that in the limit of
strong spin-orbit coupling the shape of the density of states in the Lifshitz
model changes dramatically, namely, it develops singularities at certain
energies.

##  II A single short-range impurity in the presence of spin-orbit coupling

The density of states of free electrons with a quadratic Hamiltonian
${\hat{H}}_{0}=\frac{\hbar^{2}{\bm{k}}^{2}}{2m}$, where $m$ is the effective
mass and ${\bm{k}}$ is the wave vector, is energy-independent. Upon adding a
spin-orbit term

${\hat{H}}_{SO}=\alpha[{\bm{k}}\times{\hat{\bm{\sigma}}}]\cdot{\bm{n}},$ (1)

where ${\bm{\sigma}}$ is the vector of the Pauli matrices, $\alpha$ is the
spin-orbit constant, and ${\bm{n}}$ is the normal to the 2D plane, to
${\hat{H}}_{0}$ modifies the density of states qualitatively.Galstyan ;
Chaplik  Indeed, the spectrum of the full Hamiltonian
${\hat{H}}_{0}+{\hat{H}}_{SO}$ consists of two branches

$E_{1,2}({\bm{k}})=\frac{\hbar^{2}k^{2}}{2m}\mp\alpha|{\bm{k}}|.$ (2)

The corresponding wave functions have the form

$\Psi_{{\bm{k}}}^{(1,2)}({\bm{r}})=e^{i{\bm{k}r}}\chi_{{\bm{k}}}^{(1,2)},$ (3)

where the spinors $\chi_{{\bm{k}}}^{(1,2)}$ are defined as

$\chi_{{\bm{k}}}^{(1)}=\frac{1}{\sqrt{2}}\begin{pmatrix}e^{i\Phi_{{\bm{k}}}}\\\
-1\end{pmatrix},\hskip
8.53581pt\chi_{{\bm{k}}}^{(2)}=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
e^{-i\Phi_{{\bm{k}}}}\end{pmatrix}.$

Here $\Phi_{{\bm{k}}}$ is the azimuthal angle of the vector ${\bm{k}}$. It is
crucial that the lower branch of the spectrum Eq. (2) has a minimum at

$k=k_{0}=\frac{m\alpha}{\hbar^{2}}.$ (4)

Near this minimum the spectrum can be simplified as

$E_{1}({\bm{k}})=-\Delta+\frac{\hbar^{2}(k-k_{0})^{2}}{2m},$ (5)

with the depth of the minimum
$\Delta=\frac{m\alpha^{2}}{2\hbar^{2}}=\frac{\hbar^{2}k_{0}^{2}}{2m}$. As a
result, the wave function of a localized state created by a point-like
impurity is composed of the free-electron state with energies close to
$-\Delta$. The density of these states behaves as $(E+\Delta)^{-1/2}$, i.e. it
has a one-dimensional character.Galstyan

Calculation of the level position in this setting was carried out in Ref.
Chaplik, . We reproduce it below for pedagogical reasons and in order to
generalize it later to the impurity pairs.

Denote with $U({\bm{r}})$ the impurity potential. The solution of the
Schrödinger equation

$\Bigl{[}{\hat{H}_{0}}+{\hat{H}_{SO}}+U({\bm{r}})\Bigr{]}\Psi({\bm{r}})=E_{0}\Psi({\bm{r}})$
(6)

for the level, $E_{0}$, close to $-\Delta$ can be searched in the form of the
combination of the states of the lower branch only

$\Psi({\bm{r}})=\int
d{\bm{k}}A({\bm{k}})\chi_{{\bm{k}}}^{(1)}e^{-i{\bm{k}r}}.$ (7)

Substitution of this form into Eq. (6) yields

$\int
d{\bm{k}}A({\bm{k}})\chi_{{\bm{k}}}^{(1)}e^{-i{\bm{k}r}}\Bigl{[}E_{1}({\bm{k}})-E_{0}\Bigr{]}=-U({\bm{r}})\Psi(0).$
(8)

In Eq. (8) we took into account that $U({\bm{r}})$ is short-ranged.
Multiplying this equation by $\chi_{{\bm{q}}}^{(1)*}e^{i{\bm{q}r}}$ and
integrating $d{\bm{r}}$, we get

$(2\pi)^{2}A({\bm{q}})\Bigl{(}E_{1}({\bm{q}})-E_{0}\Bigr{)}=-\Bigl{(}\Psi(0)\chi_{{\bm{q}}}^{(1)*}\Bigr{)}\int
d{\bm{r}}U({\bm{r}})e^{i{\bm{q}r}}.$ (9)

Expressing $A({\bm{q}})$, substituting it into Eq. (7) and setting
${\bm{r}}=0$, we arrive to the self-consistency equation

$\Psi(0)=-\frac{1}{(2\pi)^{2}}\int
d{\bm{q}}\frac{\chi_{{\bm{q}}}^{(1)}\Bigl{(}\Psi(0)\chi_{{\bm{q}}}^{(1)*}\Bigr{)}}{E_{1}({\bm{q}})-E_{0}}\int
d{\bm{r}}U({\bm{r}})e^{i{\bm{q}r}}.$ (10)

This equation defines the position of the level, $E_{0}$. For isotropic
potential, the dependence on the direction of ${\bm{q}}$ disappears from the
integral $d{\bm{r}}$, yielding $J_{0}(qr)$, where $J_{0}$ is the Bessel
function. Then the angular integration over ${\bm{q}}$ can be easily
performed. Concerning the integral over $|{\bm{q}}|$, it comes from the domain
$\left(|{\bm{q}}|-k_{0}\right)\ll k_{0}$. Finally, the solution of Eq. (10)
takes the formChaplik

$\varepsilon_{0}=-\Delta-
E_{0}=\frac{\pi^{2}mk_{0}^{2}}{2\hbar^{2}}\Biggl{[}\int\limits_{0}^{\infty}drrU(r)J_{0}(k_{0}r)\Biggr{]}^{2}.$
(11)

The solution is doubly degenerate with respect to the components of respective
spinors. The meaning of $\varepsilon_{0}$ is the binding energy measured from
the minimum of the spectrum of the lower branch.

To estimate the binding energy predicted by Eq. (11), we assume that the
magnitude of the potential is $U_{0}$, while the radius of potential is $a$.
To insure that the wave function does not change within the interval $r<a$ the
condition $k_{0}a<<1$ should be met, which is equivalent to the replacement of
the Bessel function in the integrand by $1$. Then, within a numerical factor,
we get

$\varepsilon_{0}\sim(k_{0}a)^{2}\frac{U_{0}^{2}}{\frac{\hbar^{2}}{ma^{2}}}.$
(12)

To test the assumptions made in course of solving of Eq. (10) this binding
energy must be much smaller than the depth $\frac{\hbar^{2}k_{0}^{2}}{2m}$ of
the minimum in the spectrum (to justify the integration over ${\bm{q}}$).
Also, this binding energy should be much smaller than $U_{0}$ (to replace
$\Psi({\bm{r}})$ by $\Psi(0)$. The first requirement leads to the usual
condition $U_{0}\ll\frac{\hbar^{2}}{ma^{2}}$. The second requirement yields a
complimentary condition $(k_{0}a)^{2}U_{0}\ll\frac{\hbar^{2}}{ma^{2}}$, which
is weaker.

The form of the wave function at distances $r\gg a$ is established from Eqs.
(7) and (9). Introducing a wave vector

$k_{c}=\left(\frac{2m\varepsilon_{0}}{\hbar^{2}}\right)^{1/2}$ (13)

and the dimensionless variable $z$ defined as

$z=\frac{q-k_{0}}{k_{c}},$ (14)

and performing the angular integration, we get for a nonzero component of a
spinor

$\Psi(r)\propto\int\limits_{0}^{\infty}\frac{J_{0}(k_{0}r+k_{c}rz)}{z^{2}+1}.$
(15)

Since characteristic $z$ is $\sim 1$, the localization length of $\Psi(r)$ is
$r_{c}\sim k_{c}^{-1}$. For $r<r_{c}$, the $z$-dependence in the argument of
the Bessel function can be neglected, so that $\Psi(r)\propto J_{0}(k_{0}r)$.
For $r>r_{c}$ the Bessel function can be replaced by a large-argument
asymptote: $J_{0}(z)\approx\left(\frac{2}{\pi
z}\right)^{1/2}\cos\left(z-\frac{\pi}{4}\right)$. This yields

$\Psi(r)\propto\frac{\cos\left(k_{0}r-\frac{\pi}{4}\right)}{r^{1/2}}\exp\left(-k_{c}r\right).$
(16)

We see that, unlike the conventional impurities, the wave function, in
addition to the exponential decay, contains an oscillating factor with a
period $2\pi/k_{0}$. This behavior is illustrated in the figure. In the next
section we recalculate this oscillations into the splitting of the levels of
two impurities.

##  III Two impurities

Let the impurities be located at $\pm\frac{1}{2}{\bm{\rho}}$, so that the net
potential has the form

${\tilde{U}}\left({\bm{r}}\right)=U\left({\bm{r}}-\frac{\bm{\rho}}{2}\right)+U\left({\bm{r}}+\frac{\bm{\rho}}{2}\right).$
(17)

It is straightforward to generalize Eq. (9) to the case of two impurities

$\displaystyle(2\pi)^{2}A({\bm{q}})\Bigl{(}E_{1}({\bm{q}})-E_{0}\Bigr{)}=-\int
d{\bm{r}}U({\bm{r}})e^{i{\bm{q}r}}\times$
$\displaystyle\Biggl{\\{}e^{i\frac{{\bm{q}}{\bm{\rho}}}{2}}\Bigl{(}\Psi\left(\frac{\bm{\rho}}{2}\right)\chi_{{\bm{q}}}^{(1)*}\Bigr{)}+e^{-i\frac{{\bm{q}}{\bm{\rho}}}{2}}\Bigl{(}\Psi\left(-\frac{\bm{\rho}}{2}\right)\chi_{{\bm{q}}}^{(1)*}\Bigr{)}\Biggr{\\}}.$
(18)

Expressing $A(\bm{q})$ from Eq. (III), substituting it into Eq. (7) and
setting $\bm{r}=\frac{\bm{\rho}}{2}$ and $\bm{r}=-\frac{\bm{\rho}}{2}$, we
arrive at the system of equations for $\Psi\left(\frac{\bm{\rho}}{2}\right)$
and $\Psi\left(-\frac{\bm{\rho}}{2}\right)$.

To cast this system in a concise form we take the advantage of the fact that
the distance, $\rho$, between the impurities is much bigger than $a$, so that
the splitting of the levels due to the overlap of the impurity wave functions
is smaller than $\varepsilon_{0}$ given by Eq. (11). We introduce the
following notations for the elements of the spinors
$\Psi\left(\frac{\bm{\rho}}{2}\right)$ and
$\Psi\left(-\frac{\bm{\rho}}{2}\right)$

$\Psi\left(\frac{\bm{\rho}}{2}\right)=\begin{pmatrix}a_{1}\\\
b_{1}\end{pmatrix},\hskip
8.53581pt\Psi\left(-\frac{\bm{\rho}}{2}\right)=\begin{pmatrix}a_{2}\\\
b_{2}\end{pmatrix}.$

Then the generalization of Eq. (10) to the case of two impurities takes the
form

$\displaystyle\kappa
a_{1}=\bigl{[}J_{0}(k_{0}\rho)a_{2}-iJ_{1}(k_{0}\rho)b_{2}\bigr{]}\exp\left(-k_{c}\rho\right),$
$\displaystyle\kappa
b_{1}=\bigl{[}iJ_{1}(k_{0}\rho)a_{2}+J_{0}(k_{0}\rho)b_{2}\bigr{]}\exp\left(-k_{c}\rho\right),$
$\displaystyle\kappa
a_{2}=\bigl{[}J_{0}(k_{0}\rho)a_{1}-iJ_{1}(k_{0}\rho)b_{1}\bigr{]}\exp\left(-k_{c}\rho\right),$
$\displaystyle\kappa
b_{2}=\bigl{[}iJ_{1}(k_{0}\rho)a_{1}+J_{0}(k_{0}\rho)b_{1}\bigr{]}\exp\left(-k_{c}\rho\right),$
(19)

where the parameter $\kappa$ is related to the energy, $\varepsilon$, and to
the single-impurity binding energy $\varepsilon_{0}$ as

$\kappa=\Bigl{(}\frac{\varepsilon}{\varepsilon_{0}}\Bigr{)}^{1/2}-1.$ (20)

The four solutions of the system Eq. (III) can be easily found, namely

$\displaystyle\kappa^{+}(\rho)=\pm\bigl{[}J_{0}(k_{0}\rho)+J_{1}(k_{0}\rho)\bigr{]}\exp\left(-k_{c}\rho\right),$
$\displaystyle\kappa^{-}(\rho)=\pm\bigl{[}J_{0}(k_{0}\rho)-J_{1}(k_{0}\rho)\bigr{]}\exp\left(-k_{c}\rho\right).$
(21)

Concerning the eigenvectors, their structure can be illustrated e.g. for
positive $\kappa^{+}$ when it is symmetric, namely, $a_{2}=a_{1}$ and
$b_{2}=b_{1}=ia_{1}$.

Eq. (III) allows to trace how the overlap-induced splitting of single-impurity
levels evolves with the increase of the distance between the impurities. For
$k_{0}\rho\ll 1$ the Bessel function $J_{0}(k_{0}\rho)$ should be replaced by
$1$, while $J_{1}(k_{0}\rho)$ should be replaced by $0$. Then the splitting is
of 1D-type and does not depend on the spin-orbit parameter, $k_{0}$. In the
opposite limit $k_{0}\rho\gg 1$, substituting the large-argument asymptotes of
the Bessel functions into Eq. (III) yields

$\displaystyle\kappa^{+}(\rho)=\pm\frac{2\sin(k_{0}\rho)}{(k_{0}\rho)^{1/2}}\exp\left(-k_{c}\rho\right),$
$\displaystyle\kappa^{-}(\rho)=\pm\frac{2\cos(k_{0}\rho)}{(k_{0}\rho)^{1/2}}\exp\left(-k_{c}\rho\right).$
(22)

The $1/\rho^{1/2}$ behavior of the prefactor confirms that, at large $\rho$,
the splitting is of the 2D-type and oscillates with distance due to the spin-
orbit coupling. Since $k_{c}\gg k_{0}$, the splitting of the impurity levels
strongly oscillates with distance, $\rho$, between the impurities.

##  IV Density of states

Assume now that the density, $n$, of the short-range impurities is finite. As
a result of the overlap of single-impurity wave functions, a level with
binding energy $\varepsilon_{0}$ gets smeared into a band.

A basic assumption of the conventional Lifshitz model,Lifshitz  is that the
width of the band is much smaller than $\varepsilon_{0}$. This amounts to the
smallness of the parameter $nk_{c}^{-2}\ll 1$. For such low densities the
overlap-induced splitting Eq. (III) of two levels at a typical distance
$\rho\sim n^{-1/2}$ contains an exponential factor,
$\exp\left[-k_{c}n^{-1/2}\right]$, which is small by virtue of this parameter.

A simplified argument which allows to establish the shape of the density of
states, $g(\varepsilon)$, in the vicinity of $\varepsilon=\varepsilon_{0}$
goes as follows.book  A shift of a given single-impurity level will be
anomalously small, if it has no neighboring impurities within a circle of a
certain radius, $\rho_{\varepsilon}$. This radius is anomalously large
compared to the typical distance between the impurities, i.e.
$\rho_{\varepsilon}\gg n^{-1/2}$. With sub-logarithmic accuracy, the condition
for $\rho_{\varepsilon}$ reads

$\exp\left(-k_{c}\rho_{\varepsilon}\right)=\frac{|\varepsilon-\varepsilon_{0}|}{\varepsilon_{0}},$
(23)

where the right-hand side is $\kappa$ defined by Eq. (20) and taken at
$|\varepsilon-\varepsilon_{0}|\ll\varepsilon_{0}$. Probability that the
neighboring impurities are absent in the circle with a radius
$\rho_{\varepsilon}$ is equal to $\exp\left(-\pi
n\rho_{\varepsilon}^{2}\right)$. Substitution of $\rho_{\varepsilon}$ found
from Eq. (23) into this probability, yields the exponent in the density of
states

$g(\varepsilon)\propto\exp\Bigg{[}-\frac{\pi
n}{k_{c}^{2}}\ln^{2}\Big{(}\frac{\varepsilon_{0}}{|\varepsilon-\varepsilon_{0}|}\Big{)}\Bigg{]}.$
(24)

Since the ratio $\frac{n}{k_{c}^{2}}$ is small, Eq. (24) describes a sharp
minimum in a parametrically narrow domain

$|\varepsilon-\varepsilon_{0}|\sim\varepsilon_{0}\exp\Big{[}-\frac{k_{c}}{(\pi
n)^{1/2}}\Big{]}.$ (25)

To incorporate the spin-orbit coupling, we use the modified splitting Eq.
(III) and present the density of states as a sum
$g(\varepsilon)=g^{+}(\varepsilon)+g^{-}(\varepsilon)$, where $g^{+}$ and
$g^{-}$ are defined as

$\displaystyle
g^{+}(\varepsilon)=\int\limits_{0}^{\infty}d{\rho}F(\rho)\delta\Bigg{[}\frac{|\varepsilon-\varepsilon_{0}|}{2\varepsilon_{0}}-\frac{4|\sin
k_{0}\rho|}{(k_{0}\rho)^{1/2}}e^{-(k_{c}\rho)}\Bigg{]},$ (26) $\displaystyle
g^{-}(\varepsilon)=\int\limits_{0}^{\infty}d{\rho}F(\rho)\delta\Bigg{[}\frac{|\varepsilon-\varepsilon_{0}|}{2\varepsilon_{0}}-\frac{4|\cos
k_{0}\rho|}{(k_{0}\rho)^{1/2}}e^{-(k_{c}\rho)}\Bigg{]}.$ (27)

Here $F(\rho)$ is the nearest neighbor distribution

$F(\rho)=2\pi n\rho\exp\left(-\pi n\rho^{2}\right).$ (28)

Our key point is that, due to the spin-orbit coupling, the level splitting
defined by $\kappa^{+}$ turns to zero at certain distances,
$\rho_{N}=\frac{\pi N}{k_{0}}$, while $\kappa^{-}$ turns to zero at
$\rho_{N+\frac{1}{2}}=\frac{\pi(N+\frac{1}{2})}{k_{0}}$. Substituting

$\rho=\rho_{N}+\delta\rho_{N}$ (29)

into the first equation of Eq. (26), we present $g^{+}(\varepsilon)$ as a sum

$\displaystyle g^{+}(\varepsilon)=\sum_{N}\frac{(\pi
N)^{1/2}e^{k_{c}\rho_{N}}}{4}F(\rho_{N})$ $\displaystyle\times\int
d\delta\rho_{N}\delta\Bigg{[}\frac{|\varepsilon-\varepsilon_{0}|}{8\varepsilon_{0}}(\pi
N)^{1/2}e^{k_{c}\rho_{N}}-|\sin(k_{0}\delta\rho_{N})|\Bigg{]}.$ (30)

In deriving Eq. (IV) we made two assumptions: $\delta\rho_{N}\ll\rho_{N}$ and
$k_{c}\delta\rho_{N}\ll 1$. As it is seen from Eq. (IV), characteristic
$\delta\rho_{N}$ is $\sim k_{0}^{-1}$, so that the first assumption is valid
for $N\gg 1$, while the second condition is ensured by the relation $k_{c}\ll
k_{0}$, which was assumed above.

The ”memory” about the spin-orbit coupling in Eq. (IV) is encoded in the
combination $\sin(k_{0}\delta\rho_{N})$. Taking the limit $k_{0}\rightarrow 0$
and replacing the summation over $N$ by integration, one recovers the result
Eq. (24).

To reveal the role of spin-orbit coupling we perform the integration in Eq.
(IV) with a help of the $\delta$-function, we find

$g^{+}(\varepsilon)=\frac{1}{4}\sum_{N}\frac{F(\rho_{N})}{\Big{[}\frac{\exp{(-2k_{c}\rho_{N})}}{\pi
N}-\Big{(}\frac{\varepsilon-\varepsilon_{0}}{8\varepsilon_{0}}\Big{)}^{2}\Big{]}^{1/2}}.$
(31)

The expression for $g^{-}(\varepsilon)$ has a similar form with replacement of
$N$ by $N+\frac{1}{2}$. Note that each term in Eq. (31) exhibits a square-root
singularity near energies

$\varepsilon_{N}^{\pm}=\varepsilon_{0}\Bigg{[}1\pm\frac{8}{(\pi
N)^{1/2}}\exp{\Big{(}-\frac{2\pi k_{c}}{k_{0}}N\Big{)}}\Bigg{]}.$ (32)

On the other hand, by virtue of the small parameter $\frac{k_{c}}{k_{0}}$, the
intervals $\\{\varepsilon_{N}^{+},\varepsilon_{N}^{-}\\}$ for different $N$
overlap. Thus, the overall shape of the density of states is smooth. To
uncover the role of the discreteness of $N$, we transform Eq. (31) using the
Poisson summation

$\sum_{N}S(N)=\int dxS(x)+2\sum_{l}\int dxS(x)\cos(2\pi lx).$ (33)

In our case, the function $S(x)$ has a form

$S(x)=\frac{1}{4}\frac{F\left(\frac{\pi
x}{k_{0}}\right)}{\Big{[}\frac{\exp{\left(-\frac{2\pi
k_{c}x}{k_{0}}\right)}}{\pi
x}-\Big{(}\frac{\varepsilon-\varepsilon_{0}}{8\varepsilon_{0}}\Big{)}^{2}\Big{]}^{1/2}}.$
(34)

The denominator in Eq. (34) turns to zero at

$x=x_{\varepsilon}\approx\frac{k_{0}}{2\pi
k_{c}}\ln\Bigl{(}\frac{8\varepsilon_{0}}{\varepsilon-\varepsilon_{0}}\Bigr{)}^{2}.$
(35)

The first term in Eq. (33) reproduces the standard result Eq. (24) for the
density of states in the absence of spin-orbit coupling. To evaluate the terms
with $l\geq 1$ one should expand the integrand around $x_{\varepsilon}$. This
leads to the oscillating spin-orbit component in the density of states
$\propto\cos(2\pi lx_{\varepsilon})$.

##  V Concluding remarks

The prime observation of the present paper is that, in the presence of spin-
orbit coupling, the pairs of impurities located at certain distances from each
other do not hybridize their levels. This leads to singularities in the in the
density of states at certain energies.

There is a certain similarity between the impurity state in a 2D electron gas
with strong spin-orbit coupling and the in-gap state created by a magnetic
impurity in a superconductor.Yu ; Shiba ; Rusinov ; Glazman  Similarly to Eq.
(16), the wave function of the in-gap state oscillates with distance as $\cos
k_{F}r$, where $k_{F}$ is the Fermi momentum. However, the limit of low
density of magnetic impurities considered in the present paper is not relevant
for superconductors. Rather, in the relevant regime, there are many impurities
within the localization radius determined by the coherence length. In this
limit, the self-consistent Born approximation applies, and the two-peak
structure of the density of the impurity states transforms into a
semicircle.Glazman

## References

  * (1) I. M. Lifshitz, ”Energy spectrum structure and quantum states of disordered condensed systems,” Usp. Fiz. Nauk 83, 617 (1964) [English transl.: Sov. Phys.-Usp. 7, 549 (1964-1965) ].
  * (2) I. M. Lifshitz, S. A. Gredeskul, and L. A. Pastur, ”Introduction to the Theory of Disordered Systems,” (Willey, New York, 1988).
  * (3) B. I. Shklovskii and A. L. Efros, ”Electronic Properties of Doped Semiconductors,” (Springer, Berlin, 1984).
  * (4) Yu. A. Bychkov and E. I. Rashba, ”Properties of a 2D electron gas with lifted spectral degeneracy,” Pis’ma Zh. Eksp. Teor. Fiz. 39, 64 (1984) [JETP Lett. 39, 78 (1984)].
  * (5) A. G. Galstyan and M. E. Raikh, ”Disorder-induced broadening of the density of states for two-dimensional electrons with strong spin-orbit coupling,” Phys. Rev. B 58, 6736 (1998).
  * (6) A. V. Chaplik and L. I. Magarill, ”Bound States in a Two-Dimensional Short Range Potential Induced by the Spin-Orbit Interaction,” Phys. Rev. Lett. 96, 126402 (2006).
  * (7) V. V. Mkhitaryan and M. E. Raikh, ”Disorder-induced tail states in gapped bilayer graphene,” Phys. Rev. B 78, 195409 (2008).
  * (8) B. Skinner, B. I. Shklovskii, and M. B. Voloshin, ”Bound state energy of a Coulomb impurity in gapped bilayer graphene,” Phys. Rev. B 89, 041405(R) (2014).
  * (9) J. Hutchinson and J. Maciejko, ”Unconventional transport in low-density two-dimensional Rashba systems,” Phys. Rev. B 98, 195305 (2018).
  * (10) L. Yu, ”Bound state in superconductors with paramagnetic impurities,” Acta Phys. Sin. 21, 75 (1965).
  * (11) H. Shiba, ”Classical Spins in Superconductors,” Progress of Theoretical Physics, 40, 435 (1968).
  * (12) A. I. Rusinov, ”Superconductivity near a paramagnetic impurity,” Zh. Eksp. Teor. Fiz. 56, 2047 (1969) [Sov. Phys. JETP 29, 1101 (1969)].
  * (13) Ya. V. Fominov, M. Houzet, and L. I. Glazman, ”Surface impedance of superconductors with weak magnetic impurities,” Phys. Rev. B 84, 224517 (2011).