The Use of Score Tests for Inference on Variance Components–A.1 December

10.1111/j.1541-0420.2005.00454.x

# The Use of Score Tests for Inference on Variance Components

Geert<EMAIL_ADDRESS>
Biostatistical Centre  Catholic University of Leuven  Leuven  Belgium  Geert
<EMAIL_ADDRESS>
Center for Statistics  Hasselt University  B-3590 Diepenbeek  Belgium

(Received April 2007\. Revised April 2007\. Accepted April 2007.; 2007)

###### Abstract

Whenever inference for variance components is required, the choice between
one-sided and two-sided tests is crucial. This choice is usually driven by
whether or not negative variance components are permitted. For two-sided
tests, classical inferential procedures can be followed, based on likelihood
ratios, score statistics, or Wald statistics. For one-sided tests, however,
one-sided test statistics need to be developed, and their null distribution
derived. While this has received considerable attention in the context of the
likelihood ratio test, there appears to be much confusion about the related
problem for the score test. The aim of this paper is to illustrate that
classical (two-sided) score test statistics, frequently advocated in practice,
cannot be used in this context, but that well-chosen one-sided counterparts
could be used instead. The relation with likelihood ratio tests will be
established, and all results are illustrated in an analysis of continuous
longitudinal data using linear mixed models.

###### keywords:

Boundary condition; Likelihood ratio test; Linear mixed model; One-sided test;
Score test; Variance component.

††volume: 63

##  1 Introduction

In a variety of applied statistical problems, there is a need for inference on
variance components. This includes a variety of applied fields, for example,
random-effects ANOVA models (Nelder 1954), linear mixed models (Verbeke and
Molenberghs 2000), generalized linear and non-linear (mixed) models (Jacqmin-
Gadda and Commenges 1995), overdispersion (Cox 1983, Smith and Heitjan 1993,
Hines 1997, Lu 1997), clustering (Britton 1997) and homogeneity in stratified
analyses (Liang 1987).

To fix ideas, we will focus on the setting of a relatively simple linear mixed
model, the so-called random-intercepts model:

$Y_{ij}=\mbox{\boldmath$x$}_{ij}^{\prime}\mbox{\boldmath$\beta$}+b_{i}+\varepsilon_{ij},$
(1)

where $Y_{ij}$ is the response for member $j=1,\dots,n_{i}$ of cluster
$i=1,\dots,N$, $\mbox{\boldmath$x$}_{ij}$ is a vector of known covariate
values, $\beta$ is a vector of unknown regression coefficients, $b_{i}\sim
N(0,\tau^{2})$ is a cluster-specific random effect, assumed to be
independently distributed from the residual error components
$\varepsilon_{ij}\sim N(0,\sigma^{2})$. Classical inferential procedures are
based on the likelihood of the marginal model, obtained by integrating (1)
over the random effects. Grouping the $Y_{ij}$ into a vector
$\mbox{\boldmath$Y$}_{i}$ and assembling the rows
$\mbox{\boldmath$x$}_{ij}^{\prime}$ into a matrix $X_{i}$, this marginal
distribution is

$\mbox{\boldmath$Y$}_{i}\sim
N(X_{i}\mbox{\boldmath$\beta$},\tau^{2}J_{n_{i}}+\sigma^{2}I_{n_{i}}),$ (2)

in which $I_{n_{i}}$ denotes the identity matrix of dimension $n_{i}$, and
where $J_{n_{i}}$ equals the $n_{i}\times n_{i}$ matrix containing only ones.

Regarding the variance component $\tau^{2}$ in the above model, one can take
two views. In the first view, where the focus is entirely on the resulting
marginal model (2), negative values for $\tau^{2}$ are perfectly acceptable
(Nelder 1954, Verbeke and Molenberghs 2000, Sec. 5.6.2), since this merely
corresponds to the occurrence of negative within-cluster correlation
$\rho=\tau^{2}/(\tau^{2}+\sigma^{2})$. This might occur, for example, in a
context of competition such as when littermates compete for the same food
resources. In such a case, the only requirement is that
$\tau^{2}+\sigma^{2}>0$, for $V_{i}=\tau^{2}J_{n_{i}}+\sigma^{2}I_{n_{i}}$ to
be a positive definite, marginal covariance matrix. Further discussions on
negative variance components can be found in Thompson (1962) and Searle,
Casella and McCulloch (1992). In the second view, when the link between the
marginal model (2) and its generating hierarchical model (1) is preserved,
thereby including the concept of random effects $b_{i}$ and perhaps even
requiring inference for them, it is imperative to restrict $\tau^{2}$ to
nonnegative values.

The first situation, which we will term the unconstrained case, is standard
regarding inference for the variance component $\tau^{2}$. In the second
situation (the constrained case), however, one typically needs one-sided tests
of the null-hypothesis

$H_{0}:\tau^{2}=0\hskip 28.45274pt\mbox{\em versus}\hskip
28.45274ptH_{A1}:\tau^{2}>0.$ (3)

As the null-hypothesis is now on the boundary of the parameter space,
classical inference no longer holds, appropriate tailored test statistics need
to be developed, and the corresponding (asymptotic) null distributions
derived.

While this has received considerable attention in the case of the likelihood
ratio test (Self and Liang 1987, Stram and Lee 1994, 1995), there is still
much confusion about the related problem for the score test. For some uses of
a score test in boundary situations, see Liang (1987), Lin (1997), Gray
(1995), Jacqmin-Gadda and Commenges (1995), Dean and Lawless (1989), Dean
(1992), Dean, Ugarte and Militino (2001), Gueorguieva (2001), Militino, Ugarte
and Dean (2001), Smith and Heitjan (1993). Some authors implicitly take the
unconstrained, two-sided view, with a few happy exceptions who explicitly
adopt such a view (Paul and Islam 1995). Jacqmin-Gadda and Commenges (1995),
Lin (1997), le Cessie and van Houwelingen (1995), and Dean, Ugarte, and
Militino (2001) do not explicitly specify the alternative model, thereby
implicitly assuming two-sided alternatives, while clearly being in a one-sided
setting. We hope to illustrate that, when required by the scientific problem,
a fully one-sided approach is both feasible and more appropriate. Silvapulle
and Silvapulle (1995) have shown how a one-sided score test can be defined,
both in the scalar as well as in the vector parameter case. Important related
work is given in Hall and Præstgaard (2001). While these authors also focus on
the restricted score tests in the context of mixed models, there are three
important differences with our take on the problem. First, Hall and Præstgaard
(2001) explicitly advocate the use of restricted score tests, thereby
improving upon earlier work (Lin 1997) in terms of efficiency. We point out
that the choice between a constrained/unconstrained setting should be tightly
linked to a constrained/unconstrained alternative space. Second, since our
score test statistics follow from the work of Silvapulle and Silvapulle
(1995), their analytic forms are slightly different from those of Hall and
Præstgaard (2001). Indeed, based on the results of Silvapulle and Silvapulle
(1995), who showed that the asymptotic equivalence of the likelihood ratio and
score tests holds, also in the constrained case, the null distribution of the
one-sided score tests will be derived. Finally, we put a lot of emphasis on
the extension of the well-known asymptotic equivalence of the likelihood ratio
and score tests to the constrained case, as follows from Silvapulle and
Silvapulle (1995). We will argue that based on this equivalence, the
researcher has full choice between both testing procedures, and moreover,
opting for a constrained likelihood ratio tests has many computational
advantages in practice. Emphasis will be on intuitive explanation of the
theoretical results, rather than on mathematical details.

In Section 2, we continue with our initial model (1), and we will show how
one-sided likelihood ratio and score tests can be constructed, and the
corresponding asymptotic null distribution will be derived heuristically.
Afterwards, in Sections 3 and 4, more general results will be discussed for
the likelihood ratio test and for the score test, respectively. Note that our
aim is not to argue for or against score tests, but rather show how to
properly use one-sided score tests for variance components. Computational
issues are discussed in Section 5. In Section 6, the results will be
illustrated in an analysis of continuous longitudinal measurements, using
linear mixed models, where the need for random effects is to be tested.
Finally, Section 7 summarizes the main results.

##  2 The Random-intercepts Model

To introduce our ideas in a simple but generic setting, we continue the
discussion of the random-intercepts model (1). Under the unconstrained
parameterization, i.e., the model under which negative values for $\tau^{2}$
are allowed, classical inferential tools are available for testing the general
two-sided hypothesis

$\displaystyle H_{0}:\tau^{2}=0\hskip 28.45274pt$ versus $\displaystyle\hskip
28.45274ptH_{A2}:\tau^{2}\neq 0.$

Wald, likelihood ratio, and score tests are then asymptotically equivalent,
and the asymptotic null distribution is well known to be $\chi^{2}_{1}$ (Cox
and Hinkley 1990). Under the constrained model, i.e., the model where
$\tau^{2}$ is restricted to the non-negative real numbers, the one-sided
hypothesis (3) is the only meaningful one.

Appropriate test statistics can now be obtained as follows. Suppressing
dependence on the other parameters, let $\ell(\tau^{2})$ denote the log-
likelihood, as a function of the random-intercepts variance $\tau^{2}$.
Further, let $\widehat{\tau}^{2}$ denote the maximum likelihood estimate of
$\tau^{2}$ under the unconstrained parameterization. We first consider the
likelihood ratio test, with statistic:

$\displaystyle T_{LR}$ $\displaystyle=$ $\displaystyle
2\ln\left[\frac{\max_{H_{1A}}\ell(\tau^{2})}{\max_{H_{0}}\ell(\tau^{2})}\right].$

Two cases, graphically represented in Figure 1, can now be distinguished.
Under Case A, $\widehat{\tau}^{2}$ is positive, and the likelihood ratio test
statistic is identical to the one that would be obtained under the
unconstrained parameter space for $\tau^{2}$. Hence, conditionally on
$\widehat{\tau}^{2}\geq 0$, $T_{LR}$ has asymptotic null distribution equal to
the classical $\chi^{2}_{1}$.

Figure 1: Graphical representation of two different situations, when
developing one-sided tests for the variance $\tau^{2}$ of the random
intercepts $b_{i}$ in model (1).

Under Case B however, we have that, under $H_{1A}$ as well as under $H_{0}$,
$\ell(\tau^{2})$ is maximized at $\tau^{2}=0$ yielding $T_{LR}=0$. Further
note that, under $H_{0}$, both cases occur with 50% probability. Hence, the
asymptotic null distribution of $T_{LR}$ is obtained from

$\displaystyle P(T_{LR}>c|H_{0})$ $\displaystyle=$ $\displaystyle
P(T_{LR}>c|H_{0},\widehat{\tau}^{2}\geq 0)P(\widehat{\tau}^{2}\geq 0|H_{0})$
$\displaystyle+P(T_{LR}>c|H_{0},\widehat{\tau}^{2}<0)P(\widehat{\tau}^{2}<0|H_{0})$
$\displaystyle=$
$\displaystyle\frac{1}{2}P(\chi^{2}_{1}>c)+\frac{1}{2}P(\chi_{0}^{2}>c),$

where $\chi_{0}^{2}$ denotes the distribution with all probability mass at
$0$. Hence, the asymptotic null distribution of the one-sided likelihood ratio
test statistic is a mixture of two chi-squared distributions, with degrees of
freedom $0$ and $1$, and with equal mixing proportions $1/2$. This was one of
Stram and Lee’s (1994, 1995) special cases. Note that, whenever
$\widehat{\tau}^{2}\geq 0$, the observed likelihood ratio test statistic is
equal to the one under the unconstrained model, but the $p$-value is half the
size of the one obtained from the classical $\chi_{1}^{2}$ approximation to
the null distribution.

We now consider the score test. The usual form of the test statistic is given
by

$\displaystyle T_{S}$ $\displaystyle=$
$\displaystyle\left[\left.\frac{\partial\ell(\tau^{2})}{\partial\tau^{2}}\right|_{\tau^{2}=0}\right]^{2}\left[\left.-\frac{\partial^{2}\ell(\tau^{2})}{\partial\tau^{2}\partial\tau^{2}}\right|_{\tau^{2}=0}\right]^{-1}.$
(4)

Nuisance parameters are suppressed from notation, and replaced by their MLE’s.
In the special case of (2) with $n_{i}\equiv n$, straightforward algebra
produces:

$T_{S}=\frac{Nn}{2}\,\frac{(C-1)^{2}}{2nC-1},$

with

