# Diagnostic tools for a multivariate negative binomial model for fitting
correlated data with overdispersion

Lizandra Castilho Fabio
Department of Statistics
Federal University of Bahia
Brazil  Cristian Villegas
Department of Exact Sciences
University of São Paulo
Brazil  Jalmar M. F. Carrasco111Corresponding author. Email:
<EMAIL_ADDRESS>
Department of Statistics
Federal University of Bahia
Brazil  Mário de Castro
Instituto de Ciências Matemáticas e de Computação
Universidade de São Paulo
Brazil

###### Abstract

We focus on the development of diagnostic tools and an R package called MNB
for a multivariate negative binomial (MNB) regression model for detecting
atypical and influential subjects. The MNB model is deduced from a Poisson
mixed model in which the random intercept follows the generalized log-gamma
(GLG) distribution. The MNB model for correlated count data leads to an MNB
regression model that inherits the features of a hierarchical model to
accommodate the intraclass correlation and the occurrence of overdispersion
simultaneously. The asymptotic consistency of the dispersion parameter
estimator depends on the asymmetry of the GLG distribution. Inferential
procedures for the MNB regression model are simple, although it can provide
inconsistent estimates of the asymptotic variance when the correlation
structure is misspecified. We propose the randomized quantile residual for
checking the adequacy of the multivariate model, and derive global and local
influence measures from the multivariate model to assess influential subjects.
Finally, two applications are presented in the data analysis section. The code
for installing the MNB package and the code used in the two examples is
exhibited in the Appendix.

keywords: Count data; Overdispersion; Multivariate negative binomial
distribution; MNB package.

##  1 Introduction

Hierarchical models have been suggested for analyzing correlated count data
due to their feature of allowing one to model the intraclass correlation and
accommodate the occurrence of overdispersion simultaneously. This is handled
through the inclusion of random effects in the systematic component of
generalized linear models. Using this approach, it is possible to relax the
assumption about the distribution of the random effects, taking into account
the empirical distributions of the data or individual profiles (Lee et al.,
2006, Molenberghs et al., 2007 and Fabio et al., 2012). Several methods have
been proposed for inferential procedures in hierarchical models due to the
intractable integrals involved in inference functions. The latter fact has
been an obstacle in the development of diagnostic tools, leading to completely
numerical procedures being used. In such cases, a multivariate model deduced
from the random effects approach (Molenberghs and Verbeke, 2010) can be an
alternative to a hierarchical model for fitting correlated count data with
extra variability. Fabio et al. (2012) proposed the random intercept Poisson
mixed regression model by assuming that the random effects follow a
generalized log-gamma (GLG) distribution (Lawless, 1987). This distribution
can be skewed to the right or skewed to the left, with the normal distribution
as a particular case. Thus, the random intercept Poisson-GLG model reduces to
a multivariate negative binomial (MNB) model when it is assumed that the scale
and shape parameters of the GLG distribution are equal. The random intercept
of a Poisson-GLG model is able to accommodate intraclass correlation and
handle overdispersion, due to the several degrees of asymmetry that the GLG
distribution can assume. The MNB inherits these features, with the
overdispersion parameter $\phi^{-1}=\lambda^{2}$ and simpler correlation
structures. In the literature, the MNB regression (MNBR) model has been used
for modeling correlated count data from several areas, for instance, health
(Solis-Trapala and Farewell, 2005), spatial (Moller and Rubak, 2010) and
economics (Sung and Lee, 2018) data.

We focus on developing diagnostic tools and an R (R Core Team, 2017) package
for the MNBR model under overdispersion, which are essential for detecting
outliers and for checking the adequacy of the model (Cook and Weisberg, 1983).
Simulation studies are also conducted to evaluate the asymptotic properties of
the maximum likelihood (ML) estimator and to analyze how these properties’ are
impacted on the random effect distribution misspecification. As noted in
Solis-Trapala and Farewell (2005), the MNBR model can provide inconsistent
estimates of the asymptotic variance of the regression coefficient when the
covariance matrix of the random effect is misspecified. We use randomized
quantile residuals (Dunn and Smyth, 1996) for checking the adequacy of the
MNBR model. Following Cook (1977, 1986) approach, global and local influence
measures are developed from the MNBR model to detect influential subjects.
These methodologies are useful to assess the impact on the estimation
procedure by removing the subject from the data set and by proposing
perturbation models, respectively. The total local influence measure suggested
by Lesaffre and Verbeke (1998) is also extended for the MNBR model. These
measures are helpful for identifying outlying subjects that matter in the MNBR
model and interpret them according to the application. The inferences and
diagnostic analysis from MNBR are performed by using the MNB package developed
by the authors.

The paper is organized as follows: In Section 2, we provide some background
about the MNBR model approach by Fabio et al. (2012). In Section 3, we discuss
diagnostic analysis methodologies for the MNBR model. In Section 4, we perform
a simulation study to evaluate the asymptotic behavior of the ML estimator. In
Section 5, the diagnostic analysis is applied to two real data sets, using the
MNB package. The code for installing the MNB package is presented in the
Appendix. Finally, in Section 6, we provide some discussions.

##  2 MNBR model

Let $y_{ij}$ denote the $j$th measurement taken on the $i$th subject or
cluster, for $j=1,\ldots,m_{i}$ and $i=1,\ldots,n$. Further, let $b_{i}$ be
random effects that follow the GLG distribution (Lawless, 2002). Assuming that
$y_{ij}|b_{i}$ are independent outcomes with probability mass function
represented by a Poisson distribution, Fabio et al. (2012) proposed a random
intercept Poisson-GLG model with the following hierarchical structure: $(i)$
$y_{ij}|b_{i}$ $\stackrel{{\scriptstyle\rm ind}}{{\sim}}$ ${\rm
Poisson}(u_{ij})$, $(ii)$ $u_{ij}=\mu_{ij}\exp(b_{i}),$ and $(iii)$
$b_{i}\stackrel{{\scriptstyle\rm iid}}{{\sim}}{\rm GLG}(0,\sigma,\lambda),$
where $\mu_{ij}=\exp{(\bm{x}^{\top}_{ij}\bm{\beta}})$, with
$\bm{x}_{ij}=(\bm{x}_{ij1},\ldots,\bm{x}_{ijp})^{\top}$ containing the values
of explanatory variables, $\bm{\beta}=(\beta_{1},\ldots,\beta_{p})^{\top}$ the
vector of regression coefficients, and $\sigma>0$ and $\lambda\in\mathbb{R}$
the scale and shape parameters of the GLG distribution. In general, in the
random effects approach, the marginal distribution of the $i$th subject does
not have an explicit form (Molenberghs and Verbeke, 2010 pp. 259-266). The ML
estimates are obtained by integrating out the random effect and maximizing the
log-likelihood function. Fabio et al. (2012) showed that the integral can be
solved analytically when the scale and shape parameters of the GLG
distribution are equal ($\sigma=\lambda$, with $\lambda\in\mathbb{R}^{+}$).
When $\phi=\lambda^{-2}$, the MNB distribution is deduced from the random
intercept Poisson-GLG model of the following form:

$f(\bm{y}_{i};\bm{\beta},\phi)=\frac{\Gamma(\phi+y_{i+})\phi^{\phi}}{\left(\prod_{j=1}^{m_{i}}y_{ij}!\right)\Gamma(\phi)}\frac{\exp\left(\sum_{j=1}^{m_{i}}y_{ij}\log(\mu_{ij})\right)}{(\phi+\mu_{i+})^{\phi+y_{i+}}},$
(1)

where $\bm{y}_{i}=(y_{i1},\ldots,y_{im_{i}})^{\top}$ is the vector of
measurements available for the $i$th subject, $\phi$ is the dispersion
parameter, $\Gamma(\cdot)$ is the gamma function,
$y_{i+}=\sum_{j=1}^{m_{i}}y_{ij}$, and $\mu_{i+}=\sum_{j=1}^{m_{i}}\mu_{ij}$.
According to Johnson et al. (1997), the MNB distribution belongs to the
discrete multivariate exponential family of distributions, and its marginals
$f(y_{ij})$ are negative binomial distributions with
$\textrm{E}(y_{ij})=\mu_{ij}$ and
$\textrm{Var}(y_{ij})=\mu_{ij}+\mu_{ij}^{2}/\phi$, for $i=1,\ldots,n$ and
$j=1,\ldots,m_{i}.$ The covariances
$\textrm{Cov}(y_{ij},y_{ij^{{}^{\prime}}})=\mu_{ij}\mu_{ij^{{}^{\prime}}}/\phi$
and intraclass correlations
$\textrm{Corr}(y_{ij},y_{ij^{{}^{\prime}}})=\sqrt{\mu_{ij}\mu_{ij^{{}^{\prime}}}}/\left(\sqrt{\phi+\mu_{ij}}\sqrt{\phi+\mu_{ij^{{}^{\prime}}}}\right)$
for $j\neq j^{{}^{\prime}}$ are always positive. When $\phi^{-1}$ indicates
the amount of excess correlation in the data (Hilbe, 2011), the parameter
$\phi$ is also called an overdispersion parameter. For large values of $\phi$,
the marginals of the MNB distribution behave approximately as independent
Poisson distributions with mean $\mu_{ij}$. By
$\bm{y}_{i}\sim\textrm{MNB}(\bm{\mu}_{i},\phi),$ we denote independent vectors
of random outcomes that follow the probability function given in (1), with
$\bm{\mu}_{i}=(\mu_{i1},\ldots,\mu_{im_{i}})^{\top}$ and $\phi>0$.

The MNBR model is defined by assuming $(i)$
$\bm{y}_{i}\stackrel{{\scriptstyle\rm ind}}{{\sim}}{\rm
MNB}(\bm{\mu_{i}},\phi)$ and $(ii)$
$\log(\mu_{ij})=\bm{x}_{ij}^{\top}\bm{\beta}$. Letting
$\bm{y}=(\bm{y}^{\top}_{1},\ldots,\bm{y}^{\top}_{n})^{\top}$ be the vector
containing all the measured outcomes for the $i$th subject, the log-likelihood
function is given by

$\displaystyle\ell(\bm{\theta})$
$\displaystyle=\sum_{i=1}^{n}\log\left\\{\frac{\Gamma(\phi+y_{i+})}{\Gamma(\phi)}\right\\}-\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}\log\left(y_{ij}!\right)+n\phi\log(\phi)$
$\displaystyle\quad-\phi\sum_{i=1}^{n}\log(\phi+\mu_{i+})+\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}y_{ij}\log\left\\{\frac{\mu_{ij}}{\phi+\mu_{i+}}\right\\},$
(2)

where $\bm{\theta}=(\bm{\beta}^{\top},\phi)^{\top}$. The ML estimates
$\widehat{\bm{\theta}}$ of $\bm{\theta}$ are computed by using the quasi-
Newton (BFGS) method. Alternatively, we can solve the nonlinear equation
obtained by setting the components of the score vector equal to zero, that is,
$\bm{U}_{\theta}=(\bm{U}_{\beta}^{\top},U_{\phi})^{\top}=\bm{0}$ (details
about the calculations are in the Appendix). For interval estimation and
hypothesis tests on the model parameters, the expected or observed Fisher
information is required. Under standard regularity conditions,
$\widehat{\bm{\theta}}-\bm{\theta}$ is asymptotically distributed as a
multivariate normal distribution with mean $\bm{0}$ and covariance matrix
equal to the inverse of the Fisher information matrix of the MNBR model (see
Fabio et al., 2012).

##  3 Diagnostic analysis

This section is devoted to diagnostic tools for the MNBR model. In what
follows, we will discuss the residual analysis, and the global and local
influence methodologies.

###  3.1 Residual analysis

Let $\bm{y}_{i}=(y_{i1},y_{i2},\ldots,y_{im_{i}})^{\top}$ be a random vector
that follows a MNB distribution. According to Tsui (1986), the distribution of
$y_{i+}=\sum_{j=1}^{m_{i}}y_{ij}$ is negative binomial with probability
distribution function given by

$\displaystyle
f(y_{i+})=\frac{\Gamma(y_{i+}+\phi)}{(y_{i+})!\Gamma(\phi)}q^{\phi}(1-q)^{y_{i+}},\
\ \ y_{i+}=0,1,2,...,$

where $q=\phi/(\phi+\mu_{i+})$. Then, the randomized quantile residuals (Dunn
and Smyth, 1996), which follow a standard normal distribution, can be used to
assess departures from the MNBR model. If $F(y_{i+};\mu,\phi)$ is the
cumulative distribution of $f(y_{i+})$, $a_{i}=\lim_{y\uparrow
y_{i+}}F(y;\widehat{\mu_{i}},\widehat{\phi})$, and
$b_{i}=F(y_{i+};\widehat{\mu}_{i},\widehat{\phi})$, then the randomized
quantile residuals for $\bm{y}_{i}$ are given by $r_{q,i}=\Phi^{-1}(u_{i})$,
where $\Phi(\cdot)$ is the cumulative distribution function of the standard
normal and $u_{i}$ is a uniform random variable on the interval
$(a_{i},b_{i}]$.

###  3.2 Global influence

Based on the case-deletion approach (Cook, 1977), we propose the generalized
Cook’s distance as a global influence measure to assess the impact on the ML
estimates of the MNBR model when the $i$th subject (or cluster) is removed
from the data set. Let $\bm{y}_{(i)}$ be a vector of random outcomes after
deleting the $i$th subject. Let $\widehat{\bm{\theta}}_{(i)}$ be the ML
estimate of $\bm{\theta}$ computed from $\ell_{(i)}(\bm{\theta})$, where
$\ell_{(i)}(\bm{\theta})$ is obtained from (2) after removing the $i$th
subject. The generalized Cook’s distance is defined as the standardized norm
of the distance between $\widehat{\bm{\theta}}_{(i)}$ and
$\widehat{\bm{\theta}}$, given by the expression

