%\documentstyle[righttag,ctagsplt,11pt]{amsart}  
\documentstyle{article}
%\documentstyle[amstex,righttag,ctagsplt]{amsart}
\renewcommand{\baselinestretch}{2.4}
\normalsize

\textwidth 6.75in
\textheight 8.50in
\oddsidemargin 0.0in
\evensidemargin 0.0in
\topmargin 0.0in


\hyphenation{Bio-met-rik-a}


\newcommand\bxi{ x \sb i}
\newcommand\bbeta{\beta}
\newcommand\bhbeta{ {\hat \beta} }
\newcommand\btbeta{ {\tilde \beta} }
\newcommand\btbeti{ {\tilde \beta \sb i} }
\newcommand\btau{\tau}

\pagestyle{plain}


\begin{document}


\begin{center}
\bf \Large
Dispersion Tests and  Adjustments 
%\end{center}
%\begin{center}
%\bf \Large
for Survival and Case-Control Studies
\end{center}




\begin{center}
Philip J.~Smith
\end{center}
\begin{center}
Centers for Disease Control and Prevention,
Division of Diabetes Translation
\end{center}
\begin{center}
4770 Buford Highway, N.E.,
Mail Stop K-10,
Atlanta, GA 30341-3724
\end{center}

\begin{center}
{January 16, 1993}
\end{center}

\section*{ \bf Summary}

Cox and Snell \cite{bib:CoxSnel}  have pointed out
that when the empirical variance of a measured response $y$
is quite different from
its nominal variance under an assumed probability model,
standard errors that are not
adjusted for 
the discrepancy 
may yield misleading inferences.
This article describes a
%multivariate conditional score 
test
for assessing 
departures from nominal dispersion 
in survival analysis
and
case control studies.
Also,
a robust procedure for adjusting 
standard errors for these models is given.

\noindent
KEY WORDS:
Conditional logistic regression;
overdispersion;
Poisson regression;
proportional hazards;
underdispersion.

\medskip

\section{\bf Introduction }


In many applications interest focuses on
the relation between
parameter estimates and their standard errors.
For example, in case control studies
an estimated log-odds which exceeds
its standard error by more than a factor of two
suggests a strong association between
a disease
and a suspected risk factor.
Implicit in this  and 
other analyses 
is the assumption that the standard errors are correct.
However, as Cox and Snell \cite{bib:CoxSnel}  point out,
when overdispersion is present, unadjusted standard errors 
are too small, and the confidence of conclusions may be overstated.
On the other hand,
if underdispersion is present unadjusted standard errors
will be too large resulting in reduced power.


We present a
%multivariate score 
test for assessing 
departures from nominal dispersion
in survival analysis and
in case-control studies.
Our methods 
are based on 
%the optimal score tests
%for mixture models 
those developed by 
White \cite{bib:White},
Neyman and Scott \cite{bib:NeySco},
and Lindsay \cite{bib:Linds1,bib:Linds2,bib:Linds3}
and extended to generalized linear models by 
Smith and Heitjan \cite{bib:SmiHeit}.
We use a test statistic that measures the discrepancy between the empirical and
nominal information for each estimated regression coefficient.
Also, we give
a method for adjusting the standard errors of 
estimated parameters when 
a departure from nominal dispersion is suspected.
%This procedure is based on a robust procedure
%described by Royal \cite{bib:Roy}.
In related literature Lin and Wei \cite{bib:LW} have
investigated the use of robust standard errors for nonparametric
survival analyses.


In Section~2, we present methodology for assessing and correcting for 
departures from nominal dispersion.
In Section~3, we apply 
the methodology to a parametric proportional hazards model
in an analysis of lung cancer data discussed by Prentice \cite{bib:Prent}.

Because the likelihoods
for the proportional hazards model (not making any parametric assumptions)
and conditional logistic regression are identical, 
an assessment of a departure from nominal dispersion 
which is based on the likelihood
would be the same for either model.
In Section~4 
we show how such an assessment is obtained
using data from a case control study of women with endometrial cancer
conducted by Mack et al. \cite{bib:Mack}
and discussed by Breslow and Day \cite{bib:BresDay}.

\section{\bf Covariance Adjustments and Specification Tests  }  

Assume that 
data $y \sb 1 , \ldots , y \sb n$
are generated from the true but unknown distribution function $G,$
but that a family of distribution functions
$F(y, \beta)$ with densities $f(y,\beta)=dF(y,\beta)/dy$
is used to model the data.
Here $\beta = ( \beta \sb 1 , \ldots , \beta \sb p) \sp T$ is a vector of $p$ unknown parameters.

This situation is common in many areas of applied statistics
and corresponds to the situation in which the chosen $F$ 
misspecifies some  distributional
aspect of the data. 
%The misspecification  with which 
%we will be concerned  is departures from nominal dispersion.
For example,
if $y$ is assumed to have a Poisson distribution 
with parameter $\mu$
but the empirical variance 
$\sum \sb {i=1} \sp n ( y \sb i - \bar y ) \sp 2 / n$
greatly exceeds the estimated nominal Poisson variance
$\bar y ,$
$y$ is said to be overdispersed.
In this case the second moment of 
$y \sb 1 , \ldots , y \sb n$
is misspecified by the Poisson model.
%
Alternatively, some other aspect of the true distribution of 
$y \sb 1 , \ldots , y \sb n$ may be misspecified. 
However, the general symptom of model misspecification is
a departure from the nominal dispersion expected under the assumed model.
%
McCullagh and Nelder \cite{bib:McNel}
maintain that overdispersion ``is
the norm in practice and nominal dispersion the exception''
and advise statisticians ``to assume that overdispersion is present
to some extent unless and until it is shown to be absent.''
%To account for overdispersion
%in generalized linear models
%McCullagh and Nelder recommend inflating the nominal
%covariance matrix of estimated regression coefficients
%$(X \sp T W X) \sp {-1} $ by a factor $\tilde \sigma \sp 2 $
%which can be estimated from the data.