$C=\frac{1}{\sigma^{2}}\,\frac{1}{Nn}\sum_{i=1}^{N}\left(\sum_{j=1}^{n}y_{ij}\right)^{2},$

in which $\sigma^{2}$ is replaced by its maximum likelihood estimate under the
null hypothesis, $\widehat{\sigma}^{2}$ say. Without loss of generality, it is
assumed that the fixed-effects parameters are zero. In the reverse case,
$y_{ij}$ needs to be replaced by appropriate residuals.

Now, score test (4) implicitly assumes a two-sided alternative. Hence, the
test statistic itself need to be redefined appropriately in order to be able
to discriminate between positive and negative alternative values for
$\tau^{2}$. The same two cases as for the likelihood ratio test can be
considered (see Figure 1). Under Case A, $\widehat{\tau}^{2}$ is positive, and
the positive score $\partial\ell(\tau^{2})/\partial\tau^{2}$ at zero is
evidence against $H_{0}$ in favor of the one-sided alternative $H_{A1}$.
Hence, (4) can be used as test statistic, provided that
$\widehat{\tau}^{2}\geq 0$. This implies that, conditionally on
$\widehat{\tau}^{2}\geq 0$ and under $H_{0}$, our test statistic
asymptotically follows the classical $\chi^{2}_{1}$ distribution. Under Case
B, however, the score at $\tau^{2}=0$ is negative, and can therefore clearly
not be used as evidence against $H_{0}$ in favor of $H_{A1}$. Hence, whenever
$\widehat{\tau}^{2}$ is negative, (4) is no longer meaningful as test
statistic. Considering that a negative score at zero supports the null
hypothesis, a meaningful test statistic is obtained from replacing (4) by

$\displaystyle T_{S}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\left[\left.\frac{\partial\ell(\tau^{2})}{\partial\tau^{2}}\right|_{\tau^{2}=0}\right]^{2}\left[\left.-\frac{\partial^{2}\ell(\tau^{2})}{\partial\tau^{2}\partial\tau^{2}}\right|_{\tau^{2}=0}\right]^{-1}&\mbox{if
}\widehat{\tau}^{2}\geq 0\\\ \\\ 0&\mbox{if
}\widehat{\tau}^{2}<0.\end{array}\right.$ (8)

The corresponding asymptotic null distribution is now obtained from

$\displaystyle P(T_{S}>c|H_{0})$ $\displaystyle=$ $\displaystyle
P(T_{S}>c|H_{0},\widehat{\tau}^{2}\geq 0)P(\widehat{\tau}^{2}\geq 0|H_{0})$
$\displaystyle+P(T_{S}>c|H_{0},\widehat{\tau}^{2}<0)P(\widehat{\tau}^{2}<0|H_{0})$
$\displaystyle=$
$\displaystyle\frac{1}{2}P(\chi^{2}_{1}>c)+\frac{1}{2}P(\chi_{0}^{2}>c),$

which is identical to the null distribution derived earlier for the likelihood
ratio test. This heuristic but insightful argument will be formalized and
generalized in Section 4.

Note that in this scalar case, an equivalent test consists of appropriately
standardizing the score function rather than embedding it in a quadratic form.
The choice between cases A and B then merely becomes a choice between
classical two-sided versus one-sided $Z$-type test procedures.

##  3  Likelihood Ratio Tests

When the use of the likelihood-ratio test is envisaged, it is now well known
that hypotheses such as (3) pose non-standard testing problems (Verbeke and
Molenberghs 2000, pp. 64–73). Such problems have been known for a long time
(Nelder 1954, Chernoff 1954). Using results of Self and Liang (1987) on
nonstandard testing situations, Stram and Lee (1994, 1995) have been able to
show that the asymptotic null distribution for the likelihood ratio test
statistic for testing hypotheses of the type (3) is often a mixture of chi-
squared distributions rather than the classical single chi-squared
distribution. For ANOVA models with independent random effects, this was
already briefly discussed by Miller (1977). This complication cannot be
relieved by considering alternative parameterizations for the variance
components, contradicting the once popular but false belief that replacing
covariance matrices by their Cholesky decomposition was able to turn the
problem in a standard one. Indeed, under the null hypothesis, the Cholesky
decomposition does not map 1 to 1 onto the original parameterization.

Stram and Lee (1994, 1995) discuss likelihood ratio tests for variance
components in linear mixed models, which are generalizations of (1) to models
with multiple, possibly correlated, random effects. These models typically
appear in the analysis of continuous longitudinal data. Let $\bm{Y_{i}}$
denote the $n_{i}$-dimensional vector of measurements available for subject
$i$, $i=1,\ldots,N$. A general linear mixed model then assumes that
$\bm{Y_{i}}$ satisfies

$\displaystyle\bm{Y_{i}}$ $\displaystyle=$ $\displaystyle
X_{i}\bm\beta+Z_{i}\bm{b_{i}}+\bm{\varepsilon_{i}},$ (9)

in which $\bm\beta$ is a vector of population-average regression coefficients
called fixed effects, and where $\bm{b_{i}}$ is a vector of subject-specific
regression coefficients. The $\bm{b_{i}}$ describe how the evolution of the
$i$th subject deviates from the average evolution in the population. The
matrices $X_{i}$ and $Z_{i}$ are $(n_{i}\times p)$ and $(n_{i}\times q)$
matrices of known covariates. The random effects $\bm{b_{i}}$ and residual
components $\bm{\varepsilon_{i}}$ are assumed to be independent with
distributions $N(\bm{0},D)$, and $N(\bm{0},\sigma^{2}I_{n_{i}})$,
respectively. Inference for linear mixed models is based on maximum likelihood
or restricted maximum likelihood estimation under the marginal model for
$\bm{Y_{i}}$, i.e., the multivariate normal model with mean $X_{i}\bm\beta$,
and covariance $V_{i}=Z_{i}DZ_{i}^{\prime}+\sigma^{2}I_{n_{i}}$ (Laird and
Ware 1982, Verbeke and Molenberghs 2000).

Similar to our simpler model (1), the marginal model does not require $D$ to
be positive definite, while a random-effects interpretation of the model does,
corresponding to the unconstrained and constrained parameterizations,
respectively. As before, inference under the unconstrained model for the
variance components in $D$ can be based on the classical chi-squared
approximation to the null distribution for the likelihood ratio test
statistic. Under the constrained model, Stram and Lee (1994, 1995) have shown
that the asymptotic null distribution for the likelihood ratio test statistic
for testing a null hypothesis which allows for $k$ correlated random effects
versus an alternative of $k+1$ correlated random effects (with positive semi-
definite covariance matrix $D_{k+1}$), is a mixture of a $\chi^{2}_{k}$ and a
$\chi^{2}_{k+1}$, with equal probability $1/2$. For more general settings,
e.g., comparing models with $k$ and $k+k^{\prime}$ ($k^{\prime}>1$) random
effects, the null distribution is a mixture of $\chi^{2}$ random variables
(Shapiro 1988, Raubertas, Lee, and Nordheim 1986), the weights of which can
only be calculated analytically in a number of special cases.

##  4 Score Tests

Our heuristic arguments in Section 2 have suggested that employment of score
tests for testing variance components under the constrained parameterization
requires replacing the classical score test statistic by an appropriate one-
sided version. This is where the general theory of Silvapulle and Silvapulle
(1995) on one-sided score tests proves very useful. They consider models
parameterized through a vector
$\mbox{\boldmath$\theta$}=(\mbox{\boldmath$\lambda$}^{\prime},\mbox{\boldmath$\psi$}^{\prime})^{\prime}$,
where testing a general hypothesis of the form

$\displaystyle H_{0}:\mbox{\boldmath$\psi$}={\bm{0}}\hskip 28.45274pt\mbox{\em
versus}\hskip 28.45274ptH_{A}:\mbox{\boldmath$\psi$}\in{\cal C}$

is of interest. In our context, the alternative parameter space ${\cal C}$
equals the nonnegative real numbers (e.g., when testing (3), or the set of
positive semi-definite covariance matrices $D$ (e.g., when testing for the
need of variance components in linear mixed models, Section 3). In general,
Silvapulle and Silvapulle (1995) allow ${\cal C}$ to be a closed and convex
cone in Euclidean space, with vertex at the origin. The advantage of such a
general definition is that one-sided, two-sided, and combinations of one-sided
and two-sided hypotheses are included.

Silvapulle and Silvapulle (1995) consider a general class of score-type test
statistics. Applying their theory to our situation yields the following
results. Let the log-likelihood function be denoted by
$\ell(\mbox{\boldmath$\theta$})$. The associated score function equals

$\displaystyle\mbox{\boldmath$S$}_{N}(\mbox{\boldmath$\theta$})$
$\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\mbox{\boldmath$\theta$}}.$

Assume the existence of a nonsingular matrix $H(\mbox{\boldmath$\theta$})$
such that, for $N\rightarrow\infty$,

$\displaystyle(A1):N^{-1/2}\mbox{\boldmath$S$}_{N}(\mbox{\boldmath$\theta$})$
$\displaystyle\stackrel{{\scriptstyle d}}{{\rightarrow}}$ $\displaystyle
N({\bm{0}},H(\mbox{\boldmath$\theta$}))$

and, for all $a\geq 0$,

$\displaystyle(A2):\sup_{||\mbox{\boldmath$h$}||\leq a}$
$\displaystyle\left[N^{-1/2}\left\\{\mbox{\boldmath$S$}_{N}(\mbox{\boldmath$\theta$}+N^{-1/2}\mbox{\boldmath$h$})-\mbox{\boldmath$S$}_{N}(\mbox{\boldmath$\theta$})\right\\}\right.$
$\displaystyle\left.+H(\mbox{\boldmath$\theta$})\mbox{\boldmath$h$}\right]=o_{p}(1).$

Further, decompose $\mbox{\boldmath$S$}_{N}$ as
$\mbox{\boldmath$S$}_{N}=(\mbox{\boldmath$S$}_{N\lambda}^{\prime},\mbox{\boldmath$S$}_{N\psi}^{\prime})^{\prime}$,
let $H_{\lambda\lambda}(\mbox{\boldmath$\theta$})$,
$H_{\lambda\psi}(\mbox{\boldmath$\theta$})$ and
$H_{\psi\psi}(\mbox{\boldmath$\theta$})$ be the corresponding blocks in
$H(\mbox{\boldmath$\theta$})$, and define
$\mbox{\boldmath$\theta$}_{H}=(\mbox{\boldmath$\lambda$}^{\prime},{\bm{0}}^{\prime})^{\prime}$.
$\mbox{\boldmath$\theta$}_{H}$ can be estimated by
$\widehat{\mbox{\boldmath$\theta$}}_{H}=(\widehat{\mbox{\boldmath$\lambda$}}^{\prime},{\bm{0}}^{\prime})^{\prime}$,
in which $\widehat{\mbox{\boldmath$\lambda$}}$ is the maximum likelihood
estimate of $\lambda$, under $H_{0}$. Finally, let $\bm{Z}_{N}$ be equal to
$\bm{Z}_{N}=N^{-1/2}\mbox{\boldmath$S$}_{N\psi}(\widehat{\mbox{\boldmath$\theta$}}_{H})$.
A one-sided score statistic can now be defined as

$\displaystyle T_{S}$ $\displaystyle:=$
$\displaystyle\bm{Z}^{\prime}_{N}H_{\psi\psi}^{-1}(\widehat{\mbox{\boldmath$\theta$}}_{H})\bm{Z}_{N}$
(10) $\displaystyle-\
\inf\left\\{(\bm{Z}_{N}-\bm{b})^{\prime}H_{\psi\psi}^{-1}(\widehat{\mbox{\boldmath$\theta$}}_{H})(\bm{Z}_{N}-\bm{b})|\bm{b}\in\mbox{${\cal
C}$}\right\\}.$

Note that our score statistic (8) derived heuristically for the random-
intercepts model is a special case of the general test statistic (10). Indeed,
when $\widehat{\tau}^{2}$ is positive, the score at zero is positive, and
therefore in ${\cal C}$, such that the infimum in (10) becomes zero. For
$\widehat{\tau}^{2}$ negative, the score at zero is negative as well and the
infimum in (10) is attained for $\bm{b}=\bm{0}$, resulting in $T_{S}=0$.

