# On variable selection in joint modeling of mean and dispersion

E. R.Edmilson Rodrigues<EMAIL_ADDRESS>[ L. A.Leandro
Alves Pereiralabel=e2 [<EMAIL_ADDRESS>Federal University of
Uberlândia, Uberlândia, Brazil,

###### Abstract

The joint modeling of mean and dispersion (JMMD) provides an efficient method
to obtain useful models for the mean and dispersion, especially in problems of
robust design experiments. However, in the literature on JMMD there are few
works dedicated to variable selection and this theme is still a challenge. In
this article, we propose a procedure for selecting variables in JMMD, based on
hypothesis testing and the quality of the model’s fit. A criterion for
checking the goodness of fit is used, in each iteration of the selection
process, as a filter for choosing the terms that will be evaluated by a
hypothesis test. Three types of criteria were considered for checking the
quality of the model fit in our variable selection procedure. The criteria
used were: the extended Akaike information criterion, the corrected Akaike
information criterion and a specific criterion for the JMMD, proposed by us, a
type of extended adjusted coefficient of determination. Simulation studies
were carried out to verify the efficiency of our variable selection procedure.
In all situations considered, the proposed procedure proved to be effective
and quite satisfactory. The variable selection process was applied to a real
example from an industrial experiment.

Model selection,

location-dispersion analysis,

coefficient of determination,

robust design experiment,

generalized linear models,

extended quasi-likelihood,

###### keywords:

and

t1The authors thank Fapemig for financial support.

##  1 Introduction

The variable selection process that will be considered in this article
concerns the joint modeling of the mean and dispersion (JMMD), which is an
extension of the generalized linear models (GLMs), introduced by (1972).

GLMs assume an exponential family distribution for the response variable and
are more general than normal linear methods in that a mean-variance
relationship appropriate for the data can be accommodated and an appropriate
scale can be chosen for modeling the mean on which the action of the
covariates is approximately linear. In generalized linear models the focus is
on modeling and estimating the mean structure of the data while treating the
dispersion parameter as a constant since GLMs automatically allow for
dependence of the variance on the mean through the distributional assumptions.
However, there are situations in which the observed data may exhibit greater
variability than the one which is implied by the mean-variance relationship
and thus the loss of efficiency in estimating the mean parameters, using
constant dispersion models when the dispersion is varying, may be substantial
(, 2016).

The joint modeling of the mean and dispersion is a general method that
overcomes this difficulty, allowing the modeling of dispersion as a function
of covariates, by the use of two interlinked generalized linear models. JMMD
has enormous potential for application in science and industry and therefore
has attracted the attention of many researchers. In the following, we present
some relevant works in this area. (1987) considered a log-linear regression
model for the variance in Gaussian models; (1987) introduced the extended
quasi-likelihood (EQL); (1987) used the pseudo-likelihood method (PL), see
(1988), p. 71; (1988) showed the asymptotic equivalence between EQL and PL and
that EQL can be affected by an asymptotic bias when the underlying
distribution is asymmetric and outside the exponential family. (1989) proposed
the class of double generalized linear models (DGLMs) for modeling of mean and
dispersion; further developments for EQL were provided by (1989) and (1989);
(1991) introduced the JMMD, using EQL, for robust parameter design, see
(1986); (1992) compared EQL and PL in simulation studies with finite samples
and showed that the maximum extended quasi-likelihood estimator is usually
superior than the pseudo-likelihood estimator in minimizing the mean-squared
error, see also (2000); (1993) used restricted maximum likelihood (REML) for
normal heteroscedastic models; (1996) introduced the class of hierarchical
generalized linear models (HGLMs) as a synthesis of JMMD, generalizing the
ideas of quasi-likelihood and EQL, see (2001) and (2006); (1998) introduced
the adjustment for the REML in the JMMD proposed by (1991); (1999a) and
(1999b) considered extensions of REML estimation of variance parameters for
DGLMs; (2000) showed that the double exponential family, see (1986), and the
EQL generate identical inferences; (2001) presented a Bayesian version of JMMD
for normal heteroscedastic linear models; (2016) considered estimation and
variable selection for JMMD in proper dispersion models, see (1987) and
(1997); (2016) proposed the class of multivariate covariance generalized
linear models (MCGLMs).

In this article, the approach considered for the JMMD uses the adjusted EQL,
as proposed by (1998). Although the JMMD can be applied in several areas, one
of its main applications and, therefore, the focus of this article, is in
robust design experiments (RDE), a very useful tool for improving quality in
industrial experiments, see (1989) and (1992). According to (1992), for data
sets in which the dispersion must be modeled with reasonably small amounts of
data, as in RDE, the best optimization criterion for model fitting in the JMMD
is the EQL.

In RDE, modeling and controlling the variance is of paramount importance,
however this task is hampered by the fact that the experimental data are, in
general, obtained from unreplicated fractional factorial experiments with few
experimental runs. The JMMD allows the dispersion to be modeled very
effectively in RDE problems, even without replication and the variance model
can be obtained from the dispersion model in a very simple way. Given the
importance of modeling the mean and the variability in industrial processes,
if we are considering the JMMD, the question that arises is how to find the
best model for mean and the best model for dispersion.

The methods used to select variables in the JMMD applied to RDE are not very
clear and are based mainly on graphical analysis, see, for example, (1996). On
the other hand, to carry out a formal process of selection of variables, we
would have to consider a large number of possible models, which makes the
process quite laborious. Thus, a clear, reliable and efficient method of
selecting variables is of fundamental importance.

The aim of this article is to propose a consistent methodology that allows the
selection of variables for both the mean and dispersion models in the JMMD. We
propose a procedure for selecting variables, based on hypothesis testing and
the quality of the model’s fit. A criterion for checking the goodness of fit
is used, in each iteration of the selection process, as a filter for choosing
the terms that will be evaluated by a hypothesis test. Three types of criteria
were considered for checking the quality of the model fit in our variable
selection procedure. The criteria used were: the extended Akaike information
criterion, the corrected Akaike information criterion and a specific criterion
for the JMMD, proposed by us, a type of extended adjusted coefficient of
determination.

The selection of variables for joint models of the mean and dispersion was
considered by (2009), where was proposed an extension of the Akaike
information criterion based on EQL; (2012), for the normal inverse
distribution; (2012), for lognormal distribution; (2014), for multilevel
models; (2015), for HGLMs; (2016), for proper dispersion models; (2017), for
beta regression models; Wu et al. (2020), for mixture of regression models.

The article is organized as follows. In Section 2, the joint modeling of mean
and dispersion is briefly described. In Section 3, the proposed procedure and
criterion are presented. In Section 4, simulation studies are carried out to
evaluate the proposed procedure and criterion. An application to real data
from an industrial experiment is presented in Section 5. Finally, a conclusion
is given in Section 6.

##  2 Joint modeling of mean and dispersion

According to (1991), the method of joint modeling of mean and dispersion
consists of finding joint models for the mean and dispersion. In their
approach, using the extended quasi-likelihood, two interlinked generalized
linear models are needed, one for the mean and the other for the dispersion.

Let $Y_{1},\ldots,Y_{n}$ be $n$ independent random variables with the same
probability distribution, representing the dependent response variables, whose
observed values are given by $y_{1},\ldots,y_{n}$. For the $i$th response
$Y_{i}$ it is assumed to be known that $E(Y_{i})=\mu_{i}$ and
$Var(Y_{i})=\phi_{i}V(\mu_{i})$, where $\phi_{i}$ is the dispersion parameter
and $V(.)$ is the variance function in GLMs. The mean and dispersion models
are built as follows.

Suppose $\mbox{\boldmath{$x$}}^{t}=(x_{1},\ldots,x_{s})$ and
$\mbox{\boldmath{$z$}}^{t}=(z_{1},\ldots,z_{r})$ are the vectors of the
independent variables that affect the mean and dispersion models,
respectively. The vectors $x$ and $z$ may contain, in addition to design
factors that are controlled by the experimenter, noise factors that cannot be
controlled and are considered as random variables.

Let $\varphi$ be a link function for the mean model, i.e., for the $i$th
response, the linear predictor for the mean model is given by
$\eta_{i}=\varphi(\mu_{i})=\mbox{\boldmath{$f$}}^{t}(\mbox{\boldmath{$x$}}_{i})\mbox{\boldmath{$\beta$}}$
with
$\mbox{\boldmath{$f$}}^{t}(\mbox{\boldmath{$x$}}_{i})=(f_{1}(\mbox{\boldmath{$x$}}_{i}),\ldots,f_{p}(\mbox{\boldmath{$x$}}_{i}))$
where $f_{j}(\mbox{\boldmath{$x$}}_{i})$, for $j=1,\ldots,p$, is a known
function of $\mbox{\boldmath{$x$}}_{i}$ and $\beta$ is a $p\times 1$ vector of
unknown parameters.

For the dispersion model is used as response variable the deviance component,
which, for each observation $y_{i}$, is given by

$d_{i}=2\int_{\mu_{i}}^{y_{i}}\frac{(y_{i}-t)}{V(t)}dt,$ (2.1)

see (1989), p. 360.

Following (1998), for the dispersion model we are assuming a Gamma model with
a log link function, i.e., for the $i$th response, the linear predictor for
the dispersion model is given by
$\xi_{i}=\log(\phi_{i})=\mbox{\boldmath{$g$}}^{t}(\mbox{\boldmath{$z$}}_{i})\mbox{\boldmath{$\gamma$}}$,
with
$\mbox{\boldmath{$g$}}^{t}(\mbox{\boldmath{$z$}}_{i})=(g_{1}(\mbox{\boldmath{$z$}}_{i}),\ldots,g_{q}(\mbox{\boldmath{$z$}}_{i}))$,
where $g_{j}(\mbox{\boldmath{$z$}}_{i})$, for $j=1,\ldots,q$, is a known
function of $\mbox{\boldmath{$z$}}_{i}$ and $\gamma$ is a $q\times 1$ vector
of unknown parameters. We also define
$\mbox{\boldmath{$X$}}=[\mbox{\boldmath{$f$}}(\mbox{\boldmath{$x$}}_{1}),\ldots,\mbox{\boldmath{$f$}}(\mbox{\boldmath{$x$}}_{n})]^{t}$
the $n\times p$ experimental matrix for the mean model and
$\mbox{\boldmath{$Z$}}=[\mbox{\boldmath{$g$}}(\mbox{\boldmath{$z$}}_{1}),\ldots,\mbox{\boldmath{$g$}}(\mbox{\boldmath{$z$}}_{n})]^{t}$
the $n\times q$ experimental matrix for the dispersion model.

The fitting for the JMMD uses as an optimizing criterion the extended quasi-
likelihood, introduced by (1987), see (1989), p. 349. In this work we use the
adjusted extended quasi-likelihood, introduced by (1998) and given by

$Q^{+}(\mbox{\boldmath{$\mu$}},\mbox{\boldmath{$\phi$}};\mbox{\boldmath{$y$}})=\sum_{i=1}^{n}-\frac{1}{2}\left(\frac{d_{i}^{*}}{\phi_{i}}+\log\\{2\pi\phi_{i}V(y_{i})\\}\right)$
(2.2)

where $d_{i}^{*}=\frac{d_{i}}{1-h_{i}}$ is the standardized deviance component
and $h_{i}$ is the $i$th element of the diagonal of the matrix
$\mbox{\boldmath{$H$}}=\mbox{\boldmath{$W$}}^{\frac{1}{2}}\mbox{\boldmath{$X$}}(\mbox{\boldmath{$X$}}^{t}\mbox{\boldmath{$WX$}})^{-1}\mbox{\boldmath{$X$}}^{t}\mbox{\boldmath{$W$}}^{\frac{1}{2}}$,
being $W$, the $n\times n$ weight matrix for the GLMs, a diagonal matrix with
elements given by
$w_{i}=\left(\frac{\partial\mu_{i}}{\partial\eta_{i}}\right)^{2}\frac{1}{\phi_{i}V(\mu_{i})}$
with $i=1,\ldots,n$. The Table 1 gives a resume of the joint modeling of mean
and dispersion. From Table 1, we can observe that the standardized deviance
component from the model for the mean becomes the response for the dispersion
model, and the inverse of fitted values for the dispersion model provides the
prior weights for the mean model.

Table 1: Summary of the JMMD for the ith response Component |  Mean model |  Dispersion model†
---|---|---
Response variable |  $y_{i}$ |  $d_{i}^{*}$
Mean |  $\mu_{i}$ |  $\phi_{i}$
Variance |  $\phi_{i}V(\mu_{i})$ |  $2\phi_{i}^{2}$
Link function |  $\eta_{i}=\varphi(\mu_{i})$ |  $\xi_{i}=\log(\phi_{i})$
Linear predictor |  $\eta_{i}=\mbox{\boldmath{$f$}}^{t}(\mbox{\boldmath{$x$}}_{i})\mbox{\boldmath{$\beta$}}$ |  $\xi_{i}=\mbox{\boldmath{$g$}}^{t}(\mbox{\boldmath{$z$}}_{i})\mbox{\boldmath{$\gamma$}}$
Deviance component |  $d_{i}=2\int_{\mu_{i}}^{y_{i}}\frac{y_{i}-t}{V(t)}dt$ |  $2\left\\{-\log\left(\frac{d_{i}^{*}}{\phi_{i}}\right)+\frac{(d_{i}^{*}-\phi_{i})}{\phi_{i}}\right\\}$
Prior weight |  $\frac{1}{\phi_{i}}$ |  $(1-h_{i})/2$
†For the dispersion model we are assuming a Gamma model with logarithmic link
function

The algorithm for parameter estimation is an extension of the standard GLMs
algorithm, in which the model for the mean is fitted assuming that the fitted
values for the dispersion are known and that the model for dispersion is
fitted using the fitted values for the mean. The fitting alternates between
the mean and dispersion models until convergence is achieved.

In the next section, we present a procedure for variable selection in JMMD,
considering two different criteria to verify the goodness of the model fit.

##  3 Procedure for variable selection in JMMD

For our procedure of variable selection in the joint modeling of the mean and
dispersion we consider two alternative criteria to verify the goodness of the
model fit. The first criterion, proposed by us, is a combination of the
criteria proposed by (2008) and Zhang (2017) for selection of variables in
generalized linear models. (2008) propose a criterion of goodness of fit given
by

$\footnotesize{R^{2}_{HS}=1-\frac{\sum_{i=1}^{n}(y_{i}-\widehat{\mu}_{i})^{2}/(n-\lambda_{n}r)}{\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}/(n-1)}},$
(3.1)

where $\overline{y}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$,
$\widehat{\mu}_{i}=\varphi^{-1}(\widehat{\eta}_{i})$, $n$ is the number of
observations, $r$ is the number of estimated parameters of the model and
$\lambda_{n}$ is some function of $n$, satisfying
$\lambda_{n}\rightarrow\infty$ and $\lambda_{n}/n\rightarrow 0$ as the sample
size $n\rightarrow\infty$. The $\lambda_{n}$ functions considered by them in
their simulation studies were $\sqrt{n}$ and $\log n$.

The adjusted criterion proposed by Zhang (2017) to verify the goodness of the
model fit is given by

$\footnotesize{R^{2}_{Z}=1-\frac{\sum_{i=1}^{n}d_{V}(y_{i},\widehat{\mu}_{i})/(n-r)}{\sum_{i=1}^{n}d_{V}(y_{i},\overline{y})/(n-1)}},$
(3.2)

where $d_{V}(a,b)=\left[\int_{a}^{b}\sqrt{1+[V^{\prime}(t)]^{2}}dt\right]^{2}$
with $V^{\prime}(t)=\frac{dV(t)}{dt}$ and $V(t)$ is a continuous and derivable
function in $(a,b)$. The function $d_{V}(a,b)$ is the squared arc length of
the variance function between $V(a)$ and $V(b)$. In equation (3.2),
$\sum_{i=1}^{n}d_{V}(y_{i},\widehat{\mu}_{i})$ represents the variation in the
response unexplained by the model and
$\sum_{i=1}^{n}d_{V}(y_{i},\overline{y})$ represents the total variation in
the response.

According to Zhang (2017), $d_{V}(a,b)$ can differ dramatically from the
Euclidean distance $(a-b)^{2}$ when the underlying variance function is
nonlinear. For example, in a particular but rather general case of quadratic
variance function, i.e., $V(\mu)=c_{2}\mu^{2}+c_{1}\mu+c_{0}$, where $c_{0}$,
$c_{1}$ and $c_{2}$ are known constants and $c_{2}\neq 0$,
$d_{V}(a,b)=\frac{1}{16c_{2}^{2}}\left\\{\log\left(\frac{\nu(b)}{\nu(a)}\right)+\nu(b)-\nu(a)\right\\}^{2}$,
where $\nu(t)=V^{\prime}(t)+\sqrt{1+[V^{\prime}(t)]^{2}}$. If $c_{2}=0$, the
variance function is linear and $d_{V}(a,b)=(1+c_{1}^{2})(a-b)^{2}$. Hence,
for the Normal model, where $V(\mu)=1$, $d_{V}(a,b)=(a-b)^{2}$; for the
Poisson model, where $V(\mu)=\mu$, $d_{V}(a,b)=2(a-b)^{2}$.

The criterion proposed by Zhang (2017) is well defined, as long as the
underlying model, like quasi-models, specifies the mean and variance
functions. In this way, $R_{Z}^{2}$ is applicable to quasi-likelihood models,
because it does not depend on the complete specification of the likelihood
function, requiring only to know the mean and the variance function.

In this way, the extension of the criterion proposed by Zhang (2017) for joint
models of the mean and dispersion is immediate, however, in the mean model,
the dispersion parameter must be taken into account for the calculation of the
criterion. In fact, as in the iterative process for the JMMD the variance for
the ith observation is given by $Var(Y_{i})=\phi_{i}V(\mu_{i})$, where
$\phi_{i}$ is known, given from the previous iteration. We can think of
$\phi_{i}$ as a measure weight for $V(\mu_{i})$.

Therefore, using the idea considered by Zhang (2017) to measure the change of
variation of the response variable, associated with the penalty function
proposed by (2008) and considering the dispersion parameter as weight for the
variance function, we propose, for the process of joint modeling of the mean
and dispersion, the following criteria of goodness of fit. For the mean model
we propose the criterion given by

$\footnotesize{\tilde{R}_{m}^{2}=1-\frac{\sum_{i=1}^{n}\tilde{d}_{V}(y_{i},\widehat{\mu}_{i})/(n-\lambda_{n}p)}{\sum_{i=1}^{n}\tilde{d}_{V}(y_{i},\overline{y})/(n-1)}}$
(3.3)

and for the dispersion model the proposed criterion is

$\footnotesize{\tilde{R}_{d}^{2}=1-\frac{\sum_{i=1}^{n}d_{V}(d_{i}^{*},\widehat{\phi}_{i})/(n-\lambda_{n}q)}{\sum_{i=1}^{n}d_{V}(d_{i}^{*},\overline{d^{*}})/(n-1)}},$
(3.4)

where
$\tilde{d}_{V}(a,b)=\left[\int_{a}^{b}\sqrt{1+\phi^{2}[V^{\prime}(t)]^{2}}dt\right]^{2}$
is calculated for the mean model as proposed by Zhang (2017), but considering
the dispersion parameter as weight for the variance function; $d_{V}(a,b)$,
calculated for the dispersion model, is the same as that given in equation
(3.2); $n$ is the number of observations, $\lambda_{n}$ is a function of $n$,
as proposed by (2008), i.e., $\lambda_{n}\rightarrow\infty$ and
$\lambda_{n}/n\rightarrow 0$ as the sample size $n\rightarrow\infty$, $p$ is
the number of parameters in the mean model, $q$ is the number of parameters in
the dispersion model, $\hat{\mu}_{i}=\varphi^{-1}(\hat{\eta}_{i})$,
$\hat{\phi}_{i}=\exp(\hat{\xi}_{i})$, $\overline{y}=(1/n)\sum_{i=1}^{n}y_{i}$
and $\overline{d^{*}}=(1/n)\sum_{i=1}^{n}d_{i}^{*}$, with $y_{i}$ and
$d_{i}^{*}$ given in Table 1. In equation (3.4) we can notice that if
$\lambda_{n}=1$ we have the criterion proposed by Zhang (2017).

Note that, as in the criterion proposed by Zhang (2017),
$\sum_{i=1}^{n}\tilde{d}_{V}(y_{i},\widehat{\mu}_{i})$ and
$\sum_{i=1}^{n}d_{V}(d_{i}^{*},\widehat{\phi}_{i})$ represent the variation in
the response unexplained by the model and
$\sum_{i=1}^{n}\tilde{d}_{V}(y_{i},\overline{y})$ and
$\sum_{i=1}^{n}d_{V}(d_{i}^{*},\overline{d^{*}})$ represent the total
variation in the response. Thus, for both $\tilde{R}_{m}^{2}$ and
$\tilde{R}_{d}^{2}$ criteria, the higher the value, the better. It is also
worth mentioning that, depending on the models considered, including their
number of parameters, the sample size and the choice of the $\lambda_{n}$
function, the $\tilde{R}_{m}^{2}$ and $\tilde{R}_{d}^{2}$ criteria may have
values outside the range $[0,1]$.

In equation (3.4), since we are considering the Gamma model for the
dispersion,
$d_{V}(a,b)=\frac{1}{16}\left\\{\log\left(\frac{2b+\sqrt{1+4b^{2}}}{2a+\sqrt{1+4a^{2}}}\right)+2b\sqrt{1+4b^{2}}-2a\sqrt{1+4a^{2}}\right\\}^{2}$.
For the mean models that will be considered in this article the values of
$\tilde{d}_{V}(a,b)$, used in equation (3.3), are calculated by:
$\tilde{d}_{V}(a,b)=(b-a)^{2}$ for the Normal type model;
$\tilde{d}_{V}(a,b)=(1+\phi^{2})(b-a)^{2}$ for the Poisson type model and for
the Binomial type model,
$\tilde{d}_{V}(a,b)=\frac{1}{16\phi^{2}}\left\\{\log\left(\frac{\beta+\sqrt{1+\beta^{2}}}{\alpha+\sqrt{1+\alpha^{2}}}\right)+\beta\sqrt{1+\beta^{2}}-\alpha\sqrt{1+\alpha^{2}}\right\\}^{2}$,
where $\alpha=\phi(1-2a)$ and $\beta=\phi(1-2b)$. Note that in the case of the
Normal type model our criterion is the same criterion proposed by (2008).

The second criterion, alternative to the first, considered for the mean model
in our procedure for selection of variables in JMMD, is the extended Akaike
information criterion ($EAIC$), proposed by (2009) and given by

$EAIC=-2Q^{+}(\mbox{\boldmath{$\mu$}},\mbox{\boldmath{$\phi$}};\mbox{\boldmath{$y$}})+F(\kappa,n),$
(3.5)

where $\kappa=p+q$ is the sum of the number of parameters of the mean and
dispersion models and $F(\kappa,n)$ is a penalty function that depends on
$\kappa$ and $n$. In our simulation studies we consider
$F(\kappa,n)=\frac{2\kappa n}{n-\kappa-1}$. Such penalty function, originally
proposed by (1978) as a bias correction to the Akaike information criterion in
the context of linear regression models, has been extended to nonlinear
regression and autoregressive time series models by (1989). For the dispersion
model, as we are considering the known likelihood, the criterion used will be
the corrected Akaike information criterion ($AIC_{c}$), given by
$AIC_{c}=-2\ln\hat{L}+\frac{2qn}{n-q-1}$, where $\hat{L}$ is the maximum value
of the likelihood function for the Gamma dispersion model. For both $EAIC$ and
$AIC_{c}$ criteria, the lower the value, the better.

The Algorithm 1 shows how the terms of the model considered (mean or
dispersion) are selected in order to find the best model.

Algorithm 1: Pseudo algorithm for variable selection in mean (dispersion)
model

In the JMMD, for a given and fixed dispersion (mean) model, choose the terms
of the mean (dispersion) model, considering the following sequence of steps.

  1. 1.

Let $\cal{V}$ be the set of terms that will be used in the variable selection
process and consider $\vartheta=\theta_{0}$ the linear predictor of the
initial model.

  2. 2.

Fit models with linear predictors $\vartheta_{j}=\vartheta+\theta_{j}v_{j}$,
for all $v_{j}\in\cal{V}$.

  3. 3.

For each model fitted in step 2, calculate the appropriate selection measure
($\tilde{R}_{m}^{2}$, $\tilde{R}_{d}^{2}$, $EAIC$, or $AIC_{c}$) to check the
quality of fit and choose the model with linear predictor, say
$\vartheta_{k}=\vartheta+\theta_{k}v_{k}$, that has the best fit. Evaluate the
value of the selection measure found in the current iteration with the value
of the selection measure found in the previous iteration. If the value of the
current selection measure is better than the previous value, go to step 4,
otherwise, go to step 5.

  4. 4.

For models with nested linear predictors $\vartheta$ and $\vartheta_{k}$,
apply an appropriate test to verify the significance of the addition of
$v_{k}$ into the linear predictor $\vartheta$. If $v_{k}$ is significant,
remove the term $v_{k}$ from $\cal{V}$, do $\vartheta=\vartheta_{k}$ and
return to step 2. Otherwise, the procedure ends and $\vartheta$ is the final
model chosen.

  5. 5.

Evaluate, in the same way as in step 4, the significance of the addition of
$v_{k}$ into the linear predictor $\vartheta$. If $v_{k}$ is significant, do
$\vartheta=\vartheta_{k}$. Stop the procedure and $\vartheta$ is the final
model chosen.

In step 4 of Algorithm 1, in the case of the dispersion model, for which we
are assuming a GLM with Gamma distribution, the test for comparing two nested
models is the usual analysis of the deviance for GLMs. In this case, if
$H_{a}$ and $H_{b}$ are two nested hypothesis of dimension $a<b$, that is
$\xi_{a}\subset\xi_{b}$, then, under $H_{a}$, the change in the deviance,
given by
$\frac{1}{2}\\{D(\mbox{\boldmath{$d$}}^{*},\mbox{\boldmath{$\phi$}}_{a})-D(\mbox{\boldmath{$d$}}^{*},\mbox{\boldmath{$\phi$}}_{b})\\}=\frac{1}{2}\left\\{\sum_{i=1}^{n}d_{d_{i}}^{a}-\sum_{i=1}^{n}d_{d_{i}}^{b}\right\\}$,
has an asymptotic $\chi^{2}_{b-a}$, where
$D(\mbox{\boldmath{$d$}}^{*},\mbox{\boldmath{$\phi$}}_{k})\\}$ is the deviance
for the model with mean $\phi_{k}$ and
$d_{d_{i}}^{k}=2\left\\{-\log\left(\frac{d_{i}^{*}}{\phi_{k_{i}}}\right)+\frac{(d_{i}^{*}-\phi_{k_{i}})}{\phi_{k_{i}}}\right\\}$
is the deviance component, as shown in Table 1. In the case of the mean model,
where $\phi$ is given, the test for comparing two nested models is obtained by
the quasi-likelihood ratio test applied to the extended quasi-likelihood. The
justification for using the quasi-likelihood ratio test in the mean model is
given by the fact that the extended quasi-likelihood is a saddle point
approximation for the log likelihood, especially when there exists a
distribution of the exponential family with a given variance function, see
(1991) and (1989), p. 350. Thus, in the same way as we considered previously
in the case of the dispersion model, if $H_{c}$ and $H_{d}$ are two nested
hypothesis of dimension $c<d$, that is $\eta_{c}\subset\eta_{d}$, then, under
$H_{c}$, the change in the extended quasi-deviance, given by
$-2\left\\{Q^{+}(\mbox{\boldmath{$\mu$}}_{c},\mbox{\boldmath{$\phi$}};\mbox{\boldmath{$y$}})-Q^{+}(\mbox{\boldmath{$\mu$}}_{d},\mbox{\boldmath{$\phi$}};\mbox{\boldmath{$y$}})\right\\}=\sum_{i=1}^{n}\frac{1}{\phi_{i}}d^{*}_{\mu_{c_{i}}}-\sum_{i=1}^{n}\frac{1}{\phi_{i}}d^{*}_{\mu_{d_{i}}}$,
has an asymptotic $\chi^{2}_{d-c}$, where $d^{*}_{\mu_{k_{i}}}$ is the
standardized deviance component in relation to $\mu_{k_{i}}$, given in Table
1. However, this result is true when $\phi_{i}$ is the correct value, which is
not the case in our iterative learning process for constructing joint mean and
dispersion models. Therefore, for our procedure, due to uncertainty in the
estimation of $\phi_{i}$, we assume that $\phi_{i}$ is replaced by
$\tau\phi_{i}$, where $\tau$ is an unknown constant. It can be said that
$\tau$ is a compensation for not being able to estimate $\phi_{i}$ exactly.
Consequently, the hypothesis test considered in our mean model selection
procedure, which overcomes the dependency on the unknown parameter $\tau$, is
the F-test and it is not difficult to verify that