%In particular, if  the  variance of 
%the data
%$y \sb 1 , \ldots , y \sb n$
%enerated from the true but unknown distribution $G$ is much greater than would be expected from
%the statistical model $F,$
%the data are said to be overdispersed with respect to $F.$

Without knowing the true distribution function, 
an estimate of $\beta$ may be obtained by maximizing 
\begin{eqnarray}
\label{eq:quasiLL}
L  (\beta) =  \prod \sb {i=1} \sp n f( y \sb i , \beta)
\end{eqnarray}
with respect to $\beta.$
White  calls the resulting estimate $\hat \beta$
a ``quasi-maximum likelihood estimate'' and
Akaike \cite{bib:Aka} notes that this estimator minimizes the
Kullback-Leibler information criterion.
Appendix A 
gives a simple illustration that 
the standard errors of $\hat \beta$ are wrong 
when $y \sb 1 , \ldots , y \sb n$
%are not nominally dispersed according to the distribution function $F.$
do not
have the variance  one would expect under the distribution function $F.$

\subsection{\bf Detecting departures from nominal dispersion}
To detect departures from nominal dispersion, let
\begin{eqnarray*}
U \sb {j} (y  ,  \beta) = 
%n \sp {-1}
%\sum \sb i \sp n
\left(
\frac
{
\partial \log f(y  ,  \beta )
}
{
\partial  \beta \sb j
}
\right) \sp 2
+
\frac
{
\partial \sp 2 \log f(y  ,  \beta )
}
{
\partial  \beta \sb j \sp 2
},
\end{eqnarray*}
$d  (y  ,  \beta)= ( U \sb 1 (y  ,  \beta) , \ldots U \sb p (y  ,  \beta) ) \sp T ,$
and 
\begin{eqnarray}
\label{eq:testvec}
D  (   \beta )
=
%n  \sp {-1}
\sum \sb {i=1} \sp n
d  (y  \sb i ,   \beta).
\end{eqnarray}
For generalized linear models
%with estimated regression parameters $\hat \beta$
%and 
with covariate values $x \sb {i1}, \ldots , x \sb {ip}$ for the $i$th observation,
Smith and Heitjan show that
the $j$th element of the vector
(\ref{eq:testvec})
%$D  (  \beta )$
is 
\begin{eqnarray}
\label{eq:teststat}
D \sb { j} (  \beta )
=
%n \sp {-1}
\sum \sb {i=1} \sp n
\left[
(y \sb i - \mu \sb i) \sp 2
-
V \sb i
\right] 
x \sb {ij} \sp 2,
\end{eqnarray}
where, 
$\mu \sb i = E \sb F (y \sb i)$ 
and
$V \sb i = Var \sb F (y \sb i)$
are the ``nominal'' mean and variance under the assumed model $F.$
Therefore, the statistic (\ref{eq:teststat})
may be regarded as a function of
the difference between 
the empirical variance $(y \sb i - \mu \sb i) \sp 2$
and nominal variance $V \sb i$
of $y \sb i.$
If $y \sb i$ is overdispersed
the empirical variance will be much larger than the nominal variance
and $D \sb { j} (  \beta )$ will be  positive
and large relative to its standard error.
Alternately, if $D \sb { j} (  \beta )$
is negative and large relative to its standard error, 
the empirical variance is much smaller
than the nominal variance and the data may be said to be underdispersed.
However, if the model $F$ is correctly specified and there is no overdispersion,
the expectation of $D ( \beta)$ is zero. 

In view of these appealing properties, 
$D  ( \hat \beta ) $ may be used as a test statistic 
to detect departures from nominal dispersion. 
In particular, a consistent estimate of the asymptotic
covariance of $ D  (\hat \beta )$ is 
\begin{eqnarray*}
V  (\hat \beta ) =
%n \sp {-1} 
\sum \sb {i=1} \sp n 
\left[
d( y \sb i , \hat \beta ) 
-
\nabla D  ( \hat \beta )
A  ( \hat \beta ) \sp {-1}
\nabla  \log f   ( y\sb i ,\hat \beta )
\right]
\cdot
\left[
d( y \sb i , \hat \beta ) 
-
\nabla D  ( \hat \beta )
A  ( \hat \beta ) \sp {-1}
\nabla \log f ( y \sb i , \hat \beta )
\right]
\sp T
\end{eqnarray*}
where
$\nabla D  ( \hat \beta ) =
\{
%n \sp {-1}
\sum \sb {i=1} \sp n
\partial d( y \sb i , \hat \beta ) / \partial \hat \beta \sb k
\}
$
and
$\nabla  \log f ( y \sb i , \hat \beta ) =
\{
\partial \log f( y \sb i , \hat \beta ) / \partial \hat \beta \sb k
\}.
$
White shows that 
\begin{eqnarray}
\label{eq:testat}
\chi \sp 2 = 
%n 
D  ( \hat \beta ) \sp T V  ( \hat \beta ) \sp {-1} D  ( \hat \beta )
\end{eqnarray}
is distributed asymptotically as 
a chi--squared variate with $p$ degrees of freedom
provided that the data $y \sb 1 , \ldots , y \sb n$ are  nominally dispersed
with respect to the model $F.$