It follows from Silvapulle and Silvapulle (1995) that, provided regularity
conditions (A1) and (A2) hold, for $N\rightarrow\infty$, the likelihood ratio
test statistic $T_{LR}$ satisfies $T_{LR}=T_{S}+o_{p}(1)$. Further, if the
observed $T_{s}$ equals $t_{s}$, then the large-sample $p$-value equals

$p=\xi\left\\{t_{s},H_{\psi\psi}(\mbox{\boldmath$\theta$}_{H}),{\cal
C}\right\\},$ (11)

where

$\displaystyle\xi(t,B,{\cal C})$ $\displaystyle=$
$\displaystyle\mbox{P}\left[\mbox{\boldmath$Z$}^{\prime}B^{-1}\mbox{\boldmath$Z$}\right.$
$\displaystyle\phantom{\mbox{P}}\left.-\inf\left\\{(\mbox{\boldmath$Z$}-\bm{b})^{\prime}B^{-1}(\mbox{\boldmath$Z$}-\bm{b})|\bm{b}\in{\cal
C}\right\\}\geq t\right]$

and $\mbox{\boldmath$Z$}\sim N(0,B)$. Shapiro (1988, Eqs. 3.1 and 3.2) has
shown that $1-\xi(t,B,{\cal C})$ equals a weighted sum of chi-squared
probabilities. The results obtained by Stram and Lee (1994, 1995) for the
linear mixed model are included in Shapiro’s results. There are a few
additional results available. For example, if the null hypothesis allows for
$k$ uncorrelated random effects (with a diagonal covariance matrix $D_{k}$)
versus the alternative of $k+k^{\prime}$ uncorrelated random effects (with
diagonal covariance matrix $D_{k+k^{\prime}}$), the null distribution is a
mixture of the form

$\sum_{m=0}^{k^{\prime}}2^{-k^{\prime}}{k^{\prime}\choose m}\chi^{2}_{m}.$

Shapiro (1988) shows that, for a broad number of cases, determining the
mixture’s weights is a complex and perhaps numerical task.

The above results show that the equivalence of the score and likelihood ratio
tests not only holds in the two-sided but also in the one-sided cases. At the
same time, an appropriate definition of the one-sided score statistic is
produced, and its asymptotic null distribution has been derived. In some
cases, the analytic null distribution is easily obtained through the results
of Stram and Lee (1994, 1995), summarized in Section 3, and the equivalence
between the likelihood ratio and score tests.

Finally, it should be emphasized that all above results are valid provided
that the conditions $(A1)$ and $(A2)$ are satisfied. In particular, $(A2)$
requires that the score $\mbox{\boldmath$S$}_{N}$ exists in a sufficiently
small neighborhood around $H_{0}$. For example, in our random-intercepts
example (Section 2), it is crucial to have valid models for sufficiently small
but negative values of $\tau^{2}$, even in the constrained setting. As a
counterexample, if we were to test $H_{0}:\sigma^{2}=0$ versus the one-sided
alternative $H_{A}:\sigma^{2}>0$ for the variance $\sigma^{2}$ in a univariate
normal $N(0,\sigma^{2})$ sample of size $N$, the score equals

$\displaystyle\mbox{\boldmath$S$}_{N}(\sigma^{2})$ $\displaystyle=$
$\displaystyle\frac{\partial\ell}{\partial\sigma^{2}}\ =\
-\frac{N}{2\sigma^{2}}\ +\ \frac{1}{2\sigma^{4}}\sum_{i=1}^{N}y_{i}^{2},$

which, evaluated at $\sigma^{2}=0$, yields $+\infty$. Then, the above theory
does not apply here as no negative values for $\sigma^{2}$ can ever yield a
valid statistical model for our sample. Hence, in this example, condition
$(A2)$ is no longer satisfied.

##  5 Computational Issues

In this section, computations for likelihood ratio and score test statistics
will be discussed, for the classical unconstrained as well as the constrained
cases. The relative complexity of all four cases will be addressed. Here, we
focus on the general concepts, while a particular example using the SAS
procedures MIXED and NLMIXED is deferred to the Appendix.

In principle, calculation of the unconstrained likelihood ratio test statistic
does not pose any specific complications, provided both the null and
alternative models can be fitted with standard software, and both log-
likelihood values at maximum are returned, minus twice the difference of which
is then referred to the appropriate chi-squared distribution, to yield the
$p$-value.

Even in the unconstrained case, the score test calculations are more involved
than their likelihood ratio counterparts, since the first and second order
derivatives of the alternative log-likelihood function are required, evaluated
under the null hypothesis. These cannot easily be obtained in many standard
packages without additional programming. Once these derivatives have been
obtained, they are the straightforward building blocks for the calculation of
the test statistic, while the $p$-value is obtained as in the previous case.

The constrained likelihood ratio test statistic can be obtained in the same
way as in the unconstrained case, provided the constraints are properly
imposed onto the alternative model. In many practical situations, this comes
down to maximizing the likelihood under a positive definiteness constraint on
a covariance matrix $D$ as, for example, in the linear mixed model setting.
Several routes can be followed. Replacing $D$ by its Cholesky decomposition
($D=L^{\prime}L$) and maximizing over $L$ rather than $D$ turns the
constrained optimization into an unconstrained one. This route has been
proposed by Lindstrom and Bates (1988). Note that the constrained testing
problem is not turned into an unconstrained one, because the Cholesky
decomposition does not map one-to-one onto the original parameterization, thus
maintaining the need for an appropriate testing theory as developed in this
article and by Hall and Præstgaerd (2001). Alternatively, a so-called barrier
type approach can be followed, for example, by adding a penalty
$a\log\\{\det(D)\\}$ to the log-likelihood function, for some pre-specified
constant $a$. While a careful consideration of the relative merits of these
and other approaches is important and interesting in its own right, it is
beyond the scope of this paper.

The constrained score test statistic (10) is composed of two parts. The first
term is identical to the unconstrained counterpart, while the second term
involves a constrained minimization of the quadratic form
$(\bm{Z}_{N}-\bm{b})^{\prime}H_{\psi\psi}^{-1}(\widehat{\mbox{\boldmath$\theta$}}_{H})(\bm{Z}_{N}-\bm{b})$
which cannot always be done analytically. In such cases, additional software
code needs to be written, invoking numerical constrained optimization
routines.

In both constrained cases, $p$-value computation is given by (11) which is a
weighted sum of chi-squared probabilities, the weights of which are known
analytically in special (but important) cases only.

##  6 Application: The Rat Data

Using a simple case study and a selected set of nested models, we illustrate
(1) likelihood ratio as well as score tests, (2) under both one-sided and two-
sided alternatives, and (3) in cases where a boundary estimate does and does
not occur.

The data considered to this end are from a randomized longitudinal experiment,
previously described and analyzed by Verdonck et al. (1998), in which 50 male
Wistar rats were randomized to either a control group or one of the two
treatment groups where treatment consisted of a low or high dose of the drug
Decapeptyl, which is an inhibitor for testosterone production in rats. The
primary aim of the study was to investigate the effect of the inhibition of
the production of testosterone in male Wistar rats on their craniofacial
growth. The treatment started at the age of 45 days, and measurements were
taken every 10 days, with the first observation taken at the age of 50 days.
One of the responses of interest was the height of the skull, measured as the
distance (in pixels) between two well-defined points on X-ray pictures of the
skull, taken after the rat has been anesthetized. The individual profiles are
shown in Figure 2. Although rats were scheduled to be followed up to the age
of 110 days, some drop out prematurely because they do not survive
anaesthesia. In fact, while 50 rats have been randomized at the start of the
experiment, only 22 of them survived the first 6 measurements, so measurements
on only 22 rats are available in the way anticipated at the design stage.

Figure 2: Rat Data. Individual profiles for each of the treatment groups in
the rat experiment separately.

As before, let $Y_{ij}$ denote the response taken at time $t_{j}$, for rat
$i$. Verbeke and Lesaffre (1999) and Verbeke and Molenberghs (2000) proposed
to model the subject-specific profiles shown in Figure 2 as linear functions
of $t=\ln(1+(\mbox{Age}-45)/10)$. More specifically, their model is of the
form

$\displaystyle Y_{ij}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l}\beta_{0}+b_{1i}+(\beta_{1}+b_{2i})t_{ij}+\varepsilon_{ij},\hskip
14.22636pt\mbox{if low dose,}\\\
\beta_{0}+b_{1i}+(\beta_{2}+b_{2i})t_{ij}+\varepsilon_{ij},\hskip
14.22636pt\mbox{if high dose,}\\\
\beta_{0}+b_{1i}+(\beta_{3}+b_{2i})t_{ij}+\varepsilon_{ij},\hskip
14.22636pt\mbox{if control.}\end{array}\right.$ (15)

Here, $\beta_{0}$ is the average response at the time of randomization, while
$\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ are the average slopes in the three
different treatment groups. Further, the $b_{1i}$ and $b_{2i}$ are rat-
specific intercepts and slopes, representing the natural heterogeneity between
rats with respect to baseline values and with respect to evolutions over time,
respectively. The above model is an example of linear mixed model (9). As in
our introductory example we have that, strictly speaking, the marginal model
does not require $D$ to be positive definite, as long as the resulting
marginal covariance $V_{i}$ is. Hence, when testing for elements in $D$, two-
sided tests can be employed. However, they then no longer allow the
hierarchical interpretation of the model, i.e., the interpretation in which
the variability in the data is believed to be generated from an underlying
random-effects model as in (15). If underlying random effects are believed to
be latently present, one-sided tests are required.

Table 1: Rat data. Summary of the results of one- as well as two-sided, likelihood ratio as well as score tests for the comparison of a series of linear mixed models fitted to the rat data. The estimate $\widehat{D}$ denotes the unconstrained maximum likelihood estimate for the matrix $D$ in the linear mixed model. Every entry consists of (1) the test statistic; (2) the null distribution; (3) the $p$-value. |  |  | LR test | Score test
---|---|---|---|---
|  | Ref. | One-sided | Two-sided | One-sided | Two-sided
1: | $D=0$ |  |  |  |  |
2: | $D=\left(\begin{array}[]{c}d_{11}\end{array}\right)$ | 1 | $1102.7-928.7=174$ | $1102.7-928.7=174$ | $27.15-0.0=27.15$ | $27.15$
|  |  | $\frac{1}{2}\chi^{2}_{0}+\frac{1}{2}\chi^{2}_{1}$ | $\chi^{2}_{1}$ | $\frac{1}{2}\chi^{2}_{0}+\frac{1}{2}\chi^{2}_{1}$ | $\chi^{2}_{1}$
| $\widehat{D}=\left(\begin{array}[]{c}3.44\end{array}\right)$ |  | $p<0.0001$ | $p<0.0001$ | $p<0.0001$ | $p<0.0001$
3: | $D=\left(\begin{array}[]{cc}d_{11}&0\\\ 0&d_{22}\end{array}\right)$ | 2 | $928.7-928.7=0.0$ | $928.7-927.4=1.3$ | $1.67-1.67=0$ | $1.67$
|  |  | $\frac{1}{2}\chi^{2}_{0}+\frac{1}{2}\chi^{2}_{1}$ | $\chi^{2}_{1}$ | $\frac{1}{2}\chi^{2}_{0}+\frac{1}{2}\chi^{2}_{1}$ | $\chi^{2}_{1}$
| $\widehat{D}=\left(\begin{array}[]{cc}3.77&0\\\ 0&-0.17\end{array}\right)$ |  | $p=1.0000$ | $p=0.2542$ | $p=1.0000$ | $p=0.1963$
4: | $D=\left(\begin{array}[]{cc}d_{11}&d_{12}\\\ d_{12}&d_{22}\end{array}\right)$ | 2 | $928.7-928.6=0.1$ | $928.7-925.8=2.9$ | $2.03-1.93=0.1$ | $2.03$
|  |  | $\frac{1}{2}\chi^{2}_{1}+\frac{1}{2}\chi^{2}_{2}$ | $\chi^{2}_{2}$ | $\frac{1}{2}\chi^{2}_{1}+\frac{1}{2}\chi^{2}_{2}$ | $\chi^{2}_{2}$
| $\widehat{D}=\left(\begin{array}[]{cc}2.83&0.48\\\ 0.48&-0.33\end{array}\right)$ |  | $p=0.8515$ | $p=0.2346$ | $p=0.8515$ | $p=0.3624$