$\frac{\left(\sum_{i=1}^{n}\frac{1}{\phi_{i}}d^{*}_{\mu_{c_{i}}}-\sum_{i=1}^{n}\frac{1}{\phi_{i}}d^{*}_{\mu_{d_{i}}}\right)\Big{/}(d-c)}{\left(\sum_{i=1}^{n}\frac{1}{\phi_{i}}d^{*}_{\mu_{d_{i}}}\right)\Big{/}(n-d)}$
(3.6)

has an asymptotic $F_{d-c,n-d}$. For a comprehensive understanding of the
relationship between the F-test and the likelihood ratio test, considering the
saddle point approximation, see (1992), p. 89.

Our selection strategy, shown in Algorithm 2, consists of a three-step scheme
in which the selected dispersion model is used to select the best mean model
and vice versa. In our selection scheme, we use a recursive procedure that is
only finalized when the value used as a measure of goodness of fit for the
mean model reaches an optimal value, for example, when $\tilde{R}_{m}^{2}$
stops increasing or when $EAIC$ stops decreasing.

Algorithm 2: Pseudo algorithm for variable selection in JMMD

  1. 1.

Assuming constant dispersion, use Algorithm 1 to find the terms of the best
current mean model and for this model calculate the value of
$\tilde{R}_{m}^{2}$ ($EAIC$).

  2. 2.