Letting
$C \sb {jj} $
denote the $j$th diagonal element in the covariance matrix $V  (\hat \beta ) ,$
to investigate   specific departures from nominal
dispersion the standardized statistics
\begin{eqnarray}
\label{eq:standard}
T \sb { j}
= 
%\sqrt n
D \sb { j} (\hat \beta) / C \sb {jj} \sp {1/2} ,
\qquad
j=1, \ldots , p,
\end{eqnarray}
may be used. 



\subsection{\bf Adjusting for departures from nominal dispersion}
Assuming that  $y \sb 1 , \ldots , y \sb n$ are distributed according
to the distribution function $F$
the nominal variance-covariance matrix of $\hat \beta$
is customarily obtained by 
evaluating the inverse of the observed information matrix
\begin{eqnarray*}
A  ( \beta) =
\{
-\sum \sb {i=1} \sp n
\partial \sp 2 \log f( y\sb i,  \beta )
/
\partial  \beta \sb j \partial  \beta \sb k
\}.\end{eqnarray*}
%
at $\beta = \hat \beta .$
However, letting
\begin{eqnarray*}
B  ( \beta)=
\{
\sum \sb {i=1} \sp n
\partial \log f( y\sb i,  \beta )
/
\partial  \beta \sb j 
\cdot
\partial  \log f( y\sb i,  \beta )
/
\partial  \beta \sb k
\},
\end{eqnarray*}
White shows that 
a consistent estimate of the asymptotic variance-covariance matrix of 
%the quasi-likelihood estimate 
$\hat \beta$
is
\begin{eqnarray}
\label{eq:nvaf}
C  ( \hat \beta ) = 
A  (\hat \beta) \sp {-1} 
B  (\hat \beta) 
A  (\hat \beta) \sp {-1} 
\end{eqnarray}
%evaluated at $\hat \beta $
whether or not the data
$ y \sb 1 , \ldots , y \sb n $
have the variance that one would expect from the distribution function $F.$
Royal \cite{bib:Roy}
refers to the the resulting variance estimates
as being ``robust'' since they compensate for model misspecification
of this type.
Letting $\mu \sb i$ and $V \sb i$
denote the mean and nominal variance of 
$y \sb i$ under $F,$ Smith and Heitjan show that 
$A  (\hat \beta) \sp {-1} 
B  (\hat \beta) $
may be regarded as a multivariate version of
the variance adjustment factor 
\begin{eqnarray*}
%\label{eq:stdscore}
\tilde \sigma \sp 2 = \sum \sb {i=1} \sp n
\frac
{
( y \sb i - \mu \sb i) \sp 2 
}
{V \sb i}
/ n-p
\end{eqnarray*}
which is used in
generalized linear models to
adjust the nominal variance-covariance matrix
of estimated regression 
coefficients  
%$A ( \hat \beta )\sp {-1} = 
%(X \sp T W X) \sp {-1} $
for departures from nominal dispersion.
In particular,
Smith and Heitjan show that
$\tilde \sigma \sp 2$ and its multivariate analog 
$A  (\hat \beta) \sp {-1} B  (\hat \beta) $
may be interpreted as the average ratio of empirical variance  $( y \sb i - \mu \sb i) \sp 2  $
to nominal variance  $V \sb i .$ 
By premultiplying  
 $A  (\hat \beta) \sp {-1}$
by $A  (\hat \beta) \sp {-1} B  (\hat \beta) ,$
each element of the nominal variance-covariance matrix
$A  (\hat \beta) \sp {-1}$
is adjusted differentially
for departures from nominal dispersion. 










\section{\bf A Parametric Survival Analysis}

Aitkin and Clayton \cite{bib:AitCla} have described how complex censored
survival data may be fitted within the framework of generalized linear models.
Let $t \sb i , \ldots , t \sb n$
be the survival times of $n$ individuals, some of which are right
censored and let $\bxi $ be the $i$th individual's corresponding explanatory
variables. Assuming the hazard function $h(t)$ to be 
associated with the explanatory
variables through a log-linear (proportional hazards) model
\begin{eqnarray}
\label{eq:prophaz}
h(t \sb i) = 
\lambda ( t \sb i )
\exp
\left(
\bbeta \sp T \bxi
\right)
\end{eqnarray}
and letting 
$y \sb i$ be an indicator variable taking the value 1 for uncensored
and 0 for censored observations, 
the likelihood function (\ref{eq:quasiLL}) becomes
\begin{eqnarray}
\label{eq:poislik0}
L(\beta) =
\prod \sb {i=1} \sp n
\left[
\mu \sb i \sp {y \sb i}  
e \sp {- \mu \sb i}
\right]
\left[
\lambda \sb i 
/
\Lambda \sb i
\right]
\sp {y \sb i}
\end{eqnarray}
where 
$\mu \sb i = \Lambda (t \sb i) \exp ( \bbeta \sp T \bxi )$
and
$\Lambda (t \sb i) = \int \sb {- \infty} \sp {t \sb i} 
\lambda ( u ) \, du$
is the cumulative hazard function.
%The first term in the likelihood is the kernel of the likelihood of $n$
%independent Poisson variates  $y \sb i $ with mean $\mu \sb i .$
%
Assuming
%$t \ge 0$ 
survival times to be exponentially distributed,
$\Lambda (t \sb i) = t$ and
%$\lambda \sb i / \Lambda \sb i = t \sp {-1}$ is free of unknown 
%parameters\footnote 
%{
%This assumption is equivalent to presuming that
%the survival times are exponentially distributed with mean
%$ \exp ( \bbeta \sp T \bxi ) .$
%For methods for maximizing 
%(\ref{eq:poislik0}) under alternative parametric distributional assumptions
%see Aitkin and Clayton
%\cite{bib:AitCla}.
%}
the likelihood function (\ref{eq:poislik0}) is
\begin{eqnarray}
\label{eq:poislik}
L(\beta) \propto
\prod \sb {i=1} \sp n
\mu \sb i \sp {y \sb i}
e \sp {- \mu \sb i}
\end{eqnarray}
\noindent
%where $\mu \sb i$ is a function of $\beta .$
which is the product of Poisson densities.
Maximizing the likelihood (\ref{eq:poislik}) with respect to $\beta$
is equvalent to fitting  a generalized linear model to
the Poisson random component 
$y \sb i$
with a logarithmic link and 
linear predictor
%\begin{eqnarray*}
$
\ln \mu \sb i = \ln t \sb i  + \bbeta \sp T \bxi 
$
%end{eqnarray*}
where 
$\ln t \sb i$ 
is an offset.
%\noindent
In terms of the notation in Section~2,
independent Poisson kernels
$f (y \sb i , \beta ) = 
	\mu \sb i \sp {y \sb i}
	e \sp {- \mu \sb i}