Several models can now be fitted and compared with one another. Table 1
summarizes some of the results obtained from fitting and comparing a series of
models to the rat data. Model 1 assume independent repeated measures and does
not include any random effects; its only variance component is the common
variance $\sigma^{2}$. Model 2 includes random intercepts only and is
therefore an example of (1), assuming all measurements $Y_{ij}$ within subject
$i$ exhibit equal correlation and common variance. Finally, Models 3 and 4
include random linear time-effects $b_{2i}$ as well, which may (Model 3) or
may not (Model 4) be correlated with the random intercepts $b_{1i}$.

Table 1 shows the results of one- as well as two-sided, likelihood ratio as
well as score tests, for model comparisons 2–1, 3–2, and 4–2. Comparison 2–1
is standard in the sense that, since the unconstrained estimate of the random-
intercepts variance $d_{11}$ under Model 2 is positive, the one- and two-sided
test statistics are identical, and only the null distribution is different.
For comparison 3–2, the unconstrained estimate for the random-slopes variance
$d_{22}$ is negative, yielding different one- and two-sided test statistics.
For the score test, for example, the infimum in (10) is attained for
$\bm{b}=0$, yielding zero as observed value for the one-sided test statistic.
For comparison 4–2, the infimum in (10) needs to be calculated numerically,
and was found to be equal to $1.93$, which is attained for
$\bm{b}=(-0.801,0.187)^{\prime}$.

In the likelihood ratio case, one might be tempted to combine a constrained
calculation of the test statistic with reference to the classical $\chi^{2}$
null distribution. This, however, would lead to $p$-values that are too large.
Therefore, the null hypothesis would be accepted too often, resulting in
incorrectly simplifying the covariance structure of the model, which may
seriously invalidate inferences, as shown by Altham (1984).

##  7 Concluding Remarks

Whenever inference for variance components is of interest, the choice between
one-sided or two-sided tests is crucial, depending on whether negative
variance components are deemed meaningful or not. For two-sided tests,
classical inferential procedures can be followed, based on likelihood ratios,
score statistics, or Wald statistics, which have the same asymptotic null
distributions. For one-sided tests, however, one-sided test statistics need to
be developed, and their null distribution derived. In contrast to the case of
likelihood ratio tests, this has thus far not received much attention in the
score-test case. Moreover, there seems to be a lot of confusion as to whether
or not classical score tests are applicable in this setting. Using heuristic
arguments in the context of a simple linear random-effects model, we have
shown why those test statistics are not appropriate for testing one-sided
hypotheses, and how one-sided versions can be obtained. Then, the general
theory of Silvapulle and Silvapulle (1995) has been invoked to derive general
one-sided score tests for variance components. Further, the well-known
equivalence between two-sided score and likelihood ratio tests is shown to
hold true for the one-sided counterparts as well.

In general, likelihood ratio tests as well as score tests are available for
testing hypotheses about variance components, and both procedures are
asymptotically equivalent, for one-sided as well as two-sided tests. The
choice between one-sided and two-sided tests should be entirely driven by the
scientific question, the data analyzed, the models fitted, and the
interpretation of the parameters in those models. A frequently quoted
justification for the use of score tests is that they do not require fitting
the alternative model. However, currently available software easily allows
pracitising statisticians to fit and compare a variety of models, containing
many variance components. Moreover, whenever one-sided tests are of interest,
the score test may require employing numerical optimization techniques for the
calculation of the infimum in (10). Therefore, it cannot be our intention to
advocate the broad use of score tests for the inference on variance
components. Instead, the aim of this paper has been to enhance insight into
the score test and to illustrate the use of score tests in this context.

We hope to have indicated that one either can take an unconstrained view and
then no additional action is needed in case variance components are negative,
or one takes a constrained view and then the inferential procedures should be
such that proper constraints are imposed. In this sense, the statement made by
Brown and Prescott (1999, p. 237): “The usual action when a negative variance
component estimate is obtained for a random coefficient would be to refit the
model with the random coefficient removed (…)”, overlooks important issues and
is therefore misleading.

## Acknowledgements

We gratefully acknowledge support from FWO-Vlaanderen Research Project
“Sensitivity Analysis for Incomplete and Coarse Data” and Belgian IUAP/PAI
network “Statistical Techniques and Modeling for Complex Substantive Questions
with Complex Data”. We are very grateful to David Cox for sharing his insights
with us and for very helpful comments on an earlier version of this
manuscript.

## References

  * Not imporant (2506) Altham, P.M.E. (1984). Improving the precision of estimation by fitting a model. Journal of the Royal Statistical Society, Series B 46, 118–119.
  * Not imporant (2506) Britton, T. (1997). Tests to detect clustering of infected individuals within families. Biometrics 53, 98–109.
  * Not imporant (2506) Brown, H. and Prescott, R. (1999). Applied Mixed Models in Medicine. Chichester: John Wiley.
  * Not imporant (2506) Chernoff, H. (1954). On the distribution of the likelihood ratio. Annals of Matematical Statistics 25, 573–578.
  * Not imporant (2506) Cox, D.R. (1983). Some remarks on overdispersion. Biometrika 70, 269–274.
  * Not imporant (2506) Cox, D.R. and Hinkley, D.V. (1990). Theoretical Statistics. London: Chapman & Hall.
  * Not imporant (2506) Dean, C. (1992). Testing for overdispersion in Poisson and binomial regression models. Journal of the American Statistical Association 87, 451–457.
  * Not imporant (2506) Dean, C. and Lawless, J.F. (1989). Tests for detecting overdispersion in Poisson regression models. Journal of the American Statistical Association 84, 467–472.
  * Not imporant (2506) Dean, C.B., Ugarte, M.D., and Militino, A.F. (2001). Detecting interaction between random region and fixed age effects in disease mapping. Biometrics 57, 192–202.
  * Not imporant (2506) Gueorguieva, R. (2001). A multivariate generalized linear mixed model for joint modelling of clustered outcomes in the exponential family. Statistical Modelling 1, 177–193.
  * Not imporant (2506) Gray, R.J. (1995). Tests for variation over groups in survival data. Journal of the American Statistical Association 90, 198–203.
  * Not imporant (2506) Hall, D.B. and Præstgaard, J.T. (2001). Order-restricted score tests for homogeneity in generalised linear and nonlinear mixed models. Biometrika 88, 739–751.
  * Not imporant (2506) Hines, R.J.O. (1997). A comparison of tests for overdispersion in generalized linear models. Journal of Statistical Computation and Simulation 58, 323–342.
  * Not imporant (2506) Jacqmin-Gadda, H. and Commenges, D. (1995). Tests of homogeneity for generalised linear models. Journal of the American Statistical Association 90, 1237–1246.
  * Not imporant (2506) Laird, N.M. and Ware, J.H. (1982). Random effects models for longitudinal data. Biometrics 38, 963–974.
  * Not imporant (2506) le Cessie, S. and van Houwelingen, J.C. (1995). Testing the fit of a regression model via score tests in random effects models. Biometrics 51, 600–614.
  * Not imporant (2506) Liang, K.–Y. (1987). A locally most powerful test for homogeneity with many strata. Biometrika 74, 259–264.
  * Not imporant (2506) Lin, X. (1997). Variance component testing in generalised linear models with random effects. Biometrika 84, 309–326.
  * Not imporant (2506) Lindstrom, M.J. and Bates, D.M. (1988). Newton-Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association 83, 1014–1022.
  * Not imporant (2506) Lu, W.S. (1997). Score tests for overdispersion in Poisson regression models. Journal of Statistical Computation and Simulation 56, 213–228.
  * Not imporant (2506) Militino, A.F., Ugarte, M.D., and Dean, C.B. (2001). The use of mixture models for identifying high risks in disease mapping. Statistics in Medicine 20, 2035–2049.
  * Not imporant (2506) Miller, J.J. (1977). Asymptotic properties of maximum likelihood estimates in the mixed model of the analysis of variance. The Annals of Statistics 5, 746–762.
  * Not imporant (2506) Nelder, J.A. (1954). The interpretation of negative components of variance. Biometrika 41, 544–548.
  * Not imporant (2506) Paul, S.R. and Islam, A.S. (1995). Analysis of proportions in the presence of over- and under-dispersion. Biometrics 51, 1400–1410.
  * Not imporant (2506) Raubertas, R.F., Lee, C.I.C., and Nordheim, E.V. (1986). Hypothesis tests for normal means constrained by linear inequalities. Communications in Statistics–Theory and Methods 15, 2809–2833.
  * Not imporant (2506) Searle, S.R., Casella, G., and McCulloch, C.E. (1992). Variance Components. New York: John Wiley & Sons.
  * Not imporant (2506) Self, S.G. and Liang, K.Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association 82, 605–610.
  * Not imporant (2506) Shapiro, A. (1988). Towards a unified theory of inequality constrained testing in multivariate analysis. International Statistical Review 56, 49–62.
  * Not imporant (2506) Silvapulle, M.J. and Silvapulle, P. (1995). A score test against one-sided alternatives. Journal of the American Statistical Association 90, 342–349.
  * Not imporant (2506) Smith, P.J. and Heitjan, D.F. (1993). Testing and adjusting for departures from nominal disperson in generalized linear models. Applied Statistics 41, 31–41.
  * Not imporant (2506) Stram, D.O. and Lee, J.W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics 50, 1171–1177.
  * Not imporant (2506) Stram, D.A. and Lee, J.W. (1995). Correction to: Variance components testing in the longitudinal mixed effects model. Biometrics 51, 1196.
  * Not imporant (2506) Thompson, W.A., Jr. (1962). The problem of negative estimates of variance components. Annals of Mathematical Statistics 33, 273–289.
  * Not imporant (2506) Verbeke, G. and Lesaffre, E. (1999). The effect of drop-out on the efficiency of longitudinal experiments. Applied Statistics 48, 363–375.
  * Not imporant (2506) Verbeke, G. and Molenberghs, G. (2000). Linear mixed models for longitudinal data, Springer Series in Statistics, Springer-Verlag, New-York.
  * Not imporant (2506) Verdonck, A., De Ridder, L., Verbeke, G., Bourguignon, J.P., Carels, C., Kuhn, E.R., Darras, V., and de Zegher, F. (1998). Comparative effects of neonatal and prepubertal castration on craniofacial growth in rats. Archives of Oral Biology 43, 861–871.

## Supporting Information

Web Appendix A, referenced in Section 5, is available with this paper at the
Biometrics website on Wiley Online Library.

## Appendix A

###  A.1 Example SAS Code

Within the context of the case study, we will show how all four tests
discussed in Section 5 can be implemented using the SAS procedures MIXED and
NLMIXED (Version 8.2). In particular, we will focus on the comparison between
Models 2 and 4.

Let us consider the likelihood ratio tests first. Fitting Model 4, with random
intercepts and random slopes is typically done using the code (Verbeke and
Molenberghs 2000, Ch. 8):

    proc mixed data = example method = ml;
    class treatment rat;
    model y = treatment*time / solution;
    random intercept time / type = un subject = rat;
    run;

However, it should be noted the SAS procedure MIXED takes by default a
constrained view-point which differs not only from the unconstrained
situation, but is also different from the constrained formulation considered
in this paper. Indeed, the only constraints imposed by SAS are for the
variances to be nonnegative; the resulting covariance matrix can still fail to
be positive semi-definite. A completely unconstrained model can be fitted by
adding the ‘nobound’ option to the PROC MIXED statement. The likelihood ratio
test statistic then follows from fitting Model 2 as well, which is obtained
from deleting the variable time from the RANDOM statement.

Fitting both models under our constraints is done by replacing the ‘type=un’
option in the RANDOM statement of the initial program by ‘type=fa0(2)’ in
Model 4 and ‘type=fa0(1)’ in Model 2.