Assuming that the selected current mean model is adequate, use Algorithm 1 to
find the terms of the best current dispersion model.

  3. 3.

Assuming that the selected current dispersion model is adequate, use Algorithm
1 to find the terms of the best current mean model and for this model
calculate a new value of $\tilde{R}_{m}^{2}$ ($EAIC$). If the updated value of
$\tilde{R}_{m}^{2}$ decreases ($EAIC$ increases) or is equal to the value
previously obtained, stop and the models previously found for mean and
dispersion are chosen. Otherwise, return to step 2.

The selection scheme proposed in this article consists of a method that, in
addition to providing a clear and objective process for selecting variables,
aims to reduce computational costs inherent in the great complexity involved
in variable selection for joint mean and dispersion models. In fact, in our
variable selection process, if we assume that $m$ regressors are considered
for each model in Algorithm 1, we have $\sum_{j=0}^{k-1}\binom{m-j}{1}+1$
$=\frac{k}{2}(2m-k+1)+1$ possibilities of joint models to be evaluated, where
$k\leq m$ represents the number of iterations of the algorithm. In our
process, governed by Algorithm 2, Algorithm 1 is used twice, except in the
first iteration, where it is used three times. Therefore, if we consider $l$
iterations of Algorithm 2, we would have to check
$(2l+1)(\frac{k}{2}(2m-k+1)+1)$ joint models. Without using our process, we
would have to evaluate
$(\sum_{j=0}^{m}\binom{m}{j})\times(\sum_{j=0}^{m}\binom{m}{j})$ $=2^{m}\times
2^{m}=2^{2m}$ different joint models. As an illustration of this result,
consider $l=m=k=10$. In our procedure, we should evaluate $21\times 56=1176$
joint models. Without using our procedure, we would have to evaluate
$2^{20}=1048576$ joint models.