%	/ y\sb i !
$
are used to model the data
$ y \sb 1 , \ldots , y \sb n$
and
the nominal mean and variance of $y \sb i$ 
under this model is $\mu \sb i .$
The test statistics (\ref{eq:teststat})
for assessing departures from nominal dispersion are
\begin{eqnarray}
\label{eq:survtest}
D \sb {j}(\hat \beta) =
%\frac 1 n
\sum \sb {i=1} \sp {n}
\left\{
\left(y \sb {i} - \mu \sb {i} \right) \sp 2 
-
\mu \sb {i} 
\right\}x \sb {ij} \sp 2
\end{eqnarray}
and are functions of the difference between the
empirical variance
$(y \sb {i} - \mu \sb {i} ) \sp 2$
and
nominal variance
$\mu \sb i $ under the Poisson model.
Algebraic expressions for 
$A(\hat \beta),$
$B(\hat \beta),$
$\nabla D(\hat \beta),$
and
$\nabla \log f (y \sb i , \hat \beta ),$
that are 
needed to test and adjust for departures from
nominal dispersion are given in Appendix B.



Because (\ref{eq:survtest}) is a function of
the difference between a measure of  empirical variance 
$\left(y \sb {i} - \mu \sb {i} \right) \sp 2$
and the nominal variance $\mu \sb i$ for a Poisson variate,
if $D \sb { j}(\hat \beta) $ is positive and large relative to its standard error,
overdispersion is indicated.
However, 
if $D \sb { j}(\hat \beta)$ is negative and large relative to its standard error, 
underdispersion is indicated.


To illustrate these methods,
we consider an example of medical survival data adapted from Prentice 
\cite{bib:Prent}. 
The data consist of survival times, T, in days of 137 lung cancer patients from a
Veterans' Administration Lung Cancer trial.
The data set also includes the following explanatory variables: 
performance STATUS, 
an ordinal variable describing ability of the individual to provide self-care;
AGE in years;
MFD, the time in months from diagnosis to starting on the study;
an indicator of PRIOR therapy;
an indicator of standard or test TREATment;
TYPE of tumor;
and a censoring indicator, Y.  

Assuming that the proportional hazards model (\ref{eq:prophaz}) is valid,
Aitkin, Anderson, Francis, and Hinde (\cite{bib:Aitetal} pp.300--302.)
conducted an analysis of these data
using GLIM \cite{bib:GLIM}
and concluded  that the  assumption of exponential survival times is
appropriate compared to a Weibull alternative. This assessment agrees
with that of Kalbfleish and Prentice 
(\cite{bib:KalPre} p.62)
who also provide a graphical justification for the exponential assumption.

Both 
Aitkin {\sl et al.}
and
Kalbfleish and Prentice
conclude 
that the test treatment appears to have no significantly different
effect from the standard treatment. 
Aitkin {\sl et al.}
show that the best fitting Poisson regression 
(having a deviance of 132.24 with 130 degrees of freedom)
suggests 
that the hazard declines with increasing performance STATUS, but this decline is
more rapid for those patients with prior therapy,
as signified by a significant
interaction between STATUS and PRIOR therapy. Also, TYPE 2 and 3 tumors have
significantly higher hazard than TYPEs 1 and 4. 
These results are concordant with those of
Kalbfleish and Prentice
who conducted an analysis stratified by PRIOR therapy status.

To test for a departure from nominal dispersion, 
the  statistic (\ref{eq:testat}) yields
$\chi \sp 2 = \mbox{16.04}$ with 7 degrees of freedom
(6 d.f. correspond to the 6 explanatory variables listed above and 1 d.f. corresponds to the ``reference cell'' intercept term)
corresponding to a p-value of $p=\mbox{0.024}.$
Thus, there is  information supporting a departure from nominal
dispersion and that standard errors of estimated regression coefficients may be
misleading.
%
Table \ref{tab:surv} lists the estimated regression coefficients\footnote
{
Results are in GLIM parameterization. 
For example, TYPE(3) is the estimated {\sl difference} of the effect of
type at level 3 versus 1 and its standard error is that of the difference. }
of the hazard
function along with their nominal standard errors, robustly estimated standard
errors (\ref{eq:nvaf}), and standardized  statistics (\ref{eq:standard}).