$GD_{i}(\bm{\theta})=(\widehat{\bm{\theta}}_{(i)}-\widehat{\bm{\theta}})^{\top}[\ddot{\ell}_{\bm{\theta}\bm{\theta}}]^{-1}(\widehat{\bm{\theta}}_{(i)}-\widehat{\bm{\theta}}),$

where $\ddot{\ell}_{\bm{\theta}\bm{\theta}}$ (see the Appendix) is the Fisher
information matrix of the MNBR model. Large values of $GD_{i}$ indicate that
the ML estimates are strongly influenced by deleting the $i$th subject.
Another popular measure of the difference between
$\widehat{\bm{\theta}}_{(i)}$ and $\widehat{\bm{\theta}}$ is the likelihood
displacement given by
$LD_{(i)}(\bm{\theta})=2\\{\ell(\widehat{\bm{\theta}})-\ell(\widehat{\bm{\theta}}_{(i)})\\}$.

###  3.3 Local influence

The local influence methodology was proposed by Cook (1986) for assessing the
sensitivity of the parameters estimated when small perturbations are
introduced in the model. Consider the perturbation vector
$\bm{\omega}=(\omega_{1},\ldots,\omega_{v})^{\top}$, varying in some open
subset $\Omega\subset\mathbb{R}^{v}$. Let $\ell(\bm{\theta}|\bm{\omega})$
denote the log-likelihood function of the perturbed model. It is assumed there
exists $\bm{\omega}_{0}\in\Omega$ such that
$\ell(\bm{\theta}|\bm{\omega}_{0})=\ell(\bm{\theta})$ for all $\bm{\theta}.$
The influence of the minor perturbation on the ML estimate
$\widehat{\bm{\theta}}$ may be assessed by the likelihood displacement
$LD(\bm{\omega})=2\\{\ell(\widehat{\bm{\theta}})-\ell(\widehat{\bm{\theta}}_{\bm{\omega}})\\}$,
where $\widehat{\bm{\theta}}_{\bm{\omega}}$ denotes the ML estimate from the
perturbed model. A plot of $LD(\bm{\omega})$ versus $\bm{\omega}$ contains
essential information about the influence of a perturbation scheme. Cook’s
idea consists of selecting a unit direction $\bm{d}$ and evaluating the plot
$LD(\omega_{0}+a\bm{d})$, where $a\in\mathbb{R}$. This plot is called lifted
line, and each fitted line can be obtained by considering the normal curvature
defined by
$C_{\bm{d}}(\bm{\theta})=2|\bm{d}^{\top}\bm{\Delta}^{\top}\ddot{\ell}(\bm{\theta})^{-1}\bm{\Delta}\bm{d}|$
around $a=0$, where
$\ddot{\ell}(\bm{\theta})=\partial^{2}\ell(\bm{\theta})/\partial\bm{\theta}\partial{\bm{\theta}}^{\top}$
is evaluated at $\widehat{\bm{\theta}}$ and
$\bm{\Delta}=\partial^{2}\ell(\bm{\theta}|\bm{\omega})/\partial\bm{\theta}\partial\bm{\omega}^{\top}$
is evaluated at $\widehat{\bm{\theta}}$ and $\bm{\omega}_{0}.$ Large values of
$C_{\bm{d}}(\bm{\theta})$ indicate the sensitivity induced by perturbation
schemes in direction $\bm{d}.$ Cook (1986) suggests the local influence
measure $C_{\bm{d}_{max}}(\bm{\theta})$, evaluated in the direction
$\bm{d}_{max}$ corresponding to the eigenvector with maximal normal curvature.
If the $i$th component of $\bm{d}_{max}$ is relatively large, this indicates
that perturbations may lead to substantial changes in the ML estimates. Based
on this approach, Lesaffre and Verbeke (1998) suggested the total local
curvature corresponding to the $i$th element. This local influence measure is
obtained by taking the direction $\bm{d}_{i}$, an $n\times 1$ vector of zeros
with a 1 in the $i$th position. The total local curvature in direction
$\bm{d}_{i}$ assumes the form
$C_{i}(\bm{\theta})=2|\bm{\Delta}_{i}^{\top}\ddot{\ell}(\bm{\theta})^{-1}\bm{\Delta}_{i}|$,
where
$\bm{\Delta}_{i}=\partial^{2}\ell_{i}(\bm{\theta}|\bm{\omega})/\partial\bm{\theta}\partial\bm{\omega}^{\top}$
denotes the $i$th row of $\bm{\Delta}$, such that
$\ell_{i}(\bm{\theta}|\bm{\omega})$ is the contribution of the $i$th
individual to the log-likelihood function. $\bm{\Delta}_{i}$ is evaluated at
$\widehat{\bm{\theta}}$ and $\bm{\omega}_{0}$, and it is recommended to look
at the index plot of $C_{i}(\bm{\theta})$ to assess influential subjects.
Usually, in the literature, three types of perturbation schemes are considered
for count data: case weights, explanatory variable, and dispersion parameter
perturbation schemes.

### Case weights perturbation for the $i$th subject

The log-likelihood function of the perturbed model takes the form

$\displaystyle\ell(\bm{\theta}|\bm{\omega})$
$\displaystyle=\sum_{i=1}^{n}\omega_{i}\Big{\\{}\log\Big{(}\Gamma(\phi+y_{i+})\Big{)}-\log\Big{(}\Gamma(\phi)\Big{)}-\sum_{j=1}^{m_{i}}\log(y_{ij}!)+\phi\log(\phi)$
$\displaystyle\quad-\phi\log(\phi+\mu_{i+})+\sum_{j=1}^{m_{i}}y_{ij}\log(\mu_{ij})-\sum_{j=1}^{m_{i}}y_{ij}\log(\phi+\mu_{i+})\Big{\\}}.$

The $\bm{\Delta}$ matrix in $C_{\bm{d}}(\bm{\theta})$ is given by
$\bm{\Delta}^{\top}=(\bm{\Delta}^{\top}_{\bm{\beta}},\bm{\Delta}^{\top}_{\phi})^{\top}$,
with
$\bm{\Delta}^{\top}_{\bm{\beta}}=(\bm{\Delta}^{\top}_{\beta_{1}},\bm{\Delta}^{\top}_{\beta_{2}},\ldots,\bm{\Delta}^{\top}_{\beta_{p}})$
and
$\bm{\Delta}_{\beta_{k}}=(\Delta_{k1},\Delta_{k2},\ldots,\Delta_{ki},\ldots,\Delta_{kn}),$
where

$\displaystyle\Delta_{ki}$
$\displaystyle=\sum_{j=1}^{m_{i}}\left\\{y_{ij}-\frac{(\widehat{\phi}+y_{i+})}{(\widehat{\phi}+\widehat{\mu}_{i+})}\exp(\bm{x}^{\top}_{ij}\widehat{\bm{\beta}})\right\\}\bm{x}_{ijk}=\sum_{j=1}^{m_{i}}(y_{ij}-\widehat{a}_{i}\widehat{\mu}_{ij})\bm{x}_{ijk},$

and $\bm{\Delta}^{\top}_{\phi}=(\Delta_{\phi 1},\Delta_{\phi
2},\ldots,\Delta_{\phi i},\ldots,\Delta_{\phi n})$, with

$\displaystyle\Delta_{\phi
i}=\psi(\widehat{\phi}+y_{i+})-\psi(\widehat{\phi})+\log\left(\frac{\widehat{\phi}}{\widehat{\phi}+\widehat{\mu}_{i+}}\right)+\left(\frac{\widehat{\mu}_{i+}-y_{i+}}{\widehat{\phi}+\widehat{\mu}_{i+}}\right),$

where $\psi(\cdot)$ is the digamma function,
$\widehat{\mu}_{i+}=\sum_{j=1}^{m_{i}}\widehat{\mu}_{ij}$, and
$\widehat{\mu}_{ij}=\exp(\bm{x}_{ij}^{\top}\widehat{\bm{\beta}})$, for
$j=1,\ldots,m_{i}$ and $i=1,\ldots,n$.

### Case weights perturbation for the $j$th measurement of the $i$th subject

For this perturbation scheme, the log-likelihood function can be expressed as

$\displaystyle\ell(\bm{\theta}|\bm{\omega})$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}\omega_{ij}\Big{\\{}\frac{1}{m_{i}}\log\Gamma(\phi+y_{i+})-\frac{1}{m_{i}}\log\Gamma(\phi)-\log(y_{ij}!)$
$\displaystyle\quad{+}\quad\frac{1}{m_{i}}\phi\log(\phi)-\frac{1}{m_{i}}\phi\log(\phi+\mu_{i+})+y_{ij}\log(\mu_{ij})-y_{ij}\log(\phi+\mu_{i+})\Big{\\}}.$

The matrix
$\bm{\Delta}=(\bm{\Delta}^{\top}_{\bm{\beta}},\bm{\Delta}_{\phi}^{\top})^{\top},$
in which element $ij$ of $\bm{\Delta}^{\top}_{\bm{\beta}}$ and
$\bm{\Delta}^{\top}_{\phi}$, respectively, is given by

$\displaystyle\Delta_{k{ij}}$
$\displaystyle=y_{ij}x_{ijk}-\frac{\widehat{\phi}/m_{i}+y_{ij}}{\widehat{\phi}+\widehat{\mu}_{i+}}\sum_{j=1}^{m_{i}}\widehat{\mu}_{ij}x_{ijk},k=1,\ldots,p\quad{\rm
and}$ $\displaystyle\Delta_{\phi{ij}}$
$\displaystyle=\frac{1}{m_{i}}(\psi(\widehat{\phi}+y_{i+})-\psi(\widehat{\phi}))+\frac{1}{m_{i}}(1+\log(\widehat{\phi})-\log(\widehat{\phi}+\widehat{\mu}_{i+}))$
$\displaystyle\quad{-}\quad\frac{\widehat{\phi}/m_{i}+y_{ij}}{\widehat{\phi}+\widehat{\mu}_{i+}}.$

### Explanatory variable perturbation

We now consider an additive perturbation on a particular continuous
explanatory variable,
$\bm{x}_{ijk},j=1,\ldots,m_{i},i=1,\ldots,n,k=1,\ldots,p$, by setting
$\bm{x}^{*}_{ijk}=\bm{x}_{ijk}+\omega_{i}S_{x}$, where $S_{x}$ is a scale
factor and $\omega_{i}\in\mathbb{R}$. This perturbation scheme leads to the
following expression for the log-likelihood function:

$\displaystyle\ell(\bm{\theta}|\bm{\omega})$
$\displaystyle=\sum_{i=1}^{n}\Bigg{\\{}\log\Big{(}\Gamma(\phi+y_{i+})\Big{)}-\log\Big{(}\Gamma(\phi)\Big{)}-\sum_{j=1}^{m_{i}}\log(y_{ij}!)+\phi\log(\phi)$
$\displaystyle\quad-\phi\log(\phi+\mu_{i+}^{*})+\sum_{j=1}^{m_{i}}y_{ij}\log(\mu_{ij}^{*})-\sum_{j=1}^{m_{i}}y_{ij}\log(\phi+\mu_{i+}^{*})\Bigg{\\}},$

where $\mu_{i+}^{*}=\sum_{j=1}^{m_{i}}\mu_{ij}^{*}$,
$\mu_{ij}^{*}=\exp(\bm{x}_{ij}^{*\top}\bm{\beta})$,
$\bm{x}_{ij}^{*\top}\bm{\beta}=\beta_{1}+\beta_{2}x_{ij2}+\ldots+\beta_{k}(x_{ijk}+\omega_{i}S_{x})+\ldots+\beta_{p}x_{ijp}$,
and $\bm{\omega}_{0}=(0,\ldots,0)^{\top}$. For $k=1,\ldots,p$, if $t=k$,
$\bm{\Delta}_{t}$ is given by

$\displaystyle\Delta_{ti}$
$\displaystyle=\frac{\widehat{\beta}_{t}S_{x}\widehat{\phi}\exp(\bm{x}^{\top}_{ij}\widehat{\bm{\beta}})}{(\widehat{\phi}+\widehat{\mu}_{i+})}\Bigg{(}1-\frac{\widehat{\mu}_{i+}}{\widehat{\phi}+\widehat{\mu}_{i+}}\Bigg{)}x_{ijt}+\frac{S_{x}\exp(\bm{x}_{ij}^{\top}\widehat{\bm{\beta}})}{\widehat{\phi}+\widehat{\mu}_{i+}}$
$\displaystyle\quad+\sum_{j=1}^{m_{i}}y_{ij}\Big{[}S_{x}+\Big{(}\widehat{\beta}_{t}S_{x}-\frac{1}{\widehat{\mu}_{ij}}\Big{)}x_{ijt}\Big{]}+\sum_{j=1}^{m_{i}}\frac{y_{ij}\exp(\bm{x}_{ij}^{\top}\widehat{\bm{\beta}})}{\widehat{\phi}+\widehat{\mu}_{i+}}$
$\displaystyle\quad\times\Big{[}S_{x}+(\widehat{\beta}_{t}S_{x}-1)x_{ijt}\Big{]}.$

If $t\neq k$,

$\displaystyle\Delta_{ti}$
$\displaystyle=\frac{\widehat{\beta}_{t}S_{x}\widehat{\phi}\exp(\bm{x}^{\top}_{ij}\widehat{\bm{\beta}})}{(\widehat{\phi}+\widehat{\mu}_{i+})}\Bigg{(}1-\frac{\widehat{\mu}_{i+}}{\widehat{\phi}+\widehat{\mu}_{i+}}\Bigg{)}x_{ijt}+\sum_{j=1}^{m_{i}}y_{ij}x_{ijt}$
$\displaystyle\quad\times\Big{(}\widehat{\beta}_{t}S_{x}-\frac{1}{\widehat{\mu}_{ij}}\Big{)}+\sum_{j=1}^{m_{i}}\frac{(\widehat{\beta}_{t}S_{x}-1)}{(\widehat{\phi}+\widehat{\mu}_{i+})}\exp(\bm{x}_{ij}^{\top}\widehat{\bm{\beta}})y_{ij}x_{ijt},$