Recall that, for the unconstrained score test, the first and second order
derivates of the alternative log-likelihood need to be evaluated at the MLE’s
under the null model. We therefore first fit the null model using the SAS
procedure MIXED, then store the parameter estimates in appropriately formatted
SAS data sets, whereafter the derivatives under the alternative model are
obtained using the NLMIXED procedure. The NLMIXED is particularly convenient
since it easily produces the first and second order derivatives at the maximum
likelihood estimates. To this end, the MLE’s under the null are supplied to
NLMIXED as starting values and the procedure is prevented from iterating.
Practically, the results are obtained using the following code.

    proc mixed data = example method = ml nobound;
    class treatmet rat;
    model y = treatment*time / solution;
    random intercept / subject = rat;
    ods output covparms = cov;
    ods output solutionf = fixed;
    run;

    data cov; set cov; keep estimate; run;
    data cov1; input parameter $; cards;
    d11
    sigma2
    ;
    run;
    data cov; merge cov1 cov; run;
    data fixed; set fixed; keep estimate; run;
    data fixed1; input parameter $; cards;
    beta0
    beta1
    beta2
    beta3
    ;
    run;
    data fixed; merge fixed1 fixed; run;
    data starting; set fixed cov; run;

    proc nlmixed data = test start hess qpoints = 50
                 maxiter = 0;
    parms d12 = 0  d22 = 1E-10 / data = starting;
    mean = beta0 + b1 + b2*time + beta1*tcon
           + beta2*thig + beta3*tlow;
    model y ~ normal(mean,sigma2);
    random b1 b2 ~ normal([0, 0],[d11, d12, d22])
           subject = rat;
    run;

For the constrained case, the above program is maintained for the calculation
of the first term in (10), except for the removal of the ‘nobound’ option from
the PROC MIXED statement. To execute the constrained minimization of
$(\bm{Z}_{N}-\bm{b})^{\prime}H_{\psi\psi}^{-1}(\widehat{\mbox{\boldmath$\theta$}}_{H})(\bm{Z}_{N}-\bm{b})$,
numerical optimization routines are often needed. While several routes could
be followed, we have chosen to use the constrained optimization features of
GAUSS. To this end, the constraints need to be written in the form of a
procedure, with a second procedure containing the objective function. In our
example, $\mbox{\boldmath$Z$}_{N}=(0.10219,-6.619967)^{\prime}$ and

$H_{\psi\psi}(\widehat{\mbox{\boldmath$\theta$}}_{H})=\left(\begin{array}[]{cc}16.5355&-9.1388\\\
-9.1388&26.2642\end{array}\right).$

The code then becomes:

    library co;

    proc ineqp(b);
    local d11, d12, d22, d;
    d11=3.4374;
    d12=b[1,.];
    d22=b[2,.];
    d=(d11~d12)|(d12~d22);
    retp(minc(eigh(d))-1e-5);
    endp;

    _co_IneqProc=&ineqp;

    proc factor2(b);
    local zn, hmat;
    zn={0.10219,  -6.619967};
    hmat={16.5355 -9.1388, -9.1388 26.2642};
    retp((zn-b)’*inv(hmat)*(zn-b));
    endp;

    start={-1,2};
    _co_Algorithm=5;
    {x,f,g,ret}=co(&functie,start);


# Product-free sets in approximate subgroups of distal groups