##  4 Numerical evaluation

To evaluate the accuracy of the variable selection procedure, Monte Carlo
simulations were performed in order to verify the percentage of hits. The main
objective is to evaluate how the variable selection scheme behaves for
different sample sizes and how the criteria considered affect the proposed
scheme.

In the evaluation process of the selected models, three possible situations
were considered: category 1 (Type 1 model): when the selected model does not
contain all the terms of the simulated model; category 2 (Type 2 model): when
the selected model contains all the terms of the simulated model plus other
terms and correct model (Optimal): when the selected model contains exactly
the terms of the simulated model. All simulations were performed using the
statistical software R (, 2019).

The simulations were performed for mean models in which the response variable
was considered to follow the Normal, Binomial or Poisson distribution. For the
dispersion model, the response variable was considered to have a Gamma
distribution. The general linear predictors, chosen for the data generation
process, were:

$\eta=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{3}$ (4.1)

and

$\xi=\gamma_{0}+\gamma_{1}z_{1}+\gamma_{2}z_{2}+\gamma_{3}z_{3},$ (4.2)

where (4.1) was the linear predictor used for the mean model and (4.2) was the
linear predictor used for the dispersion model. For the simulations, $16$
different scenarios were considered, which originated from the combination of
four possibilities for the goodness of fit criteria with four different sample
sizes. The sample sizes used in the simulations were $25$, $50$, $100$ and
$150$. The criteria considered were the extended Akaike information criterion,
the corrected Akaike information criterion and three variants of our proposed
criterion, one using $\lambda_{n}=1$, another with $\lambda_{n}=\log n$ and
another with $\lambda_{n}=\sqrt{n}$. The extended Akaike information criterion
was used in the mean model and the corrected Akaike information criterion was
used in the dispersion model. For each scenario considered 1000 Monte Carlo
replications were performed. The values of
$\mbox{\boldmath{$x$}}^{t}=(x_{1},x_{2},x_{3})$ and
$\mbox{\boldmath{$z$}}^{t}=(z_{1},z_{2},z_{3})$ were randomly obtained from a
Uniform distribution between $-1$ and $1$.

Regarding the level of significance, considered for the hypothesis tests in
step 4 of Algorithm 1, we understand that a level of 0.05 is very strict and
could exclude important variables from the model. (1997) examined the question
of the level of significance in stepwise procedures for the case of logistic
regression and suggested that the level of significance adopted should be
between 0.15 and 0.20. Based on the discussions presented by (1997) and in the
references considered there, we decided to use an intermediate value between
0.05 and 0.15, that is, a significance level of $0.10$. In all simulated
studies performed, $i=1,\ldots,n$ refers to the $i$th response observed and
$n$ is the sample size. The simulation studies are shown below.

###  4.1 Simulation study for Normal distribution

For the mean model with Normal probability distribution or Normal type, i.e.,
with $V(\mu_{i})=1$, the data simulation process must ensure that
$E(Y_{i})=\mu_{i}$ and $Var(Y_{i})=\phi_{i}$. Such results can be directly
obtained by taking $Y_{i}\sim Normal(\mu_{i},\phi_{i})$ or considering
$S_{i}\sim\phi_{i}\chi^{2}_{1}$, where $\chi^{2}_{1}$ represents the Chi-
square distribution with 1 degree of freedom, and taking
$Y_{i}=\pm\sqrt{S_{i}}+\mu_{i}$.

