# A general approach to maximise information density in neutron reflectometry
analysis

Andrew R. McCluskey<EMAIL_ADDRESS><EMAIL_ADDRESS>Diamond Light Source, Rutherford Appleton Laboratory, Harwell Science and
Innovation Campus, Didcot, OX11 0DE, UK  Department of Chemistry, University
of Bath, Claverton Down, Bath BA2 7AY, UK  Joshaniel F. K. Cooper  ISIS Pulsed
Neutron and Muon Source, Science and Technology Facilities Council, Rutherford
Appleton Laboratory, Harwell Science and Innovation Campus, Didcot,
Oxfordshire, OX11 0QX, UK  Tom Arnold  European Spallation Source, SE-211 00
Lund, Sweden  Department of Chemistry, University of Bath, Claverton Down,
Bath BA2 7AY, UK  ISIS Pulsed Neutron and Muon Source, Science and Technology
Facilities Council, Rutherford Appleton Laboratory, Harwell Science and
Innovation Campus, Didcot, Oxfordshire, OX11 0QX, UK  Tim Snow  Diamond Light
Source, Rutherford Appleton Laboratory, Harwell Science and Innovation Campus,
Didcot, OX11 0DE, UK  School of Chemistry, University of Bristol, Bristol, BS8
1TS, UK

###### Abstract

Neutron and X-ray reflectometry are powerful techniques facilitating the study
of the structure of interfacial materials. The analysis of these techniques is
ill-posed in nature requiring the application of model-dependent methods. This
can lead to the over- and under- analysis of experimental data when too many
or too few parameters are allowed to vary in the model. In this work, we
outline a robust and generic framework for the determination of the set of
free parameters that are capable of maximising the in-formation density of the
model. This framework involves the determination of the Bayesian evidence for
each permutation of free parameters; and is applied to a simple phospholipid
monolayer. We believe this framework should become an important component in
reflectometry data analysis, and hope others more regularly consider the
relative evidence for their analytical models.

††preprint: APS/123-QED

##  I  Introduction

Neutron and x-ray reflectometry methods offer an important insight into the
structure of interfacial materials; such as surfactant systems McCluskey  _et
al._ (2019a), biological membranes Clifton  _et al._ (2019), and polymeric
photovoltaics Pérez  _et al._ (2019). However, the analysis of reflectometry
data is a non-trivial task, due in part to the phase problem. This is where a
direct inversion of the reflectometry profile, to recover the scattering
length density, $\rho(z)$ profile is not possible since the phase of the
reflected radiation is lost Majkrzak and Berk (1995a). While there has been
much work to retrieve the phase information from a reflectometry profile Kirby
_et al._ (2012); from using reference layers Majkrzak and Berk (1995a); de
Haan  _et al._ (1995); Majkrzak  _et al._ (2003); Nikova  _et al._ (2019),
polarized neutrons Leeb  _et al._ (1998), a dwell time approach Fiedeldey  _et
al._ (1992), and from the investigation of multiple contrasts Majkrzak and
Berk (1995b); Majkrzak  _et al._ (2000, 2003); Koutsioubas (2019), these
methods are non-trivial and often require a specific experimental setup.

Where it is not possible to retrieve the phases, a model-dependent analysis
method can be used to determine the structure of the experimental system
Pedersen and Hamley (1994); de Haan and Drijkoningen (1994); Nelson (2006);
van der Lee  _et al._ (2007); Gerelli (2016a, b); Nelson and Prescott (2019).
This is where a model, typically a series of contiguous layers, is created and
the model reflectometry is determined using the Abelès matrix formalism for
stratified media Abelès (1948) or Parratt recursive method Parratt (1954).
However, the use of model-dependent methods is ill-posed in nature, where the
solution may not necessarily be unique Klibanov and Sacks (1992, 1994). This
necessitates the integration of _a priori_ information about the underlying
physics and chemistry of the experimental system, such as assumptions about
molecular volume Waldie  _et al._ (2018); Campbell  _et al._ (2018); McCluskey
_et al._ (2019a). This is typically achieved through the re-parameterisation
of the model, to describe the system in terms of physiochemical descriptors
from which a set of layers can be determined Schalke and Lösche (2000);
Heinrich and Lösche (2014); Heinrich (2016). The inclusion of this prior
information acts to reduce the dimensionality of the optimisation space and
therefore reduces the number of possible structures that would offer
satisfactory agreement with the data.

The inclusion of prior information is often achieved by taking literature
values from other techniques and applying them inflexibly to constrain values
in the re-parameterised model, such as molecular volumes from computational
simulation or solution scattering Sun  _et al._ (1994); Armen  _et al._
(1998). However, as was shown in the work of Campbell _et al._ Campbell  _et
al._ (2018) and McCluskey _et al._ McCluskey  _et al._ (2019a), these values
are not always strictly representative when the physical or chemical
environment is varied. This has lead to the use of Markov-chain Monte Carlo
(MCMC) methods for the analysis of reflectometry data, resulting in the
inclusion of these methods in many common reflectometry analysis packages
Hughes (2017); Kienzle  _et al._ (2017); Nelson and Prescott (2019). The
inclusion of these methods enables the prior probabilities for a series of
variables to be outlined, which are considered in a Bayesian fashion during
the sampling of the parameter space Bayes (1763); McCluskey  _et al._ (2019a).
This allows for more rational inclusion of the prior information about the
experimental system in the analysis, in addition to, the ability to obtain
information about the sampled posterior probabilities that are allowed from
the experimental data.

