数理解析研究所講究録
第 1910 巻 2014 年 20-28

20Model Selection Based on Quasi-likelihood with
Application to Overdispersed Data
Yiping Tang and Jinfang Wang
Graduate School of Science, Chiba University
1-33 Yayoi-cho, Inage-ku, Chiba 263-8522 Japan
E-mail: tou.kihei@gmail.com
Abstract

In analyzing complicated data, we are often unwilling or not confident to impose
a parametric model for the data-generating structure. One important example is data
analysis for proportional or count data with overdispersion. The obvious advantage
of assuming full parametric models is that one can resort to likelihood analyses, forinstance, to use AIC or BIC to choose the most appropriate regression models. For
overdispersed proportional data, possible parametric models include the Beta-binomial
models, the double exponential models, etc. In this paper, we extend the generalized
linear models by replacing the full parametric models with a finite number of moment
restrictions on both the data and the structural parameters. For such semiparametric statistical models, we propose a method for selecting the best possible regression
model in the semiparametric model class. We will apply the proposed model selection
technique to overdispersed data. We will demonstrate the use of the proposed semiparametric information criterion using the well-known data on germination of 0robanche.
Key Words: Akaike information criterion, Generalized linear models, KullbackLeibler information, Overdispersion, Quasi-likelihood

1

Introduction

Let $y_{i}(i=1, \ldots, n)$ be independent response variables which may be proportions or counts. As in the usual generalized hnear models (GLM; McCulwith mean and variancelagh and Nelder, 1989), we assume in addition that $E[y_{i}]=\mu_{i}=t^{-1}(\eta_{i})$ , whereis a vector of the unis a vector of covariates,is the inverse link function,and does not depend on . Let $X$ be theknown parameters
and $q<$ n.In
matrix of covariates with the ith row being
many situations, it is diﬃcult to impose a specific distributional assumption on as in theclassical GLM. In such cases we may use the quasi-likelihood method (Wedderburn, 1974)
to estimate the regression parameters . Under weak regularity conditions, it can be shownthat the quasi-likelihood estimator is consistent and asymptotically normally distributed$\mu_{i}$$\phi V(\mu_{i})$$\eta_{i}=x_{i}’\beta,$$t^{-1}(\cdot)$$\beta$$q$$x_{i}$$\beta’=(\beta_{1}, \ldots, \beta_{q})$$q$$n\cross q$$\beta$$\phi$$X_{(i)},$$\theta=(\beta’, \phi)’\in\Theta\subset \mathbb{R}^{q+1}$$y_{\iota’}$$\beta$$\hat{\theta}$(Moore, 1986).

21
However, works on model selection based on quasi-likelihoods are rather limited. There
are mainly two approaches for existing works in this area. The first approach is to restrict
the choice of the variance function
so that one can integrate the quasi-score function to
obtain the scalar quasi-likelihood function and formally use AIC with the likelihood function being replaced by the quasi-likelihood function (e.g., Lebreton, et al. (1992), Anderson,
Burnham and White (1994), Qian, Gabor and Gupta (1996), Pan (2001)). The second approach makes no restriction on the choice of the variance functions, instead the quasi-score
vector is projected on a subspace of estimating functions so that the projection becomesintegrable (e.g., McLeish and Small (1992), Li (1993), Hanfelt and Liang (1995)), consequently one may use the first approach to formally construct the information criterion (e.g.,
Lin (2011)). In this paper we shall propose a diﬀerent approach by projecting the true distribution function onto the subspace of probability distribution functions satisfying the first
and the second moment assumptions. This projection will be used as the semiparametricpredictive distribution for the underlying statistical model. Using the obtained semiparametric predictive distribution, we then formally extend the idea underlying the derivationof the information criterion AIC to construct a new semiparametric information criterion
SIC. We shall demonstrate the use of the proposed semiparametric information criterion SIC
using the well-known data on germination of 0robanche.
$V$