\begin{table}[b]
\renewcommand{\baselinestretch}{1}
\normalsize
\caption
{\label{tab:surv}
Dispersion statistics for the lung cancer data
}
\begin{center}
\begin{tabular}
{ l c c c c  }
parameter & MLE & nominal s.e. & robust s.e. & standardized statistics \\
\hline 
intercept    &-3.75 & 0.39 & 0.38 & 0.01 \\
STAT         &-0.24 & 0.06 & 0.05 & -0.71 \\
PRIOR(2)     &~1.78 & 0.65 & 0.60 & -0.24 \\
STAT.PRIOR(2)&-0.29 & 0.11 & 0.09 & -0.30 \\
TYPE(2)      &~0.65 & 0.24 & 0.27 &  0.46 \\
TYPE(3)      &~1.04 & 0.27 & 0.22 & -2.84 \\
TYPE(4)      &~0.26 & 0.27 & 0.23 & -1.40 \\
\hline
\end{tabular}
\end{center}
\end{table}
%\noindent
Standardized statistics in this table 
indicate that the departure from nominal dispersion
suggested by the  statistic \ref{eq:testat} is due to the TYPE variable:
the standardized statistics
is large and negative
for TYPE(3)
suggesting underdispersion. 
In this case, nominal standard errors are expected to be 
too large
and the robustly estimated standard error for
this level adjusts its value downwards appropriately from 0.27 to 0.23.

One referee wondered whether the  departure from nominal
dispersion in the lung cancer data 
could be attributed to a specific source of model misspecification
such as a departure from the proportional hazards assumption.   
%
Figure~1 gives a plot of Schoenfeld \cite{bib:Sch} residuals versus
survival time along with a lowess smooth \cite{bib:Cle}.
If the proportional hazards assumption is valid then
the plot of the residuals versus survival time should show no apparent
pattern with increasing time and be centered around zero. 
However, as suggested by the lowess smoother
the Schoenfeld residuals may increase with survival
time between zero and approximately 200 days.
Harrell's \cite{bib:Har}  test for 
nonproportional odds yields a p--value of 0.085 giving moderate
support  for the null hypothesis that the proportional
odds assumption is correct.

\section{\bf Case Control Studies}

Breslow and Day \cite{bib:BresDay}
and Pregibon \cite{bib:Preg}
describe the use of conditional logistic regression
in case control studies
in which one case is matched to several controls.


Let
$x \sb {0k} = ( x \sb {10k}, \ldots , x \sb {p0k} ) \sp T$ 
denote the exposure covariates of the $k$th case
and 
$x \sb {rk} = ( x \sb {1rk}, \ldots , x \sb {prk} ) \sp T $ 
denote the 
exposure covariates of the $r$th control
in the $k$th matched set,
$r=1, \ldots R \sb k ,$  
$k=1, \ldots , K.$ 
Also, let
$y \sb {rk}$ be an indicator variable equal to 1 if $r=0$ to denote
the case in the $k$th matched set
and equal to 0 if $r=1 , \ldots, R \sb k$ to denote the
controls.
Liddell, McDonald and Thomas \cite{bib:LMT}
show the likelihood of the conditional logistic model is
\begin{eqnarray}
L( \beta ) &=&
\prod \sb {k=1} \sp K
\prod \sb {r=1} \sp{R \sb k}
\left[
\frac
{
\exp ( \bbeta \sp T x \sb {rk} )
}
{
\sum \sb {r=0} \sp {R \sb k}
\exp ( \bbeta \sp T x \sb {rk} )
}
\right]
\sp{y \sb {rk}}
\nonumber
\\
&=&
\prod \sb {k=1} \sp K
\frac
{
\exp ( \bbeta \sp T x \sb {0k} )
}
{
\sum \sb {r=0} \sp {R \sb k}
\exp ( \bbeta \sp T x \sb {rk} )
}.
\label{eq:condlik}
\end{eqnarray}
%Whitehead \cite{bib:Whitcc} has shown that this likelihood
%is identical to that of a Poisson for variables
%$y \sb {ij}$
%with mean 
%$\mu \sb {ij} = \exp ( \alpha \sb i + \bbeta \sp T x \sb {ij} )$
%and with realized values
%$y \sb {i0} = 1$
%for cases
%and
%$y \sb {ij} = 0$ 
%for controls.
\noindent
The expression (\ref{eq:condlik})
is identical to the partial likelihood of the Cox
proportional hazards model. In this regard tests and adjustments for
departures from nominal dispersion based on the (\ref{eq:condlik})
are the same for case--control
studies and survival analyses using the nonparametric Cox model. 
%Lin and Wei
%\cite{bib:LW} have discussed the use of robust standard errors for
%regression coefficients in this context but do not give a test for a
%eparture from nominal dispersion.