The ability to include ranges of values, in the form of uninformative uniform
prior probabilities, or normally distributed priors centred on a literature
value in the analysis of an experimental dataset is an important development.
However, “with great power, there must also come – great responsibility!” Lee
and Ditko (1962) and the ability to sample or fit all of the parameters in an
experimental analysis can lead to a high risk of over-fitting Mayer  _et al._
(2010). As a result, it is necessary to investigate methods to quantify the
information that is available from a given experimental dataset and to enable
the comparison of different analytical models; including the number of free
parameters that may be considered. A specific approach to this problem is
outlined in early works from Sivia and co-workers Sivia  _et al._ (1991);
Geoghagan  _et al._ (1996); Sivia and Webster (1998), where a two-parameter
model was compared with a three-parameter model; with the additional parameter
being related to the scaling of the model reflectometry. Sivia and co-workers
assumed that the likelihood function for both of these model would have “only
one significant maximum”, allowing numerical simplifications to be performed
and therefore the problem to be solved analytically. This may not be possible
for all models investigated by reflectometry, due again to the ill-posed
nature of the analysis problem. Furthermore, the work of Hughes and co-workers
Hughes  _et al._ (2019) more recently applied the use of Bayesian model
comparison to determine the mixing of lipids in a model bacterial membrane.

In this work, we outline a general Bayesian framework for the comparison of
different analytical models for the rationalisation of reflectometry data.
This method can allow the maximum information to be obtained from the analysis
of a given dataset, without risking over-fitting. In the supplementary
information for this paper, we provide detailed code showing how this
framework may be applied to a problem leveraging exclusively open-source
software packages.

##  II  Methods

###  II.1  Bayesian model selection framework

Bayesian model selection involves the comparison of the posterior
probabilities $p(H|\mathbf{D})$ 111In this work, a reduced notation has been
used, where the background information parameter, $I$ is implied., for two
models, $H_{x}$ and $H_{y}$, given some dataset, $\bm{D}$. The posterior
probability for a given model, $H$, can be determined from Bayes theorem Bayes
(1763),

$p(H_{x}|\bm{D})=\frac{p(\bm{D}|H_{x})p(H_{x})}{p(\bm{D})},$ (1)

where, $p(\bm{D}|H)$ is the evidence for the model given the data, $p(H)$ is
the prior probability for the model and $p(\bm{D})$ is the probability
associated with the measured data, which is equal to the sum of the evidences
for every potential model (which is potentially large if not infinite),

$p(\bm{D})=\sum_{x=1}{p(\bm{D}|H_{x})}.$ (2)

Therefore, in order to describe the evidence as a truely conditional
probability, it is necessary to evaluate very possible model that maybe
applied to the experimental data. This is not feasible, therefore in this work
we will use the notation, $p(\bm{D}|H_{x})$ for the evidence of a given model,
however this is not a strict conditional probability. When comparing the
models, concerning the same data, the following relationship can be applied
Pullen and Morris (2014),

$\frac{p(H_{x}|\bm{D})}{p(H_{y}|\bm{D})}=\frac{p(\bm{D}|H_{x})}{p(\bm{D}|H_{y})}\times\frac{p(H_{x})}{p(H_{y})},$
(3)

where, the probabilities associated with the data only cancel out, and the
relative posterior probabilities depend only on the evidence and prior
probabilities for each model. Throughout this work, we have used the
uninformative assumption that the prior probability for all models is $1$.
Therefore, it is only the evidence for each model, for the given data, that
influences the posterior probability.

Nested sampling is a method for Bayesian evidence determination that was
developed by Skilling Skilling (2006) and is implemented in the Python package
dynesty Speagle (2019). This sampling method allows for the evidence
determination through the estimation of the, possibly multi-dimensional,
integral under the likelihood, $\mathcal{L}(\bm{X}|H)$, and parameter specific
prior $p(\bm{X}|H)$ Sivia and Skilling (2005),

$p(\mathbf{D}|H)=\iint_{\bm{R}}{\mathcal{L}(\bm{X}|H)p(\bm{X}|H)}\;\text{d}^{M}\bm{X},$
(4)

where, $\bm{X}$ is a vector of length $M$ of the varying parameters in the
model, $\bm{R}$ is a $2\times M$ matrix that describes the range over which
the integral should be evaluated and $p(\bm{X}|H)$ is the parameter specific
prior. Nested sampling involves “intelligently” sampling this integral until
some stopping criterion is achieved. This stopping criterion is when the
following is less than a given value,

$\displaystyle\Delta\ln\\{\hat{p}(\bm{D}|H)_{i}\\}=\ln\\{\hat{p}(\bm{D}|H)_{i}-\Delta\hat{p}(\bm{D}|H)_{i}\\}$
(5) $\displaystyle-\ln\hat{p}(\bm{D}|H)_{i},$

where, $\hat{p}(\bm{D}|H)_{i}$ is the current estimated evidence and
$\Delta\hat{p}(\bm{D}|H)_{i}$ can be estimated from the prior volume of the
previous point in the sampling. In this work, a stopping criteria of
$\Delta\ln\\{\hat{p}(\bm{D}|H)\\}=0.5$ was used.

The likelihood is the difference between the data and model, and is defined in
this work as it is implemented in the refnx reflectometry package Nelson and
Prescott (2019); Nelson  _et al._ (2019),