2Semiparametric Model Selection Criteria

Akaike (1973) proposed a model selection criterion called AIC which has been widely usedin many areas of applications. Derivation of AIC is based upon minimizing the KullbackLeibler information between the predictive distribution and the true distribution generating
the data. In this section we will extend this idea to the semiparametric setting with thefull distributional assumption replaced by the assumptions on the mean and variance of the
response variables.
We propose semiparametric model selection criteria SIC and $SIC_{T}$ . They apply generally
to data no matter whether they are discrete or continuous. For convenience, we suppose thatthe data in the subsequent paragraphs are continuous. If they are discrete then integrationis replaced by summation.

2.1

Semiparametric models

We assume that the true probability density function of $Y$ is $h(y)$ and the true distributionfunction is $H(y)$ . We will relax the assumption of a complete specification of the probability
model, only using some general form of the variance function about the mean as used in
GLM. Now suppose that the statistician makes the assumption that the random variables$y_{i}(i=1, \ldots, n)$ have distribution
density functionmean
and variance. The first and the second moments can be combined in the following way:
$G$

$g$

$\mu_{i}(\beta)$$\phi V(\mu_{i}(\beta))$$m(y_{i}, \theta)=m(y_{i}, \beta, \phi)=((y_{i}-\mu_{i}(\beta))^{2}-\phi V(\mu_{i}(\beta))y_{i}-\mu_{i}(\beta))=(\begin{array}{l}m_{1i}m_{2i}\end{array})$22
By assumption, we have that
$E_{G}[m(y_{i}, \theta)]=\int m(y_{i}, \theta)dG(y)=0.$

Let $M$ denote the set of all possible distribution functions, let$\mathcal{G}(\theta)=\{G\in M:\int m(y, \theta)dG(y)=0\}.$Defineas the set of all possible distribution functions that are restrictedby the moments conditions. The class is the statistical model we shall consider in thispaper.
$\mathcal{G}=\bigcup_{(\theta\in\Theta)}\mathcal{G}(\theta)$$\mathcal{G}$2.2

Kullback-Leibler information

Extending the idea of AIC, we shall generalize the data-based model selection approachbased also on the Kullback-Leibler information. We use the Kullback-Leibler informationLet
to measure the closeness between $H(y)$ and the assumed semiparametric model
vary in the statistical model class . We consider the constrained minimization problem
$G$

$g$

$\mathcal{G}$(2.1)$v( \beta, \phi)= \inf_{g\in \mathcal{G}}\int\log(\frac{h(y)}{g(y)})h(y)dy$

subject to$\int m(y, \beta, \phi)g(y)dy=0$

and$\int g(y)dy=1.$

Using the results in convex analysis (e.g., Kitamura, 2006), the solution to infinite dimencan be equivalently written as the solution to the finite
sional optimization problem
as follows:
dimensional optimization problem$v(\beta, \phi)$$v^{*}(\beta, \phi)$$v^{*}( \beta, \phi)=\sup_{\lambda\in R,r\in R^{2}}[1+\lambda+\int\log(-\lambda-r’m(y,\beta, \phi))h(y)dy]$

.

(2.2)

with respect to and , we canSince we are only interested in the maximization ofignore other terms and proceed to define the semiparametric information of the statisticalmodel as follows:$v^{*}(\beta, \phi)$$\beta$$\phi$$\mathcal{G}$DEFINITION 2.1 The semiparametric Kullback-Leibler information for quasi-likelihood function is defined as follows:$SKL(\theta)=\sup_{r\in R^{2}}E_{H}[\rho(y, \theta, r)]=\sup_{r\in R^{2}}\int\rho(y, \theta, r)h(y)dy$