In terms of the notation in Section 2,
the $j$th element of the 
vector of test statistics (\ref{eq:testvec}) is
\begin{eqnarray}
\label{eq:survscore}
D \sb { j} ( \beta) =
% \frac 1 n
\sum \sb {k=1} \sp K
\left\{
\left(x \sb { j0k} - \mu \sb {j \cdot k} \right) \sp 2 
-
\sum \sb { c=0} \sp {R \sb k}
[
x \sb { jck} - \mu \sb {j \cdot k}
] \sp 2
\pi \sb {ck}
\right\},
\end{eqnarray}
$j=1, \ldots, p,$
where
$\pi \sb {ck}
=
\exp ( x \sb {ck} \sp T  \beta  )
/
\sum \sb {c=0} \sp {R \sb k}
\exp ( x \sb {ck} \sp T  \beta  )
$
is the conditional probability of risk factors  $ x \sb {ck}$ 
given that case $c$ in matched set $k$ has disease,
and
$\mu \sb {j \cdot k}
=
\sum \sb {c=0} \sp {R \sb k}
x \sb {jck}
\pi \sb {ck}
$
is the expected value of the $j$th risk factor 
for the case
in  matched set $k.$
It is noteworthy that (\ref{eq:survscore}) is a function of
the difference between the empirical variance 
$\left(x \sb { j0k} - \mu \sb {j \cdot k} \right) \sp 2$
and the nominal variance $\sum \sb { c=0} \sp {R \sb k}
[
x \sb { jck} - \mu \sb {j \cdot k}
] \sp 2
\pi \sb {jk}
$
of the $j$th risk factor
for the case in matched set $k.$
Thus, if $D \sb { j}(\hat \beta)$ 
is positive and large relative to its standard error,
overdispersion is indicated.
Algebraic expressions for the statistics 
that are
needed to test and adjust for departures from
nominal dispersion are given in Appendix C.



To illustrate these methods,
we consider data 
from the Los Angeles study of endometrial cancer
\cite{bib:Mack}.
The purpose of the study was to obtain
an overall relative
risk for estrogens,
to detect a possible interaction with age for the risk associated with
gall--bladder disease, 
and to examine the joint effects of gall bladder disease and hypertension.
Breslow and Day \cite{bib:BresDay}
present a comprehensive analysis of this data that includes
63 women with endometrial cancer
each of whom is
matched to 4 controls.
In our analysis 
we assess the impact of possible departures from nominal dispersion
on the model investigating how the risk of
endometrial cancer depends upon
history of estrogen use (EST)
and gall bladder disease (GALL).


Table \ref{tab:case} lists the estimated coefficients
of the conditional logistic regression model 
using main effects EST and GALL
along with 
their nominal standard errors, 
robustly estimated standard errors (\ref{eq:nvaf}), 
and standardized statistics (\ref{eq:standard}).
The standardized statistics in this table 
indicate small departures from nominal dispersion:
the standardized statistic for EST is small and negative
indicating  underdispersion
and its robust standard error is adjusted downwards. 
Also, 
the standardized statistic for GALL is small and positive
indicating  overdispersion
and its robust standard error is adjusted upwards. 
The minor degree of adjustment is confirmed by
the test for a departure from nominal dispersion, 
the  statistic (\ref{eq:testat}) yields
$\chi \sp 2 = \mbox{1.74}$ with 2 degrees of freedom,
corresponding to a p-value of $p=0.41$.
Thus, there is little evidence indicating a departure from nominal
dispersion for the model.


\begin{table}[ht]
\renewcommand{\baselinestretch}{1}
\normalsize
\caption
{\label{tab:case}
Dispersion statistics for the endometrial cancer data
}
\begin{center}
\begin{tabular}
{ l c c c c  }
parameter & MLE & nominal s.e. & robust s.e. & standardized statistics \\
\hline 
EST           & 2.11 & 0.43 & 0.39 & -0.35  \\
GALL          & 1.27 & 0.41 & 0.44 &  1.01 \\
\hline
\end{tabular}
\end{center}
\end{table}
%\noindent

\section*{\bf Acknowledgments}

Thanks are given to Frank Harrell and Terry Therneau for
recommending methods for assessing the proportional hazards
assumption and Mitchell Gail for bringing the work of Lin and Wei to the author's attention. Also, the author thanks Ross Brechner and two
anonymous referees for valuable comments that aided in improving the
manuscript.

\section*{\bf Appendix A}

Let $y = (y \sb 1 , \ldots , y \sb n) \sp T$ be independent random variables 
generated from the true but unknown distribution function $G$
and let $V \sb G ( y ) = \mbox{diag}\{ v \sb i  \sp {\prime} \}$
denote the true  variance-covariance matrix of $y.$
Also, let $V \sb F( y ) = \mbox{diag}\{ v \sb i   \}$
denote the nominal variance-covariance matrix of $y$ under the model $F.$
Letting $\hat \beta \sb j = h \sb j (y)$ denote the estimate of $\beta \sb j$
obtained by maximizing
(\ref{eq:quasiLL}),
a Taylor's series approximation yields
\begin{eqnarray*}
V \sb F (\hat \beta \sb j)
\simeq
\sum \sb {i=1} \sp n
\left( \frac {\partial h \sb j (y) }{ \partial y \sb i} \right) \sp 2 
v \sb i
\end{eqnarray*}
as
the variance  of  $\hat \beta \sb j = h\sb j (y)$ under the nominal model $F.$
If the data $y$ are overdispersed then 
$v \sb i   < v \sb i \sp {\prime} ,$
for $i=1 , \ldots , n$
and
\begin{eqnarray*}
V \sb F ( \hat \beta \sb j )
\simeq
\sum \sb {i=1} \sp n
\left( \frac {\partial h \sb j (y) }{  \partial y \sb i} \right) \sp 2 
v \sb i
<
\sum \sb {i=1} \sp n
\left( \frac {\partial h \sb j (y) }{ \partial y \sb i} \right) \sp 2 
v \sb i \sp {\prime}
\simeq
V \sb G ( \hat \beta \sb j )
\end{eqnarray*}
for
$j=1 , \ldots , p.$ 