$\displaystyle\ln\mathcal{L}=-\frac{1}{2}\sum_{i=1}^{N}\Bigg{(}\bigg{\\{}\frac{R(q_{i})-R(q_{i})_{m}}{\delta
R(q_{i})}\bigg{\\}}^{2}$ (6) $\displaystyle+\ln[2\pi\\{\delta
R(q_{i})\\}^{2}]\Bigg{)},$

where, $N$ is the number of measured $q$-vectors, $q_{i}$ is the $i$th
$q$-vector, $R(q_{i})$ is the experimental reflected intensity at the $i$th
$q$-vector, $\delta R(q_{i})$ is the uncertainty in the experimental reflected
intensity at the $i$th $q$-vector, and $R(q_{i})_{m}$ is the model reflected
intensity at the $i$th $q$-vector determined from the implementation of the
Abelès matrix formalism for stratified media Abelès (1948); Parratt (1954) in
refnx.

The evidence allows for the comparison of two models, where the data are the
same. It is only possible to compare a given pair of models for the same data,
as in Equation 3 the data-only dependent probabilities were cancelled out as
they were equal. The Bayes factor, $B_{xy}$, is the name given to the ratio
between the evidence for two models,

$\displaystyle\ln{B_{xy}}$
$\displaystyle=\ln{\Bigg{[}\frac{p(\bm{D}|H_{x})}{p(\bm{D}|H_{y})}\Bigg{]}}$
(7) $\displaystyle=\ln[p(\bm{D}|H_{x})]-\ln[p(\bm{D}|H_{y})].$

The interpretation of the Bayes factor was outlined by Kass and Raftery Kass
and Raftery (1995), which is given in Table 1.

Table 1:  Interpretation of the Bayes factor, $B_{xy}$ between models $x$ and $y$. $2\ln{B_{xy}}$ | Evidence against $H_{y}$
---|---
$[0,2)$ | Not worth more than a bare mention
$[2,6)$ | Positive
$[6,10)$ | Strong
$[10,\infty)$ | Very strong

This method is easily generalised to any model-dependent experiment, and has
been applied in areas such as systems biology, astronomy, computational
chemistry, and indeed, in a less general fashion than the current work, to
neutron reflectometry Pullen and Morris (2014); Cornish and Littenberg (2007);
Ensign and Pande (2010); Sivia and Webster (1998).

The Bayesian approach allows for the comparison of models where there are
different numbers of free parameters, without the risk of overfitting and
“fitting an elephant” Mayer  _et al._ (2010). This is due to the Bayesian
evidence being derived from an integral in parameter space, and therefore
scaling with the number of parameters. This means that the addition of an
additional free parameter to the model must offer a significant improvement to
the likelihood determined. However, it is important also to note that the
accuracy of the determined evidence depends on the prior probabilities chosen
for each of the free parameters and therefore care should be taken to ensure
that these are meaningful.

###  II.2  Neutron reflectometry measurements

The neutron reflectometry measurements investigated in this work were
previously published by Hollinshead _et al._ Hollinshead  _et al._ (2009) and
subsequently re-analysed by McCluskey _et al._ McCluskey  _et al._ (2019b).
The measurements concern the study of a monolayer of 1,2-distearoyl-_sn_
-glycero-3-phosphocholine (DSPC) at the air/water interface and were conducted
on seven different isotopic contrasts of the phospholipid and water. These
contrasts were made up from four phospholipid types; fully-hydrogenated
(h-DSPC), head-deuterated (d_13-DSPC), tail-deuterated (d_70-DSPC), and fully
deuterated (d_83-DSPC). These were paired with two water contrasts; D2O and
air-contrast matched water (ACMW, where D2O and H2O are mixed such that the
resulting scattering length density is zero). The pairing of the fully-
hydrogenated phospholipid with ACMW was not performed, most likely due to the
lack of the scattering available from such a system. Table 2 gives the
shorthands that are used in this work to refer to the different contrast
pairings investigated.

Table 2:  The shorthands used for the different combinations of phospholipid monolayer and water investigated. Shorthand | Phospholipid contrast | Water contrast
---|---|---
h-D2O | h-DSPC |  D2O
d_13-D2O |  d_13-DSPC |  D2O
d_13-ACMW |  d_13-DSPC |  ACMW
d_70-D2O |  d_70-DSPC |  D2O
d_70-ACMW |  d_70-DSPC |  ACMW
d_83-D2O |  d_83-DSPC |  D2O
d_83-ACMW |  d_83-DSPC |  ACMW

The data analysed in this work was taken at a surface pressure of
$30\text{\,}\mathrm{m}\mathrm{N}\mathrm{m}^{-1}$. Additional details of the
data collection may be found in the original publication Hollinshead  _et al._
(2009).

###  II.3  Analytical model

In the model-dependent analysis of neutron reflectometry data, it is necessary
to have a description of the experimental system to which the Abelès matrix
formalism may be applied. Campbell _et al._ investigated different models for
the study of phospholipid monolayers and found a two-layer description that
differentiates between the phospholipid heads, $h$, and tails, $t$, to be the
best at reproducing the reflected intensity Campbell  _et al._ (2018). The
head layer is adjacent to the solvent, while the tail layer is adjacent to the
air. The thickness of the tail layer is defined as $d_{t}$ and the head layer
is $d_{h}$. At each of the three interfaces, there is some interfacial
roughness, $\sigma$, which is modelled with an error function Névot and Croce
(1980). In this work, the roughness was taken to be conformal, which is to say
that it does not vary between interfaces, in agreement with the work of
Campbell _et al._ Campbell  _et al._ (2018). This final assumption is
reasonable given the fact that this system consists of a monolayer of a single
lipid type, however, this may not be reasonable in more complex systems.