and the $i$th element of $\bm{\Delta}^{\top}_{\phi}$ is

$\Delta_{\phi_{i}}=\frac{\widehat{\beta}_{t}S_{x}\widehat{\mu}_{i+}}{(\widehat{\phi}+\widehat{\mu}_{i+})}\Big{(}\frac{\widehat{\phi}}{\widehat{\phi}+\widehat{\mu}_{i+}}-1\Big{)}+\widehat{\beta}_{t}S_{x}\sum_{j=1}^{m_{i}}\frac{y_{ij}\widehat{\mu}_{i+}}{(\widehat{\phi}+\widehat{\mu}_{i+})^{2}}.$

### Dispersion parameter perturbation

Let $\phi^{*}_{i}=\omega_{i}\times\phi$ be the dispersion parameter
perturbation, where $\omega_{i}\in\mathbb{R}^{+}$. The log-likelihood function
of the perturbed model takes the form

$\displaystyle\ell(\bm{\theta}|\bm{\omega})$
$\displaystyle=\sum_{i=1}^{n}\Bigg{\\{}\log\Big{(}\Gamma(\phi^{*}_{i}+y_{i+})\Big{)}-\log\Big{(}\Gamma(\phi^{*}_{i})\Big{)}-\sum_{j=1}^{m_{i}}\log(y_{ij}!)+\phi^{*}_{i}\log(\phi^{*}_{i})$
$\displaystyle\quad{-}\phi^{*}_{i}\log(\phi^{*}_{i}+\mu_{i+})+\sum_{j=1}^{m_{i}}y_{ij}\log(\mu_{ij})-\sum_{j=1}^{m_{i}}y_{ij}\log(\phi^{*}_{i}+\mu_{i+})\Bigg{\\}}.$

$\bm{\Delta}=(\bm{\Delta}_{\bm{\beta}}^{\top},\bm{\Delta}^{\top}_{\phi})^{\top}$,
where the $i$-th elements of $\bm{\Delta}^{\top}_{\bm{\beta}}$ and
$\bm{\Delta}^{\top}_{\phi}$, respectively, are given by

$\displaystyle\Delta_{k_{i}}$
$\displaystyle=\frac{\widehat{\phi}^{*}_{i}\widehat{\mu}_{i+}(y_{i+}-\widehat{\mu}_{i+})}{(\widehat{\phi}^{*}_{i}+\widehat{\mu}_{i+})^{2}}x_{ijk}~{}~{}{\rm
and}$ $\displaystyle\Delta_{\phi_{i}}$
$\displaystyle=\psi(\widehat{\phi}^{*}_{i}+y_{i+})-\psi(\widehat{\phi}^{*}_{i})+\widehat{\phi}^{*}_{i}[\psi^{\prime}(\widehat{\phi}^{*}_{i}+y_{i+})-\psi^{\prime}(\widehat{\phi}^{*}_{i})]+1$
$\displaystyle\quad{+}\log\left(\frac{\widehat{\phi}^{*}_{i}}{\widehat{\phi}^{*}_{i}+\widehat{\mu}_{i+}}\right)+\frac{\widehat{\mu}_{i+}}{\widehat{\phi}^{*}_{i}+\widehat{\mu}_{i+}}+\frac{y_{i+}\widehat{\mu}_{i+}}{(\widehat{\phi}^{*}_{i}+\widehat{\mu}_{i+})^{2}}+\frac{\widehat{\phi}^{*}_{i}(2\widehat{\mu}_{i+}+\widehat{\phi}^{*}_{i})}{(\widehat{\phi}^{*}_{i}+\widehat{\mu}_{i+})^{2}}.$

For all the perturbation schemes, the $\ddot{\ell}(\bm{\theta})$ matrix in
$C_{\bm{d}}(\bm{\theta})$ is described in the Appendix.

##  4 Numerical results

First, a simulation study is conducted to evaluate the asymptotic behavior of
the ML estimator, $\widehat{\bm{\theta}}$. The vector
$\bm{y}=(\bm{y}^{\top}_{1},\ldots,\bm{y}^{\top}_{n})^{\top},$ where
$\bm{y}_{i}\sim\text{MNB}(\bm{\mu}_{i},\phi)$, is generated from a random
intercept Poisson-GLG model with $(i)$ $y_{ij}|b_{i}$
$\stackrel{{\scriptstyle\rm ind}}{{\sim}}$ ${\rm Poisson}(u_{ij})$, $(ii)$
$\log(u_{ij})=\beta_{0}+\beta_{1}x_{1ij}+\beta_{2}x_{2ij}+b_{i}$, and $(iii)$
$b_{i}\stackrel{{\scriptstyle\rm iid}}{{\sim}}{\rm GLG}(0,\lambda,\lambda),$
where $x_{1ij}\sim N(0,1)$ and $x_{2ij}$ is a dummy variable with two levels,
for $i=1,\ldots,n$ and $j=1,2,3$. Further, it is assumed that
$\bm{\beta}=(1.5,1.0,0.0)^{\top}$ and $\lambda=\phi^{-1/2}$ in the GLG
distribution. The Bias, root mean square error (RMSE) and coverage
probabilities for the $95\%$ confidence level are obtained from $R=10,000$
Monte Carlo replications performed for a sample size of $n=50$, $100$, $150$,
and $200$ and for three different values of the dispersion parameter, namely,
$\phi=3$, $5$ and $7$. The estimates of the Bias and RMSE for $\bm{\theta}$
are obtained from: ${\rm Bias}(\theta_{s})=\overline{\theta}_{s}-\theta_{s}$
and ${\rm
RMSE}(\theta_{s})=\sqrt{\sum_{r=1}^{R}(\widehat{\theta}^{(r)}_{s}-\theta_{s})^{2}/R}$,
respectively, where
$\overline{\theta}_{s}=\sum_{r=1}^{R}\widehat{\theta}_{s}^{(r)}/R$ and
$\widehat{\theta}^{(r)}_{s}$ the estimate from $s$th parameter obtained in
$r$th Monte Carlo replication. These simulation results are presented in Table
1. As expected, the Bias and RMSE values of the regression coefficients
estimates decrease as the sample size and values of the $\phi$ parameter
increase. The coverage probabilities are close to 0.95. Moreover, the results
show that the dispersion parameter exhibits a desirable asymptotic behavior
when $\phi$ assumes small values. Since $\phi=\lambda^{-2}$, it is possible to
affirm that the asymptotic properties of the dispersion parameter estimator
are associated with the asymmetry of the GLG distribution.

Table 1: Bias, root mean square error (RMSE), and coverage (%) of 95% confidence intervals. $\phi(\lambda)$ | Measure | $n$ | $\phi$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$
---|---|---|---|---|---|---
3.0(0.58) | Bias | 50 | 0.3931 | -0.0078 | -0.0014 | -0.0029
100 | 0.1778 | -0.0041 | 0.0003 | 0.0000
150 | 0.0995 | -0.0034 | 0.0001 | 0.0000
200 | 0.0735 | -0.0031 | -0.0002 | 0.0000
RMSE | 50 | 1.0345 | 0.1298 | 0.0516 | 0.1821
100 | 0.5951 | 0.0920 | 0.0335 | 0.1280
150 | 0.4387 | 0.0748 | 0.0256 | 0.1031
200 | 0.3793 | 0.0642 | 0.0251 | 0.0886
Coverage | 50 | 96.72 | 93.92 | 94.61 | 93.73
100 | 96.00 | 94.33 | 95.17 | 94.23
150 | 95.75 | 94.80 | 95.15 | 94.73
200 | 95.28 | 94.93 | 95.22 | 94.93
5.0(0.44) | Bias | 50 | 0.8473 | -0.0063 | 0.0002 | 0.0000
100 | 0.3264 | -0.0020 | -0.0003 | 0.0000
150 | 0.2179 | -0.0021 | 0.0001 | 0.0000
200 | 0.1583 | -0.0009 | 0.0001 | 0.0000
RMSE | 50 | 2.1448 | 0.1078 | 0.0492 | 0.1505
100 | 1.1021 | 0.0741 | 0.0299 | 0.1038
150 | 0.8531 | 0.0613 | 0.0264 | 0.0829
200 | 0.7190 | 0.0534 | 0.0228 | 0.0724
Coverage | 50 | 97.22 | 93.93 | 94.94 | 93.51
100 | 95.94 | 94.59 | 95.02 | 94.26
150 | 95.32 | 94.59 | 94.87 | 95.21
200 | 95.18 | 94.64 | 95.24 | 94.87
7.0(0.38) | Bias | 50 | 1.3985 | -0.0057 | 0.0000 | 0.0000
100 | 0.5524 | -0.0031 | 0.0000 | 0.0000
150 | 0.3459 | -0.0019 | 0.0000 | 0.0000
200 | 0.2485 | -0.0014 | 0.0000 | 0.0000
RMSE | 50 | 3.6728 | 0.0976 | 0.0470 | 0.1343
100 | 1.7556 | 0.0668 | 0.0304 | 0.0901
150 | 1.3277 | 0.0540 | 0.0256 | 0.0731
200 | 1.1060 | 0.0473 | 0.0219 | 0.0632
Coverage | 50 | 97.13 | 93.97 | 95.01 | 93.36
100 | 96.92 | 94.54 | 94.94 | 94.83
150 | 95.83 | 94.80 | 94.85 | 95.05
200 | 95.81 | 94.76 | 94.76 | 95.07

A second simulation study is performed to evaluate the impact of misspecifying
the random effect distribution on the ML estimates of the MNBR model. It is
assumed that the random intercept is normally distributed and negatively
correlated. For sample size $n=50,100,$ and $150$, the vector
$\bm{y}=(\bm{y}^{\top}_{1},\ldots,\bm{y}^{\top}_{n})^{\top}$ is generated from
a random intercept Poisson-Normal distribution with the following hierarchical
structure: $(i)$ $y_{ij}|b_{i}$ $\stackrel{{\scriptstyle\rm ind}}{{\sim}}$
${\rm Poisson}(u_{ij})$, $(ii)$
$\log(u_{ij})=\beta_{0}+\beta_{1}x_{1ij}+b_{i}$ and $(iii)$
$b_{i}\stackrel{{\scriptstyle\rm iid}}{{\sim}}{\rm N}(0,\sigma^{2})$,
$\sigma^{2}=\phi^{-1}$ and $(i)$ $y_{ij}|b_{ij}$ $\stackrel{{\scriptstyle\rm
ind}}{{\sim}}$ ${\rm Poisson}(u_{ij})$, $(ii)$
$\log(u_{ij})=\beta_{0}+\beta_{1}x_{1ij}+z_{1ij}b_{ij}$, and $(iii)$
$\bm{b}{\sim}{\rm N}_{3}(\bm{0},\bm{\Sigma}),$ where
$(\beta_{0},\beta_{1})^{\top}=(1.5,1.0)$, $x_{1ij}\sim U(0,1)$,

$\displaystyle\bm{\Sigma}=\Bigg{[}\begin{tabular}[]{ccc}$\phi^{-1}$&-0.5&-0.1\\\
-0.5&$\phi^{-1}$&-1.0\\\ -0.1&-1.0&$\phi^{-1}$\\\ \end{tabular}\Bigg{]},$

and $\phi^{-1}=4.0$. In Table 2, we present the Bias and root mean square
error (RMSE), calculated by fitting the MNBR model on 1,000 Monte Carlo
replications for different sample sizes. Based on the assumption that
$b_{i}\stackrel{{\scriptstyle\rm iid}}{{\sim}}{\rm N}(0,\sigma^{2})$, we
observe that the parameter estimates of $\phi$ and $\beta_{0}$ are biased. To
summarize, we conclude that the MNBR model provides inconsistent estimates of
the asymptotic variance of the ML estimators when the covariance matrix of the
random effect is misspecified.

Table 2: Bias and root mean square error (RMSE) when the covariance matrix of the random response variable is misspecified. Variance | $n$ | Measure | $\phi$ | $\beta_{0}$ | $\beta_{1}$
---|---|---|---|---|---
$\sigma^{2}$ | 50 | Bias | -0.4285 | -1.8818 | 0.0561
RMSE | 0.4584 | 1.9805 | 1.0931
100 | Bias | -0.3799 | -1.9408 | 0.0497
RMSE | 0.3968 | 2.0068 | 0.9257
150 | Bias | -0.3633 | -1.9282 | -0.0129
RMSE | 0.3750 | 1.9754 | 0.6798
$\bm{\Sigma}$ | 50 | Bias | -0.4965 | -1.8623 | 0.0068
RMSE | 0.5301 | 1.9669 | 1.0656
100 | Bias | -0.4375 | -1.9014 | -0.0315
RMSE | 0.4560 | 1.9601 | 0.8784
150 | Bias | -0.4196 | -1.9228 | -0.0372
RMSE | 0.4353 | 1.9953 | 0.8282

A third simulation study is conducted to evaluate the variance-to-mean ratio,
${\rm VMR}={\rm Var}(\bm{y})/{\rm E}(\bm{y})$. Considering three different
values for the dispersion parameter, $\phi=0.1,1.0$ and $50$, random samples
are generated and the VMR statistic are computed. The ranges to VMR, for each
value of $\phi$, are $(122.39,220.24)$, $(19.38,29.35)$ and $(5.85,16.87)$,
respectively. In the Table 3, we reported the bias, RMSE and coverage of the
confidence intervals. As expected, Table 3 shows that the ML estimates present
a desirable asymptotic behavior when the sample size increases. In addition,
we observe that when $\phi$ increases the VMR statistic decreases and assumes
values more than one. Thus, it is possible to conclude that the MNBR model is
an overdispersion model. When the dispersion parameter is more than 50 similar
results are obtained.