Atticus Stonestrom (Email<EMAIL_ADDRESS>
Abstract: Recall that a subset $X$ of a group $G$ is ‘product-free’ if
$X^{2}\cap X=\varnothing$, ie if $xy\notin X$ for all $x,y\in X$. Let $G$ be a
group definable in a distal structure. We prove there are constants $c>0$ and
$\delta\in(0,1)$ such that every finite subset $X\subseteq G$ distinct from
$\\{1\\}$ contains a product-free subset of size at least
$\delta|X|^{c+1}/|X^{2}|^{c}$. In particular, every finite $k$-approximate
subgroup of $G$ distinct from $\\{1\\}$ contains a product-free subset of
density at least $\delta/k^{c}$.

The proof is short, and follows quickly from Ruzsa calculus and an iterated
application of Chernikov and Starchenko’s distal regularity lemma.

Acknowledgements: I would like to thank Artem Chernikov, Ben Green, Zach
Hunter, and Sergei Starchenko for helpful discussion, and I would especially
like to thank Jonathan Pila and Udi Hrushovski for much encouragement and
feedback.

Notation: Throughout we will use the standard terminology and notation of
‘Ruzsa calculus’; thus a ‘multiplicative set’ is simply a finite subset of a
group. Given multiplicative sets $X,Y$ in the same ambient group, we have
$XY=\\{xy:x\in X,y\in Y\\}$ and $X^{-1}=\\{x^{-1}:x\in X\\}$ and, for any
natural number $n$, $X^{n}=\\{x_{1}\cdots x_{n}:x_{i}\in X\\}$. Also, we take
$\log$ to mean the logarithm base $2$.

##  1 Introduction

Product-free sets – especially in the abelian context, where they are called
‘sum-free’ – are a rich topic of study in additive combinatorics; one of the
earliest results is Erdős’ 1965 proof that every finite set of non-zero
integers contains a sum-free subset of density at least $1/3$. Given a
multiplicative set $X$, two natural questions of interest include the
following: what is the largest density of a product-free subset of $X$, and
how many product-free subsets does $X$ contain? There is an extensive
literature on these and related topics, and we refer the reader to the survey
[20] for a picture of some of this landscape in the abelian case. The starting
point of our motivation here is a question raised by Babai and Sós in [4]: is
there a constant $\delta>0$ such that every finite non-trivial group contains
a product-free subset of density at least $\delta$?

In [10], Gowers gave a negative answer to this question. Indeed, let $H$ be a
non-trivial finite group and $d$ the minimal dimension of a non-trivial
complex representation of $H$. Gowers then shows that $H$ does not contain any
product-free set of size greater than $|H|/d^{1/3}$; by a classical result of
Frobenius, this is enough to conclude that, for every prime power $q$,
$\mathrm{PSL}_{2}(q)$ contains no product-free set of size greater than
$2|\mathrm{PSL}_{2}(q)|^{8/9}$.

Gowers also obtained a bound in the reverse direction, showing that $H$ does
contain a product-free set of size at least $|H|/2000^{d}$. In [14], Nikolov
and Pyber improved the exponential dependence on $d$ to a polynomial one,
showing the existence of an absolute constant $\delta\in(0,1)$ such that $H$
contains a product-free set of size at least $\delta|H|/d$.

These results (and others – see eg [12] for a survey, and see [15] for pseudo-
finite analogues due to Palacín) give perspective on the maximal densities of
product-free subsets of finite groups. In this paper we are motivated by the
natural generalization of this line of inquiry to the class of finite
approximate groups. We show that, for fixed $k$, a positive answer to the
analogue of Babai and Sós’ question for $k$-approximate groups can be obtained
in a certain model-theoretic context: namely, in the context of a group
interpretable in a distal structure. Specifically, given such a group $G$, we
show the existence of constants $c>0$ and $\delta\in(0,1)$ such that every
finite subset $X\subseteq G$ distinct from $\\{1\\}$ contains a product-free
subset of density at least $\delta/k^{c}$, where $k=|X^{2}|/|X|$. In the case
$G=\mathrm{GL}_{n}(\mathbb{C})$, we note that, modulo a weighted version of an
argument of Alon and Kleitman, this also follows from Breuillard, Green, and
Tao’s classification of the approximate subgroups of
$\mathrm{GL}_{n}(\mathbb{C})$ (again with polynomial dependence on $k$).

##  2 Preliminaries

###  2.1 Approximate groups

Recall that the ‘doubling’ and ‘tripling’ constants of a multiplicative set
$X$ are the ratios $|X^{2}|/|X|$ and $|X^{3}|/|X|$, respectively. For
$k\geqslant 1$, $X$ is said to be a $k$-approximate group if $1\in X=X^{-1}$
and $X^{2}$ is covered by at most $k$-many translates of $X$; note that the
doubling constant of a $k$-approximate group is at most $k$. There is a vast
literature on $k$-approximate groups, and we refer the reader to the survey
[5] for context on the area. For our main result, we will need only two basic
facts. The first is a consequence of Ruzsa’s ‘triangle inequality’:

######  Fact 2.1.

If $X$ is a multiplicative set with $|X^{2}|\leqslant k|X|$, then
$|XX^{-1}|\leqslant k^{2}|X|$.

The second is a consequence of a theorem of Petridis; see Theorem 1.5 in [16]:

######  Fact 2.2.

If $X$ is a multiplicative set with $|X^{2}|\leqslant k|X|$, then there is a
subset $Y\subseteq X$ such that $|Y|\geqslant|X|/k$ and $|Y^{3}|\leqslant
k^{3}|Y|$.

In some follow-up discussion, it will also be convenient to cite the following
results; the first is a consequence of Theorem 4.6 in [19].

######  Fact 2.3.

Let $Y,Z$ be multiplicative sets of size $n$, and suppose $|YZ|\leqslant kn$.
Then, in the ambient group, there is an $O(k^{O(1)})$-approximate group $W$
and an element $u$ such that $|W|\leqslant O(k^{O(1)})n$ and $|Y\cap
uW|\geqslant|Y|/O(k^{O(1)})$.

The second is a consequence of the celebrated Balog-Szemerédi-Gowers theorem;
see for example Lemma 5.1 in [19].

######  Fact 2.4.

Let $X$ be a multiplicative set, and suppose there are at least $|X|^{2}/k$
many pairs $(x,y)\in X\times X$ satisfying $xy\in X$. Then there are subsets
$Y,Z\subseteq X$, each of size at least $|X|/O(k^{O(1)})$, such that
$|YZ|\leqslant O(k^{O(1)})|X|$.

###  2.2 Distal structures

The model-theoretic context for this paper is the notion of distality, a
property of first-order structures introduced by Simon in [17]. We will not
need the definition of a distal structure here, and refer the reader to the
introduction of [7] for an overview. Instead, we point out two motivating
examples: (i) the field $(\mathbb{R},\cdot,+)$, or more generally any
o-minimal structure, and (ii) the field $(\mathbb{Q}_{p},\cdot,+)$ for any
prime $p$.

In this note we are concerned with groups $(G,\cdot)$ that are definable in
some distal structure. One non-trivial fact is that these are precisely the
groups interpretable in a distal structure; see [13] or [3] for detailed
proofs of this. (As an aside, note that it also follows directly from Fact 2.7
below alone that every relation interpretable in a distal structure has the
strong Erdős-Hajnal property.)

As a concrete example, the reader may throughout consider $G$ to be the
$\mathbb{R}$-points of any semialgebraic group.

The main tool we employ comes from the paper [7], wherein Chernikov and
Starchenko showed that relations definable in distal structures enjoy strong
combinatorial regularity properties. To state their theorem, we first recall
the notion of a ‘homogeneous’ tuple of sets in an $n$-partite hypergraph; let
$R\subseteq X_{1}\times\dots\times X_{n}$ be an $n$-ary relation on sets
$X_{i}$.

######  Definition 2.5.

A tuple $(U_{1},\dots,U_{n})$ of subsets $U_{i}\subseteq X_{i}$ is
‘$R$-homogeneous’ if $U_{1}\times\dots\times U_{n}$ is either a subset of $R$
or disjoint from $R$.

######  Definition 2.6.

$R$ has the ‘strong Erdős-Hajnal property’ if there exists a constant
$\delta\in(0,1)$ with the following property: for any finite subsets
$W_{i}\subseteq X_{i}$, there are subsets $U_{i}\subseteq W_{i}$, with each
$U_{i}$ of size at least $\delta|W_{i}|$, such that $(U_{1},\dots,U_{n})$ is
$R$-homogeneous.

The relevant result from [7] is then the following, which is an immediate
consequence of Corollary 4.5 there:

######  Fact 2.7.

Every relation definable in a distal structure has the strong Erdős-Hajnal
property.

For some examples of the varied interactions between distality and additive
combinatorics, we refer the reader to the papers [2], [8], and [9].

##  3 Main Result

We now have all the necessary prerequisites. Our result holds in any group in
which the $6$-ary relation $R(\overline{x},\overline{y})\equiv
x_{1}x_{2}x_{3}=y_{1}y_{2}y_{3}$ has the strong Erdős-Hajnal property; by Fact
2.7, it applies in particular to any group definable in a distal structure.
For the rest of this section thus fix a group $G$ in which $R$ has the strong
Erdős-Hajnal property.

######  Lemma 3.1.

There is a constant $c_{0}>0$ such that, for any $\alpha\in(0,1)$, there is
some $\varepsilon\in(0,1)$ with the following property: every finite subset
$Y\subseteq G$ of size at least $8/\alpha$ contains subsets $U,V,W\subseteq
Y$, all of size at least $\varepsilon|Y|^{c_{0}+1}/|Y^{3}|^{c_{0}}$, such that
$|UVW|\leqslant\alpha|Y|$.

###### Proof.

Let $\delta\in(0,1)$ be a constant witnessing the strong Erdős-Hajnal property
for $R$; ie any finite subsets $U_{1},U_{2},U_{3},V_{1},V_{2},V_{3}\subseteq
G$ contain respective subsets
$U^{\prime}_{1},U^{\prime}_{2},U^{\prime}_{3},V^{\prime}_{1},V^{\prime}_{2},V^{\prime}_{3}$,
all of density at least $\delta$, such that
$(U^{\prime}_{1},\dots,V^{\prime}_{3})$ is $R$-homogeneous. I claim we may
take $c_{0}=\log(1/\delta)$; to see this, fix $\alpha\in(0,1)$ and let
$\varepsilon=\delta\alpha^{c_{0}}$.

First note that if $U_{1},U_{2},U_{3},V_{1},V_{2},V_{3}$ are subsets of $G$ of
size at least $2$, then we cannot have $u_{1}u_{2}u_{3}=v_{1}v_{2}v_{3}$ for
all $u_{i}\in U_{i}$ and $v_{i}\in V_{i}$. Thus if $(U_{1},\dots,V_{3})$ is
$R$-homogeneous, then in fact $U_{1}U_{2}U_{3}$ and $V_{1}V_{2}V_{3}$ are
disjoint, and so in particular one of $U_{1}U_{2}U_{3}$ and $V_{1}V_{2}V_{3}$
has size at most $|W_{1}W_{2}W_{3}|/2$ for any sets $W_{i}\supseteq U_{i}\cup
V_{i}$. The strong Erdős-Hajnal property hence tells us that any finite
subsets $W_{1},W_{2},W_{3}\subseteq G$, all of size at least $2$, contain
respective subsets $U_{1},U_{2},U_{3}$ with $|U_{i}|\geqslant\delta|W_{i}|$
for each $i$ and with $|U_{1}U_{2}U_{3}|\leqslant|W_{1}W_{2}W_{3}|/2$.

Now we simply iterate this fact; fix a finite subset $Y\subseteq G$ with
$|Y|\geqslant 8/\alpha$ and let $U_{1}=V_{1}=W_{1}=Y$. Applying the previous
remark, we inductively define sets $U_{n},V_{n},W_{n}$ such that (i)
$U_{n+1},V_{n+1},W_{n+1}$ are contained in $U_{n},V_{n},W_{n}$, respectively,
(ii) $U_{n+1},V_{n+1},W_{n+1}$ have size at least
$\delta|U_{n}|,\delta|V_{n}|,\delta|W_{n}|$, respectively, and (iii)
$|U_{n+1}V_{n+1}W_{n+1}|\leqslant|U_{n}V_{n}W_{n}|/2$; the $(n+1)$-th tuple of
sets in this sequence can be constructed as long as all of $U_{n},V_{n},W_{n}$
have size at least $2$, and hence as long as $\delta^{n-1}|Y|\geqslant 2$.

Now let $n=\lceil\log(|Y^{3}|/\alpha|Y|)\rceil+1$. Then

$\delta^{n-2}\geqslant\delta^{\log(|Y^{3}|/\alpha|Y|)}=(\alpha|Y|/|Y^{3}|)^{\log(1/\delta)}=(\alpha|Y|/|Y^{3}|)^{c_{0}},$

whence $\delta^{n-1}|Y|\geqslant\varepsilon|Y|^{c_{0}+1}/|Y^{3}|^{c_{0}}$, and
so it suffices to find $U,V,W\subseteq Y$ of size at least $\delta^{n-1}|Y|$
and with $|UVW|\leqslant\alpha|Y|$. If $\delta^{n-2}|Y|\geqslant 2$, then the
sets $U_{n},V_{n},W_{n}$ exist and are of appropriate size, and further
satisfy $|U_{n}V_{n}W_{n}|\leqslant|Y^{3}|/2^{n-1}\leqslant\alpha|Y|$, as
needed. If instead $2>\delta^{n-2}|Y|$, then in particular
$2>\delta^{n-1}|Y|$, and so we may take $U,V,W$ to be any subsets of $Y$ of
size at least $2$; these will satisfy $|UVW|\leqslant 8\leqslant\alpha|Y|$,
again as needed. ∎

######  Corollary 3.2.

There are constants $c_{1}>0$ and $\varepsilon_{1}\in(0,1)$ such that every
finite subset $Y\subseteq G$ of size at least $16$ contains a subset
$Z\subseteq Y$, of size at least
$\varepsilon_{1}|Y|^{c_{1}+1}/|Y^{3}|^{c_{1}}$, such that
$|Z^{-1}ZZ^{-1}|\leqslant|Y|/2$.

###### Proof.

Let $c_{0}$ be given by Lemma 3.1, and let $\varepsilon_{0}\in(0,1)$ witness
the lemma for the case $\alpha=1/2$; I claim we may take $c_{1}=3c_{0}$ and
$\varepsilon_{1}=4\varepsilon_{0}^{3}$. To see this, fix a finite subset
$Y\subseteq G$ of size at least $16$, and for notational convenience let
$k=|Y^{3}|/|Y|$. By Lemma 3.1 we can find subsets $U,V,W\subseteq Y$, each of
size at least $\varepsilon_{0}|Y|/k^{c_{0}}$, such that $|UVW|\leqslant|Y|/2$.

Now, for each $g,h\in G$, let $Z_{g,h}=U^{-1}g\cap V\cap hW^{-1}$. The map
$(g,h,z)\mapsto(gz^{-1},z,z^{-1}h)$ gives a bijection from $\\{(g,h,z):g,h\in
G,z\in Z_{g,h}\\}$ to $U\times V\times W$, with inverse given by
$(u,v,w)\mapsto(uv,vw,v)$, whence $\sum_{g,h\in G}|Z_{g,h}|=|U||V||W|$. On the
other hand, $Z_{g,h}\neq\varnothing$ only if $g\in UV$ and $h\in VW$. The size
of these sets is bounded above by $|UVW|\leqslant|Y|/2$, so we have
$|U||V||W|\leqslant(|Y|^{2}/4)\sup_{g,h\in G}|Z_{g,h}|$, and there thus exist
some $g,h\in G$ such that

$|Z_{g,h}|\geqslant 4|U||V||W|/|Y|^{2}\geqslant
4\varepsilon_{0}^{3}|Y|/k^{3c_{0}}=\varepsilon_{1}|Y|^{c_{1}+1}/|Y^{3}|^{c_{1}}.$

On the other hand, $Z_{g,h}^{-1}Z_{g,h}Z_{g,h}^{-1}\subseteq g^{-1}UVWh^{-1}$,
whence $|Z_{g,h}^{-1}Z_{g,h}Z_{g,h}^{-1}|\leqslant|UVW|\leqslant|Y|/2$. So
taking $Z=Z_{g,h}$ gives the desired result. ∎

######  Theorem 3.3.

There are constants $c_{2}>0$ and $\varepsilon_{2}\in(0,1)$ such that the
following holds: if $X\subseteq G$ is finite and distinct from $\\{1\\}$, and
$k=|X^{2}|/|X|$, then $X$ contains a product-free subset of size at least
$\varepsilon_{2}|X|/k^{c_{2}}$.

###### Proof.

Let $c_{1}$ and $\varepsilon_{1}$ be given by Corollary 3.2; I claim we may
take $c_{2}=3c_{1}+4$ and $\varepsilon_{2}=\min\\{\varepsilon_{1}/2,1/16\\}$.
Thus fix a finite subset $X\subseteq G$ distinct from $\\{1\\}$. If $|X|<16k$,
then $\varepsilon_{2}|X|/k^{c_{2}}<16k/16k^{c_{2}}\leqslant 1$, so we may pick
any $z\in X\setminus\\{1\\}$ and then $\\{z\\}$ will be a product-free subset
of $X$ of appropriate size. Hence we may assume $|X|\geqslant 16k$. By Fact
2.2, there is a subset $Y\subseteq X$ with $|Y|\geqslant|X|/k$, hence
$|Y|\geqslant 16$, and $|Y^{3}|\leqslant k^{3}|Y|$. By Corollary 3.2, there is
then a subset $Z\subseteq Y$ with $|Z|\geqslant\varepsilon_{1}|Y|/k^{3c_{1}}$
and $|Z^{-1}ZZ^{-1}|\leqslant|Y|/2$.

For each $g\in G$, let $Y_{g}=gZ\cap Y$. The map $(g,y)\mapsto(y,g^{-1}y)$
gives a bijection from $\\{(g,y):g\in G,y\in Y_{g}\\}$ to $Y\times Z$, whence
$\sum_{g\in G}|Y_{g}|=|Y||Z|$. Moreover, $|Y_{g}|\leqslant|gZ|=|Z|$ for all
$g\in G$, and $|Z^{-1}ZZ^{-1}|\leqslant|Y|/2$, so $\sum_{g\in
Z^{-1}ZZ^{-1}}|Y_{g}|\leqslant|Y||Z|/2$; combining these bounds gives
$\sum_{g\notin Z^{-1}ZZ^{-1}}|Y_{g}|\geqslant|Y||Z|/2$.

On the other hand, $Y_{g}\neq\varnothing$ only if $g\in YZ^{-1}$. Since
$Y,Z\subseteq X$ and $|X^{2}|\leqslant k|X|$, by Fact 2.1 we have
$|YZ^{-1}|\leqslant|XX^{-1}|\leqslant k^{2}|X|\leqslant k^{3}|Y|$. So in
particular

$|Y||Z|/2\leqslant k^{3}|Y|\sup_{g\notin Z^{-1}ZZ^{-1}}|Y_{g}|,$

and there thus exists some $g\in G\setminus Z^{-1}ZZ^{-1}$ such that
$|Y_{g}|\geqslant|Z|/2k^{3}$.

By definition, $Y_{g}\subseteq Y\subseteq X$. We further have
$|Y_{g}|\geqslant\varepsilon_{1}|Y|/2k^{3c_{1}+3}$ since
$|Z|\geqslant\varepsilon_{1}|Y|/k^{3c_{1}}$, and therefore

$|Y_{g}|\geqslant\varepsilon_{1}|X|/2k^{3c_{1}+4}\geqslant\varepsilon_{2}|X|/k^{c_{2}}$

since $|Y|\geqslant|X|/k$. So we need only show that $Y_{g}$ is product-free
and we will be done; for this, since $Y_{g}\subseteq gZ$, it suffices to show
that $gZ$ is product-free. Thus suppose otherwise; then there are
$z_{1},z_{2},z_{3}\in Z$ such that $gz_{1}gz_{2}=gz_{3}$. Rearranging this
gives $g=z_{1}^{-1}z_{3}z_{2}^{-1}$, contradicting that $g\notin
Z^{-1}ZZ^{-1}$, and so we are done. ∎

This concludes the main result. For completeness, we would like to note an
alternative proof of Theorem 3.3 in the case that
$G=\mathrm{GL}_{n}(\mathbb{C})$. The $k$-approximate subgroups of $G$ were
classified qualitatively by Hrushovski in [11], and then with polynomial
dependence on $k$ by Breuillard, Green, and Tao in [6]; the following is an
immediate consequence of Theorem 2.5 in [6].

######  Fact 3.4.

For every $n$, there are constants $c>0$ and $\varepsilon\in(0,1)$ such that
the following holds: for any finite $k$-approximate group
$X\subseteq\mathrm{GL}_{n}(\mathbb{C})$, there is a nilpotent group
$H\leqslant\mathrm{GL}_{n}(\mathbb{C})$ of step at most $n-1$ and an element
$g\in\mathrm{GL}_{n}(\mathbb{C})$ such that $|X\cap
gH|\geqslant\varepsilon|X|/k^{c}$.

By Fact 2.3 and Corollary A.3, and noting that any non-trivial coset of a
subgroup is product-free, Fact 3.4 implies Theorem 3.3 in the case
$G=\mathrm{GL}_{n}(\mathbb{C})$, with the same polynomial dependence on $k$.
It is tempting to ask whether the dependence on $k$ can be removed entirely,
either in this case or in the more general distal context:

######  Question 3.5.

Suppose $G$ is a group definable in a distal structure. Is there a constant
$\delta>0$ such that every finite subset $X\subseteq G$ distinct from
$\\{1\\}$ contains a product-free subset of size at least $\delta|X|$? How
about in the case $G=\mathrm{GL}_{n}(\mathbb{C})$?

For completeness we include some remarks on this question in the next section.

##  4 Remarks on Question 3.5

In this section we will briefly discuss Question 3.5; these are standard
observations but perhaps worth noting nonetheless. Let $G$ be a group
definable in a distal structure and let $X\subseteq G$ be a finite subset
distinct from $\\{1\\}$; we wish to find a product-free subset of $X$ of fixed
positive density. On the one hand, if there are ‘very few’ pairs $(x,y)\in
X\times X$ with $xy\in X$, say at most $k|X|$-many for some fixed $k\geqslant
1$, then one can check that $X$ contains a product-free subset of density at
least $1/O(k)$; indeed this is true without any hypotheses on $G$ whatsoever.

On the other hand, if ‘most’ of the pairs $(x,y)\in X\times X$ satisfy $xy\in
X$, then we can couple Theorem 3.3 with the Balog-Szemerédi-Gowers theorem to
obtain the desired set. Specifically, if there are at least $|X|^{2}/k$ such
pairs, then $X$ contains a product-free subset of density at least
$1/O(k^{O(1)})$; let us quickly sketch the proof of this.

The first relevant observation is that every dense subset of a ‘coset’ of an
approximate group contains a dense product-free subset. More precisely,
suppose that $W\subseteq G$ is a finite $k$-approximate group, that $u\in G$
is arbitrary, and that $V$ is a subset of $uW$ distinct from $\\{1\\}$ and of
size $|W|/O(k^{O(1)})$. We have two cases; if $u\notin W^{-1}WW^{-1}$, then
$uW$ and hence $V$ as well are already product-free, so there is nothing to
show. If instead $u$ does lie in $W^{-1}WW^{-1}=W^{3}$, then $(uW)^{2}$ is
contained in $W^{8}$, whence $uW$ has doubling at most $k^{7}$. Since
$|uW|/|V|=O(k^{O(1)})$, thus $V$ has doubling $O(k^{O(1)})$, and now by
Theorem 3.3 $V$ contains a product-free subset of density $1/O(k^{O(1)})$, as
needed.

Now fix $k\geqslant 1$, and suppose there are at least $|X|^{2}/k$ pairs
$(x,y)\in X\times X$ with $xy\in X$. By the Fact 2.4, there are then subsets
$Y,Z\subseteq X$, each of size $|X|/O(k^{O(1)})$, such that
$|YZ|=O(k^{O(1)})|X|$. By Fact 2.3, there is then an $O(k^{O(1)})$-approximate
subgroup $W\subseteq G$ of size at most $O(k^{O(1)})|X|$ and an element $u\in
G$ such that $|Y\cap uW|$ has size at least $|X|/O(k^{O(1)})$; now one applies
the previous paragraph to $V=Y\cap uW$ to obtain a product-free subset of $Y$
of size $|X|/O(k^{O(1)})$, giving the desired result.

So, we are able to handle the two edge cases, where either ‘very many’ or
‘very few’ pairs of elements of $X$ have product landing in $X$. However, it
is entirely unclear to us how to handle the more general question, and we thus
close our remarks here.

##  Appendix A Product-free sets in solvable groups

In [1], Alon and Kleitman showed that every finite set of non-zero elements of
an abelian group contains a sum-free subset of density at least $2/7$. In this
section we point out an analogous fact for solvable groups, which follows from
a weighted version of Alon and Kleitman’s theorem. The proof is a completely
standard adaptation of their argument.

First we need the following easy fact, which is observed in [1]; let
$\alpha=1/4$.

######  Fact A.1.

For every non-zero finite cyclic group $(H,+)$, there is a sum-free subset
$I\subseteq H$ such that $|K\cap I|\geqslant\alpha|K|$ for every non-zero
subgroup $K\leqslant H$. (Indeed when $H=\mathbb{Z}/n\mathbb{Z}$ one may take
$I$ to be the image of $\\{\lfloor n/3\rfloor+1,\dots,\lfloor 2n/3\rfloor\\}$
in $H$.)

Now the following is the weighted version of Alon and Kleitman’s result, using
an identical argument as theirs:

######  Lemma A.2.

Let $(G,+,0)$ be an abelian group and $B\subseteq G\setminus\\{0\\}$ a finite
subset, and for each $b\in B$ let $w_{b}$ be a positive integer. Then there is
a sum-free subset $A\subseteq B$ with $\sum_{a\in
A}w_{a}\geqslant\alpha\sum_{b\in B}w_{b}$.

###### Proof.

Every finite subset of an abelian group is Freiman isomorphic to a subset of
$(\mathbb{Z}/n\mathbb{Z})^{s}$ for some $n,s$, so we may assume without loss
of generality that $G=(\mathbb{Z}/n\mathbb{Z})^{s}$. Let
$I\subseteq\mathbb{Z}/n\mathbb{Z}$ be given for $H=\mathbb{Z}/n\mathbb{Z}$ by
Fact A.1.

Now, for each $b=(b_{1},\dots,b_{s})\in B$, consider the homomorphism
$f_{b}:G\to\mathbb{Z}/n\mathbb{Z}$ defined by
$(c_{1},\dots,c_{s})\mapsto\sum_{i\in[s]}c_{i}b_{i}$. Note that, since
$(0,\dots,0)\notin B$, $\operatorname{im}(f_{b})$ is a non-zero subgroup of
$\mathbb{Z}/n\mathbb{Z}$ for every $b\in B$. Now, pick $c\in G$ uniformly at
random. Consider the random quantity $W_{c}=\sum_{b\in B}w_{b}1_{f_{b}(c)\in
I}$. Since $f_{b}$ is a group homomorphism, its fibers have constant size, and
so $\mathbb{P}[f_{b}(c)\in I]=|\operatorname{im}(f_{b})\cap
I|/|\operatorname{im}(f_{b})|$ for each $b\in B$. Since also
$\operatorname{im}(f_{b})$ is non-zero, by Fact A.1 $\mathbb{P}[f_{b}(c)\in
I]\geqslant\alpha$ for each $b\in B$, and so in particular
$\mathbb{E}[W_{c}]\geqslant\alpha\sum_{b\in B}w_{b}$.

Thus we may find $c_{0}\in G$ with $W_{c_{0}}\geqslant\alpha\sum_{b\in
B}w_{b}$. Now, let $A=\\{a\in B:f_{a}(c_{0})\in I\\}$. Then $\sum_{a\in
A}w_{a}=W_{c_{0}}\geqslant\alpha\sum_{b\in B}w_{b}$, and I claim that $A$ is
sum-free. Indeed, suppose there are $a,a^{\prime}\in A$ with $a+a^{\prime}\in
A$. Note that $f_{a}(c_{0})+f_{a^{\prime}}(c_{0})=f_{a+a^{\prime}}(c_{0})$; by
definition of $A$ all three of these quantities lie in $I$, contradicting that
$I$ is sum-free. ∎

######  Corollary A.3.

Let $(G,\cdot,1)$ be a solvable group admitting a length-$(n+1)$ subnormal
series with abelian factor groups. Then every finite subset $C\subseteq
G\setminus\\{1\\}$ contains a product-free subset of size at least
$\alpha|C|/2^{n}$.

###### Proof.

By induction on $n$. When $n=0$, $G$ is abelian, and so the base case is
handled. For the inductive step, suppose $n>0$. Then there is a normal
subgroup $H\triangleleft G$ such that $G/H$ is abelian and $H$ admits a
length-$n$ subnormal series with abelian factor groups. Let $C\subseteq G$ be
any finite subset. If $|C\cap H|\geqslant|C|/2$, then we are done by the
inductive hypothesis. Thus suppose $|C\setminus H|\geqslant|C|/2$. Let
$\pi:G\to G/H$ denote the projection map, and let $B=\pi(C\setminus H)$. For
each $b\in B$, let $w_{b}=|\pi^{-1}(b)\cap C|$; by Lemma A.2, there is a
product-free subset $A\subseteq B$ with $\sum_{a\in
A}w_{a}\geqslant\alpha\sum_{b\in B}w_{b}$. But $\sum_{a\in
A}w_{a}=|\pi^{-1}(A)\cap C|$, and $\sum_{b\in B}w_{b}=|C\setminus
H|\geqslant|C|/2$. So $|\pi^{-1}(A)\cap C|\geqslant\alpha|C|/2$. The preimage
of a product-free set under a group homomorphism is product-free, and so
taking $\pi^{-1}(A)\cap C$ gives the desired product-free set. ∎

## References

  * [1] Noga Alon and Daniel Kleitman. Sum-free subsets, in: A Tribute to Paul Erdős, Cambridge University Press, 1990.
  * [2] Aaron Anderson. Combinatorial bounds in distal structures. Preprint. 2021.
  * [3] Matthias Aschenbrenner, Artem Chernikov, Allen Gehret, Martin Ziegler. Distality in valued fields and related structures. Transactions of the AMS, 2022.
  * [4] László Babai, Vera Sós. Sidon sets in groups and induced subgraphs of Cayley graphs. European Journal of Combinatorics, 1985.
  * [5] Emmanuel Breuillard. A brief introduction to approximate groups, in: Thin Groups and Superstrong Approximation, Cambridge University Press, 2013.
  * [6] Emmanuel Breuillard, Ben Green, Terence Tao. Approximate subgroups of linear groups. Geometric and Functional Analysis, Springer, 2011.
  * [7] Artem Chernikov and Sergei Starchenko. Regularity lemma for distal structures. Journal of the European Mathematical Society, 2018.
  * [8] Artem Chernikov, David Galvin, Sergei Starchenko. Cutting lemma and Zarankiewicz’s problem in distal structures. Selecta Mathematica, 2020.
  * [9] Gabe Conant, Anand Pillay, Caroline Terry. Structure and regularity for subsets of groups with finite VC-dimension. Journal of the EMS, 2022.
  * [10] Timothy Gowers. Quasirandom groups. Combinatorics, Probability, and Computing, 2008.
  * [11] Ehud Hrushovski. Stable group theory and approximate subgroups. Journal of the American Mathematical Society, 2011.
  * [12] Kiran Kedlaya. Product-free subsets of groups, then and now. Contemporary Mathematics, 2009.
  * [13] Krzysztof Krupiński and Adrián Portillo. On stable quotients. Notre Dame Journal of Formal Logic, 2022.
  * [14] Nikolay Nikolov, László Pyber. Product decompositions of quasirandom groups and a Jordan type theorem. Journal of the EMS, 2011.
  * [15] Daniel Palacín. On compactifications and product-free sets. Journal of the LMS, 2019.
  * [16] Giorgis Petridis. New proofs of Plünnecke-type estimates for product sets in groups. Combinatorica, 2012.
  * [17] Pierre Simon. Distal and non-distal NIP theories. Annals of Pure and Applied Logic, 2013.
  * [18] Pierre Simon. A note on “Regularity lemma for distal structures”. Proceedings of the American Mathematical Society, 2016.
  * [19] Terence Tao. Product set estimates for non-commutative groups. Combinatorica, 2008.
  * [20] Terence Tao, Van Vu. Sum-free sets in groups: a survey. Journal of Combinatorics, 2017.


#  LABORATORI NAZIONALI DI FRASCATI SIS-Pubblicazioni
LNF–08/27(P) November 29, 2008
SENSITIVITY AND ENVIRONMENTAL RESPONSE OF THE CMS RPC GAS GAIN MONITORING
SYSTEM

L. Benussi1, S. Bianco1, S. Colafranceschi${}^{1},2,3$, F. L. Fabbri1, M.
Giardoni1
B. Ortenzi1, A. Paolozzi,${}^{1},2$ L. Passamonti1, D. Pierluigi1
B. Ponzio1, A. Russo1, A. Colaleo4, F. Loddo4, M. Maggi4
A. Ranieri4, M. Abbrescia${}^{4},5$, G. Iaselli${}^{4},5$, B.
Marangelli${}^{4},5$, S. Natali${}^{4},5$
S. Nuzzo${}^{4},5$, G.Pugliese${}^{4},5$, F. Romano${}^{4},5$, G.
Roselli${}^{4},5$, R. Trentadue${}^{4},5$
S. Tupputi${}^{4},5$, R. Guida3, G. Polese${}^{3},6$, A. Sharma3, A.
Cimmino${}^{7},8$
D. Lomidze8, D. Paolucci8, P. Piccolo8, P. Baesso9, M. Necchi9
D. Pagano9, S. P. Ratti9, P. Vitulo9, C. Viviani9
1) INFN Laboratori Nazionali di Frascati, Via E. Fermi 40, I-00044 Frascati,
Italy.
2) Università degli Studi di Roma ”La Sapienza”, Piazzale A. Moro.
3) CERN CH-1211 Genéve 23 F-01631 Switzerland.
4) INFN Sezione di Bari, Via Amendola, 173I-70126 Bari, Italy.
5) Dipartimento Interateneo di Fisica, Via Amendola, 173I-70126 Bari, Italy.
6) Lappeenranta University of Technology, P.O. Box 20 FI-538 1 Lappeenranta,
Finland.
7) INFN Sezione di Napoli, Complesso Universitario di Monte Sant’Angelo,
edificio 6, 80126 Napoli, Italy.
8) Università di Napoli Federico II, Complesso Universitario di Monte
Sant’Angelo, edificio 6, 80126 Napoli, Italy.
9) INFN Sezione di Pavia, Via Bassi 6, 27100 Pavia, Italy and Università degli
studi di Pavia, Via Bassi 6, 27100 Pavia, Italy.