The scattering length density, $\rho_{i}$, for a layer $i$, from which the
model reflectometry profile is obtained, is calculated from the volume
fraction of the chemical component, $\phi_{i}$, the scattering length of the
component, $b_{i}$ (found in Table 3), and its molecular volume, $V_{i}$,

$\rho_{i}=\frac{b_{i}}{V_{i}}\phi_{i}+\rho_{\ce{H2O}}(1-\phi_{i}).$ (8)

where, $\rho_{\ce{H2O}}$ is the scattering length density of water, which is
$0\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-2}$ for ACMW and
$6.35\text{\times}{10}^{-6}\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-2}$ for
D2O. To preserve the chemical bonding between the heads and tails of the
phospholipid, the following constraint is included in the model,

$\phi_{h}=\frac{d_{t}V_{h}\phi_{t}}{V_{t}d_{h}}.$ (9)

This ensures that the number of head groups is equal to the number of pairs of
tail groups.

Table 3:  The scattering lengths, $b_{i}$, for each of the DSPC phospholipid components in this work; where the subscript indicates either the head, $h$, or tail, $t$. Phospholipid contrast |  $b_{h}$/$\mathrm{fm}$ |  $b_{t}$/$\mathrm{fm}$
---|---|---
h-DSPC | $6.01$ | $-3.58$
d_13-DSPC | $19.54$ | $-3.58$
d_70-DSPC | $11.21$ | $69.32$
d_83-DSPC | $24.75$ | $69.32$

The model outlined above has up to six possible free parameters; the two-layer
thicknesses, $d_{t}$ and $d_{h}$, the two molecular volumes, $V_{t}$ and
$V_{h}$, the volume fraction of the tail layer, $\phi_{t}$, and the conformal
interfacial roughness, $\sigma$. In this work, the evidence as a function of
these parameters is of interest as it is possible to use the evidence and
Bayes factor to determine the best set of free parameters for a given set of
data. Therefore, in Table 4 the uniform prior probability range and value when
constrained is given for each of the parameters.

Table 4:  The constrained and uniform prior probability range for each of the free parameters in the analytical model. Parameter | Constrained Value | Prior Range
---|---|---
$d_{h}$/$\mathrm{\SIUnitSymbolAngstrom}$ |  $10.0$222The work of Campbell _et al._ uses a head layer thickness of $10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ for the phosphatidylcholine head layer based on “molecular dimensions” Campbell  _et al._ (2018). | $[8.0,16.0)$
$V_{h}$/$\mathrm{\SIUnitSymbolAngstrom}^{3}$ |  $339.5$333From the work of Nučerka _et al._ Kučerka  _et al._ (2004) and Balgavý Balgavý  _et al._ (2001). | $[300.0,380.0)$
$d_{t}$/$\mathrm{\SIUnitSymbolAngstrom}$ |  $21.0$444The maximum tail length for DSPC is found as $24.3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ from the Tanford equation Tanford (1980), the value of $21\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ comes from assuming a chain tilt of $30\text{\,}\mathrm{\SIUnitSymbolDegree}$. | $[10.0,26.0)$
$\phi_{t}$ |  $1.0$555Based on the hydrophobic interaction of the tail groups, it is unlikely that water would intercalate. |  $[0.5,1.0)$666A maximum of $1$ was set as a volume fraction greater than $1$ would be unphysical.
$V_{t}$/$\mathrm{\SIUnitSymbolAngstrom}^{3}$ |  $850.4$777From the previous analysis of this dataset McCluskey  _et al._ (2019b), which agrees with the expected compression of $\sim 12\text{\,}\mathrm{\char 37\relax}$ discussed by Campbell _et al._ Campbell  _et al._ (2018). | $[800.0,1000.0)$
$\sigma$/$\mathrm{\SIUnitSymbolAngstrom}$ |  $2.9$888From the work of Campbell _et al._ Campbell  _et al._ (2018). |  $[2.9,5.0)$999The minimum of $2.9\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ is based on that of water Braslau  _et al._ (1985); Sinha  _et al._ (1988), and the fact that surfactants have been shown to increase surface roughness Tikhonov  _et al._ (2000).

These six parameters can give rise to $63$ permutations where there is at
least one free parameter. These permutations constitute every possible
combination of the six parameters, either held constant or being allowed to
varies within the prior probability. This allows for a complete understanding
of which combinations of constrained and free parameters give rise to the
highest evidence, and therefore represents the maximum information density
that may be obtained from the given data. We note again, that both the chosen
prior and the constrained values are not infallible and will influence the
result.

All seven isotopic neutron contrasts were analysed with a single structural
model. In doing so we make the common assumption that there is no effect on
the structure as a result of the contrast variation method Nelson (2006).
Additionally, the reflected intensity from the model was scaled and a uniform
background was added based on values obtained in the previous analysis of this
dataset McCluskey  _et al._ (2019b). We accept that these parameters may
influence the resulting evidence, as was shown in the work of Sivia and
Webster Sivia and Webster (1998).

##  III  Results & Discussion

###  III.1 All parameters

Figure 1 shows the variation of the estimated evidence as a function of
iterations, for the model where all six parameters were unconstrained.

Figure 1:  The estimation of evidence as a function of number of samples; for
the model where all parameters were unconstrained.