Table 3: Bias, root mean square error (RMSE), coverage (%) of 95% confidence intervals. $\phi$ | Measure | $n$ | $\phi$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | VMR
---|---|---|---|---|---|---|---
0.1 | Bias | 50 | 0.0093 | -0.1769 | -0.0006 | 0.0575 | (122.39,220.24)
100 | 0.0051 | -0.1201 | -0.0004 | 0.0347
150 | 0.0027 | -0.0838 | -0.0003 | 0.0310
200 | 0.0020 | -0.0331 | -0.0001 | 0.0226
RMSE | 50 | 0.0299 | 0.7048 | 0.0707 | 0.9794
100 | 0.0198 | 0.5105 | 0.0411 | 0.6742
150 | 0.0151 | 0.3891 | 0.0298 | 0.5531
200 | 0.0130 | 0.3203 | 0.0253 | 0.4468
Coverage | 50 | 95.4 | 91.3 | 95.0 | 91.1
100 | 95.0 | 92.8 | 94.2 | 92.8
150 | 95.3 | 93.7 | 94.9 | 92.9
200 | 95.5 | 94.5 | 95.1 | 94.7
1.0 | Bias | 50 | 0.0979 | -0.0201 | 0.0017 | -0.0026 | (19.38,29.35)
100 | 0.0390 | -0.0059 | 0.0000 | -0.0189
150 | 0.0317 | -0.0045 | 0.0000 | -0.0160
200 | 0.0164 | -0.0039 | 0.0000 | -0.0006
RMSE | 50 | 0.2706 | 0.2188 | 0.0552 | 0.3046
100 | 0.1627 | 0.1550 | 0.0380 | 0.2132
150 | 0.1351 | 0.1249 | 0.0291 | 0.1712
200 | 0.1065 | 0.1069 | 0.0244 | 0.1477
Coverage | 50 | 95.7 | 93.9 | 94.7 | 92.2
100 | 95.4 | 94.3 | 95.6 | 93.5
150 | 95.3 | 94.9 | 95.2 | 94.1
200 | 95.1 | 95.2 | 94.9 | 94.7
50 | Bias | 50 | 36.1329 | -0.0034 | 0.0000 | 0.0022 | (5.85,16.87)
100 | 17.2463 | -0.0030 | 0.0000 | 0.0015
150 | 11.0381 | -0.0009 | 0.0000 | 0.0011
200 | 9.1020 | -0.0008 | 0.0000 | 0.0021
RMSE | 50 | 73.4328 | 0.0621 | 0.0332 | 0.0739
100 | 45.1901 | 0.0433 | 0.0272 | 0.0519
150 | 30.6868 | 0.0367 | 0.0213 | 0.0454
200 | 25.7042 | 0.0312 | 0.0206 | 0.0374
Coverage | 50 | 95.0 | 95.0 | 93.2 | 95.3
100 | 94.2 | 95.5 | 94.1 | 95.4
150 | 94.0 | 93.9 | 94.0 | 94.1
200 | 95.7 | 93.4 | 94.1 | 94.5

##  5 Data analysis

In this section, we present two examples illustrating the diagnostic tools for
the MNBR model. The data sets and programs can be found in the MNB package
available on the R platform.

#### Seizures data

The data set described in Diggle et al. (2013) refers to an experiment in
which 59 epileptic patients were randomly assigned to one of two treatment
groups: treatment (progabide drug) and placebo groups. The number of seizures
experienced by each patient during the baseline period (week eight) and the
four consecutive periods (every two weeks) was recorded. The main goal of this
application is to analyze the drug effect with respect to the placebo. Two
dummy covariates are considered in this study; Group which assumes values
equal to 1 if the patient belongs to treatment group and 0 otherwise, and
Period which assumes values equal to 1 if the number of seizures are recorded
during the treatment and 0 if are measured in the baseline period. Taking into
account the irregular measurement of rate seizures during the time, the
variable Time is considered as an offset for fitting the data, where Time
assumes values equal 8 if the number of seizures is observed in the baseline
period and 2 otherwise. The individual profiles of the patients belonging to
the placebo and progabide groups are shown in Figure 1 (a) and (b),
respectively. The atypical individual profiles correspond to patient $\sharp
25(18,24,76,25)$ of the placebo group, who presented a high number of seizures
in the third visit compared to other clinic visits, and patient $\sharp
49(102,65,72,63)$ of the progabide group, who suffered a high number of
seizures in every clinic visit, indicating the ineffectiveness of the drug in
patients with complex seizures. In both groups, it is possible to see the
right skewness of the empirical distributions of individual profiles.

(a)

(b)

(c)

Figure 1: Individual profiles of the patients in (a) the placebo group, (b)
the progabide group, and (c) the simulated envelope plot of the Pearson
residuals.

Figure 1(c) shows a normal probability plot with simulated envelope for the
Pearson residual (Faraway, 2016, p.135) computed by fitting the Poisson
regression model to seizures data. We confirm the occurrence of the
overdispersion phenomenon. Thus, based on Figure 1(a) - (c), which shows the
asymmetric behavior of the empirical distribution of individual profiles, the
MNBR model is proposed for modeling this behavior according to the following
structure: $(i)$ $\bm{y}_{i}\stackrel{{\scriptstyle\rm ind}}{{\sim}}{\rm
MNB}(\bm{\mu}_{i},\phi)$ and $(ii)$ $\log(\mu_{ij})=\beta_{0}+\beta_{1}{\rm
Group}_{i}+\beta_{2}{\rm Period}_{ij}+\beta_{3}({\rm Group}_{i}\times{\rm
Period}_{ij})+\log({\rm Time}_{ij}),$ where
$\bm{y}_{i}=(y_{i1},y_{i2},y_{i3},y_{i4})^{\top},$
$\bm{\mu}_{i}=(\mu_{i1},\mu_{i2},\mu_{i3},\mu_{i4})^{\top},$ for
$i=1,\ldots,59$, $\beta_{1}$ is the logarithm of the ratio of the average rate
of the treatment group to the placebo group at baseline, $\beta_{2}$ is the
logarithm of the ratio of the seizure mean after treatment period to before
treatment period for the placebo group, and $\exp(\beta_{3})$ is the treatment
effect, and it is the ratio of post- to pre-treatment mean seizure ratios
between treatment and placebo groups. The parameter estimates obtained by
using the fit.MNB function from the MNB package are shown in Table 4. It is
not observed evidence of the treatment effect. The shape parameter estimate
indicates that the variability of individual profiles with respect their
average is asymmetric to the left with dispersion equal to 1.607
($\lambda=1/\sqrt{\phi}=0.789$).

Table 4: Parameter estimates with their respective approximate standard errors (Std. error), z-values, and $p$-values for the MNBR model fitted to the epileptic data. Parameter | Estimate | Std. error | z-value |  $p$-value
---|---|---|---|---
$\phi$ | 1.607 | 0.278 | $-$ | $-$
$\beta_{0}$ | 1.348 | 0.153 | 8.813 | $<0.001$
$\beta_{1}$ | 0.028 | 0.211 | 0.131 | 0.896
$\beta_{2}$ | 0.112 | 0.047 | 2.386 | 0.017
$\beta_{3}$ | -0.105 | 0.065 | -1.610 | 0.107
$\lambda$ | 0.789 |  |  |

The randomized quantile residual presented in Section 3.1 is used to
investigate the presence of outliers or any indication of lack of fit. Figure
2(a) and (b) shows the absence of extra variability and evidence that patient
$\sharp 49$ is atypical. The qMNB and envelope.MNB functions from the MNB
package is used to display the graphics of the randomized quantile residual.

(a)

(b)

Figure 2: Index plot of the randomized quantile residuals (a) and (b) the
simulated envelope plot of the randomized quantile residuals.

(a)

(b)

Figure 3: Index plot of (a) the generalized Cook distance and (b) the
likelihood displacement.

(a)

(b)

Figure 4: Index plot of $|\bm{d}_{max}|$ for the case weight perturbation
scheme: (a) local influence and (b) total local influence on the $i$th
patient.

(a)

(b)

Figure 5: Index plot of $|\bm{d}_{max}|$ for case weight perturbation scheme:
(a) local influence and (b) total influence of the $j$th measurement of each
patient.

The global influential graphics presented in Figure 3(a) and (b) reveal that
patients $\sharp 49$, and $\sharp 25$ can impact the ML estimates of the MNBR
model when they are removed from the model. The local and total local
influential graphics exhibited in Figure 4(a) and (b) indicate the sensitivity
of the estimate associated with patient $\sharp 49$ when the minor
perturbation is induced in the directions of $\bm{d}_{max}$ and $\bm{d}_{i},$
respectively. According to Figure 5(a) and (b), the number of epileptic
seizures of patient $\sharp 49$ recorded on the first and last clinic visits
of the trial are more influential. The dispersion perturbation scheme did not
show evidence of any possible influential observations.

Finally, the percentage relative deviations
$PRD$=$[(\widehat{\theta}-\widehat{\theta}^{*})/\widehat{\theta}]\times
100\%$, where $\widehat{\theta}^{*}$ is the estimator of $\theta$ obtained
after deleting one or more atypical subjects, are calculated and the results
are presented in Table 5. We can observe significatively changes are associate
with the estimates of $\beta_{1}$ and $\beta_{3}$. Noting that the signal of
the estimate $\beta_{1}$ is indicating that the treatment decreases the number
of seizures, its descriptive-level keeps nonsignificant though. For the
parameter $\beta_{3}$ we observe that the interaction effect incresed 50% and
that its descriptive-level becomes significant.

Table 5: Parameter estimates, standard errors (Std. error), $p$-values, and percentage relative deviations ($PRD$). Dropping | Parameter | Estimate | Std. error |  $p$-value |  $PRD($%$)$
---|---|---|---|---|---
| $\phi$ | 2.060 | 0.371 | $-$ | -28.21
| $\beta_{0}$ | 1.348 | 0.136 | $<0.000$ | 0.00
$\sharp 49$ | $\beta_{1}$ | -0.107 | 0.189 | 0.573 | 487.95
| $\beta_{2}$ | 0.112 | 0.047 | 0.017 | 0.00
| $\beta_{3}$ | -0.302 | 0.070 | $<0.000$ | -188.73

#### Accident data

Karlis (2003) provides a data set which refers to the number of car accidents
in 24 central roads in Athens, for a time period of 1987 to 1991. In this
paper, we are interested in modeling the number of car accidents per road
length during the period of 5 years. Thus, the dummy covariate years is
considered for fitting the data including the road length as an offset. Figure
6(a) presents the individual profiles of the central roads for the number of
car accidents recorded on each road over time. In general, the highest number
of car accidents occurred in 1989 and the smallest in 1991. The empirical
distribution of individual profiles is skewed to the right. Figure 6(b)
reveals that $Patision$ and $Peiraios$ roads are considered as outliers and
evidence of variability between the individual profiles and across the years
is suggested in Figure 6(c).

(a)

(b)

(c)

Figure 6: (a) Individual profiles, (b) box-plots of the number of car
accidents on central Athens roads during the 5 years, and (c) simulated
envelope plot of Pearson residuals.

In the fitting the Poisson regression model to the accidents data, Figure
6(c), we observe the occurrence of the overdispersion phenomenon. Based on
Figure 6(a) - (b) showing the asymmetric behavior of the empirical
distribution of individual profiles, we propose the MNBR model with the
following structure for modeling the data: $(i)$
$\bm{y}_{i}\stackrel{{\scriptstyle\rm ind}}{{\sim}}{\rm
MNB}(\bm{\mu}_{i},\phi)$ and $(ii)$ $\log(\mu_{ij})=\beta_{0}+\beta_{t}{\rm
year}_{j}+\log(\text{length}_{i})$, where
$\bm{y}_{i}=(y_{i1},y_{i2},y_{i3},y_{i4},y_{i5})^{\top}$,
$\bm{\mu}_{i}=(\mu_{i1},\mu_{i2},\mu_{i3},\mu_{i4},\mu_{i5})^{\top}$,
$t=1,2,3,4$, $\beta_{t}$ is the logarithm of the ratio of the average car
accidents per road length for year $t$ to year 1987. The ML estimates for the
MNBR model are presented in Table 6. We observe the estimates are significant
and that the small difference in the logarithm of the rate average of the
number of car accidents per length of the roads occurred between the years the
1987 and 1991.

Table 6: Parameter estimates, standard errors (Std. error), z-values, and $p$-values for the MNBR model fitted to the accidents data. Parameter | Estimate | Std. error | z-value |  $p$-value
---|---|---|---|---
$\phi$ | 3.364 | 0.954 | 3.526 | $-$
$\beta_{0}$ | 2.250 | 0.119 | 18.899 | $<0.001$
$\beta_{1}$ | 0.366 | 0.053 | 6.890 | $<0.001$
$\beta_{2}$ | 0.586 | 0.051 | 11.490 | $<0.001$
$\beta_{3}$ | 0.474 | 0.052 | 9.102 | $<0.001$
$\beta_{4}$ | -0.322 | 0.063 | -5.112 | $<0.001$

The randomized quantile residuals in Figure 7(a) shows the absence of extra
variability and atypical road, and Figure 7(b) shows a normal probability plot
with simulated envelope for the randomized quantile residual.

(a)

(b)

Figure 7: Index plot of the randomized quantile residuals (a) and (b) the
simulated envelope plot of the randomized quantile residuals.

The global influential graphic in Figure 8(a) indicates that
$Leof.~{}Kavalas$, $Leof.~{}Athinon$ and $Panepistimiou$ can impact the ML
estimates when they are removed. $Leof.~{}Kavalas$ (4,6,12,9,3) and
$Panepistimiou$ (24, 58, 40, 36, 05) are small roads, 2.0 and 1.1 kilometers
long, respectively. $Leof.~{}Athinon$ (15,11,16,21,28) is a large road, 6.1
kilometers long. Furthermore, two roads that can impact the ML estimates are
showed in Figure 8(b), $Patision$ (80,108,114,113,86) and
$Katehaki$(2,3,27,24,7), which are 4.1 and 1.4 kilometers long, respectively.
In Figure 9(a) and (b) are showed influential roads, $Peraios$ (86 89 109 90
49) and $Aharnon$ (44, 79, 91, 88, 33), which are 8.0 and 5.5 kilometers long,
respectively.