Therefore, if the data are overdispersed
the standard errors of $\{\hat \beta \sb j \}$
will be too small.


\section*{\bf Appendix B}

For the parametric survival model in Section 3,
letting $\mu \sb i = t \sb i \exp ( \beta \sp T x \sb i ),$
\begin{eqnarray*}
A ( \hat \beta ) =
\{
\sum \sb {i=1} \sp n
\mu \sb i x \sb {ij} x \sb {ik}
\},
\end{eqnarray*}
\begin{eqnarray*}
B ( \hat \beta ) =
\{
\sum \sb {i=1} \sp n
( y \sb i - \mu \sb i ) \sp 2 x \sb {ij} x \sb {ik}
\},
\end{eqnarray*}
\begin{eqnarray*}
\nabla D ( \hat \beta ) =
\{
\sum \sb {i=1} \sp n
-
\left[2( y \sb i - \mu \sb i ) - 1 \right] \mu \sb i x \sb {ij} \sp 2 x \sb {ik}
\},
\end{eqnarray*}
\noindent and
\begin{eqnarray*}
\nabla \log f ( y \sb i , \hat \beta ) =
\{
( y \sb i - \mu \sb i ) x \sb {ik} \}.
\end{eqnarray*}


\section*{\bf Appendix C}

For the conditional logistic regression model in Section 4,
\begin{eqnarray*}
A ( \hat \beta ) =
\{
\sum \sb {k=1} \sp K
\sum \sb { p=0} \sp {R \sb k}
[
x \sb { vpk} - \mu \sb {v \cdot k}
] 
[
x \sb { wpk} - \mu \sb {w \cdot k}
] 
\pi \sb {pk}
\},
\end{eqnarray*}
\begin{eqnarray*}
B ( \hat \beta ) =
\{
\sum \sb {k=1} \sp K
[
x \sb { v0k} - \mu \sb {v \cdot k}
] 
[
x \sb { w0k} - \mu \sb {w \cdot k}
] 
\},
\end{eqnarray*}
\begin{eqnarray*}
\nabla D (\hat \beta )
&=&
-2 \sum \sb {k=1} \sp K
(x \sb { v0k} - \mu \sb {v \cdot k})
\sum \sb { p=0} \sp {R \sb k}
[
x \sb { vpk} - \mu \sb {v \cdot k}
] 
[
x \sb { wpk} - \mu \sb {w \cdot k}
] 
\pi \sb {pk}
\\
\qquad& -&
 \sum \sb {k=1} \sp K
\sum \sb { p=0} \sp {R \sb k}
[
x \sb { vpk} - \mu \sb {v \cdot k}
] \sp 2
[
x \sb { wpk} - \mu \sb {w \cdot k}
] 
\pi \sb {pk},
\end{eqnarray*}
\noindent and
\begin{eqnarray*}
\nabla f ( y \sb i , \hat \beta ) =
\{
x \sb { v0k} - \mu \sb {v \cdot k}
\}.
\end{eqnarray*}






\begin{thebibliography}{99}



\bibitem{bib:CoxSnel}
Cox, D.R. and Snell, E.J.
{\sl Analysis of Binary Data}, second edition.
Chapman and Hall, New York.
1989.

\bibitem{bib:White}
White, H. 
Maximum likelihood estimation of misspecified models.
{\sl Econometrika,}
{\bf 50}, 1-25. 1982.


\bibitem{bib:NeySco}
Neyman, J. and Scott, E.
On the use of $C(\alpha)$ optimal tests of composite hypotheses.
{\sl Bulletin of the Institute of International Statistics,}
{\bf 41}, 477--497. 1966.

\bibitem{bib:Linds1}
Lindsay, B.G. 
Nuisance parameters, mixture models, and the efficiency
of partial likelihood estimators.
{\sl Philosophical Transactions of the Royal Society
of London,}
{\bf 296}(1427), 639-665. 1980.

\bibitem{bib:Linds2}
Lindsay, B.G. 
Conditional score function: some optimality results.
{\sl Biometrika,}
{\bf 69}(3), 503-512. 1982.

\bibitem{bib:Linds3}
Lindsay, B.G. 
Efficiency of the conditional score in a mixture setting.
{\sl The Annals of Statistics,}
{\bf 11}(2), 486-497. 1983.


\bibitem{bib:SmiHeit}
Smith, P.J. and Heitjan, D.F.
Testing and adjusting for overdispersion in generalized linear models.
{\sl Journal of the Royal Statistical Society, Series C: Applied Statistics.}
{\bf 42}, 31-42.

\bibitem{bib:LW}
Lin, D.Y. and Wei, L.J.
The robust inference for the Cox proportional hazards model.
{\sl Journal of the American Statistical Association,}
{\bf 84}(408), 1074--1078. 1989.



\bibitem{bib:Prent}
Prentice, R.L.
Exponential survivals with censoring and explanatory variables.
{\sl Biometrika,}
{\bf 60}, 279-288. 1973.

\bibitem{bib:Mack}
Mack, T.M., Pike, M.C., Henderson, B.E., Pfeffer, R.I., Gerkins, V.R., Arthur,
B.S., and Brown, S.E.
Estrogens and endometrial cancer in a retirement community.
{\sl New England Journal of Medicine.}
{\bf 294}, 1262-1267. 1976.