This plot shows clearly the ability for the nested sampling method to
determine a good estimate of the evidence, with the plot plateauing at
$6859.9\pm 0.3$ after more than $10\,000$ samples. This was the slowest model
to converge, as it was sampling the high-dimensionality parameter space, on a
workstation computer this nested sampling process took around
$20\text{\,}\mathrm{min}$ to converge leading to the results provided here.

In addition to the ability for the nested sampling method to determine the
evidence for a given model, it also allows for the posterior probabilities of
the analysis to be evaluated. Figure 2 gives the posterior probabilities and
cross-correlations between each of the six parameters.

Figure 2:  The posterior probabilities and correlation plots for each of the
six variables in the model.

It is clear that the volume fraction of the tail layer, $\phi_{t}$, and the
interfacial roughness, $\sigma$, both are not complete represented by the
defined prior probabilities. Hence, there lacks a full Gaussian distribution
and the extremely skewed probability distribution. However, maximum for
$\phi_{t}$ and minimum for $\sigma$ both come from physical limitations on the
system; the former being that the volume fraction cannot be greater than 1,
and the latter being that a surfactant monolayer will increase the surface
roughness compared to water Tikhonov  _et al._ (2000), which is $\sim
2.9\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ Braslau  _et al._ (1985); Sinha
_et al._ (1988). The skew-ness of $p(\phi_{t})$ is to be expected given the
hydrophobic interaction between the phospholipid tails which create the
situation where it is highly unlikely for water to be intercalated into the
tail structure. The nature of $p(\sigma)$ is similar in that the expected
value of $2.9\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ is the minimum value
available to the distribution. However, it may be the case that by fitting all
possible parameters, the data is being overfitted. Luckily, the Bayesian model
comparison method allows for the evidence comparison to fit the model which
offers the maximum information to be obtained, without over-fitting to the
data.

###  III.2 Greatest evidence

The evidence for all of the $63$ possible permutations of the free parameters
was determined. This allows for the comparison of these to determine which
model (set of varying parameters) offers the most information from the
dataset, without the over-fitting to the data that may be occurring when all
six parameters are varied. This comparison plotted in Figure 3 101010These
values are tabulated in Table S1..

Figure 3:  A comparison of the evidence for each of the $63$ permutations of
free parameters available to the analytical model. The dashed line indicated
the maximum evidence observed.

The dotted black line indicates the set of free parameters that offer the
greatest evidence for the co-refinement of the seven contrasts of neutron
reflectometry data, which in this case is $d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ ,
which has an evidence of $6872.1\pm 0.2$ . The set of parameters with the next
greatest evidence is $d_{h}$/$d_{t}$/$V_{t}$ , however, the evidence for this
set of parameters is $6866.8\pm 0.2$ . The interpretation of the Bayes factor
between these two sets of parameters gives $2\ln{B_{xy}}=$$10.7\pm 0.7$ where,
$d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ is $x$ and $d_{h}$/$d_{t}$/$V_{t}$ is $y$.
This indicates that there strong or very strong evidence for the more complex
model, and that the maximum information density that can be obtained from the
analysis at hand, without over-fitting, is the fitting of
$d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ .

Figure 4 shows the posterior probability distributions that were obtained from
the $d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ set of parameters.

Figure 4:  The posterior probabilities and correlation plots when
$d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ were the variable parameters.

The median values for each of the probability distributions and their
$95\text{\,}\mathrm{\char 37\relax}$ confidence intervals are given in Table
III.2.

Table 5:  The median values and $95\text{\,}\mathrm{\char 37\relax}$
confidence intervals for each of the free parameters in the model with the
greatest evidence.

Free Parameter | Value
---|---
$d_{h}$/$\mathrm{\SIUnitSymbolAngstrom}$ |

$V_{h}$/$\mathrm{\SIUnitSymbolAngstrom}^{3}$ & $347.0^{+3.6}_{-3.6}$
$d_{t}$/$\mathrm{\SIUnitSymbolAngstrom}$ $17.6^{+0.1}_{-0.1}$
$V_{t}$/$\mathrm{\SIUnitSymbolAngstrom}^{3}$ $828.0^{+4.3}_{-4.2}$ The value
found here for the phosphatidylcholine is slightly larger than the
$10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ used by Campbell _et al._ Campbell
_et al._ (2018) and the $10.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ quoted
by Li _et al._ Li  _et al._ (1998). The thickness of the tail layer that is
observed is less than that suggested previously Campbell  _et al._ (2018).
However, it agrees well with the results from the original analysis of this
data Hollinshead  _et al._ (2009), which itself was consistent with previous
x-ray reflectometry measurements DeWolf  _et al._ (1999); Brezesinski  _et
al._ (2001); Struth  _et al._ (2001). The tail molecular volume agrees well
with the compression of $\sim 12\text{\,}\mathrm{\char 37\relax}$, when
compared to the values from other experimental techniques Sun  _et al._
(1994); Armen  _et al._ (1998), suggested by Campbell _et al._.

The reflectometry profiles for the median values is compared with the
experimental data in Figure 5.

Figure 5:  The experimental (coloured lines) reflectometry and the median
values for the model with greatest evidence (black lines). The different
contrasts are offset be an order of magnitude in reflected intensity;
d_13-ACMW (blue), d_13-D2O (orange), h-D2O (green), d_70-ACMW (red), d_70-D2O
(purple), d_83-ACMW (brown), and d_83-D2O (pink).