(a)

(b)

Figure 8: Index plot of (a) the generalized Cook distance and (b) the
likelihood displacement for the number of car accidents on central Athenian
roads in 1989-1991.

(a)

(b)

Figure 9: Index plot of $|\bm{d}_{max}|$ for the case weight perturbation
scheme: (a) local influence and (b) total local influence for the number of
car accidents on central Athenian roads in 1989-1991.

(a)

(b)

Figure 10: Index plot of $|\bm{d}_{max}|$ for the case weight perturbation
scheme: (a) local influence and (b) total local influence on $j$th measurement
of the number of car accidents for each atypical central Athenian road.

(a)

(b)

Figure 11: Index plot of $|\bm{d}_{max}|$ for the dispersion perturbation
scheme, (a) local influence and (b) total local influence for the number of
car accident on central Athenian roads of 1989-1991.

Finally, the percentage relative deviations ($PRD$) are calculated and we
verify that the most significant changes in the ML estimates are associated
with the removal of the influential roads $Patision$ and $Peiraios$. The
results in Table 7 show that the estimates are affected when $Patision$ and
$Peiraios$ are removed simultaneously.

Table 7: Parameter estimates, standard errors (Std. Error), z-values, $p$-values, and $PRD(\%)$ when $Patision$ road is deleted from the accidents data. Dropping | Parameter | Estimate | Std. Error | z-value |  $p$-value | $PRD(\%)$
---|---|---|---|---|---|---
$Patision$ | $\phi$ | 3.507 | 1.021 | $-$ | $-$ | $-$
| $\beta_{0}$ | 2.199 | 0.120 | 18.291 | $<0.001$ | 2.27
| $\beta_{1}$ | 0.376 | 0.057 | 6.596 | $<0.001$ | -2.65
| $\beta_{2}$ | 0.617 | 0.054 | 11.331 | $<0.001$ | -5.36
| $\beta_{3}$ | 0.492 | 0.055 | 8.831 | $<0.001$ | -3.86
| $\beta_{4}$ | -0.399 | 0.069 | -5.765 | $<0.001$ | -23.97
$Peiraios$ | $\phi$ | 3.235 | 0.934 | $-$ | $-$ | $-$
| $\beta_{0}$ | 2.215 | 0.124 | 17.795 | $<0.001$ | 1.54
| $\beta_{1}$ | 0.413 | 0.057 | 7.249 | $<0.001$ | -12.63
| $\beta_{2}$ | 0.634 | 0.056 | 11.605 | $<0.001$ | -8.23
| $\beta_{3}$ | 0.531 | 0.055 | 9.534 | $<0.001$ | -11.98
| $\beta_{4}$ | -0.287 | 0.067 | -4.257 | $<0.001$ | 10.91
$Peiraios~{}and~{}Patision$ | $\phi$ | 3.366 | 0.934 | $-$ | $-$ | $-$
| $\beta_{0}$ | 2.154 | 0.126 | 17.070 | $<0.001$ | 4.24
| $\beta_{1}$ | 0.432 | 0.062 | 7.001 | $<0.001$ | -17.95
| $\beta_{2}$ | 0.678 | 0.059 | 11.491 | $<0.001$ | -15.76
| $\beta_{3}$ | 0.561 | 0.060 | 9.319 | $<0.001$ | -18.47
| $\beta_{4}$ | -0.370 | 0.075 | -4.924 | $<0.001$ | -14.91

##  6 Discussion

In this study, we developed diagnostic tools for the MNBR model derived from
the Poisson mixed model where the GLG distribution is assumed for the random
effect. Simulation results exhibited the features of the MNBR model and its
association with the hierarchical model, which was essential for explaining
that the asymptotic consistency of the estimators of the regression
coefficients and the dispersion parameter depends on the asymmetry of the GLG
distribution. As expected, the MNBR model provides inconsistent estimates of
the asymptotic variance of the ML estimators when the covariance matrix of the
response variable is misspecified. The randomized quantile residuals can be
used to assess possible departures of the data from the MNBR model
assumptions. Following the approaches of Cook (1977, 1986) and Lesaffre and
Verbeke (1998), global and local measures for the MNBR model were derived and
implemented in the authors’ MNB package. The application of the MNB package to
two data sets was presented and it was shown that the asymmetric behavior of
the empirical distributions of the individual profiles can indicate the need
to use the multivariate model to better fit these profiles. The proposed
methodology was helpful for identifying outlying subjects that matter in the
MNBR model over time and handle the overdispersion phenomenon. The code for
installing the MNB package is presented in the Appendix.

## Acknowledgements

The work of the fourth author is partially funded by CNPq, Brazil. We also
thank anonymous referees for constructive comments and suggestions.

## References

  * Cook (1977) Cook, R. D. (1977).  Detection of influential observations in linear regression.  Technometrics, 19, 15–18.
  * Cook (1986) Cook, R. D. (1986).  Assessment of local influence (with discussion).  Journal of the Royal Statistical Society Series B, 48, 133–169.
  * Cook and Weisberg (1983) Cook, R. D. and Weisberg, S. (1983).  Residuals and Influence in Regresion.  Chapman and Hall, New York.
  * Diggle et al. (2013) Diggle, P. J., Liang, K. Y., and Zeger, S. L. (2013).  Analysis of Longitudinal Data.  Oxford University Press, N.Y., 2 edition.
  * Dunn and Smyth (1996) Dunn, P. K. and Smyth, G. K. (1996).  Randomized quantile residuals.  Journal of Computational and Graphical Statistics, 5, 236–244.
  * Fabio et al. (2012) Fabio, L., Paula, G. A., and de Castro, M. (2012).  A Poisson mixed model with nonormal random effect distribution.  Computational Statistics and Data Analysis, 56, 1499–1510.
  * Faraway (2016) Faraway, F. (2016).  Extending the Linear Model with R: Generalized Linear, Mixed Effects and nonparametric regression models.  Chapman & Hall/CRC, New York.
  * Hilbe (2011) Hilbe, J. M. (2011).  Negative Binomial Regression.  Cambridge, United Kingdom.
  * Johnson et al. (1997) Johnson, N., Kotz, S., and Balakrishnan, N. (1997).  Discrete Multivariate Distributions.  Wiley, New York.
  * Karlis (2003) Karlis, D. (2003).  An EM algorithm for multivariate Poisson distribution and related models.  Journal of Applied Statistics, 30, 63–77.
  * Lawless (1987) Lawless, J. (1987).  Negative binomial and mixed Poisson regression.  The Canadian Journal of Statistics, 15, 209–225.
  * Lawless (2002) Lawless, J. F. (2002).  Statistical Models and Methods for Lifetime Data.  Wiley-Interscience, New York.
  * Lee et al. (2006) Lee, Y., Nelder, J. A., and Pawitan, Y. (2006).  Generalized linear models with random effects. Unified analysis via H-likelihood.  Chapman and Hall/CRC, Boca Raton.
  * Lesaffre and Verbeke (1998) Lesaffre, E. and Verbeke, G. (1998).  Local influence in linear mixed models.  Biometrics, 54, 570–582.
  * Molenberghs and Verbeke (2010) Molenberghs, G. and Verbeke, G. (2010).  Models for Discrete Longitudinal Data.  Springer, New York.
  * Molenberghs et al. (2007) Molenberghs, G., Verbeke, G., and Demétrio, C. G. B. (2007).  An extended random-effects approach to modeling repeated, overdispersed count data.  Lifetime Data Analysis, 13, 513–531.
  * Moller and Rubak (2010) Moller, J. and Rubak, E. (2010).  A model for positively correlated count variables.  International Statistical Review, 78, 65–80.
  * R Core Team (2017) R Core Team (2017).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Solis-Trapala and Farewell (2005) Solis-Trapala, I. and Farewell, V. (2005).  Regression analysis of overdispersed correlated count data with subject specific covariates.  Statistics in Medicine, 24, 2557–2575.
  * Sung and Lee (2018) Sung, Y. and Lee, K. (2018).  Negative binomial loglinear mixed models with general random effects covariance matrix.  Communications for Statistical Applications and Methods, 25, 61–70.
  * Tsui (1986) Tsui, K.-W. (1986).  Multiparameter estimation for some multivariate discrete distributions with possibly dependent components.  Annals of the Institute of Statistical Mathematics, 38, 45–56.

## Appendix

## 1\. The score vector and the observed Fisher information matrix

The score vector $\bm{U}_{\theta}=(\bm{U}_{\beta}^{\top},U_{\phi})^{\top}$ is
obtained by deriving the log-likelihood function (2) with respect to the
parameters $\bm{\beta}$ and $\phi$, respectively. Thus,

$\displaystyle\bm{U}_{\beta}$
$\displaystyle=\sum_{i=1}^{n}\bm{X}_{i}^{\top}(\bm{y}_{i}-a_{i}\bm{\mu}_{i})\quad\text{and}$
$\displaystyle U_{\phi}$
$\displaystyle=\sum_{i=1}^{n}\bigg{\\{}\psi(\phi+y_{i+})-\psi(\phi)+\log\left(\frac{\phi}{\phi+\mu_{i+}}\right)+\left(\frac{\mu_{i+}-y_{i+}}{\phi+\mu_{i+}}\right)\bigg{\\}},$
(4)

where $a_{i}=(\phi+y_{i+})/(\phi+\mu_{i+}),$ $\psi(\cdot)$ is the digamma
function, and $\bm{X}_{i}$ is an $m_{i}\times p$ matrix with row
$\bm{x}^{\top}_{ij}$ for $i=1,\ldots,n$ and $j=1,\ldots,m_{i}.$ However, using
the fact that
$\Gamma(\phi+y_{i+})/\Gamma(\phi)=\phi(\phi+1)(\phi+2)\ldots(\phi+y_{i+}-1)$
in (4), $U_{\phi}$ takes the following form

$U_{\phi}=\sum_{i=1}^{n}\left\\{\sum_{j=0}^{y_{i+-1}}(j+\phi)^{-1}-\frac{y_{i+}}{\phi+\mu_{i+}}-\log(1+\phi^{-1}\mu_{i+})+\frac{\mu_{i+}}{(\phi+\mu_{i+})}\right\\}.$

The observed information matrix is obtained by deriving the score vector with
respect to $\bm{\beta}$ and $\phi$. Thus,

$-\ddot{\ell}(\bm{\theta})=\begin{bmatrix}-\ddot{\ell}_{\beta\beta}&-\ddot{\ell}_{\beta\phi}\\\
-\ddot{\ell}_{\phi\beta}&-\ddot{\ell}_{\phi\phi},\end{bmatrix}$

where

$\ddot{\ell}_{\beta\beta}=-\sum_{i=1}^{n}\frac{(\phi+y_{i+})}{(\phi+\mu_{i+})}\sum_{j=1}^{m_{i}}\bm{x}_{ij}\bm{x}^{\top}_{ij}\mu_{ij}+\sum_{i=1}^{n}\frac{(\phi+y_{i+})}{(\phi+\mu_{i+})^{2}}\sum_{j=1}^{m_{i}}\bm{x}_{ij}\mu_{ij}\sum_{j=1}^{m_{i}}\bm{x}^{\top}_{ij}\mu_{ij},$

$\ddot{\ell}_{\beta\phi}=\ddot{\ell}_{\phi\beta}=-\sum_{j=1}^{m_{i}}x_{ijk}\mu_{ij}\frac{(\mu_{i+}-y_{i+})}{(\phi+\mu_{i+})^{2}}\
\ {\rm and}$

$\ddot{\ell}_{\phi\phi}=-\sum_{i=1}^{n}\left\\{\sum_{s=0}^{y_{i}^{*}}\frac{1}{(s+\phi)^{2}}+\frac{\phi^{-1}\mu_{i+}}{(\phi+\mu_{i+})}-\frac{\mu_{i+}}{(\phi+\mu_{i+})^{2}}-\frac{y_{i+}}{(\phi+\mu_{i+})^{2}}\right\\}.$

## 2\. How to install and use the MNB package

We present the R platform code used to install the MNB package.

⬇

1require(devtools)

2devtools::install_github("carrascojalmar/MNB")

### 2.1. Codes for Seizures data

⬇

1require(MNB)

2data(seizures)

3

4#

5star <-list(phi=1, beta0=1, beta1=1, beta2=1, beta3=1)

6mod <- fit.MNB(formula=Y ~ trt + period +

7trt:period + offset(log(weeks)), star=star, dataSet=seizures,tab=FALSE)

8mod

9

10#

11par <- mod$par

12names(par)<-c()

13res.q <- qMNB(par=par,formula=Y ~ trt + period + trt:period +

14offset(log(weeks)),dataSet=seizures)

15

16plot(res.q,ylim=c(-3,4.5),ylab="Randomized quantile residual",

17xlab="Index",pch=15,cex.lab = 1.5, cex = 0.6, bg = 5)

18abline(h=c(-2,0,2),lty=3)

19

20#

21envelope.MNB(formula=Y ~ trt + period + trt:period +

22offset(weeks),star=star,nsim=21,n.r=6,

23dataSet=seizures,plot=TRUE)

24

25#

26global.MNB(formula=Y ~ trt + period +

27trt:period + offset(log(weeks)),star=star,dataSet=seizures,plot=TRUE)

28

29#

30local.MNB(formula=Y ~ trt + period + trt:period +
offset(log(weeks)),star=star,dataSet=seizures,

31schemes="weight",plot=TRUE)

### 2.2. Codes for Accident data

⬇

1

2require(MNB)

3accident <- read.table("...\\\accident.txt",h=TRUE)

4temp <- apply(accident[,2:6],1,t)