\bibitem{bib:BresDay}
Breslow, N.E. and Day, N.E.
{\sl Statistical Methods in Cancer Research.}
Oxford University Press. Cary, NC. 1980.

\bibitem{bib:McNel}
McCullagh, P. and Nelder, J.A.
{\sl Generalized Linear Models.}
Chapman and Hall: New York. 1989.

\bibitem{bib:Aka}
Akaike, H. 
Information theory and an extension of the likelihood theory
in
{\sl Proceedings of the Second International Symposium of Information Theory.}
B.N. Petrov and F. Csaki,
editors.
Budapest: Akademiai Kiado, 1973.

\bibitem{bib:Roy}
Royall, R.R. 
Model robust confidence intervals using maximum likelihood estimators.
{\sl International Statistical Review,}
{\bf 54}, 221--226. 1986.





\bibitem{bib:AitCla}
Aitkin, M.A. and Clayton, D.
The fitting of exponential, Weibull and extreme value distributions to complex
censored survival data using GLIM.
{\sl Journal of the Royal Statistical Society, Series C: Applied Statistics.}
{\bf 29}(2), 156-163. (1980).

\bibitem{bib:Aitetal}
Aitkin, M., Anderson, D., Francis, B., and Hinde, J. 
{\sl Statistical Modeling in GLIM.}
New York: Oxford Science Publications. 1989.



\bibitem{bib:GLIM}
Payne, C.D. 
{\sl The GLIM System Release 3.77 Manual.}
Oxford: Numerical Algorithms Group. 1987.

\bibitem{bib:KalPre}
Kalbfleish, J.D. and Prentice. RlL. (1980).
{\sl The Analysis of Failure Time Data.}
New York: John Wiley and Sons.

\bibitem{bib:Sch}
Schoenfeld, D. (1982).
Partial residuals for the proportional hazards model.
{\sl Biometrika,}
{\bf 69}, 239--241.

\bibitem{bib:Cle}
Cleveland, W.S. (1979).
Robustly local weighted regression and smoothing scatterplots.
{\sl Journal of the American Statistical Association,}
{\bf 74}, 829--836.

\bibitem{bib:Har} 
Harrell, F. Jr. (1986).
The PHGLM procedure,
In {\sl SUGI Supplemental Library User's Guide, Version 5 Edition}, pp
437-66.


\bibitem{bib:Preg}
Pregibon, D. 
Data analytic methods for matched case-control studies.
{\sl Biometrics,}
{\bf 40}, 639-651. 1984.

\bibitem{bib:LMT}
Liddell, F.D.K., McDonald, J.C., an Thomas, D.C.
Methods of cohort analysis: appraisal by application to asbestos mining.
{\sl Journal of the Royal Statistical Society, Series A,}
{\bf 140}, 469-491. 1977.



\end{thebibliography}




\end{document}
L( \beta ) &=&
\prod \sb {k=1} \sp K
\prod \sb {r=1} \sp{R \sb k}
\left[
\frac
{
\exp ( \bbeta \sp T x \sb {rk} )
}
{
\sum \sb {r=0} \sp {R \sb k}
\exp ( \bbeta \sp T x \sb {rk} )
}
\right]
\sp{y \sb {rk}}
\\
&=&
\prod \sb {k=1} \sp K
\frac
{
\exp ( \bbeta \sp T x \sb {0k} )
}
{
\sum \sb {r=0} \sp {R \sb k}
\exp ( \bbeta \sp T x \sb {rk} )
}.
\end{eqnarray*}








It is interesting to note that after deleting the 2 most extreme cases
corresponding to death times $t=991$ and $t=999$ from the data set the
test statistic (\ref{eq:testat}) summarizing evidence of a departure
from nominal dispersion is $\chi \sp 2 = 14.306$ on 7 degrees with
freedom corresponding to a p--value of 0.046. 
In this case the strength of evidence against the null hypothesis of
nominal dispersion is nearly the same as when the two extreme
death times are included (in this case $p=0.053$).
However, with the 2 extreme death times deleted from the analysis 
Harrell's test yields a p--value of 0.230 (compared to 0.085 with the 2
extreme cases included in the analysis) indicating that the null
hypothesis that the proportional hazards model is likely to be
valid.
Figure~2 gives a plot of the Schoenfeld residuals versus death time
along with a lowess smooth. This figure gives graphical evidence that
is concordant with the conclusion of Harrell's test.

In view of these results it seems that  with regard to the VA lung
cancer data, evidence supporting a departure from the proportional
hazards assumption may be sustained primarily by the information
provided by the 2 extreme death times. With these 2 observations
deleted, however, evidence supporting a departure from nominal
dispersion is as strong as when the 2 cases are included in the
analysis. This suggests that model
misspecification attributable to a departure from the proportional
hazards assumption may not be the cause of a possible departure from
nominal dispersion in the VA lung cancer data. However, this is not
to imply that a departure from the proportional hazards model
cannot cause a departure from nominal dispersion in other data.


remove(ls())
q()

#### statistics for 2 factor with interaction model ####
> $stats:
              [,1]      [,2]      [,3]       [,4] 
     est  2.699556 0.6116253 0.6360016 -0.3664990
    gall  2.893759 0.8829538 0.9041171 -0.8311788
est:gall -2.052167 0.9948858 0.9876457 -0.8138596

$chi:
        [,1] 
[1,] 1.07793

$pval:
[1] 0.7824043

$p:
[1] 3