There is clear agreement between the experiment and the model across all of
the contrasts. This indicates that it is likely that there is little variation
due to the isotopic contrast variation that was used, and the assumption that
is discussed in II.3 is valid.

###  III.3 Comparison with previous analysis

As mentioned above, this particular experimental dataset was previously
analysed by McCluskey _et al._ McCluskey  _et al._ (2019b). The analytical
model used in this previous work was similar to that used herein, however, the
tail volume of the phospholipid was constrained based on the assumption that
the measured area per molecule was accurate and the volume fraction of the
tail layer is $1$. That model had three free parameters, namely $d_{h}$,
$d_{t}$, and $\sigma$.

This previous analysis made no effort to evaluate different models of analysis
or numbers of free parameters. Therefore, it is not clear if the three
parameters chosen offered the maximum information density available. While not
completely comparable, as the constrained values in this work were different
to those used previously, the suggestion from this work that having
$d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ as the free parameters do not match with those
probed previously. In future, we hope that by outlining this general
methodology, others will consider the evidence of their particular model and
the free parameter being used during their analysis.

###  III.4 Best per number of parameters

Finally, we will briefly present the best set of parameters for each number of
parameters. These are shown in Figure 6, which has the characteristic shape
found in the previous work by Sivia and others Sivia  _et al._ (1991, 1993);
Sivia and Webster (1998); Sivia and Skilling (2005).

Figure 6:  A comparison of the evidence for the set of parameters with the
greatest evidence, per number of parameters. The evidence for a single
parameter has been cut off for clarity, however, it had a value of $6414.9\pm
0.1$ .

This plot can be seen to peak when there are four free parameters, before
dropping off slowly. This drop off indicates that when more parameters are
included there is an indication that over-fitting is occurring. This plot
shows clearly where the maximum information is available in the analysis for
this particular reflectometry data.

##  IV  Conclusions

This work outlines a general application of a Bayesian framework to the
problem of model selection in neutron reflectometry. Building on the work of
Sivia and co-workers in the 1990s Sivia and Webster (1998), we have utilised a
nested sampling method to generalise the methodology beyond the assumption of
only one significant maximum Sivia  _et al._ (1991). This offers the ability
to go beyond the three-parameter space previously investigated to consider
every possible parameter in an analysis problem. In doing so, it is possible
to compare the evidence of different models, where different numbers of
parameters are allowed to vary, to determine that which offers the maximum
information, without over-fitting.

To show this in action, we have investigated the analysis of seven isotopic
contrasts of neutron reflectometry data from a phospholipid monolayer at the
air/water interface. This particular analysis had a maximum of six variable
parameters, and in this work, the evidence of all $63$ possible permutations
of these parameters was determined. This allowed for the permutation that
offered the maximum evidence to be found, which under the circumstances of the
chosen priors was when $d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ were varied. This
particular combination of parameters had an evidence of $6872.1\pm 0.2$ ,
while the next best, which had more free parameters had a value of $6866.8\pm
0.2$ . This was interpreted in terms of the Bayes factor, and it was shown
that there was strong evidence for the $d_{h}$/$V_{h}$/$d_{t}$/$V_{t}$ model.
Finally, we have compared these free parameters with those used previously to
analyse the same dataset McCluskey  _et al._ (2019b). This comparison
indicates that this previous analysis did not necessarily maximise the
information density available to the data.

We accept that the relative simplicity of this model system allows for the use
of nested sampling on a relatively short timescale, with all of the
computation capable of being performed in a single day on a workstation
machine. However, the applicability of the methods discussed herein is
relevant to any model-dependent reflectometry analysis, but we note that as
the number of free parameters increases so does the dimensionality of the
parameter space being sampled, leading to longer computational timescales.
Therefore, we propose that for analyses with a large number of parameters the
analyst starts by reducing the number of free parameters as much as is
possible to ensure that the nested sampling will converge in a reasonable
timescale and gradually increasing the complexity. We hope that by developing
and sharing of this methodology, in future those analysing reflectometry data
will consider which of the free parameters in their models may be varied and
will try to determine the evidence for their particular model for other
options before completing their analysis.

## Reproducibility Statement

Electronic Supplementary Information (ESI) available: All analysis/plotting
scripts and data files allowing for a fully reproducible, and automated,
analysis workflow for the work presented is available at
https://github.com/arm61/model_select (DOI: 10.5281/zenodo.3820690) under a CC
BY-SA 4.0 license.

## Author Contributions

A.R.M. proposed, implemented, and applied the Bayesian model comparison
framework, with input from T.A., J.F.K.C. and T.S.; J.F.K.C. and T.S. sourced
funding; A.R.M. wrote the manuscript, with input from all authors.

## Acknowledgements

A.R.M. would like to acknowledge David J. Barlow and M. Jayne Lawrence for
kindly sharing the neutron reflectometry data used, and Simon Titmuss for
suggesting the use of the dynesty package. This work is supported by the Ada
Lovelace Centre – a joint initiative between the Science and Technology
Facilities Council (as part of UK Research and Innovation), Diamond Light
Source, and the UK Atomic Energy Authority.