Results from the gas gain monitoring (GGM) system for the muon detector using
RPC in the CMS experiment at the LHC is presented. The system is designed to
provide fast and accurate determination of any shift in the working point of
the chambers due to gas mixture changes.

PACS: 07.77.Ka; 95.55.Vj; 29.40.Cs

Presented by S. Colafranceschi at the IEEE 08 - 23 October 2008, Dresden,
Germany

##  1 Introduction

Resistive Plate Chambers (RPC) detectors are widely used in HEP experiments
for muon detection and triggering at high-energy, high-luminosity hadron
colliders, in astroparticle physics experiments for the detection of extended
air showers, as well as in medical and imaging applications. At the LHC, muon
systems of the CMS experiment rely on Drift Tubes (DT), Cathode Strip Chambers
(CSC) and RPCs for their muon trigger system, with a total gas volume of about
50 m3. Utmost attention has to be paid to the possible presence of gas
contaminants which degrade the chamber performance. The gas gain monitoring
(GGM) system monitors the gas quality online and is based on small RPC
detectors. The working point - gain and efficiency - is continuously monitored
along with environmental parameters, such as temperature, pressure and
humidity, which are important for the operation of the muon detector system.
Design parameters, construction, prototyping and preliminary commissioning
results of the CMS RPC Gas Gain Monitoring (GGM) system have been presented
previously [1],[7]. In this paper, results on the response of the GGM
detectors to environmental changes are presented.