where $\rho(y, \theta, r)=\log(1+r’m(y,$$\beta,$,(2.3)$\phi$

Since the true distribution function $H$ is unknown, we use the empirical distribution of
the data to estimate the true distribution function $H$

23
2.3

Regularity conditions

In this paper, we shall use $SKL(\theta)$ to select an appropriate semiparametric model. We willconsider both the correctly specified case and the misspecified case. A semiparametric modelis called correctly specified if the moments assumed evaluated at a particular value of theparameter are equal to the moments of the underlying true distribution. First, we state theregularity conditions which will be used in the sequel to prove the asymptotic unbiasednessof the proposed information criteria.
Cl Defineas the pseudo-true value, such that the semiparametric Kullback$SKL(\theta)$
Leibler information
is minimized at
C2
is twice continuously diﬀerentiable inandNote that C2 implies that$\theta_{*}=(\beta_{*}’, \phi_{*})’$$\theta_{*}.$$r\in \mathbb{R}^{2}.$$\theta\in\Theta$$\rho(y, \theta, r)$(2.4)$\tilde{r}(\theta)=\arg\sup_{r\in R^{2}}\int\rho(y, \theta, r)h(y)dy$

is continuous with respect tocorresponding to :$\theta$. Furthermore, from Cl we know that there is a unique$\tilde{r}_{*}$$\theta_{*}$$\tilde{r}_{*}=\arg\sup_{r\in R^{2}}\int\rho(y, \theta_{*}, r)h(y)dy$

.

(2.5)

, whereC3 The problem of
has a unique saddle point solution
is defined in (2.5). When the model is correctly specified, there exists a uniquesatisfying $SKL(\theta_{*})=0.$
C4 The following law of large numbers holds:$(\theta_{*}’,\tilde{r}_{*}’)$$\inf_{\theta\in\Theta}SKL(\theta)$$\theta\sup_{\in\Theta}\sup_{r\in R^{2}}|\frac{1}{n}\sum_{i=1}^{n}\rho(y_{i}, \theta, r)-E_{H}[\rho(y, \theta, r)]|=o_{p}(1)$2.4$\tilde{r}_{*}$$(\theta_{*}’,\tilde{r}_{*}’)$$\mathcal{G}$.

The expected information

When we consider $SKL(\theta)$ as defined in (2.3), the minimizer of $SKL(\theta)$ can be viewed as
the solution to the following saddle point problem:$\theta_{*}=\arg\min\sup E_{H}[\rho(y, \theta, r)]$

(2.6)

,

$\theta\ominus$$L(\theta)=\{r : r’m(y, \theta)\in \mathcal{I}\}$ , andwhere denotes the parameter space foris an open
interval containing zero.
Next we describe the process of solving the above saddle point problem. Firstly, theexpectation ofsatisfies the following
is maximized for given , so that thecondition:$\Theta$$\mathcal{I}$$\theta,$$\rho(y, \theta, r)$$\theta$$\theta_{*}$$E_{H}[\rho_{1}(y, \theta, r)m(y, \theta)]=0,$

where $\rho_{1}(y, \theta, r)=1/(1+r’m(y, \theta))$ is the first derivative of with respect to $r’m(y, \theta)$ .
C5 Assume that $\int\rho(y, \theta, r)h(y)dy$ and $\int\rho_{1}(y, \theta, r)h(y)dy$ are diﬀerentiable under the integral sign.
$\rho$

24
Secondly, Since
that$\theta_{*}$

is the minimizer of the quantity$E_{H}[\rho(y,$for a proper , implying$r$$\theta,$$r$$E_{H}[\rho_{1}(y, \theta, r)M(y, \theta)’r(\theta)]=0,$

where
Now denote the matrix of the first derivative of$M(y, \theta)=\partial m(y, \theta)/\partial\theta’.$$\rho(y, \theta, r)$with respect to$\theta$

and

$r$

by

$\Upsilon(\theta,r)=(\begin{array}{l}\rho_{1}(y,\theta,r)M(y,\theta)’r\rho_{1}(y,\theta,r)m(y,\theta)\end{array}),$The values$(\theta_{*},\tilde{r}(\theta_{*}))$

must satisfy the following first-order condition:$E_{H}[T(\theta_{*},\tilde{r}(\theta_{*}))]=0$

.

(2.7)

C6 Define$Q=E_{H}[\frac{\partial\Upsilon(\theta_{*},\tilde{r}_{*})}{\partial(\theta’,r’)^{J}}], S=E_{H}[\Upsilon(\theta,\tilde{r}_{*})\Upsilon(\theta,\tilde{r}_{*}$

and are finite and is invertible.
However,
is not computable because it involves the expectation with respect to thetrue density function $h(y)$ . One way is to use the value to minimize the empirical version of
the expected semiparametric Kullback-Leibler information:

Assume that

$Q$

$S$

$S$

$\theta_{*}$,$\hat{\theta}_{*}=\arg_{\theta}\min_{\in\Theta}\{\frac{1}{n}\sum_{i=1}^{n}\rho(y_{i}, \theta,\tilde{r}(\theta))\}$

whereis defined in (2.4). To approximate
function in place of the true distribution function $H$$\tilde{r}(\theta)$$\tilde{r}(\theta)$(2.8), we use the empirical distribution
to obtain$\hat{r}(\theta)=\arg_{r}\max_{\in R^{2}}\{\frac{1}{n}\sum_{i=1}^{n}\rho(y_{i},\theta, r)\}$

.

(2.9)

We can show thatare consistent and asymptotically normally distributed along
similar lines of Chen et al. (2007), we omit the proof here. The calculation ofwith being theis instable and usually diﬃcult to compute. We assume that
guasi-likelihood estimator and the estimator from the usual Pearson residual. Since theis also a consistent estimator, it is reasonable to use
instead. It is easy to showthatholds the properties of consistency and asymptotic normality too.
$(\hat{\theta}_{*}, \hat{r}(\hat{\theta}_{*}))$$(\hat{\theta}_{*}, \hat{r}(\hat{\theta}_{*}))$$\hat{\theta}=(\sqrt{}’,$$\phi\wedge$$\sqrt{}\wedge$$\hat{\phi}$$(\hat{\theta},\hat{r}(\hat{\theta}))$$\theta$$(\hat{\theta}, \hat{r}(\hat{\theta}))$DEFINITION 2.2 In the statistical model , the expected semiparametric information is de$\mathcal{G}$fined as follows:$E_{H}[ SKL(\hat{\theta})]=E_{H}[\int\rho(y,\hat{\theta},\tilde{r}(\hat{\theta}))h(y)dy]$.

(2.10)

25
For later use, we write explicitly the components ofas follows:$\Upsilon(\theta, r)$$M(y, \theta)=\frac{\partial m(y_{)}\theta)}{\partial\theta}=(\begin{array}{ll}M_{1} M_{3}M_{2} M_{4}\end{array}),$

where

2.5

$M_{1}= \frac{\partial m_{1}}{\partial\sqrt{}’}=-\frac{\partial\mu_{i}(\sqrt{})}{\partial\sqrt{}},$$M_{2}= \frac{\partial m_{2}}{\partial\beta’}=2(\mu_{i}(\beta)-y)\frac{\partial\mu_{i}}{\partial\beta}-\phi\frac{\partial V(\mu_{i}(\sqrt{}))}{\partial\mu_{i}}\frac{\partial\mu_{i}}{\partial\beta},$$M_{3}= \frac{\partial m_{1}}{\partial\phi}=0,$$M_{4}= \frac{\partial m_{2}}{\partial\phi}=-V(\mu_{i}(\beta))$.

The proposed information criteria

DEFINITION 2.3 Let$k$be the dimension of , then we define the following information$\theta$criterion:

SIC$= \sum_{i=1}^{n}\rho(y_{i}, \hat{\theta},\hat{r}(\hat{\theta}))+k$.

(2.11)We shall call this information criterion the semiparametric information criterion, $or$ SIC
for short.

The proposed criterion SIC has the following properties.
PROPOSITION 2.1 Suppose that the model is correctly specified. Further suppose that $(C1)-$(C6) hold. Then SIC is an asymptotically unbiased estimator of the expected information.
For possibly misspecified model, we propose the following criterion for model selection.

DEFINITION 2.4 Consider a possibly misspecified semiparametric modeltion 2.1. Let$\mathcal{G}$as defined in Sec-$\hat{S}=\frac{1}{n}\sum_{i=1}^{n}[\Upsilon_{i}(\hat{\theta},\hat{r}(\hat{\theta}))\Upsilon_{i}(\hat{\theta},\hat{r}(\hat{\theta}))’], \hat{Q}=\frac{1}{n}\sum_{i=1}^{n}[\frac{\partial\Upsilon_{i}(\hat{\theta},\hat{r}(\hat{\theta}))}{\partial(\theta’,r)’}]$

Then we defined the following information criterion$SIC_{T}=\sum_{i=1}^{n}\rho(y_{i},\hat{\theta},\hat{r}(\hat{\theta}))-trace(\hat{S}\hat{Q}^{-1})$.

(2.12)The proposed criterion $SIC_{T}$ has the following properties.
PROPOSITION 2.2 Under the regularity conditions (C1) $-(C6)$ listed in Section 2.3 and
2.4, $SIC_{T}$ is an asymptotically unbiased estimator of the expected information.

can be understood as a semiparametric extension of the Takeuchi information criterion
for parametric likelihood analysis (Takeuchi, 1976).
$SIC_{T}$26

3Data Analysis

Now we illustrate the use of SIC using a well-known data set Orobanche on germination of
the seed variates, which was considered in Crowder (1978, Table 3). In this experiment,
two types of seeds from Orobanche (O.aegyptiaca75 and O.aegyptiaca73) were tested for
germination when they were given with two diﬀerent root extracts (Bean and cucumber).
The goal of the experiment was to determine the eﬀect of the root extract on inhibiting the
growth of the parasitic plant. There is suspicion of overdispersion for this data set and anumber of authors have used this data set to study methodologies extending the classical
GLM to take into consideration the problem of overdispersion. For instance, this data
was analyzed by Breslow and Clayton (1993) to illustrate the use of the generalized linear
mixed models for overdispersion. In these studies, however, the authors have not considered
the aspects of model selection. Now we shall reanalyze it by applying the model selection
technique proposed in the previous section based on the quasi-likelihood.
Now let be the proportions of germinated seeds, and the seeds in the ith data set$(i=1, \ldots, 21)$ . Suppose that $E[y_{i}]=\pi_{i}$ and $var[y_{i}]=\phi v(\pi_{i})=\phi\pi_{i}(1-\pi_{i})/n_{i}$ , whereis theoverdispersion parameter. Further, we shall suppose that logit $(\pi_{i})=\log(\pi_{i}/(1-\pi_{i}))=x_{i}’\beta,$
where $x_{i}=(1,$where $x_{1i}=1$ if the root extract is cucumber, otherwise if
the root extract is Bean; and $x_{2i}=1$ if the type of seeds is O.aegyptiaca75, otherwise if
the type of seeds is O.aegyptiaca73. We shall compare the folowing five candidate models:$n_{i}$$y_{i}$$\phi$$0$$x_{1i},$ $x_{2i},$ $x_{1i}x_{2i}$$0$(1) logit$(\pi_{i})=\beta_{0}$(2) logit$(\pi_{i})=\beta_{0}+\beta_{1}x_{1i}$(3) logit$(\pi_{i})=\beta_{0}+\beta_{2}x_{2i}$(4) logit$(\pi_{i})=\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}$(5) logit$(\pi_{i})=\beta_{0}+\beta_{1}x_{1i}+\beta_{2}x_{2i}+\beta_{3}x_{1i}x_{2i}$Table 1: Model selection criteria for the Orobanche data$\overline{\frac{Models(1)(2)(3)(4)(5)}{SIC3.8212.7325.9713.8815.521}}$$SIC_{T}$$AIC_{BB}$

0.551
133.0

1.260 0.305
120.4 133.1

1.860 3.380
119.7 117.5

We computed the values of the various information criteria and the results are summarized in Table 1. The proposed SIC has the smallest value for model (2) among all candidate

27
models. SIC ranks the models in the following way: (2) $>(1)>(4)>(5)>(3)$ . This result
is consistent with the results obtained in Crawley (2005, p.255-260). Crawley (2005) usedquasi-likelihood method to analyze the Orobanche data. He applied the -test to compare
the full model (5) with other simpler models. He found that there is no compelling evidence
that the types of seeds and the interaction should be kept in the model, and the minimal
adequate model is model (2). In Table 1, we have also shown the values of $AIC_{BB}$ , the valuesof AIC based on the beta-binomial model. $AIC_{BB}$ has minimum value for model (5), themost complicated model considered here. In Table 1 the values $SIC_{T}$ are also shown; $SIC_{T}$is a modified version of SIC which is not discussed here. $SIC_{T}$ is a corresponding version of
the Takeuchi type information criteria in full likelihood analysis.
$F$

References[1] Akaike, H. (1973). Information theory and an extension of the maximum likelihood
principle. In Second International Symposium on Information Theory, B. N., Petrov
and F., Csaki (eds), 267-281, Budapest: Akademiai Kiado.

[2] Anderson, D. R., Burnham, K. P. and White, G. C. (1994). AIC model selection in
overdispersed capture-recapture data. Ecology, 75, 1780-1793.
[3] Breslow, N E. and Clayton, D. G. (1993). Approximate inference in generalized linear
mixed models. J Am. Statist. Ass., 88, 9-25.
[4] Crawley, M. J. (2005). Statistics: an introduction using R. John Wiley & Sons, Ltd.

[5] Crowder, M. J. (1978). Beta-binomial anova for proportions. Appl. Statist., 27, 34-37.
[6] Godambe, V. P. (1999). Linear Bayes and optimal estimation. Annals of the Instituteof Statistical Mathematics, 51, 201-215.
[7] Hanfelt, J. J. and Liang, K.-Y. (1995). Approximate likelihood ratios for general estimating functions. Biometrika, 82, 461-477.

[8] Lebreton, J. D., K. P. Burnham, J. Clobert., and Anderson, D. R. (1992). Modelingsurvival and testing biological hypothesis using marked animals: a unified approachwith case studies. Ecological Monographs, 62, 67-118.
[9] Li, B. (1993). A deviance function for the quasi-likelihood method. Biometrika, 80,

741-753.
[10] Lin, P.-S. (2011). Quasi-deviance functions for spatially correlated data. Statistica
Sinica, 21, 1785-1806.
[11] McCullagh, P. and Nelder, J. A. (1989). Generalized linear models. Chapman and Hall:
London.

28[12] McLeish, D. L. and Small, C. G. (1992). A projected likelihood function for semiparametric models. Biometrika, 79, 93-102.
[13] Moore, D. F. (1986). Asymptotic properties of moment estimators for overdispersed
counts and proportions. Biometrika, 73, 283-288.
[14] Pan, W. (2001). Akaike’s information criterion in generalized estimating equations.
Biometrics, 57, 120-5.

[15] Qian, G., Gabor, G. and Gupta, R. P. (1996). Generalized linear model selection by thepredictive least quasi-deviance criterion. Biometrika, 83, 41-54.
[16] Wedderburn, R. W. M. (1974). Quasi-likelihood functions, generalized linear models,
and the Gauss-Newton method. Biometrika, 61, 439-447.