## References

  * McCluskey  _et al._ (2019a) A. R. McCluskey, A. Sanchez-Fernandez, K. J. Edler, S. C. Parker, A. J. Jackson, R. A. Campbell, and T. Arnold, Phys. Chem, Chem. Phys. 21, 6133 (2019a).
  * Clifton  _et al._ (2019) L. A. Clifton, N. Paracini, A. V. Hughes, J. H. Lakey, N.-J. Steinke, J. F. K. Cooper, M. Gavutis, and M. W. A. Skoda, Langmuir 35, 13735 (2019).
  * Pérez  _et al._ (2019) G. E. Pérez, G. Bernardo, H. Gaspar, J. F. K. Cooper, F. Bastianini, A. J. Parnell, and A. D. F. Dunbar, ACS Appl. Mater. Interfaces 11, 13803 (2019).
  * Majkrzak and Berk (1995a) C. F. Majkrzak and N. F. Berk, Phys. Rev. B 52, 10827 (1995a).
  * Kirby  _et al._ (2012) B. J. Kirby, P. A. Kienzle, B. B. Maranville, N. F. Berk, J. Krycka, F. Heinrich, and C. F. Majkrzak, Curr. Opin. Colloid Interface Sci 17, 44 (2012).
  * de Haan  _et al._ (1995) V.-O. de Haan, A. A. van Well, S. Adenwalla, and G. P. Felcher, Phys. Rev. B 52, 10831 (1995).
  * Majkrzak  _et al._ (2003) C. F. Majkrzak, N. F. Berk, and U. A. Perez-Salas, Phys. Rev. B 19, 7796 (2003).
  * Nikova  _et al._ (2019) E. S. Nikova, Y. A. Salamatov, E. A. Kravtsov, and V. V. Ustinov, J. Phys.: Conf. Ser. 1389, 012153 (2019).
  * Leeb  _et al._ (1998) H. Leeb, J. Kasper, and R. Lipperheide, Phys. Lett. A 239, 147 (1998).
  * Fiedeldey  _et al._ (1992) H. Fiedeldey, R. Lipperheide, H. Leeb, and S. A. Sofianos, Phys. Lett. A 170, 347 (1992).
  * Majkrzak and Berk (1995b) C. F. Majkrzak and N. F. Berk, Phys. Rev. B 58, 15416 (1995b).
  * Majkrzak  _et al._ (2000) C. F. Majkrzak, N. F. Berk, S. Kreuger, J. A. Dura, M. Tarek, D. Tobias, V. Silin, J. Meuse, C. W. Woodward, and A. L. Plant, Biophys. J. 79, 3330 (2000).
  * Koutsioubas (2019) A. Koutsioubas, J. Appl. Crystallogr. 52, 538 (2019).
  * Pedersen and Hamley (1994) J. S. Pedersen and I. W. Hamley, Physica B 198, 16 (1994).
  * de Haan and Drijkoningen (1994) V.-O. de Haan and G. G. Drijkoningen, Physica B 198, 24 (1994).
  * Nelson (2006) A. R. J. Nelson, J. Appl. Crystallogr. 39, 273 (2006).
  * van der Lee  _et al._ (2007) A. van der Lee, F. Salah, and B. Harzallah, J. Appl. Crystallogr. 40, 820 (2007).
  * Gerelli (2016a) Y. Gerelli, J. Appl. Crystallogr. 49, 330 (2016a).
  * Gerelli (2016b) Y. Gerelli, J. Appl. Crystallogr. 49, 712 (2016b).
  * Nelson and Prescott (2019) A. R. J. Nelson and S. W. Prescott, J. Appl. Crystallogr. 52, 193 (2019).
  * Abelès (1948) F. Abelès, Ann. Phys. 12, 504 (1948).
  * Parratt (1954) L. G. Parratt, Phys. Rev. 95, 359 (1954).
  * Klibanov and Sacks (1992) M. V. Klibanov and P. E. Sacks, J. Math. Phys. 33, 3813 (1992).
  * Klibanov and Sacks (1994) M. V. Klibanov and P. E. Sacks, J. Comp. Phys. 112, 273 (1994).
  * Waldie  _et al._ (2018) S. Waldie, T. K. Lind, K. Browning, M. Moulin, M. Haertlein, V. T. Forsyth, A. Luchini, G. A. Strohmeier, H. Pichler, S. Maric, and M. Cárdenas, Langmuir 34, 472 (2018).
  * Campbell  _et al._ (2018) R. A. Campbell, Y. Saaka, Y. Shao, Y. Gerelli, R. Cubitt, E. Nazaruk, D. Matyszewska, and M. J. Lawrence, J. Colloid Interface Sci. 531, 98 (2018).
  * Schalke and Lösche (2000) M. Schalke and M. Lösche, Adv. Colloid Interface Sci. 88, 243 (2000).
  * Heinrich and Lösche (2014) F. Heinrich and M. Lösche, BBA-Biomembranes 1838, 2341 (2014).
  * Heinrich (2016) F. Heinrich, in  _Isotope Labeling of Biomolecules - Applications_, Methods in Enzymology, Vol. 566, edited by Z. Kelman (Academic Press, 2016) pp. 211–230.
  * Sun  _et al._ (1994) W.-J. Sun, R. M. Suter, M. A. Knewtson, C. R. Worthington, S. Tristram-Nagle, R. Zhang, and J. F. Nagle, Phys. Rev. Lett. 49, 4665 (1994).
  * Armen  _et al._ (1998) R. S. Armen, O. D. Uitto, and S. E. Feller, Biophys. J. 75, 734 (1998).
  * Hughes (2017) A. V. Hughes, “Rascal sourceforge,” https://sourceforge.net/projects/rscl/ (2017), accessed: 2019-10-09.
  * Kienzle  _et al._ (2017) P. A. Kienzle, B. B. Maranville, K. V. O’Donovan, J. F. Ankner, N. K. Berk, and C. F. Majkrzak, “NCNR reflectometry software,” https://www.nist.gov/ncnr/reflectometry-software (2017), accessed: 2019-09-15.
  * Bayes (1763) T. Bayes, Phil. Trans. Roy. Soc. 53, 370 (1763).
  * Lee and Ditko (1962) S. Lee and S. Ditko, Amazing Fantasy (Marvel Comics) 15 (1962).
  * Mayer  _et al._ (2010) J. Mayer, K. Khairy, and J. Howard, Am. J. Phys. 78, 648 (2010).
  * Sivia  _et al._ (1991) D. S. Sivia, W. A. Hamilton, and G. S. Smith, Physica B 173, 121 (1991).
  * Geoghagan  _et al._ (1996) M. Geoghagan, R. A. L. Jones, D. S. Sivia, J. Penfold, and A. S. Clough, Thin Solid Films 53, 825 (1996).
  * Sivia and Webster (1998) D. S. Sivia and J. R. P. Webster, Physica B 248, 327 (1998).
  * Hughes  _et al._ (2019) A. V. Hughes, D. S. Patel, G. Widmalm, J. B. Klauda, L. A. Clifton, and W. Im, Biophys. J. 116, 1095 (2019).
  * Note (1) In this work, a reduced notation has been used, where the background information parameter, $I$ is implied.
  * Pullen and Morris (2014) N. Pullen and R. J. Morris, PLOS ONE 9, e88419 (2014).
  * Skilling (2006) J. Skilling, AIP Conference Proceedings 735, 395 (2006).
  * Speagle (2019) J. S. Speagle, “dynesty: A Dynamic Nested Sampling Package for Estimating Bayesian Posteriors and Evidences,” (2019), arXiv:1904.02180 .
  * Sivia and Skilling (2005) D. S. Sivia and J. Skilling,  _Data Analysis: A Bayesian Tutorial_ , 2nd ed. (Oxford University Press, Oxford, UK, 2005).
  * Nelson  _et al._ (2019) A. Nelson, S. Prescott, A. Hannum, A. R. McCluskey, and I. Gresham, “refnx/refnx v0.1.7,”  (2019).
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, J. Am. Stat. Assoc. 90, 773 (1995).
  * Cornish and Littenberg (2007) N. J. Cornish and T. B. Littenberg, Phys. Rev. D 76, 083006 (2007).
  * Ensign and Pande (2010) D. L. Ensign and V. S. Pande, J. Phys. Chem. B 114, 280 (2010).
  * Hollinshead  _et al._ (2009) C. M. Hollinshead, R. D. Harvey, D. J. Barlow, J. R. P. Webster, A. V. Hughes, A. Weston, and M. J. Lawrence, Langmuir 25, 4070 (2009).
  * McCluskey  _et al._ (2019b) A. R. McCluskey, J. Grant, A. J. Smith, J. L. Rawle, D. J. Barlow, M. J. Lawrence, S. C. Parker, and K. J. Edler, J. Phys. Comm. 3, 075001 (2019b).
  * Névot and Croce (1980) L. Névot and P. Croce, Rev. Phys. Appl. (Paris) 15, 761 (1980).
  * Kučerka  _et al._ (2004) N. Kučerka, M. A. Kiselev, and P. Balgavý, Eur. Biophys. J. 33, 328 (2004).
  * Balgavý  _et al._ (2001) P. Balgavý, N. Kučerka, V. I. Gordeliy, and V. G. Cherezov, Acta. Phys. Slovaca 51, 53 (2001).
  * Tanford (1980) C. Tanford,  _The Hydrophobic Effect: Formation of Micelles and Biological Membranes_ , 2nd ed. (John Wiley & Sons, New York, USA, 1980).
  * Braslau  _et al._ (1985) A. Braslau, M. Deutsch, P. S. Pershan, A. H. Weiss, J. Als-Nielsen, and J. Bohr, Phys. Rev. Lett. 54, 114 (1985).
  * Sinha  _et al._ (1988) S. K. Sinha, E. B. Sirota, S. Garoff, and H. B. Stanley, Phys. Rev. B 38, 2297 (1988).
  * Tikhonov  _et al._ (2000) A. M. Tikhonov, D. M. Mitrinovic, M. Li, Z. Huang, and M. L. Schlossman, J. Phys. Chem. B 104, 6336 (2000).
  * Note (2) These values are tabulated in Table S1.
  * Li  _et al._ (1998) Z. X. Li, P. N. Thirtle, R. K. Thomas, J. Penfold, J. R. P. Webster, and A. R. Rennie, Physica B 248, 171 (1998).
  * DeWolf  _et al._ (1999) C. DeWolf, S. Leporatti, C. Kirsch, R. Klinger, and G. Brezesinski, Chem. Phys. Lipids 97, 129 (1999).
  * Brezesinski  _et al._ (2001) G. Brezesinski, H. J. Müller, J. L. Toca-Herrera, and R. Krustev, Chem. Phys. Lipids 110, 183 (2001).
  * Struth  _et al._ (2001) B. Struth, F. Rieutord, O. Konovalov, G. Brezesinski, G. Grübel, and P. Terech, Phys. Rev. Lett. 88, 025502 (2001).
  * Sivia  _et al._ (1993) D. S. Sivia, W. I. F. David, K. S. Knight, and S. F. Gull, Physica D 66, 234 (1993).