The CMS RPCs are bakelite-based double-gap RPC with strip readout (for
construction details see [2] and reference therein) operated with 96.2% C2H2F4
\- 3.5% Iso-C4H10 \- 0.3% SF6 gas mixture humidified at about 40%. The large
volume of the whole CMS RPC system and the cost of gas used make mandatory the
operation of RPC in a closed-loop gas system (for a complete description see
[3]), in which the gas fluxing the gaps is reused after being purified by a
set of filters[4].

The operation of the CMS RPC system is strictly correlated to the ratio
between the gas mixture components, and to the presence of pollution due to
contaminants that can be be produced inside the gaps during discharges (i.e.
HF produced by SF6 or C2H2F4 molecular break-up and further fluorine
recombination), accumulated in the closed-loop or by pollution that can be
present in the gas piping system (tubes, valves, filters, bubblers, etc.) and
flushed into the gaps by the gas flow. The monitoring of the presence of these
contaminants, as well as the gas mixture stability, is therefore mandatory to
avoid RPC damage and to ensure their correct functionality.

A monitoring system of the RPC working point due to changes of gas composition
and pollution must provide a faster and sensitive response than the CMS RPC
system itself in order to avoid irreversible damage of the whole system. Such
a Gas Gain Monitoring system monitors efficiency and signal charge
continuously by means of a cosmic ray telescope based on RPC detectors. In the
following will be briefly described the final setup of the GGM system, and the
first results obtained during its commissioning at the ISR test area (CERN).

##  2 The Gas Gain Monitoring System

The GGM system is composed by the same type of RPC used in the CMS detector
but of smaller size (2mm Bakelite gaps, 50$\times$50 cm2). Twelve gaps are
arranged in a stack located in the CMS gas area (SGX5 building) in the
surface, close to CMS assembly hall (LHC-P5). The choice to install the system
in the surface instead of underground allows one to profit from maximum cosmic
muon rates. In order to ensure a fast response to working point shifts with a
precision of 1%, $10^{4}$ events are are required, corresponding to about 30
minutes exposure time on surface, to be compared with a 100-fold lesser rate
underground. The trigger is provided by four out of twelve gaps of the stack,
while the remaining eight gaps are used to monitor the working point
stability.

The eight gaps are arranged in three sub-system: one sub-system (two gaps) is
fluxed with the fresh CMS mixture and its output sent to vent. The second sub-
system (three gaps) is fluxed with CMS gas coming from the closed-loop gas
system and extracted before the gas purifiers, while the third sub-system
(three gaps) is operated with CMS gas extracted from the closed-loop extracted
after the gas filters. The basic idea is to compare the operation of the three
sub-systems and, if some changes are observed, to send a warning to the
experiment. In this way, the gas going to and coming from the CMS RPC detector
is monitored by using the two gaps fluxed with the fresh mixture as reference
gaps. This setup will ensure that pressure, temperature and humidity changes
affecting the gaps behavior do cancel out by comparing the response of the
three sub-system operating in the same ambient condition.

The monitoring is performed by measuring the charge distributions of each
chamber. The eight gaps are operated at different high voltages, fixed for
each chamber, in order to monitor the total range of operating modes of the
gaps. The operation mode of the RPC changes as a function of the voltage
applied. A fraction of the eight gaps will work in pure avalanche mode, while
the remaining will be operated in avalanche+streamer mode. Comparison of
signal charge distributions and the ratio of the avalanche to streamer
components of the ADC provides a monitoring of the stability of working point
for changes due to gas mixture variations.

Details on the construction of GGM can be found in [7]. Each chamber of the
GGM system consists of a single gap with double sided pad read-out: two copper
pads are glued on the two opposite external side of the gap. The signal is
read-out by a transformer based circuit A3 (Fig.1). The circuit allows to
algebraically subtract the two signal, which have opposite polarities, and to
obtain an output signal with subtraction of the coherent noise, with an
improvement by about a factor 4 of the signal to noise ratio. The output
signals from circuit A3 are sent to a CAEN V965 ADC [6] for charge analysis.

Figure 1: The electric scheme of the read-out circuit providing the algebraic
sum of the two pad signal (PAD + and PAD -).

A typical ADC distribution of a GGM gap is shown in fig.2 for two different
effective operating voltage, defined as the high voltage set on the HV power
supply corrected for the local atmospheric pressure and temperature. Fig.2 a)
corresponding to HVeff=9.9kV shows a clean avalanche peak well separated from
the pedestal. Fig.2 b) shows the charge distribution at HVeff=10.7kV with two
signal regions corresponding to the avalanche and to avalanche+streamer mode.

Figure 2: Typical ADC charge distributions of one GGM chamber at two operating
voltages. Distribution (a) correspond to HVeff = 9.9kV while distribution (b)
to HVeff=10.7kV. In (b) is clearly visible the streamer peak around 1900 ADC
channels. The events on the left of the vertical line (1450 ADC channels in
this case) are assumed to be pure avalanche events.

Figure 3: Efficiency plot (full dots) of GGM chambers as a function of HVeff.
The efficiency is defined as the ratio between the number of ADC entries above
3$\sigma_{ped}$ and the number of acquired triggers. Open dot plots correspond
to the streamer fraction of the chamber signal as a function of HVeff.

Fig.3 shows the GGMS single gap efficiency (full dots), and the ratio between
the avalanche and the streamer component (open circles), as a function of the
effective high voltage. Each point corresponds to a total of 10000 entries in
the full ADC spectrum. The efficiency is defined as the ratio between the
number of triggers divided by the number of events above 3$\sigma_{ped}$ over
ADC pedestal, where $\sigma_{ped}$ is the pedestal width. The avalanche to
streamer ratio is defined by counting the number of entries in the avalanche
(below the ADC threshold (fig.2 b) and above the pedestal region) and dividing
it by the number of streamer events above the avalanche threshold. Both
efficiency and avalanche plateau are in good agreement with previous results
[5].

In order to determine the sensitivity of GGM gaps to working point shifts, the
avalanche to streamer transition was studied by two methods, the charge method
and the efficiency method. In the charge method, the mean value of the ADC
charge distribution in the whole ADC range is studied as a function of HVeff
(fig.4). Each point corresponds to 10000 events in the whole ADC spectrum. In
the plot three working point regions are identified

  1. 1.

inefficiency (HV${}_{eff}<$ 9.7 kV);

  2. 2.

avalanche (9.7 kV $<$ HV${}_{eff}<$ 10.6 kV;

  3. 3.

avalanche+streamer mode (HV${}_{eff}>$ 10.6 kV).

The best sensitivity to working point shifts is achieved in the
avalanche+streamer region, estimated to be about 25 ADC ch/10 V or 1.2pC/10V.

Figure 4: Avarege avalanche charge of the eight monitor chamber signal as a
function of HVeff. The slope is about 25 ADC ch/10 or 1.2pC/10V. Each point
corresponds to 10000 triggers.

In the efficiency method, the ADC avalanche event yield is studied as a
function of HVeff (5). The avalanche signal increases by increasing the HV
applied to the gap, until it reaches a maximum value after which the streamer
component starts to increase. The 9.0kV-10.0kV shows a sensitivity to work
point changes of approximately 1.3/%/10V.

Figure 5: Streamer and avalanche yields as a function of HVeff. Each point
corresponds to 10000 collected triggers. The solid line has a slope of
approximately 130 events/10 V corresponding to a sensitivity of 1.3%/10V.

##  3 Response of GGM to environmental effects

The workpoint of GGM is affected by environment. However, all environmental
effects cancel out thanks to redundancy of the system. Each environmental
effect not connected with a modification of gas mixture will be cancelled out
by a comparison between different RPC chamber flown with the same gas, which
are affected by the same environmental parameters.

An example of such cancellation is shown in Fig. 6, where the average charge
distribution (black dots) is plotted across a changeover of gas bottles. Data
show a sudden increase in the average charge distribution which may
interpreted as a shift of working point due to changes in gas mixture
composition. By weighing the average charge with a correction factor linearly
depending on atmosferic pressure, however, no significant increase is left in
the distribution of corrected average charge (green dots) which may signal an
anomalous shift due to gas mixture. The cancellation algorithm is applied by
correcting variables withing gaps belonging to the same subdetector. Fig. 7
shows the average charge for two chambers working in different regimes at
different voltages. The average charge of both chambers is completely
correlated, and very well correlated to the atmosferic pressure variations.

Figure 6: Average charge and pressure-corrected charge

Figure 7: Average charge of two chambers at different voltages as influenced
by pressure

##  4 Conclusions

Results from the Gas Gain Monitoring System for the CMS RPC Detector have been
reported on. The purpose of GGM is to monitor any shift of the working point
of the CMS RPC detector. The GGM is being commissioned at CERN and is planned
to start operation by the end of 2008. Preliminary results show good
sensitivity to working point changes. The system redundancy allows for
effectively cancelling out the environmental effects. Further tests are in
progress to determine the sensitivity to gas variations.

## References

  * [1] M. Abbrescia et al., “Gas analysis and monitoring systems for the RPC detector of CMS at LHC”, presented by S.Bianco at the IEEE 2006, San Diego (USA), arXiv:physics/0701014.
  * [2] R. Adolphi et al. [CMS Collaboration], JINST 3 (2008) S08004
  * [3] R. Guida et al., “ The CMS RPC gas system and the Closed Loop recirculation system”, presented at the RPC Conference, Mumbai, 2008.
  * [4] G. Saviano et al., “Material, Filter and Gas analysis for the CMS RPC detector in the Closed Loop test setup”, presented at the RPC Conference, Mumbai, 2008 to appear on Journal of Instrumentation
  * [5] M. Abbrescia et al., “Cosmic ray tests of double-gap resistive plate chambers for the CMS experiment”, Nucl. Instrum. Meth. A 550, 116 (2005).
  * [6] C.A.E.N. Costruzioni Apparecchiature Elettroniche Nucleari S.p.A. Via Vetraia 11 - 55049 Viareggio (Italy).
  * [7] L. Benussi et al., “The CMS RPC Gas Gain Monitoring System: an Overview and Preliminary Result”, presented in Mumbai08 to appear on Nuclear Instruments and Methods arXiv:0812.1108 [physics.ins-det].