5Y <- matrix(temp,ncol=1)

6length.x <- rep(accident$length,each=5)

7ind <- rep(1:24,each=5)

8year <- as.factor(rep(1:5,24))

9accident2 <- data.frame(Y,year,length.x,ind)

10

11#

12star=list(phi=10, beta1=2, beta2=0.4, beta3= 0.6, beta4=0.5, beta5=-0.3)

13mod <- fit.MNB(Y ~ year +
offset(log(length.x)),star=star,dataSet=accident2,tab=FALSE)

14mod

15

16#

17par <- mod$par

18names(par)<-c()

19qMNB(par=par,formula=Y ~ year + offset(log(length.x)),dataSet=accident2)

20

21#

22envelope.MNB(formula=Y ~ year +
offset(log(length.x)),star=star,nsim=21,n.r=6,

23dataSet=accident2,plot=TRUE)

24

25#

26global.MNB(Y ~ year +
offset(log(length.x)),star=star,dataSet=accident2,plot=TRUE)

27

28#

29local.MNB(Y ~ year + offset(log(length.x)),star=star,dataSet=accident2,

30schemes="weight",plot=TRUE)


# The Turán number of Berge-$K_{4}$ in triple systems

András Gyárfás Renyi Institute of Mathematics, Budapest, HungaryResearch
supported in part by NKFIH Grant No. K116769.

###### Abstract

A Berge-$K_{4}$ in a triple system is a configuration with four vertices
$v_{1},v_{2},v_{3},v_{4}$ and six distinct triples $\\{e_{ij}:1\leq i<j\leq
4\\}$ such that $\\{v_{i},v_{j}\\}\subset e_{ij}$ for every $1\leq i<j\leq 4$.
We denote by $\mathcal{B}$ the set of Berge-$K_{4}$ configurations. A triple
system is $\mathcal{B}$-free if it does not contain any member of
$\mathcal{B}$. We prove that the maximum number of triples in a
$\mathcal{B}$-free triple system on $n\geq 6$ points is obtained by the
balanced complete $3$-partite triple system: all triples $\\{abc:a\in A,b\in
B,c\in C\\}$ where $A,B,C$ is a partition of $n$ points with

$\left\lfloor{n\over 3}\right\rfloor=|A|\leq|B|\leq|C|=\left\lceil{n\over
3}\right\rceil.$

##  1 Introduction

Extending the notion of Berge paths and cycles, Gerbner and Palmer [7]
introduced the notion of Berge-$G$ hypergraphs as follows. A hypergraph
$H=(V,F)$ is called Berge-$G$ if $G=(V,E)$ and there exists a bijection
$g:E(G)\mapsto E(H)$ such that for $e\in E(G)$ we have $e\subset g(e)$. For a
fixed $G$, let $\mathcal{B}(G)$ denote the family of Berge-$G$ hypergraphs.
For a family $\mathcal{F}$ of $r$-uniform hypergraphs, the Turán number ${\rm
ex}_{r}(n,\mathcal{F})$ is the largest number of edges in an $r$-uniform
hypergraph on $n$ vertices that does not contain any member from a family
$\mathcal{F}$ as a subhypergraph.

Many asymptotic results are known for ${\rm ex}_{r}(n,\mathcal{B}(G))$, Győri
and G. Y. Katona [11], Győri, Lemons [12], Gerbner and Palmer [7], Füredi and
Özkahya [4], Gerbner, Methuku and Vizer [6], Palmer, Tait, Timmons and Wagner
[16], Grósz, Methuku, Tompkins. [9]. Results of Alon and Shikhelman [1],
Gerbner and Palmer [8] also imply bounds for ${\rm ex}_{r}(n,\mathcal{B}(G))$.

In case of $G=K_{3}$ a spectacular sharp result is that ${\rm
ex}_{3}(n,\mathcal{B}(K_{3}))=\lfloor{n^{2}\over 8}\rfloor$ for $n\geq 3$.
This was proved (among other results) independently by Győri [10] and by
Frankl, Füredi and Simonyi [3] with an elegant reduction to Mantel’s theorem.
For $t\geq 13$, Maherani and Shahsiah [13] determined exactly ${\rm
ex}_{3}(n,\mathcal{B}(K_{t}))$.

Here we focus on ${\rm ex}_{3}(n,\mathcal{B}(K_{4}))$ and from now on
$\mathcal{B}=\mathcal{B}(K_{4})$. We refer to $3$-uniform hypergraphs as
triple systems. A (special case of a) result of Mubayi [14] is that ${\rm
ex}_{3}(n,\mathcal{B})$ is asymptotic to the maximum number of triples in a
$3$-partite triple system on $n$ vertices, i.e.

${\rm ex}_{3}(n,\mathcal{B})\sim f(n)=\left\lfloor{n\over
3}\right\rfloor\left\lfloor{n+1\over 3}\right\rfloor\left\lfloor{n+2\over
3}\right\rfloor.$

Mubayi’s result was sharpened by Pikhurko [17] who proved that ${\rm
ex}_{3}(n,\mathcal{B})=f(n)$, for large enough $n$. In fact, both results were
proved in a stronger form, only one specific member of $\mathcal{B}$ was
forbidden, the expansion, where the bijection $g$ maps each edge $e\in K_{4}$
to a triple $e\cup v_{e}$ so that the vertices $v_{e}$ are all different and
different from the vertices of $K_{4}$ as well. Thus the expansion of $K_{4}$
has six triples on ten vertices, it is the largest member of $\mathcal{B}$.

We prove that ${\rm ex}_{3}(n,\mathcal{B})$ is equal to $f(n)$ for every
$n\geq 6$. For the small cases, ${\rm ex}_{3}(3,\mathcal{B})=1$, ${\rm
ex}_{3}(4,\mathcal{B})=4$ are obvious, ${\rm ex}_{3}(5,\mathcal{B})=5$ will be
proved in Theorem 10.

######  Theorem 1.

For $n\geq 6$, ${\rm ex}_{3}(n,\mathcal{B})=f(n)$.

The proof of Theorem 1 does not use that all configurations of $\mathcal{B}$
are excluded. In fact, the exclusion of the “maximal” member, the extension of
$K_{4}$, is not used. This gives the following corollary.

######  Corollary 2.

Let $B^{*}$ denote the expansion of $K_{4}$. Then ${\rm
ex}_{3}(n,\mathcal{B}\setminus B^{*})=f(n)$ for $n\geq 6$.

Note that Pikhurko’s result and Corollary 2 shows that there is no unique
minimal subset $\mathcal{B}^{\prime}\subset\mathcal{B}$ for which ${\rm
ex}_{3}(n,\mathcal{B}^{\prime})=f(n)$ for large enough $n$.

Assume that $H$ is a triple system and $A$ is a proper nonempty subset of
$V(H)$. The trace of $H$ is a multigraph $G_{H}(A)$ (multiloops allowed) on
vertex set $V(H)\setminus A$. Every triple $xyz$ of $H$ with $x,y\in A,z\in
V(H)\setminus A$ defines a loop on $z$ with label $\\{x,y\\}$ and every triple
$xyz$ of $H$ with $x\in A,y,z\in V(H)\setminus A$ defines an edge $yz$ with
label $\\{x\\}$. If $H$ is clear from the context, we use simply $G(A)$
instead of $G_{H}(A)$. The label of a loop or an edge is denoted by
$\ell(v,v),\ell(v,w)$, respectively. Let $\mu(v),\mu(e)$ denote the
multiplicity of a loop $vv$ and an edge $e$, respectively. When $|A|=1$, the
trace is a simple graph, usually called the link of $A$. In a triple system
the degree of a vertex $v$ is the number of triples containing $v$ and denoted
by $d(v)$.

Theorem 1 is proved by combining induction on $n$ with ${\rm
ex}(m,K_{4})=\lfloor{m^{2}\over 3}\rfloor$ (Turán theorem) applied to a
subgraph of $G=G(A)$ where $A$ is a triple in a suitable Berge-triangle of a
$\mathcal{B}$-free $H_{n}$. I have been recently informed by D. Gerbner that a
different approach [5] also implies Theorem 1 for $n\geq 9$.

##  2 Proof of Theorem 1 and its corollary

###### Proof.

We prove Theorem 1 by induction. The base case $n=6$ (together with ${\rm
ex}_{3}(5,\mathcal{B})=5$) is proved in Theorem 10. Assuming the theorem is
true for $n-1\geq 6$, let $H_{n}$ be a $\mathcal{B}$-free triple system with
$n$ vertices and $f(n)+1$ edges. From this we will get a contradiction.

We can apply induction if for some vertex $v$ of $H_{n}$, $d(v)\leq
f(n)-f(n-1)$. Indeed, then deleting $v$ and the triples containing $v$ we get
a $B(K_{4})$-free triple system $H^{\prime}$ with $n-1$ vertices and more than
$f(n-1)$ edges, contradiction. Thus we may assume that the minimum degree
$\delta=\delta(H_{n})$ is larger than $f(n)-f(n-1)$. The actual values of
$f(n)-f(n-1)$ are stated in the next observation.

######  Observation 3.