The values chosen for the parameter vector
$(\beta_{0},\beta_{1},\beta_{2},\beta_{3},$
$\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ used to perform simulations
were $(4,15,13,0,0.3,0,3,0)$, where the zero value means that the parameter
was not used in the simulation. For each Monte Carlo simulation, the responses
were generated considering log link function for the dispersion model and
identity link function for the mean model.

Table 2 presents the percentage of Optimal, Type 1 and Type 2 models, obtained
from 1000 Monte Carlo simulations of the proposed procedure of variable
selection, as previously established at the beginning of Section 4.

Analyzing the percentages for the mean models, referring to Optimal model in
Table 2, we can see that the simulations performed with 25 sample size already
have high success rates, but with considerable indices of Type 1 models. When
we increased the sample size, we can observe the rapid increase in percentages
for the Optimal model and a sharp drop in percentages for the Type 1 model.
For the dispersion models, there is a lower rate of Optimal model, but an
acceptable rate is achieved for samples of size 100 and 150. With respect to
the Type 2 models, we noticed that for both mean and dispersion models the
variable selection scheme presented low percentages and continued with low
percentages even increasing the size of the sample.

Table 2: Percentages of model types found in 1000 Monte Carlo simulations of the proposed variable selection procedure for the Normal distribution |  |  |  | Optimal |  |  |  | Type 2 |  |  |  | Type 1 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | Sample size |  |  |  | Sample size |  |  |  | Sample size |
Model | Criterion |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |
Mean | $\tilde{R}_{m}^{2}(\log n)$ |  |  | 71.0 | 83.5 | 92.0 | 95.1 |  |  |  | 4.9 | 5.1 | 5.5 | 4.8 |  |  |  | 24.1 | 11.4 | 2.5 | 0.1 |
$\tilde{R}_{m}^{2}(\sqrt{n})$ |  |  | 76.6 | 88.5 | 92.0 | 98.9 |  |  |  | 1.0 | 0.2 | 5.5 | 0.1 |  |  |  | 22.4 | 11.3 | 2.5 | 1.0 |
$\tilde{R}_{m}^{2}(1)$ |  |  | 69.2 | 81.9 | 90.4 | 90.4 |  |  |  | 5.8 | 7.0 | 5.9 | 8.4 |  |  |  | 25.0 | 11.1 | 3.7 | 1.2 |
$EAIC$ |  |  | 66.6 | 88.5 | 91.7 | 94.7 |  |  |  | 6.7 | 0.2 | 6.3 | 5.0 |  |  |  | 26.7 | 11.3 | 2.0 | 0.3 |
Dispersion | $\tilde{R}_{d}^{2}(\log n)$ |  |  | 30.9 | 53.4 | 62.8 | 68.8 |  |  |  | 0.8 | 2.1 | 4.5 | 4.0 |  |  |  | 68.3 | 44.5 | 32.7 | 27.2 |
$\tilde{R}_{d}^{2}(\sqrt{n})$ |  |  | 32.0 | 53.8 | 62.9 | 67.5 |  |  |  | 1.2 | 3.5 | 4.4 | 3.4 |  |  |  | 66.8 | 42.7 | 32.7 | 29.1 |
$\tilde{R}_{m}^{2}(1)$ |  |  | 30.9 | 52.1 | 63.8 | 67.9 |  |  |  | 1.8 | 3.6 | 2.7 | 1.5 |  |  |  | 67.3 | 44.3 | 33.5 | 30.6 |
$AIC_{c}$ |  |  | 36.1 | 58.1 | 71.9 | 78.2 |  |  |  | 1.6 | 5.6 | 7.0 | 4.0 |  |  |  | 62.3 | 36.3 | 21.1 | 17.8 |
Figure 1: Graphs, obtained from the simulations of the Normal distribution,
regarding the relationship between sample size and percentage found for the
Optimal and Type 1 models, for the $\tilde{R}_{m}^{2}$, $\tilde{R}_{d}^{2}$,
$EAIC$ and $AIC_{c}$ criteria, where the solid line represents the graph for
$\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$ with $\lambda_{n}=\log n$, the
dashed line represents the graph for $\tilde{R}_{m}^{2}$
$(\tilde{R}_{d}^{2})$, with $\lambda_{n}=\sqrt{n}$, the dotted line represents
the graph for $\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$, with $\lambda_{n}=1$
and the dash-dotted line represents the graph for $EAIC(AIC_{c})$.

The results presented in Table 2 seem to show that the proposed procedure is
efficient to find the correct model from 150 sample size, but also it is
acceptable with smaller sample size. When we comparing the percentages found
for the mean and dispersion models, we found that in both the mean and
dispersion models the scheme of variable selection is efficient, however for
the dispersion model the efficiency is lower than that of the mean model,
mainly for small samples. For the dispersion models, we had maximum percentage
of Optimal model equal to $78.2\%$, considering the $AIC_{c}$; while for the
mean models, the maximum percentage was $98.9\%$, considering the criterion
$\tilde{R}_{m}^{2}$ with $\lambda_{n}=\sqrt{n}$.

When comparing the criteria $\tilde{R}_{m}^{2}$ ($\tilde{R}_{d}^{2}$) and
$EAIC(AIC_{c})$, used as selection measures, we note that, for the mean
models, the criterion $\tilde{R}_{m}^{2}$ with $\lambda_{n}=\sqrt{n}$ was
better and for the dispersion models, the $AIC_{c}$ performed better. These
conclusions were observed for all sample sizes.

The graphs displayed in Figure 1, for the percentages shown in Table 2, show
only the results for the Optimal and Type 1 models and visually confirm part
of the conclusions obtained previously. The solid line show the percentage of
hits of the proposed procedure using the criterion $\tilde{R}_{m}^{2}$
($\tilde{R}_{d}^{2}$) with $\lambda_{n}=\log n$, the dashed line show the
percentage of hits using the criterion $\tilde{R}_{m}^{2}$
($\tilde{R}_{d}^{2}$) with $\lambda_{n}=\sqrt{n}$, the dotted line show the
percentage of hits using the criterion $\tilde{R}_{m}^{2}$
($\tilde{R}_{d}^{2}$) with $\lambda_{n}=1$ and the dash-dotted line presents
the percentage of hits using the $EAIC(AIC_{c})$. The graphs represent well
the efficiency of the proposed procedure, since they show the growth of the
percentages of the Optimal models and the considerable decrease of the
percentages of the Type 1 models. With respect to the comparison of the
criteria, the graphs referring to the mean model have the graph lines very
close, with an advantage for the $\tilde{R}_{m}^{2}$ with
$\lambda_{n}=\sqrt{n}$. In the graphs referring to the dispersion model, the
$AIC_{c}$ is represented by the line farthest from the others, indicating a
best efficiency than $\tilde{R}_{d}^{2}$, regardless of the $\lambda_{n}$
considered.

###  4.2 Simulation study for Binomial distribution

For the mean model with Binomial probability distribution with index $m$ or
Binomial type, i.e., with $V(\mu_{i})=\mu_{i}(1-\mu_{i}/m)$, the data
simulation process must ensure that $E(Y_{i})=\mu_{i}$ and
$Var(Y_{i})=\phi_{i}\mu_{i}(1-\mu_{i}/m)$ and thus, assuring overdispersion.
According to (1989), overdispersion occurs when the variance of the response
exceeds the nominal variance of theoretical model. For the Binomial
distribution the nominal variance is $\mu(1-\mu/m)$. Overdispersion can arise
in a number of ways, but the simplest and most common mechanism is clustering
in the population. The generation of data is adapted from (1989), p. 125.

Clusters usually vary in size, but, in the same way as considered by (1989),
we also assume that the cluster size, $k$, is fixed and that the $m$
individuals sampled come from $\ell=m/k$ clusters. Let’s consider the total
number of positive respondents, represented by the conditional random variable
$Y|\Pi=S_{1}+\ldots+S_{\ell}$, where $S_{j}$, the number of positive
respondents in the cluster $j$, with $j=1,\ldots,\ell$, are independent random
variables with the same distribution as the random variable $S\sim
Binomial(k,\pi)$. We also assume that the random variable $\Pi$ has a Beta
distribution. Note that, for given $\pi_{i}$, $Y_{i}|\pi_{i}\sim
Binomial(m,\pi_{i})$. In this way, for $\Pi_{i}\sim Beta(a_{i},b_{i})$, we can
obtain $E(Y_{i})=m\lambda_{i}$ and
$Var(Y_{i})=m\lambda_{i}(1-\lambda_{i})[1+(m-1)\delta_{i}]$, where
$\lambda_{i}=\frac{a_{i}}{a_{i}+b_{i}}$ and
$\delta_{i}=\frac{1}{a_{i}+b_{i}+1}$. The random variable $Y_{i}$ has the
Beta-Binomial distribution, see (1989), p. 140. Thus, by taking
$\mu_{i}=m\lambda_{i}$ and $\phi_{i}=1+(m-1)\delta_{i}$ we have the desired
result.

For the Binomial distribution, the values chosen for the parameter vector
$(\beta_{0},\beta_{1},\beta_{2},$
$\beta_{3},\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ used to perform
simulations were $(0.2,0.6,0,0.8,$ $1,0,0,2.5)$, where the zero value means
that the parameter was not used in the simulation. We have also considered
$k=5$ and $m=10$. For each Monte Carlo simulation, the responses were
generated considering log link function for the dispersion model and logistic
link function for the mean model.

Table 3: Percentages of model types found in 1000 Monte Carlo simulations of the proposed variable selection procedure for the Binomial distribution |  |  |  | Optimal |  |  |  | Type 2 |  |  |  | Type 1 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | Sample size |  |  |  | Sample size |  |  |  | Sample size |
Model | Criterion |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |
Mean | $\tilde{R}_{m}^{2}(\log n)$ |  |  | 33.9 | 63.1 | 66.1 | 88.8 |  |  |  | 0.4 | 2.1 | 2.0 | 1.4 |  |  |  | 65.7 | 34.8 | 31.9 | 9.8 |
$\tilde{R}_{m}^{2}(\sqrt{n})$ |  |  | 13.4 | 51.8 | 82.2 | 92.4 |  |  |  | 0.4 | 0.2 | 0.3 | 0.3 |  |  |  | 86.2 | 48.0 | 17.5 | 7.3 |
$\tilde{R}_{m}^{2}(1)$ |  |  | 35.3 | 69.0 | 80.8 | 87.7 |  |  |  | 2.9 | 3.4 | 4.7 | 4.6 |  |  |  | 61.8 | 27.6 | 14.5 | 7.7 |
$EAIC$ |  |  | 53.9 | 88.6 | 94.4 | 96.2 |  |  |  | 2.4 | 4.3 | 5.3 | 0.0 |  |  |  | 43.7 | 7.1 | 0.3 | 3.8 |
Dispersion | $\tilde{R}_{d}^{2}(\log n)$ |  |  | 13.8 | 34.6 | 38.4 | 80.8 |  |  |  | 0.2 | 1.4 | 1.3 | 6.5 |  |  |  | 86.0 | 64.0 | 60.3 | 12.7 |
$\tilde{R}_{d}^{2}(\sqrt{n})$ |  |  | 11.9 | 33.8 | 66.2 | 83.6 |  |  |  | 0.4 | 1.8 | 4.3 | 4.9 |  |  |  | 87.7 | 64.4 | 29.5 | 11.5 |
$\tilde{R}_{d}^{2}(1)$ |  |  | 11.2 | 33.4 | 65.5 | 79.4 |  |  |  | 0.8 | 2.3 | 4.1 | 6.4 |  |  |  | 88 | 64.3 | 30.4 | 14.5 |
$AIC_{c}$ |  |  | 17.3 | 43.5 | 69.7 | 83.6 |  |  |  | 0.4 | 3.1 | 4.1 | 5.5 |  |  |  | 81,9 | 53,4 | 26,2 | 10,9 |
Figure 2: Graphs, obtained from the simulations of the Binomial distribution,
regarding the relationship between sample size and percentage found for the
Optimal and Type 1 models, for the $\tilde{R}_{m}^{2}$, $\tilde{R}_{d}^{2}$,
$EAIC$ and $AIC_{c}$ criteria, where the solid line represents the graph for
$\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$ with $\lambda_{n}=\log n$, the
dashed line represents the graph for $\tilde{R}_{m}^{2}$
$(\tilde{R}_{d}^{2})$, with $\lambda_{n}=\sqrt{n}$, the dotted line represents
the graph for $\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$, with $\lambda_{n}=1$
and the dash-dotted line represents the graph for $EAIC(AIC_{c})$.

In the Binomial case, the simulations had a low success rate for the Optimal
models for 25 sample size, however there was a rapid growth of this rate,
reaching excellent values for 150 sample size. This behavior was observed for
both the mean and dispersion models. Unlike the simulations for the Normal
case, the $EAIC$ criterion showed better results for the mean model, however
it is observed that the criteria converge to a similarity when the sample size
is increased. For smaller samples, we observed a large difference between the
criteria, indicating that in this case the $EAIC$ criterion must be chosen,
but for large samples both criteria can be chosen. For the dispersion model,
the $AIC_{c}$ criterion was slightly better than the other criteria, although
in some cases the values were close. These results show that the choice of the
criterion must be made with care.

The graphs displayed in Figure 2, for the percentages shown in Table 3, show
what was discussed before, where we can see the lines, which are divergent for
small samples, converging to a common point. The dash-dotted line is generally
above the others in the Optimal models graphs and below the others in the Type
1 models graphs, indicating the best result for the $EAIC$ ($AIC_{c}$)
criterion for small samples and the approximation of lines at 150 sample size
indicating that the criteria converge to be similar.

###  4.3 Simulation study for Poisson distribution

For the mean model with Poisson probability distribution or Poisson type,
i.e., with $V(\mu_{i})=\mu_{i}$, the data simulation process must ensure that
$E(Y_{i})=\mu_{i}$ and $Var(Y_{i})=\phi_{i}\mu_{i}$. The generation of data is
adapted from (1989), p. 198. Consider the random variable
$Y|N=S_{1}+\ldots+S_{N}$, where $S_{j}$ are independent and identically
distributed random variables and $N$ is a random variable with Poisson
distribution, independent of $S_{j}$, $\forall j=1,\ldots,N$. In this way we
can obtain $E(Y_{i})=E(N_{i})E(S_{i})$ and $Var(Y_{i})=E(N_{i})E(S_{i}^{2})$.
Note that there will be overdispersion if $Var(Y_{i})>E(Y_{i})$, that is, if
$E(S_{i}^{2})>E(S_{i})$. Thus, in order to obtain the desired result, it is
sufficient to take $N_{i}\sim Poisson(\mu_{i}/\rho_{i})$ and $S_{i}\sim
Poisson(\rho_{i})$ with $\rho_{i}=\phi_{i}-1>0$.

The values chosen for the parameter vector
$(\beta_{0},\beta_{1},\beta_{2},\beta_{3},$
$\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ used to perform simulations
were $(1.5,3,2,0,0.2,0,3,0)$, where the zero value means that the parameter
was not used in the simulation. For each Monte Carlo simulation, the responses
were generated considering log link function for both mean and dispersion
models.

Table 4: Percentages of model types found in 1000 Monte Carlo simulations of the proposed variable selection procedure for the Poisson distribution |  |  |  | Optimal |  |  |  | Type 2 |  |  |  | Type 1 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | Sample size |  |  |  | Sample size |  |  |  | Sample size |
Model | Criterion |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |  |  |  | 25 | 50 | 100 | 150 |
Mean | $\tilde{R}_{m}^{2}(\log n)$ |  |  | 88.3 | 90.5 | 91.1 | 92.4 |  |  |  | 11.3 | 9.3 | 8.8 | 7.5 |  |  |  | 0.4 | 0.2 | 0.1 | 0.1 |
$\tilde{R}_{m}^{2}(\sqrt{n})$ |  |  | 92.6 | 92.6 | 93.3 | 95.3 |  |  |  | 7.00 | 6.8 | 6.5 | 4.6 |  |  |  | 0.4 | 0.6 | 0.2 | 0.1 |
$\tilde{R}_{m}^{2}(1)$ |  |  | 89.5 | 90.7 | 91.3 | 90.8 |  |  |  | 10.4 | 9.2 | 8.7 | 9.2 |  |  |  | 0.1 | 0.1 | 0.1 | 0.0 |
$EAIC$ |  |  | 85.9 | 91.8 | 91.3 | 92.4 |  |  |  | 6.7 | 6.6 | 8.6 | 7.5 |  |  |  | 7.4 | 1.6 | 0.1 | 0.1 |
Dispersion | $\tilde{R}_{d}^{2}(\log n)$ |  |  | 10.6 | 27.3 | 47.3 | 58.3 |  |  |  | 0.4 | 3.4 | 9.5 | 13.4 |  |  |  | 89.0 | 69.3 | 43.2 | 28.3 |
$\tilde{R}_{d}^{2}(\sqrt{n})$ |  |  | 13.5 | 27.2 | 51.6 | 59.8 |  |  |  | 0.3 | 2.3 | 7.7 | 13.1 |  |  |  | 86.2 | 70.5 | 40.7 | 27.1 |
$\tilde{R}_{d}^{2}(1)$ |  |  | 11.2 | 29.7 | 46.0 | 57.6 |  |  |  | 1.0 | 3.2 | 9.2 | 15.6 |  |  |  | 87.8 | 67.1 | 44.8 | 28.8 |
$AIC_{c}$ |  |  | 51.7 | 64.9 | 78.8 | 84.0 |  |  |  | 2.4 | 3.1 | 3.1 | 4.7 |  |  |  | 51,7 | 32.0 | 18.1 | 11.3 |
Figure 3: Graphs, obtained from the simulations of the Poisson distribution,
regarding the relationship between sample size and percentage found for the
Optimal and Type 1 models, for the $\tilde{R}_{m}^{2}$, $\tilde{R}_{d}^{2}$
and $EAIC$ criteria, where the solid line represents the graph for
$\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$ with $\lambda_{n}=\log n$, the
dashed line represents the graph for $\tilde{R}_{m}^{2}$
$(\tilde{R}_{d}^{2})$, with $\lambda_{n}=\sqrt{n}$, the dotted line represents
the graph for $\tilde{R}_{m}^{2}$ $(\tilde{R}_{d}^{2})$, with $\lambda_{n}=1$
and the dash-dotted line represents the graph for $EAIC(AIC_{c})$.

Table 4 shows that, for the mean models, the procedure for variable selection
remains efficient for the Poisson distribution and we even have higher
efficiency for smaller samples compared to previous distributions. The
$\tilde{R}_{m}^{2}$ with $\lambda_{n}=\sqrt{n}$ criterion had a slight
advantage over the other criteria.

For the dispersion models, the $AIC_{c}$ criterion was better, with an high
rate of Optimal models for large samples. The $\tilde{R}_{d}^{2}$ criterion
had an improvement in the Optimal models rate with the growth of the sample
size, reaching close to $60\%$. Also note that Type 1 models showed high rates
for small samples, however there was an accelerated decrease in that rate when
the sample size was increased.

In the graphs displayed in Figure 3, we can see lines that are above the
others, that is, for the mean models we observe the advantage of the
$\tilde{R}_{m}^{2}$ criterion with $\lambda_{n}=\sqrt{n}$ (dashed line) and
for the dispersion models we observe the advantage of the $AIC_{c}$ criterion
(dash-dotted line).

For type 1 models, the $AIC_{c}$ criterion was better for dispersion models,
with the dash-dotted line always below the others. For the mean models, the
criterion $\tilde{R}_{m}^{2}$, for any value of $\lambda_{n}$, was better than
$EAIC$ for small samples, however the criteria become equivalent as the sample
size increases.

##  5 Application

The proposed variable selection procedure was applied to data from an
injection molding experiment.

The experiment was performed to study the influence of seven controllable
factors and three noise factors on the mean value and the variation in the
percentage of shrinkage of products made by injection molding. Noise factors
are fixed during the experiment but are expected to vary randomly outside the
experimental context.

The aim of the experiment was to determine the process parameter settings so
that the shrinkage percentage was close to the target value and robust against
environmental variations. The responses were percentages of shrinkage of
products made by injection molding, the design for seven controllable factors
( $A=$ cycle time, $B=$ mould temperature, $C=$ cavity thickness, $D=$ holding
pressure, $E=$ injection speed, $F=$ holding time, and $G=$ gate size) was a
$2^{7-4}$ fractional factorial. At each setting of the controllable factors,
four observations were obtained from a $2^{3-1}$ fractional factorial with
three noise factors ($M=$ percentage regrind, $N=$ moisture content and $O=$
ambient temperature). The experimental data from the injection molding are
shown in Table 5.

The data set considered is well known in the literature of industrial
experiments and has been analyzed by several authors such as (1992), (1996)
and (1998). Although the articles by (1996) and (1998) have used the JMMD, the
variable selection process was not discussed. In our approach, in addition to
obtaining estimates of the parameters of the mean and dispersion models by
using JMMD, we also built a consistent investigative procedure to verify the
adequacy of the final joint model found. For the application of Algorithm 2,
we considered $V(\mu)=1$ and identity link function for the mean model. For
the dispersion model we considered the Gamma distribution with logarithmic
link function. The joint models we considered for mean and dispersion were the
same as those proposed by (1996) and (1998).

The selection measure considered in Algorithm 1 for the mean model was
$\tilde{R}_{m}^{2}$ with $\lambda_{n}=\sqrt{n}$ and for the dispersion model
the selection measure was the $AIC_{c}$. As explained in Section 4, the level
of significance for the hypothesis tests performed was $0.10$.

Table 6 presents the results of all steps of the variable selection process
(Algorithm 2) taken to obtain the final joint model for the mean and
dispersion. In each iteration of Algorithm 2, Algorithm 1 was used to find the
terms of the mean and dispersion models. The additional terms to the current
models, which appear in sequence in Table 6, are those that, for all terms not
belonging to the current model, have the highest $\tilde{R}_{m}^{2}$ value, in
case of the mean model, and the smallest $AIC_{c}$, in the case of the
dispersion model. For these terms, which were candidates to enter the current
model, Chi-square hypothesis tests, in the case of the dispersion model, and
$F$ test, in the case of the mean model, were used to verify if, in fact, the
term chosen must be added to the current model. Measurements $D$ and $D^{*}$
were used to obtain test statistics. The values of $D$, $D^{*}$, $\chi^{2}$,
$F$ and $\tilde{R}_{m}^{2}$, with $\lambda_{n}=\sqrt{n}$ and with
$\lambda_{n}=1$, are shown for each step of the Algorithm 1, where $D$ is the
deviance for the Gamma model,
$D^{*}=\sum_{i=1}^{n}\frac{1}{\phi_{i}}d_{i}^{*}$ (see Algorithm 1),
$\chi^{2}$, the test statistic for the Chi-square test in the dispersion
model, is the deviance difference for nested Gamma models and $F$ is the test
statistic, given by equation (3.6), for the $F$ test in the mean model. It is
worth mentioning that the measure $\tilde{R}_{m}^{2}$ with $\lambda_{n}=1$,
which was not used as a selection measure, was presented in Table 6 only to
show the impact on the value of $\tilde{R}_{m}^{2}$, when the number of terms
in the model increases, considering the penalty with $\lambda_{n}=\sqrt{n}$
and with $\lambda_{n}=1$. In each iteration, the final models for both mean
and dispersion are shown in bold.

Following Algorithm 2, since the value of $\tilde{R}_{m}^{2}$ for the mean
model, obtained in the third iteration, was lower than that obtained in the
second iteration, the procedure stopped. In this way, three iterations of
Algorithm 2 were performed, obtaining at the end, in the second iteration, a
model for the mean with $\tilde{R}_{m}^{2}(\sqrt{n})$ equal to 0.803 and
$\tilde{R}_{m}^{2}(1)$ equal to 0.973. All the mean models found in the three
iterations of Algorithm 2 had the same terms (see Table 6), differing only in
the value of the effect estimates. Table 7 presents the final estimates, with
their standard deviations and Wald tests. Note that all estimates were
significant at the 1% level, except the terms $E$, $N$ and $C$ in the mean
model, which entered in the model only after the variable selection procedure
ends, due to the $CN$ and $EN$ interaction terms being significant.

The diagnostic graphs for the mean model, displayed in Figure 4, show that the
model obtained is adequate, since it showed no tendency to violate the
assumptions of independence, constant variance and normality of the residuals,
there is no evidence of inadequate of the link and variance functions or
points of influence.

For the dispersion model, the models obtained in each iteration of Algorithm 2
were different. This is due to the fact that the response variable, $d^{*}$,
is not the same in each iteration of Algorithm 2. The final dispersion model
had a value of $AIC_{c}$ equal to $-74.2471$ (see Table 6) and showed highly
significant effects by the Wald test (see Table 7). Diagnostic graphs for the
dispersion model are presented in Figure 5.

Table 5: Experimental data from the injection molding |  |  |  |  |  |  |  |  |  | Noise factors
---|---|---|---|---|---|---|---|---|---|---
|  | $M$ | $-1$ | $-1$ | 1 | 1
|  | $N$ | $-1$ | 1 | $-1$ | $1$
| Controllable factors |  | $O$ | $-1$ | 1 | 1 | $-1$
run |  | $A$ | $B$ | $C$ | $D$ | $E$ | $F$ | $G$ |  |  |  |  |  |
1 |  | $-1$ | $-1$ | $-1$ | $-1$ | $-1$ | $-1$ | $-1$ |  |  | 2.2 | 2.1 | 2.3 | 2.3
2 |  | $-1$ | $-1$ | $-1$ | 1 | 1 | 1 | 1 |  |  | 2.5 | 0.3 | 2.7 | 0.3
3 |  | $-1$ | 1 | 1 | $-1$ | $-1$ | 1 | 1 |  |  | 0.5 | 3.1 | 0.4 | 2.8
4 |  | $-1$ | 1 | 1 | 1 | 1 | $-1$ | $-1$ |  |  | 2.0 | 1.9 | 1.8 | 2.0
5 |  | 1 | $-1$ | 1 | $-1$ | 1 | $-1$ | 1 |  |  | 3.0 | 3.1 | 3.0 | 3.0
6 |  | 1 | $-1$ | 1 | 1 | $-1$ | 1 | $-1$ |  |  | 2.1 | 4.2 | 1.0 | 3.1
7 |  | 1 | 1 | $-1$ | $-1$ | 1 | 1 | $-1$ |  |  | 4.0 | 1.9 | 4.6 | 2.2
8 |  | 1 | 1 | $-1$ | 1 | $-1$ | $-1$ | 1 |  |  | 2.0 | 1.9 | 1.9 | 1.8
Table 6: Steps of variable selection algorithm for JMMD applied to data from the injection molding experiment Iteration | Dispersion model | $AIC_{c}$ |  $D$ ‡ |  $\chi^{2}_{c}$ value⋆ | $Pr(>\chi^{2}_{c})$ |
---|---|---|---|---|---|---
1 | 1 | - | - | - | - |
Mean model | $\tilde{R}_{m}^{2}(\sqrt{n})$ | $\tilde{R}_{m}^{2}(1)$ |  $D^{*}$ ‡ |  $F_{c}$ value | $Pr(>F_{c})$
1 | - | - | 34.6839 | - | -
1+CN | -0.0060 | 0.3063 | 24.0587 | 13.2491 | 0.0010
1+CN+EN | 0.2232 | 0.5974 | 13.9628 | 20.9687 | 0.0001
1+CN+EN+A | 0.3234 | 0.7735 | 7.8557 | 21.7672 | 0.0001
1+CN+EN+A+D | -0.0782 | 0.8516 | 5.1467 | 14.2120 | 0.0008
Iteration | Dispersion model | $AIC_{c}$ |  $D$ ‡ |  $\chi^{2}_{c}$ value⋆ | $Pr(>\chi^{2}_{c})$ |
2 | 1 | - | 80.8899 | - | - |
1+E | -64.1778 | 72.0353 | 4.4273 | 0.0354 |
1+E+B | -68.4785 | 61.7544 | 5.1404 | 0.0234 |
1+E+B+G | -74.2471 | 50.8136 | 5.4704 | 0.0193 |
1+E+B+G+D | -77.5280 | 44.3609 | 3.2263 | 0.0725 |
Mean model | $\tilde{R}_{m}^{2}(\sqrt{n})$ | $\tilde{R}_{m}^{2}(1)$ |  $D^{*}$ ‡ |  $F_{c}$ value | $Pr(>F_{c})$
1 | - | - | 1,098.102 | - | -
1+A | 0.133 | 0.402 | 699.342 | 17.106 | 0.0003
1+A+CN | 0.378 | 0.678 | 404.023 | 21.197 | 0.0001
1+A+CN+EN | 0.841 | 0.947 | 52.454 | 187.666 | 0.0000
1+A+CN+EN+D | 0.803 | 0.973 | 26.624 | 26.195 | 0.0000
Iteration | Dispersion model | $AIC_{c}$ |  $D$ ‡ |  $\chi^{2}_{c}$ value⋆ | $Pr(>\chi^{2}_{c})$ |
3 | 1 | - | 75.4049 | - | - |
1+D | 42.5318 | 65.4232 | 4.9908 | 0.0255 |
1+D+F | 41.9160 | 61.3397 | 2.0417 | 0.1530 |
Mean model | $\tilde{R}_{m}^{2}(\sqrt{n})$ | $\tilde{R}_{m}^{2}(1)$ |  $D^{*}$ ‡ |  $F_{c}$ value | $Pr(>F_{c})$
1 | - | - | 59.6112 | - | -
1+CN | -0.0078 | 0.3051 | 41.8419 | 12.7404 | 0.0012
1+CN+EN | 0.2266 | 0.5992 | 24.3629 | 20.8059 | 0.0001
1+CN+EN+A | 0.5173 | 0.8384 | 9.4400 | 44.2627 | 0.0000
1+CN+EN+A+D | 0.2260 | 0.8935 | 6.2942 | 13.4942 | 0.0010
‡$D$ is the deviance for the Gamma model and
$D^{*}=\sum_{i=1}^{n}\frac{1}{\phi_{i}}d_{i}^{*}$.
⋆$\chi^{2}_{c}$ is the value of the Chi-square statistic, calculated from the
difference in deviations
for two nested Gamma models.
⋄$F_{c}$ is the value of the F-test statistic, calculated by the equation
(3.6).
Table 7: Regression coefficients and Wald test for the mean and dispersion models obtained from injection molding data Mean model |  |  |  |
---|---|---|---|---
Terms | Estimate | Std. Error |  $t$ value | $Pr(>|t|)$
1 | 2.24903 | 0.03322 | 67.693 | 0.0000
$A$ | 0.42802 | 0.06575 | 6.509 | 0.0000
$C$ | 0.07172 | 0.05802 | 1.236 | 0.2279
$N$ | -0.00433 | 0.05994 | -0.072 | 0.9430
$D$ | -0.28639 | 0.06407 | -4.470 | 0.0001
$E$ | 0.06528 | 0.05971 | 1.093 | 0.2851
$CN$ | 0.58684 | 0.03322 | 17.663 | 0.0000
$EN$ | -0.55727 | 0.05994 | -9.297 | 0.0000
Dispersion model |  |  |  |
Terms | Estimate | Std. Error |  $t$ value | $Pr(>|t|)$
1 | -2.2973 | 0.1754 | -13.097 | 0.0000
$E$ | -0.8670 | 0.1754 | -4.942 | 0.0000
$B$ | 0.6773 | 0.1754 | 3.861 | 0.0006
$G$ | -0.6015 | 0.1754 | -3.429 | 0.0019
Figure 4:  Diagnostic plots for mean model obtained from injection molding
data.  Figure 5: Diagnostic plots for dispersion model obtained from injection
molding data.

##  6 Conclusion

The article presents a procedure for variable selection in joint modeling of
the mean and dispersion with application in robust design experiments. The
proposed procedure provides an iterative mechanism for selecting variables in
the JMMD. The process takes advantage of the seesaw structure from the
algorithm to estimate the parameters in the JMMD, where the weights from the
dispersion model are incorporated into the mean model. For both mean and
dispersion models, a selection measure is used to choose, in each iteration of
the procedure, the variables that contribute to the best fit of the model.
Appropriate hypothesis tests are used to confirm the choice of these variables
or to terminate the procedure. As the weights from the dispersion model are
incorporated into the mean model, the adjustment of the mean model becomes
better. The final joint model obtained is the one that has the best fit for
the mean model. In this way, the dispersion model is constituted by the terms
that provide the best weights for the mean model.

For each of the mean and dispersion models, two alternative selection measures
were used in the simulation studies. The $\tilde{R}_{m}^{2}$ and $EAIC$
selection measures were used in the mean model and $\tilde{R}_{d}^{2}$ and
$AIC_{c}$ were used in the dispersion model. Each of the $\tilde{R}_{m}^{2}$
and $\tilde{R}_{d}^{2}$ criteria were used considering the variations based on
the choices of $\lambda_{n}$ equal to $\sqrt{n}$, $\log n$ and 1.

In the simulation studies, performed for the Normal, Poisson and Binomial
distributions, the $\tilde{R}_{m}^{2}$ criterion with $\lambda_{n}=\sqrt{n}$
was better for the Normal and Poisson cases, while the $EAIC$ was better for
the Binomial case. The $AIC_{c}$ was better in all cases considered.

The efficiency of the proposed procedure depends on the criterion used as a
selection measure and the sample size. As should be expected, the efficiency
grows with the increase of the sample size. The simulations showed that, for a
large sample size, the proposed procedure considerably reduces the risk of
finding poor models and ensures that satisfactory models are obtained for mean
and dispersion. However, the procedure proved to be reasonably good for small
samples.

In application, discussed in Section 5, the procedure for selecting variables
was explained step by step and the results obtained for the mean and
dispersion models were satisfactory, indicating the effectiveness of the
proposed methodology and corroborating the results found by the simulations.

[Acknowledgments] The authors would like to thank the anonymous referees and
the Editor for their constructive comments that improved the quality of this
paper.

This work was supported by Research Support Foundation of the State of Minas
Gerais (FAPEMIG) under Grant Number APQ-03365-18.

## References

  * (1) Aitkin, M. (1987). Modelling variance heterogeneity in normal regression using GLIM. Journal of the Royal Statistical Society \- Series C - Applied Statistics, 36, 332-339.
  * (2) Antoniadis, A., Gijbels, I., Lacroix, S.L. and Poggi, J.M. (2016). Joint estimation and variable selection for mean and dispersion in proper dispersion models. Electronic Journal of Statistics, 10, 1630-1676.
  * (3) Bayer, F.M. and Cribari Neto, F. (2017). Model selection criteria in beta regression with varying dispersion. Communications in Statistics - Simulation and Computation, 46, 729-746.
  * (4) Bonat, W.H. and Jorgensen, B. (2016). Multivariate covariance generalized linear models. Journal of the Royal Statistical Society \- Series C - Applied Statistics, 65, 649-675.
  * (5) Carroll, R.J., Ruppert, D. (1988). Transformation and weighting in regression. New York: Chapman and Hall.
  * (6) Charalambous, C., Pan, J. and Tranmer, M. (2014). Variable selection in joint mean and dispersion models via double penalized likelihood. Sankhyã: The Indian Journal of Statistics, 76, 276-304.
  * (7) Charalambous, C., Pan, J. and Tranmer, M. (2015). Variable selection in joint modelling of the mean and variance for hierarchical data. Statistical Modelling, 15, 24-50.
  * (8) Cuervo, E. and Gamerman, D. (2001). Bayesian modeling of variance heterogeneity in normal regression models. Brazilian Journal of Probability and Statistics, 14, 207-221.
  * (9) Davidian, M. and Carroll, R.J. (1987). Variance function estimation. Journal of the American Statistical Association, 82, 1079-1091.
  * (10) Davidian, M. and Carroll, R.J. (1988). A note on extended quasi-likelihood. Journal of the Royal Statistical Society, Series B, 50, 74-82.
  * (11) Efron, B. (1986). Double exponential families and their use in generalized linear regression. Journal of the American Statistical Association, 81, 709–721.
  * (12) Engel, J. (1992). Modeling variation in industrial experiments. Applied Statistics, 41, 579-593.
  * (13) Engel, J. and Huele, A. F. (1996). A generalized linear modeling approach to robust Design. Technometrics, 38, 365-373.
  * (14) Godambe, V.P. and Thompson, M.E. (1989). An extension of quasi-likelihood estimation. Journal of Statistical Planning and Inference, 22, 137-152.
  * (15) Hu, B. and Shao, J. (2008). Generalized linear model selection using $R^{2}$. Journal of Statistical Planning and Inference, 138, 3705-3712.
  * (16) Hurvich, C.M. and Tsai, C.L. (1989). Regression and time series model selection in small samples. Biometrika, 76, 297-307.
  * (17) Jorgensen, B. (1987). Exponential dispersion models (with discussion). Journal of the Royal Statistical Society, Series B, 49, 127–162.
  * (18) Jorgensen, B. (1992). The theory of exponential dispersion models and analysis of deviance. Rio de Janeiro: IMPA.
  * (19) Jorgensen, B. (1997). The theory of dispersion models. New York: Chapman and Hall.
  * (20) Lee, K. I., Koval, J. J. (1997). Determination of the best significance level in forward stepwise logistic regression. Communications in Statistics - Simulation and Computation, 26, 559-575.
  * (21) Lee, Y and Nelder, J.A. (1996). Hierarchical generalized linear models. Journal of the Royal Statistical Society, Series B, 58, 619-678.
  * (22) Lee, Y. and Nelder, J.A. (1998). Generalized linear models for analysis of quality improvement experiments. The Canadian Journal of Statistics, 26, 95-105.
  * (23) Lee, Y. and Nelder, J.A. (2000). The relationship between double-exponential families and extended quasi-likelihood families, with application to modelling Geissler’s human Sex ratio data. Applied Statistics, 49, 413-419.
  * (24) Lee, Y. and Nelder, J.A. (2001). Hierarchical generalised linear models: A synthesis of generalised linear models, random-effect models and structured dispersions. Biometrika, 88, 987-1006.
  * (25) Lee, Y., Nelder, J.A., Pawitan, Y. (2006). Generalized Linear Models with Random Effects: Unified Analysis via H-likelihood. Boca Raton: Chapman and Hall.
  * (26) McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models. 2nd ed., London: Chapman and Hall.
  * (27) Nair, V.N., Ed. (1992). Taguchi’s Parameter Design: A Panel Discussion. Technometrics, 34, 127-161.
  * (28) Nelder, J.A. (2000). Quasi-likelihood and pseudo-likelihood are not the same thing. Journal of Applied Statistics, 27, 1007-1011.
  * (29) Nelder, J.A. and Lee, Y. (1991). Generalized linear models for the analysis of Taguchi-type experiments. Applied Stochastic Models and Data Analysis, 7, 107-120.
  * (30) Nelder, J.A. and Lee, Y. (1992). Likelihood, quasi-likelihood and pseudolikelihood: Some Comparisons. Journal of the Royal Statistical Society, Series B, 54, 273-284.
  * (31) Nelder, J.A. and Pregibon, D. (1987). An extended quasi-likelihood function. Biometrika, 74, 221-232.
  * (32) Nelder, J.A. and Wedderburn, R.W.M. (1972). Generalized linear models. Journal of the Royal Statistical Society, Series A, 135, 370-384.
  * (33) Paul, S. and Saha, K.K. (2007). The generalized linear model and extensions: a review and some biological and environmental applications. Environmetrics, 18, 421-443.
  * (34) Phadke, M.S. (1989). Quality Engineering Using Robust Design, New York: Prentice-Hall.
  * (35) R Development Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austrian. URL http://www.R-project.org/.
  * (36) Smyth, G.K. (1989). Generalized linear models with varying dispersion. Journal of The Royal Statistical Society – Series B, 51, 47-60.
  * (37) Smyth G.K. and Verbyla A.P. (1999a). Double generalized linear models: approximate REML and diagnostics. In Statistical Modelling: Proceedings of the 14th International Workshop on Statistical Modelling (IWSM14) ( Editors, Friedl H., Berghold A. and Kauermann G. eds.) 66-88. Graz: Austrian.
  * (38) Smyth G.K. and Verbyla A.P. (1999b). Adjusted likelihood methods for modeling dispersion in generalized linear models. Environmetrics, 10, 696-709.
  * (39) Sugiura, N. (1978). Further analysis of the data by Akaike’s information criterion and the finite corrections. Communications in Statistics - Theory and Methods, 7, 13-26.
  * (40) Taguchi, G. (1986). Introduction to Quality Engineering. New York: Unipub/Kraus International Publications.
  * (41) Verbyla, A.P. (1993). Modelling variance heterogeneity: residual maximum likelihood and diagnostics. Journal of the Royal Statistical Society – Series B, 55, 493-508.
  * (42) Wang, D. and Zhang, Z. (2009). Variable selection in joint generalized linear models. Chinese Journal of Applied Probability and Statistics, 25, 245-256.
  * (43) Wu, L. and Li, H. (2012). Variable selection for joint mean and dispersion models of the inverse Gaussian distribution. Metrika, 75, 795-808.
  * Wu et al. (2020) Wu, L, Li, S. and Tao, Y. (2020). Estimation and variable selection for mixture of joint mean and variance models. Communications in Statistics - Theory and Methods, 1-18.
  * (45) Wu, L., Zhang, Z. and Xu, D. (2012). Variable selection for joint mean and dispersion models of the lognormal distribution. Hacettepe Journal of Mathematics and Statistics, 41, 307-320.
  * Zhang (2017) Zhang, D. (2017). A coefficient of determination for generalized linear models. The American Statistician, 71, 310-316.