$f(n)-f(n-1)=\left\\{\begin{array}[]{ll}k^{2}+2&\mbox{if $n=3k$}\\\
k^{2}&\mbox{if $n=3k+1$}\\\ k^{2}+k&\mbox{if $n=3k+2$}.\end{array}\right.$

Two specific members of $\mathcal{B}(K_{3})$ are $K_{4}^{3}-e$, three triples
within four vertices and the tight path, with triples $abc,bcd,cde$ on vertex
set $\\{a,b,c,d,e\\}$.

######  Proposition 4.

$H_{n}$ contains either a $K_{4}^{3}-e$ or a tight path.

###### Proof.

If the statement is not true then the trace $G=G(\\{v\\})$ of any vertex $v\in
V(H_{n})$ contains no triangles or $P_{4}$-s, thus $G$ must be a star forest.
Therefore $f(n)-f(n-1)<\delta\leq n-2$ and from Observation 3 we get
contradiction for every $n\geq 6$ except for $n=7$. However, in this case we
have $f(7)+1=13$ triples in $H_{7}$ thus there exists a vertex of degree at
least six. Applying the argument to this vertex we get a contradiction. ∎

Using Proposition 4, we can select a Berge-triangle $B$ in $H_{n}$ such that

$E(B)=\\{123,12x,23y\\}$

where $x,y\in V(H_{n})\setminus[3]$. If there exists $K_{4}^{3}-e$ in $H_{n}$
then we assume $x=y$, otherwise $x\neq y$.

Consider $G=G(B)=G(\\{1,2,3\\})$, the trace of $H_{n}$ on
$V(H_{n})\setminus[3]$. By our assumptions, there is a loop on $x$ with label
$\\{1,2\\}$ and a loop on $y$ with label $\\{2,3\\}$. If $x=y$ then there is a
double loop on $x$ with labels $\\{1,2\\},\\{2,3\\}$, respectively. Set
$Z=V(G)\setminus\\{x,y\\}$.

We partition $Z$ into eight sets $Z_{I}$ where $I\subseteq[3]$; a vertex $v$
belongs to $Z_{I}$ if the union of the labels on the edges containing $v$ is
$I$.

Vertices of $Z_{\emptyset}$ are isolated in $G$. Vertices of
$Z_{1},Z_{2},Z_{3}$ cannot be incident to loops or multiple edges of $G$.
Edges between $Z_{ij}$ and $Z_{ik}$ are not multiedges and labeled with $i$
for any choice of three different indices.

######  Proposition 5.

Assume that $z\in Z$ and $e_{1},e_{2},e_{3}$ are distinct edges of $G$
containing $z$. Then the sets $\ell(e_{1}),\ell(e_{2}),\ell(e_{3})$ have no
distinct representatives.

###### Proof.

Indeed, otherwise $\\{1,2,3,z\\}$ would span a member of $\mathcal{B}$: the
pairs of $\\{1,2,3\\}$ are covered by the three triples of the Berge-triangle
$B$ and the pairs $1z,2z,3z$ can be covered by the triples containing the
three $e_{i}$-s. This is a contradiction. ∎

We collect some consequences of Proposition 5.

######  Proposition 6.

The following properties hold.

  * •

1\. Loops and non-loop edges in $Z$ have multiplicity at most two.

  * •

2\. Double loops in $Z$ form one-vertex components in $G$ (inside $Z_{123}$).

  * •

3\. Multiple edges within $Z$ are inside some $Z_{ij}$.

  * •

4\. Assume $v_{1}v_{2},v_{2}v_{3}$ are double edges in $Z_{ij}$. Then no
further edges or loops can be on $v_{1}$ or on $v_{3}$.

###### Proof.

Properties 1-3 are immediate consequences of Proposition 5. By symmetry, it is
enough to prove property 4 for $i=1,j=2$ and for $v_{3}$. Suppose
$v_{1}v_{2},v_{2}v_{3}$ are double edges in $Z_{12}$ and $v_{3}z\in E(G)$
where $z\neq v_{2}$ and w.l.o.g. the label of $z$ is $\\{1\\}$. Then we have a
member of $B(K_{4})$ spanned by $\\{1,2,v_{2},v_{3}\\}$ shown by the
assignments

$1v_{2}\mapsto 1v_{1}v_{2},1v_{3}\mapsto 1v_{3}z,2v_{2}\mapsto
2v_{1}v_{2},2v_{3}\mapsto 2v_{2}v_{3},v_{2}v_{3}\mapsto 1v_{2}v_{3},12\mapsto
123,$

leading to contradiction. ∎

Let $G^{*}$ be the simple graph obtained from $G$ by replacing its multiple
edges by single edges and removing all loops. Obviously
$|E(G^{*})|=|E(G)|-s(G)$ where

$s(G)=\sum_{v\in V(G)}\mu(v,v)+\sum_{e\in E(G)}(\mu(e)-1).$

The surplus subgraph $S$ of $G$ has vertex set $V(G)$ and its edges are the
(possible multiple) loops and multiple edges. The connected components of $S$
are called blocks. A block $Q$ is a bad block if $s(S[Q])>|V(Q)|$, otherwise
it is a good block. For example, a block containing a single loop is good, a
block which has a double edge with a loop on both ends is bad.

We use the bad blocks to determine the bad connected components of $G$,
defined similarly: a connected component $C$ is a bad component if
$s(G[C])>|V(C)|$, otherwise it is a good component.

######  Proposition 7.

Every bad component contains at least one bad block.

###### Proof.

In a connected component $C$ the blocks $Q_{1},\dots,Q_{m}$ are vertex
disjoint. If all blocks are good, then

$s(G[C])=\sum_{i=1}^{m}s(G[Q_{i}])\leq\sum_{i=1}^{m}|V(G[Q_{i}])|\leq|V(C)|$

thus $C$ would be a good component, contradiction. ∎

######  Lemma 8.

With a suitable choice of the Berge-triangle $B$, a bad component $C$ is one
of the following.

  * •

1\. A triple loop on $x$, implying $y=x,s(C)=3,|V(C)|=1$.

  * •

2\. An $m$-star on $x$, defined as follows. There is a double loop on $x$
(labeled by $\\{1,2\\},\\{2,3\\}$ consequently $y=x$). There are $m-1\geq 2$
double edges $xz_{1},xz_{2},\dots,xz_{m-1}$. The set
$\\{z_{1},\dots,z_{m-1}\\}$ is independent in $G$, i.e. contains no edges of
$G$. Thus $s(C)=m+1,|V(C)|=m\geq 3$.

  * •

3\. A double loop, $s(C)=2,|V(C)|=1$.

  * •

4\. A dumbbell: a double edge with single loops on both ends or a double loop
on one end, $s(C)=3,|V(C)|=2$.

###### Proof.

We prove the following statements.

  * •

1\. For any choice of $B$, a bad component $C$ not containing $x$ or $y$ is a
dumbbell.

Assume $C$ is a bad component not containing $x$ or $y$. Consider any bad
block $Q$ in $C$. By proposition 6 (2 and 3), $Q$ is either a double loop
component in $Z_{123}$ or $Q$ is inside some $Z_{ij}$, say inside $Z_{12}$.
Clearly $Q$ must contain some double edges and by Proposition 6 (4) the double
edges of $Q$ form a star. If the star has one double edge then both ends must
be a loop otherwise $Q$ is not bad. We claim that this is the only possibility
for a bad block. Indeed, if the star has at least two double edges, then there
must be at least one loop at an endpoint of the star and this is excluded in
Proposition 6 (4). Thus every bad block is a dumbbell. However the dumbbells
are connected components at the same time, since any edge, say $zz_{1}$ with
label $\\{1\\}$ touching the dumbbell $z_{1}z_{2}$, would define a member of
$\mathcal{B}$ spanned by $\\{1,2,z_{1},z_{2}\\}$ with the following
assignments:

$12\mapsto 123,z_{1}z_{2}\mapsto 1z_{1}z_{2},1z_{1}\mapsto
1z_{1}z,2z_{1}\mapsto 12z_{1},1z_{2}\mapsto 12z_{2},2z_{2}\mapsto
2z_{1}z_{2}.$

  * •

2\. Assume there exists a $K_{4}^{3}$ in $H_{n}$. Then, for some choice of
$B\subset K_{4}^{3}$ and $x\in V(B)$, the bad component containing $x$ is a
triple loop.

Assume that $x$ contains a triple loop, i.e. $\\{1,2,3,x\\}$ spans $K_{4}^{3}$
(and $x=y$). If the triple loop is not a component of $G$ then there exists an
edge $xz\in E(G)$ say with label $\\{1\\}$.

We claim that no vertex of $Z$ has a loop. Indeed, suppose that $z_{1}\in Z$
has a loop ($z_{1}=z$ is possible). By symmetry, we may assume that the label
of the loop is either $\\{1,2\\}$ or $\\{2,3\\}$. In both cases
$\\{1,2,3,x\\}$ would span a member of $\mathcal{B}$ with the following
assignments.

$1x\mapsto 1xz,12\mapsto 12z_{1},13\mapsto 123,23\mapsto 23x,2x\mapsto
12x,3x\mapsto 13x,$ $1x\mapsto 1xz,12\mapsto 12x,13\mapsto 123,23\mapsto
23z_{1},2x\mapsto 23x,3x\mapsto 13x,$

proving the claim.

Consider $G^{\prime}=G(1,2,x)$, now vertex $3$ plays the role of $x$. If the
triple loop on $3$ is not a component of $G^{\prime}$, there exists
$3z^{\prime}\in E(G^{\prime})$ ($z=z^{\prime}$ is possible). The label of this
edge can be one of $\\{x,1,2\\}$. For all the three choices $\\{1,2,3,x\\}$
spans a member of $\mathcal{B}$ with the following assignments (the first
assignment shows the choices).

$3x\mapsto 3xz^{\prime},2x\mapsto 23x,1x\mapsto 1xz,23\mapsto 123,13\mapsto
13x,12\mapsto 12x;$ (1) $13\mapsto 13z^{\prime},2x\mapsto 23x,1x\mapsto
1xz,23\mapsto 123,3x\mapsto 13x,12\mapsto 12x;$ $23\mapsto
23z^{\prime},2x\mapsto 23x,1x\mapsto 1xz,3x\mapsto 13x,13\mapsto 123,12\mapsto
12x,$

leading to contradiction. Thus in $G$ or in $G^{\prime}$ the triple loop is
the (bad) component containing $x$.

  * •

3\. Assume that $H_{n}$ contains $K_{4}^{3}-e$ but does not contain
$K_{4}^{3}$. Then, for some choice of $B\subset K_{4}^{3}-e$ and $x\in V(B)$,
the bad component containing $x$ is a double loop or a dumbbell or an
$m$-star.

Starting from a $K_{4}^{3}-e$, we have a double loop on $x$, with labels
$\\{1,2\\},\\{2,3\\}$ (again, $x=y$ in this case). We determine first the
block $Q$ containing $x$. If no multiedge is incident to $x$ then $Q$ is the
double loop on $x$. Suppose we have a double edge $xz$ with label $\\{i,j\\}$
in $E(G)$.

We claim that no other loops or edges can be incident to $z$. Indeed, a loop
on $z$ must be labeled with $\\{i,j\\}$ and the label of an edge on $zz_{1}$
must be $\\{i\\}$ or $\\{j\\}$, say $\\{i\\}$, because of Proposition 5.
Therefore $\\{i,j,x,z\\}$ would span a member of $\mathcal{B}$. Indeed,
$iz\mapsto ijz$ if we have a loop on $z$, $iz\mapsto izz_{1}$ if we have an
edge $zz_{1}$. Then $jz\mapsto jxz,xz\mapsto ixz$ and the triple
$\\{ix,jx,ij\\}$ has obviously a bijection to $123,12x,23x$. Thus we get a
contradiction.

Therefore $Q$ consists of a double loop and some (at least one) double edges
on $x$. If $Q$ is a component of $G$ we have a star or a dumbbell as a bad
component. Otherwise there exists $z\in Z$ different from all vertices of the
star, such that $xz\in E(G)$, w.l.o.g. with label $\\{1\\}$.

As in the previous case, we consider $G^{\prime}=G(1,2,x)$, now vertex $3$
plays the role of $x$. With the argument of the previous paragraph we get that
a star on $3$ is the bad block $Q$ of $G^{\prime}$. If $Q$ is a component of
$G^{\prime}$ we have the star or a dumbbell as a bad component. Otherwise
there exists $z^{\prime}\in Z$ different from all vertices of the star such
that $3z^{\prime}\in E(G^{\prime})$, implying $3z^{\prime}w\in E(H_{n})$ where
$w\in\\{x,1,2\\}$. However, for all choices of $w$ we get a contradiction
exactly as in (1).

  * •

4\. Assume that $H_{n}$ does not contain $K_{4}^{3}-e$. Then a bad component
on $x$ (or on $y$) is a dumbbell on $xy$.

Suppose there is a double edge $zz^{\prime}$ in some $Z_{ij}$, say in $Z_{12}$
which is incident to a loop, say $zz$. We can define
$B=\\{12z,1zz^{\prime},2zz^{\prime}\\}$ and $G^{\prime\prime}=G(\\{1,2,z\\})$.
Now $z^{\prime}$ can be in the role of $x$, it has a double loop, thus there
is a $K_{4}^{3}-e$ in $H_{n}$, contradiction. The same argument eliminates a
double edge from $x$ or from $y$ to a loop on $z$ for $z\in Z$. Likewise, a
double edge $xy$ with label $\\{1,2\\}$ or label $\\{2,3\\}$ can be eliminated
this way.

It is impossible to have a double edge from $x$ or from $y$ to a double edge
in $Z$. Indeed, assume that there are double edges $xz,zz^{\prime}$ with
$z,z^{\prime}\in Z$. Let the label of $xz$ be w.l.o.g. $\\{1,i\\}$, the label
of $zz^{\prime}$ must be also $\\{1,i\\}$. Then $\\{1,i,x,z\\}$ spans a member
of $\mathcal{B}$ with assignment

$1i\mapsto 123,1x\mapsto 12x,ix\mapsto ixz,1z\mapsto 1zz^{\prime},iz\mapsto
izz^{\prime},xz\mapsto 1xz,$

giving contradiction.

By a similar argument, it is also impossible that $x,y$ both send a double
edge to the same vertex $z\in Z$:

$1i\mapsto 123,1x\mapsto 12x,ix\mapsto iyz,1z\mapsto 1yz,iz\mapsto
iyz,xz\mapsto 1xz.$

We conclude that either $x,y$ belong to distinct blocks (both a single loop
with some double edges), or $x,y$ are joined with a double edge with label
$\\{1,3\\}$. In the former case the blocks are good blocks, otherwise $x,y$
define a (bad) dumbbell block $Q$. We claim that in this case $Q$ is a
component.

Indeed, if there is an edge from $\\{x,y\\}$ to $Z$, say $xz$ with $z\in Z$
and with label $\\{1\\}$ or label $\\{2\\}$ then $\\{1,2,x,y\\}$ spans a
member of $\mathcal{B}$; and if the label is $\\{3\\}$ then $\\{1,3,x,y\\}$
spans a member of $\mathcal{B}$:

$1x\mapsto 1xz,2x\mapsto 12x,1y\mapsto 1xy,2y\mapsto 23y,12\mapsto
123,xy\mapsto 3xy;$ $2x\mapsto 2xz,1x\mapsto 12x,1y\mapsto 1xy,2y\mapsto
23y,12\mapsto 123,xy\mapsto 3xy;$ $3x\mapsto 3xz,1x\mapsto 12x,3y\mapsto
23y,1y\mapsto 1xy,13\mapsto 123,xy\mapsto 3xy,$

proving the claim.

This finishes the proof of Lemma 8. ∎

Assume that the connected components of $G$ are $C_{1},C_{2},\dots$ and let
$I,J$ denote the index sets of the good and bad components, respectively. Let
$U\subset V(G)$ be the set of vertices in $G$ uncovered by the bad components.
Then we have

$s(G)=\sum_{j\in J}s(G[C_{j}])+\sum_{i\in I}s(G[C_{i}])\leq\sum_{j\in
J}s(G[C_{j}])+|U|.$ (2)

Now we estimate the number of edges in $G^{*}$. Using Proposition 8, assume
that we have the following bad components: $p$ double loops, $q$ dumbbells and
an $m$-star or a triple loop on $x$. The bad components cover $m+p+2q$
vertices where the $m=1$ case occurs if $x$ is covered by a triple loop,
$m\geq 3$ if $x$ is covered by an $m$-star (the $m=2$ case is considered as a
dumbbell). The $m=0$ case is when no bad component covers $x$. Then
$|U|=n-3-(m+p+2q)$ and $s(G)=m^{\prime}+2p+3q$ where $m^{\prime}=3$ if $m=1$,
otherwise $m^{\prime}=m+1$. Therefore, by (2)

$s(G)\leq n-3-(m+p+2q)+(m^{\prime}+2p+3q)=n-3+(\rho+p+q)$

where $\rho=2$ if $m=1$ otherwise $\rho=1$.

Using $\delta(H_{n})\geq f(n)-f(n-1)+1$ for the vertices $\\{1,2,3\\}$ (and
subtracting one because the edge $\\{1,2,3\\}$ does not contribute to edges of
$G$), we have

$|E(G^{*})|=|E(G)|-s(G)\geq 3\delta(H_{n})-s(G)\geq
3(f(n)-f(n-1))-(n-3+(\rho+p+q)).$

The union of the bad blocks contain $m-1+q$ edges of $G^{*}$ ($m-1$ in the
$m$-star, $q$ in the dumbbells) thus

$|E(G^{*}[U])|\geq 3(f(n)-f(n-1))-(n-3+(m-1+\rho+p+2q))=M.$ (3)

Since $0\leq m+p+2q\leq n-3$, we can write $m+p+2q=\alpha(n-3)$ with
$0\leq\alpha\leq 1$. Then

$|V(G^{*}[U])|=(1-\alpha)(n-3),M=3(f(n)-f(n-1))-((n-3)(\alpha+1)+\rho-1).$

We use Turán’s theorem to finish the proof of Theorem 1, by showing that
$G^{*}[U]$ contains a $K_{4}$, thus $H_{n}$ contains a member of
$\mathcal{B}$, leading to contradiction. We need to show that
$M\geq{|V(G^{*}[U])|^{2}\over 3}$, i.e. (using that $\rho-1\leq 1$)

######  Proposition 9.

$9(f(n)-f(n-1))-3(1+\alpha)(n-3)-(1-\alpha)^{2}(n-3)^{2}-3\geq 0$.

###### Proof.

Using Observation 3, we have three very similar cases.

  * •

1\. $n=3k$. Then $9(f(n)-f(n-1))=9(k^{2}+2)=(3k-3)^{2}+18k+9$ thus

$18k+9+(3k-3)^{2}-3(1+\alpha)(n-3)-(1-\alpha)^{2}(n-3)^{2}-3=$
$=18k+9-3(3k-3)(1+\alpha)+(3k-3)^{2}(2\alpha-\alpha^{2})-3\geq$
$=\alpha((3k-3)^{2}-3(3k-3))+9k+15\geq 0$

if $k\geq 2$ (for the last inequality we used $\alpha\geq\alpha^{2}$).

  * •

2\. $n=3k+1$. Then $9(f(n)-f(n-1))=9k^{2}=(3k-2)^{2}+12k-4$ thus

$12k-4+(3k-2)^{2}-3(1+\alpha)(n-3)-(1-\alpha)^{2}(n-3)^{2}-3=$
$=12k-4-3(3k-2)(1+\alpha)+(3k-2)^{2}(2\alpha-\alpha^{2})-3\geq$
$=\alpha((3k-2)^{2}-3(3k-2))+3k-1\geq 0$

if $k\geq 2$ (using again $\alpha\geq\alpha^{2}$).

  * •

3\. $n=3k+2$. Then $9(f(n)-f(n-1))=9(k^{2}+k)=(3k-1)^{2}+15k-1$ thus

$15k-1+(3k-1)^{2}-3(1+\alpha)(n-3)-(1-\alpha)^{2}(n-3)^{2}-=$
$=15k-1-3(3k-1)(1+\alpha)+(3k-1)^{2}(2\alpha-\alpha^{2})-3\geq$
$=\alpha((3k-1)^{2}-3(3k-1))+6k-1\geq 0$

if $k\geq 2$ (using again $\alpha\geq\alpha^{2}$).

This proves Proposition 9 and finishes the proof of Theorem 1. ∎

∎

Proof of Corollary 2. One can easily check that the configurations of
$\mathcal{B}$ whose absence was used in the proof of Theorem 1 had at most
nine vertices. In fact, the only place where a nine-vertex configuration of
$\mathcal{B}$ could appear was in the proof of Proposition 5, where the vertex
set was $\\{1,2,3,x,y,z\\}$ and the possible three vertices of
$e_{i}\setminus\\{z\\}$ for $i=1,2,3$. In all other places we referred to
configurations of at most seven vertices. ∎

##  3 Launching the induction

It is possible that the next theorem can be proved by a computer program, but
we prove it by traditional ways. The first statement, ${\rm
ex}_{3}(5,\mathcal{B})=5$, was proved originally in [2] as a lemma to show
that the $2$-color Ramsey number of $\mathcal{B}(K_{4})$ in triple systems is
equal to six.

######  Theorem 10.

${\rm ex}_{3}(5,\mathcal{B})=5$ and ${\rm ex}_{3}(6,\mathcal{B})=8$.

###### Proof.

A pair of vertices in a triple system $H$ is uncovered if no triple of $H$
contains the pair. Let $W$ be the graph formed by the uncovered pairs.

I. $n=5$. Five triples clearly cannot form a member of $\mathcal{B}$ thus we
have to show ${\rm ex}_{3}(5,\mathcal{B})<6$. Assume $H_{5}$ is a triple
system with six triples on vertex set $[5]$ without any member of
$\mathcal{B}$. Observe that the maximum degree of $H_{5}$ is at least
$\lceil{6\times 3\over 5}\rceil=4$.

  * •

1\. Suppose that $W$ has an edge, for some $1\leq i<j\leq 5$, the pair $ij$ is
not covered by any triple of $H$. By symmetry, let $i=1,j=2$. Then $H_{5}$
either contains the six triples meeting $\\{1,2\\}$ in one vertex or one of
them, say $234$ is missing. In the first case the assignment

$e_{13}\mapsto 134,e_{14}\mapsto 145,e_{15}\mapsto 135,e_{34}\mapsto
234,e_{35}\mapsto 235,e_{45}\mapsto 245$

defines a member of $\mathcal{B}$, otherwise the assignment $e_{34}\mapsto
234$ is replaced by the assignment $e_{34}\mapsto 345$ to get a member of
$\mathcal{B}$. In both cases we have a contradiction.

  * •

2\. Every pair of $[5]$ is covered by some triple of $H_{5}$. Let $1$ be a
vertex of maximum degree. If its degree is at least five, then the trace
$G=G(\\{1\\})$ contains a $4$-cycle $23452$. Then $\\{2,3,4,5\\}$ spans a
member of $\mathcal{B}$ with the two distinct triples covering pairs $24,35$
and the four triples on $1$, contradiction. Thus the degree of $1$ is four,
and we may assume that $G$ does not contain a $4$-cycle. Thus $G$ is a
triangle, say $2342$ with a pendant edge $45$. There are two triples of
$H_{5}$ within $\\{2,3,4,5\\}$ and apart from symmetry there are four possible
choices (with respect to the four edges of $G$). In each case we can easily
find a member of $\mathcal{B}$ spanned by $\\{1,2,3,4\\}$ as shown below (the
first two assignments are the ones using the two edges of $H_{5}$ inside
$\\{2,3,4,5\\}$.

$34\mapsto 234,23\mapsto 235,12\mapsto 123,13\mapsto 134,14\mapsto
145,24\mapsto 124;$ $24\mapsto 234,34\mapsto 345,12\mapsto 124,13\mapsto
134,14\mapsto 145,23\mapsto 123;$ $23\mapsto 235,24\mapsto 245,12\mapsto
124,13\mapsto 123,14\mapsto 145,34\mapsto 134;$ $24\mapsto 245,34\mapsto
345,12\mapsto 124,13\mapsto 134,14\mapsto 145,23\mapsto 123.$

II. $n=6$. Suppose we have $9=f(6)+1$ triples in $H_{6}$. The maximum degree
is at least $\lceil{9\times 3\over 6}\rceil=5$. The minimum degree is at least
$4$ otherwise deleting the vertex we get six triples on five vertices
contradicting Case I.

This implies that $W$ has no vertex with degree at least two, thus its edges
form a matching $M$. If $|M|=3$ then some vertex on $M$ has degree at least
$5$, contradiction.

  * •

1\. Suppose $|M|=2$, say $12,34$ are the edges of $W$ and some vertex on $M$,
say $1$, has degree $5$. Then the five triples containing $1$ are determined:
$e_{1}=135,e_{2}=136,e_{3}=145,e_{4}=146,e_{5}=156$. Then, since the degree of
$2$ is at least four, one of the triples $235,236$ must be an edge of $H_{6}$.
Then $\\{1,3,5,6\\}$ spans a member of $\mathcal{B}$ in $H_{6}$, shown by one
of the assignments

$35\mapsto 235,13\mapsto 135,15\mapsto 145,16\mapsto 146,36\mapsto
136,56\mapsto 156,$ $36\mapsto 236,13\mapsto 136,15\mapsto 145,16\mapsto
146,35\mapsto 135,56\mapsto 156.$

  * •

2\. Suppose $|M|=2$, say $12,34$ are the edges of $W$ and vertices $1,2,3,4$
have degree four. Then $d(5)+d(6)=11$ implying that the trace $G(\\{5,6\\})$
of $H_{6}$ must contain a $P_{4}$, say $1,3,2,4$ with double edges, and (by
Property 4. in Proposition 6) $\\{2,3,5,6\\}$ spans a member of $\mathcal{B}$.

  * •

3\. Suppose $|M|=1$, $12$ is an edge of $W$. Now the trace $G=G(\\{1,2\\})$
has no loops. Note that $G$ cannot contain a $4$-cycle, say $3,4,5,6,3$
because then with the edges covering the pairs $35,46$ we get a member of
$\mathcal{B}$ on $\\{3,4,5,6\\}$. Since $G$ contains at least eight edges, the
trace must be a triangle with a pendant edge formed by double edges, say
$34,45,35,36$. Now $\\{1,3,4,5\\}$ spans a member of $\mathcal{B}$,
contradiction.

  * •

4\. Suppose $|M|=0$ i.e. $W$ has no edges, every pair of $[6]$ is covered by
some edge of $H_{6}$. No vertex of $H_{6}$ has degree at least 6. Indeed, if
$1$ is such a vertex, then the trace $G=G(\\{1\\})$ has at least six edges and
either contains a $4$-cycle, say $23452$ or isomorphic to two edge disjoint
triangles, say with edges $23,34,24,25,26,56$. In the first case
$\\{2,3,4,5\\}$ spans a member of $\mathcal{B}$ with assignments

$23\mapsto 123,34\mapsto 134,45\mapsto 145,25\mapsto 125$

and the pairs $24,35$ are mapped to the triples containing them. In the second
case one can easily see that if any triple in $\\{3,4,5,6\\}$ is not in
$H_{6}$ then we have a member of $\mathcal{B}$. For example, assume that $345$
is not in $H_{6}$, thus the edges $e_{35},e_{45}$ that cover the pairs $35,45$
are different and $\\{2,3,4,5\\}$ spans a member of $\mathcal{B}$ with the
assignments

$35\mapsto e_{35},45\mapsto e_{45},23\mapsto 123,24\mapsto 124,25\mapsto
125,34\mapsto 134,$

giving a contradiction.

Thus the degree sequence of $H_{6}$ is $4,4,4,5,5,5$. Counting the codegrees,
i.e. the number of times the vertex pairs of $H_{6}$ are covered by the
triples, we get that at least three pairs $x_{i}y_{i}$ have codegree one.
Indeed, otherwise the sum of the codegrees is at least $2+13\times 2>3\times
9$, contradiction.

We claim that for $i=1,2,3$, the degrees of $x_{i},y_{i}$ are equal to $4$.
Indeed, otherwise the trace $G=G(\\{x_{i},y_{i}\\})$ has at least seven edges,
implying that it must contain a pair of incident double edges plus a further
loop or edge incident to an endpoint, leading to a member of $\mathcal{B}$. In
fact, this is Property 4. in Proposition 6.

We are left with one case to consider: the pairs $x_{i}y_{i}$ form a triangle,
say $1,2,3,1$ their degrees are equal to $4$ and vertices $4,5,6$ has degree
$5$. Then $\\{1,2,3\\}$ cannot be in $E(H_{6})$, that would allow at most
eight triples. The trace $G=G(\\{1,2,3\\})$ on $\\{4,5,6\\}$ cannot contain a
triple edge, say $45$ with label $\\{1,2,3\\}$ because then $\\{1,2,3,4\\}$
would span a member of $\mathcal{B}$. Indeed,

$14\mapsto 145,24\mapsto 245,34\mapsto 345$

and the pairs $12,13,23$ are covered by the triples of $H_{6}$ containing
them. This implies that the double edges $45,56,46$ are labeled w.l.o.g. with
$\\{1,2\\},\\{1,3\\},\\{2,3\\}$, respectively, and each vertex of $T$ has one
loop. Now $\\{1,2,3,4\\}$ spans a member of $\mathcal{B}$ again,

$14\mapsto 145,24\mapsto 245,34\mapsto 346$

and the pairs $12,13,23$ are covered by the edges of $H_{6}$ containing them.
This is a contradiction, finishing the proof of Theorem 10.

∎

Acknowledgement. Thanks to Maria Axenovich, Dániel Gerbner and Zoltán Füredi
for conversations on the subject.

## References

  * [1] N. Alon, C. Shikhelman, Many $T$ copies in $H$-free graphs, Journal of Combinatorial Theory B 121 (2016) 146-172.
  * [2] M. Axenovich, A. Gyárfás, A note on Ramsey numbers of Berge-$G$ hypergraphs, in preparation.
  * [3] P.Frankl, Z. Füredi, G. Simonyi, Triangle-free triple systems, manuscript, 2003
  * [4] Z. Füredi, L. Özkahya, On $3$-uniform hypergraphs without a cycle of given length, arXiv:1412.8083v2 (2014)
  * [5] D. Gerbner, A. Methuku, C. Palmer, A general lemma for Berge-Turán hypergraph problems, in preparation.
  * [6] D. Gerbner, A. Methuku, M. Vizer, Asymptotics for the Turán number of Berge-$K_{2,t}$, arXiv:1705.04134v2 (2017)
  * [7] D. Gerbner, C. Palmer, Extremal results for Berge-hypergraphs, SIAM Journal on Discrete Mathematics 31 (2017) 2314-2327.
  * [8] D. Gerbner, C. Palmer, Counting copies of a fixed subgraph in $F$-free graphs, arXiv:1805.07520v2 (2018)
  * [9] D. Grósz, A. Methuku, C. Tompkins, Uniformity thresholds for the asymptotic size of extremal Berge-$F$-free hypergraphs, submitted
  * [10] E. Győri, Triangle-free hypergraphs. Combinatorics, Probability and Computing 15 (2006), 185–191.
  * [11] E. Győri, G. Y. Katona, N. Lemons, Hypergraph extensions of the Erdős - Gallai theorem, Electronic Notes in Discrete Mathematics 36 (2010) 665-662.
  * [12] E. Győri, N. Lemons, $3$-uniform hypergraphs avoiding a given odd cycle, Combinatorica 32 (2012), 187-203.
  * [13] L. Maherani, M. Shahsiah, Turán numbers of complete $3$-uniform Berge-hypergraphs, arXiv:1612.08856v1 (2016)
  * [14] D. Mubayi, A hypergraph extension of Turán’s theorem, Journal of Combinatorial Theory B. 96 (2006) 122-134.
  * [15] D. Mubayi, J. Verstraëte, A survey of Turán problems for expansions, arXiv:1505.08078v1 (2015)
  * [16] C. Palmer, M. Tait, C.Timmons, A.Z Wagner, Turán numbers for Berge-hypergraphs and related extremal problems, arXiv:1706.04249v1 (2017)
  * [17] O. Pikhurko, Exact computation of the hypergraph Turán function for expanded complete $2$-graphs, Journal of Combinatorial Theory B. 103 (2013) 220-225.