\documentclass[coll,linenumbers]{imsart}

\RequirePackage[OT1]{fontenc}
%\usepackage{natbib}
\usepackage{amsthm,amsmath,natbib}
\RequirePackage[colorlinks,citecolor=blue,urlcolor=blue]{hyperref}
\RequirePackage{hypernat}
\usepackage{epsfig, amsmath, amssymb, amscd}
\usepackage{color, graphics}
\usepackage{amsfonts,amsthm}
\RequirePackage{bm}
%\RequirePackage{cmmib57}

% settings
\pubyear{2009}
\volume{0}
%\volumetitle{IMS Collections}
\firstpage{66}
\lastpage{92}
\doi{0003/0000.0000}
\arxiv{math.PR/0000003}


\startlocaldefs
\numberwithin{equation}{section}
\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\endlocaldefs

\newcommand{\undertilde}[1]{\mathord{\vtop{\ialign{##\crcr
   $\hfil\displaystyle{#1}\hfil$\crcr\noalign{\kern1.5pt\nointerlineskip}
   $\hfil\tilde{}\hfil$\crcr\noalign{\kern1.5pt}}}}}

\theoremstyle{definition}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}{Remark}[section]
\numberwithin{equation}{section}

\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\T}{\overline{T}}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\op}[1]{\mathbf{#1}}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\abs}[1]{\left| #1 \right|}

\begin{document}
\date{empty}
\begin{frontmatter}
\title{Semiparametric Models and Likelihood - The Power of Ranks}%\protect\thanksref{T1}}
\runtitle{Semiparametric Models and Likelihood - The Power of Ranks}
%\thankstext{T1}{Prepared with partial support of N.S.F. Grant G.P. 2593.}

\begin{aug}
%\thanksref{a,t2}
\author{\fnms{Kjell} \snm{Doksum}\thanksref{a1,t1}\ead[label=e1]{doksum@stat.wisc.edu}}
\and \author{\fnms{Akichika} \snm{Ozeki}\thanksref{a2,t1}\ead[label=e2]{}}
\address[a1]{Department of Statistics, University of Wisconsin, Madison, WI 53706, email:
\url{doksum@stat.wisc.edu}}
\thankstext{t1}{Supported in part by NSF grant DMS-0505651.}
\address[a2]{Department of Statistics, University of Wisconsin, Madison, WI 53706,
email:
\url{ozeki@stat.wisc.edu}}
%\thankstext{t2}{Supported in part by NSF grant DMS-0505651.} 


\runauthor{Doksum and Ozeki}


\affiliation{University of Wisconsin, Madison}

%\end{center}

\contributor{K. Doksum}{University of Wisconsin, Madison}
\contributor{A. Ozeki}{University of Wisconsin, Madison}


\end{aug}

\begin{abstract}
We consider classes of models related to those introduced by \linebreak Lehmann in 1953 and Sklar in 1959.  Recently developed algorithms for finding profile NP likelihood procedures are discussed, extended and implemented for such models by combining them with the MM algorithm. In particular we consider statistical procedures for a regression model with proportional expected hazard rates, and for transformation models including the normal copula. A variety of likelihoods introduced to deal with semiparametric models are considered.  They all generate rank results,
not only tests, but also estimates, confidence regions, and optimality theory,
thereby, to paraphrase Lehmann (1953),
demonstrating ``the power of ranks''.
\end{abstract}

\begin{keyword}[class=AMS]
\kwd[Primary ]{62G05, 62G20}
%\kwd{60K35}
\kwd[; secondary ]{62N02}
\end{keyword}

\begin{keyword}
\kwd{Lehmann model}
\kwd{proportional hazard model}
\kwd{profile NP likelihood}
\kwd{nonparametric maximum likelihood}
\kwd{MM algorithm}
\kwd{copula models}
\kwd{Box-Cox models}
\end{keyword}

\tableofcontents



\end{frontmatter}

\section {Introduction}

We will focus on statistical inference for models where the distribution of the data can be expressed as a parametric function of unknown distribution functions.

\subsection{Lehmann Type Models. Cox Regression}
Suppose T is a random variable with a continuous distribution function F. For testing the null hypothesis $H_0 : F = F_0 $, 
Lehmann (1953) considered alternatives of the form 
\begin{equation}
\label{eq1.1}
F_{\theta}(\cdot)=C_{\theta}(F_0(\cdot)), 
\end{equation}
for some continuous distribution $C_{\theta}(\cdot)$ on 
$[0,1]$, 
which is known except for the parameter
$\theta$. We consider the problem of estimating $\theta$ when $F_0(\cdot)$ is an unknown baseline distribution.
In this case, if $T_1, \cdots, T_n$ are independent with $T_i \sim C_{\theta_i}(F_0(\cdot))$ and we set $U_i = F_0(T_i)$, then $U_i$ has distribution $C_{\theta_i}(\cdot)$ Moreover $R_i\equiv Rank(T_i) = Rank(U_i)$, which implies that the distribution of any statistical method based on $R_1, \cdots, R_n$ will not depend on $F_0$. 

For regression experiments with observations $(T_i, \bm{x}_i), i=1, \cdots , n$, where $T_i$ is a response and $\bm{x}_i$ is a vector of nonrandom 
covariates, Cox (1972) considered the parametrization  $\theta_i=g(\bm{\beta}^T\bm{x}_i)$ with $g(\cdot)$ a known function and $\bm{\beta}$ a vector of regression coefficients. He considered statistical inference procedures based on the Cox (1972, 1975) partial likelihood in very general frameworks.  These procedures are based on generalized ranks and show how powerful ranks are in generating statistical inference procedures.

In this paper we consider a special case of (\ref{eq1.1}) obtained from the Lehmann models 
$[F_0(t)]^{N}$ and
$1-[1-F_0(t)]^{N}$ by letting N be a zero truncated Poisson variable whose parameter depends on covariates and regression coefficients.
%%% INSERT to page 1  %%%%%%%%%%%%%%%%%%
We call this model ``SINAMI'' after 
SIbuya (1968) and NAbeya and MIura (1972).
For a subset of the parameter space, the model has proportional expected hazard rate (PEHR).  We show that semiparametric likelihood methods for the SINAMI model give more weight to intermediate survival times than the Cox proportional hazard model which heavily weights long survival times. Recently developed algorithms for finding profile nonparametric maximum likelihood estimates (profile NPMLE's) are combined with the MM algorithm to produce estimates.  In the two sample case, we carry out a
Monte Carlo comparison of the NPMLE with a parametric MLE and a method of moment (MOM) estimate of the two sample parameter.  The NPMLE is nearly unbiased but only about $70\%$ as efficient in terms of root MSE as the parametric estimate if the parametric model is true.  The MOM estimate is slightly less efficient than the NPMLE.

\subsection{Sklar Type Models. Copula Regression}
Suppose X and Y are random variables with continuous joint distribution $H(\cdot,\cdot)$ and marginals $F_1(\cdot)$ and $F_2(\cdot)$. 
Sklar (1959) considered models that include models of the form
\begin{equation}
\label{eq1.2}
H_{\theta}(\cdot, \cdot)=C_{\theta}(F_1(\cdot),F_2(\cdot)),
\end{equation}
for some continuous distribution $C_{\theta}(\cdot,\cdot)$ on $[0,1]\times [0,1]$, 
which is known except for the parameter 
$\theta$. We consider the problem of estimating $\theta$ when $F_1(\cdot)$ and $F_2(\cdot)$ are unknown baseline distributions.
If we set $U=F_1(X)$, $V=F_2(Y)$, then $(U,V)$ has distribution $C_{\theta}(\cdot,\cdot)$, and $C_{\theta}(\cdot,\cdot)$ is called a \emph{copula}. 
Note that if $(X_1,Y_1), \cdots, (X_n,Y_n)$ are independent with $(X_i,Y_i) \sim H_{\theta_i}(\cdot,\cdot)$, then 
$R_i \equiv Rank(X_i)=Rank(F_1(X_i))$ and 
$S_i \equiv Rank(Y_i)=Rank(F_2(Y_i))$,
which shows that the distribution of any statistical method based on these ranks will not depend on $(F_1,F_2)$. This model extends in the natural way to the d-dimensional case.

In this paper we consider the bivariate normal copula model where
$C_{\theta}(u,v)$ $=\Phi_{\theta}(\Phi^{-1} (u),\Phi^{-1} (v))$ with 
$C_{\theta}$ the bivariate $N(0,0,1,1,\theta)$ distribution.  We also consider 
the multivariate normal copula model and show that in regression experiments it can be 
used to construct a ``transform both sides regression'' transformation model (copula regression model.)
Klaassen and Wellner (1997) have shown that the normal scores correlation coefficient is semiparametrically efficient for the bivariate normal copula. We use simulations to compare this estimate with the profile MLE for the transform both sides Box-Cox regression model and a nonparametric estimate based on splines thereby augmenting the comparisons made by Zou and Hall (2002). The normal scores estimate is nearly as efficient as the parametric MLE for estimating median regression when the transform both sides Box-Cox model is correct. We also consider the performance of the estimates for models outside the copula regression model and find that the normal scores based estimate of median regression is remarkably robust with respect to both bias and variance.  On the other hand, the profile MLE of median regression derived from the transform both sides Box-Cox model is very sensitive to deviations from the model.  The nonparametric spline estimate is the best for extreme deviations from the copula regression model.                   

\section{Proportional Hazard and Proportional Expected Hazard Rate Models} 

Interesting special cases of (\ref{eq1.1}) are obtained by
considering the distributions of 
$T_1 =$ $\min(T_{01},$ $\cdots, T_{0k})$ and $T_2 = \max(T_{01}, \cdots, T_{0k})$ where $T_{01}, T_{02},\cdots$ are i.i.d. as $T_0 \sim F_0$.  Then, for $k\geq 1$,
\begin{equation}
\label{eq2.1}
T_1 \sim F_1(t)=1-[1-F_0(t)]^k, 
\end{equation}
and 

\begin{equation}
\label{eq2.2}
T_2 \sim F_2(t)=[F_0(t)]^k, 
\end{equation}
with $t>0, k=1,2, \cdots$. More general forms  (Lehmann (1953); Savage (1956)) are

\begin{equation}
\label{LEHMANN_MODEL}
T_1 \sim F_1(t)=1-[1-F_0(t)]^\Delta, 
\end{equation}
and 
\begin{equation}
T_2 \sim F_2(t)=[F_0(t)]^\Delta,
\end{equation}
with $t>0, \Delta>0$. Here (\ref{LEHMANN_MODEL}) can be derived by considering two-sample models where the two samples follow distributions of the form (\ref{eq2.1}) with different 
k's (Bickel and Doksum (2007),  Problem 1.1.12.) 

For $T_1$, the hazard rate is
\begin{equation}
\label{HR_COX}
\lambda(t) \equiv \ \frac{f(t)}{1-F(t)}=\triangle \frac{f_0(t)}{1-F_0(t)}  \equiv \triangle \lambda_0(t).
\end{equation}
In regression experiments, we set $\triangle_i=g(\bm{\beta}^T \bm{x}_i)$, and note that (\ref{HR_COX}) is the Cox proportional hazard (PH) 
model (Cox (1972)).

Nabeya and Miura (1972) proposed replacing k in (\ref{eq2.1}) and (\ref{eq2.2}) by a random variable.   
In particular, they considered 
$T_1 = \min(T_{01}, \cdots, T_{0N})$, where $N$ is independent of $T_{01},T_{02}, \cdots $, and has a zero truncated Poisson($\theta$) distribution with $\theta>0$. 
They also considered $T_2 =$ $\max(T_{01}, \cdots, T_{0M})$, $T_{0i} \sim F_0$ where $M$ is independent of $T_{01},T_{02}, \cdots $, and has a zero truncated Poisson($-\theta$) distribution with $\theta<0$.

Using Sibuya (1968), they found
\begin{eqnarray}
\label{SINAMI_PLUS}
 T_1 \sim F_1(t)=\frac{1-e^{-\theta F_0(t)}}{1-e^{-\theta}}, \quad \theta > 0,             
\end{eqnarray}
\begin{eqnarray}
\label{SINAMI_MINUS}
 T_2 \sim F_2(t)=\frac{1-e^{-\theta F_0(t)}}{1-e^{-\theta}}, \quad \theta < 0.            
\end{eqnarray}
Combining (\ref{SINAMI_PLUS}) and (\ref{SINAMI_MINUS}), we get
\begin{eqnarray}
\begin{aligned}
\label{SINAMI}
 T \sim F(t) &=& \frac{1-e^{-\theta F_0(t)}}{1-e^{-\theta}}, \quad \theta \not= 0,  \\
 \label{SINAMI0}
             & = & F_0(t),  \quad \quad \quad \quad \theta=0.
\end{aligned}
\end{eqnarray}
Note that model (\ref{SINAMI_PLUS}) is a mixture of proportional hazard models for individuals with the same baseline hazard rate $\lambda_0(\cdot)$ but different hazard factors 
$\triangle$ in the factorization (\ref{HR_COX}) of the hazard rate.  Let $\lambda(t;k)$ denote the hazard rate of 
$T_1$
given $N=k$; then by (\ref{HR_COX})
%% PP3
\begin{eqnarray}
\quad \quad \quad \quad  E\lambda(t;N)=\displaystyle\sum_{k=1}^{\infty}k\lambda_0(t)p_{\theta}(k)
=\tau(\theta) \lambda_0(t),\ \theta>0,
\end{eqnarray}
where $p_{\theta}(x)$ is the zero truncated Poisson($\theta$) probability and 
\begin{eqnarray}
\tau (\theta) =E(N)= \frac{\theta}{1-\exp(-\theta)}.
\end{eqnarray}
Thus (\ref{SINAMI_PLUS}) is a model with proportional expected hazard rate. Note that  (\ref{SINAMI}) does not have this property for $\theta<0$.  We will refer to
(\ref{SINAMI_PLUS}) and
(\ref{SINAMI})  
as the PEHR and SINAMI models, respectively.\\ 
\\
\noindent $\bm{Remark\ 2.1:}$ In regression experiments, the traditional frailty models
are also constructed by introducing a random element in the PH model.  However, these models are different from the PEHR and SINAMI models. To see this recall that in the frailty model the
conditional hazard rate given the covariate vector $\bm{x}$
(see Oakes (1992)) for the history and interpretation of frailty models) is of the form
\begin{eqnarray}
\label{FRAILTY}
\lambda_W(t|\bm{x})=\lambda_0 (t) W \exp[ \bm{\beta}^T \bm{x}], 
\end{eqnarray}
where W is a random effect that incorporates potential unobservable covariates that represent frailties. 
Semiparametric optimality theory for model (\ref{FRAILTY}) has been developed by Kosorok, Lee, and Fine (2004). 

Consider model (\ref{HR_COX}) with $\Delta=N$ and N a zero truncated Poisson($\theta$) random variable with 
$\theta=g(\bm{\beta}^T \bm{x})$, i.e., the conditional hazard rate given $\bm{x}$ is
\begin{equation}
\label{SINAMI_HR}
\lambda^{(N)}(t|\bm{x})=N \lambda_0(t).
\end{equation}
Here N plays the role of 
$ W \exp[ \bm{\beta}^T \bm{x}]$ in (\ref{FRAILTY}). 
However, (\ref{FRAILTY}) and (\ref{SINAMI_HR}) are different because N is an integer and $ W \exp[ \bm{\beta}^T \bm{x}]$ is not when $\bm{\beta} \not= 0$.
In model (\ref{SINAMI_HR}), 
N represents the effect of both observed covariates and frailties. In deriving the distribution function (\ref{SINAMI_PLUS}), the unobservable covariates are averaged out, that is, we compute 
$ P(T \leq t) = E[P(T \leq t|N)]$.\\
\\
\noindent $\bm{Remark\ 2.2:}$  Model  (\ref{SINAMI}) was considered by Bell and Doksum (1966), 
Example 5.2 and Table 8.1) and Ferguson (1967, p.257, Problem 5.7.7) without any of the above interpretations. Nabeya and Miura (1972) did not use any proportional hazard or frailty interpretation.  These concepts had not been invented yet.\\

Fig.\ref{FIG_HAZRATIO} gives a plot of the relative hazard rate $\lambda(t|x=1)/\lambda(t|x=0)$ with \linebreak $\theta = -3, -1.5, 1.5, 3$ for model (\ref{SINAMI}) with $\theta=\beta x$, $F_0(t)=1-\exp(-t),\ t>0$, and 
\begin{eqnarray}
\lambda(t|x)=\frac{\theta f_0(t)}{1-e^{-\theta(1-F_0(t))}}. 
\end{eqnarray}
In the PEHR and SINAMI models, the hazard ratio between two covariate values converge to unity as time increases.  This explains why the likelihoods for these models give less weight to long survival times than the likelihood for the Cox model (see Section \ref{SEC_RANK}).
The hazard rate is decreasing for the PEHR model for any continuous $F_0$.
\begin{figure}[h]
\begin{minipage}{1\hsize}
\begin{center}
\includegraphics[width=4in,]{doksum_Fig_HazRatio01_New.pdf}
\caption{ \label{FIG_HAZRATIO}\noindent SINAMI and PEHR hazard ratios for $\theta = -3, -1.5, 1.5, 3$.}
\end{center}
\end{minipage}
 \end{figure}

\section{Rank, Partial and Marginal Likelihood}

\label{SEC_RANK}

In regression experiments, we observe $(T_i,\bm{x}_i),\ i=1,\cdots,n$, where $T_1$,$\cdots,$ $T_n$ are independent responses and $\bm{x}_i$ is a nonrandom covariate vector.  In the proportional hazard model, it is customary to use model (\ref{HR_COX}) for $T_i$ with
$\triangle_i=\exp(\bm{\beta}^T \bm{x}_i)$
because $\triangle_i$ needs to be positive.  In the PEHR model, $\theta_i=\bm{\beta}^T \bm{x}_i$ is a
possible parametrization, but 
$\theta_i=\exp(\bm{\beta}^T \bm{x}_i)$ could also be used.
Let $\bm{R}=(R_1, \cdots, R_n)$ where $R_i=\mbox{Rank}(T_i)$, then
$l_{\bm{r}}(\bm{\beta}) = P(\bm{R}=\bm{r})$ is the \emph{rank likelihood}
(Hoeffding (1951)).


We first consider the one covariate case. Using the rank likelihood, the locally most powerful (LMP) rank test statistic for $H_0:\beta=0$ versus $H_1: \beta>0$ is (approximately) for the Cox model (Savage (1956), 1957), 
Cox (1964)), Oakes and Jeong (1998)): 
\begin{eqnarray}
\sum_{i=1}^n [-\log (1-\frac{R_i}{n+1})](x_i-\bar{x}) \quad \mbox{ (Savage or log rank)},
\end{eqnarray}
and for the PEHR and SINAMI models, the LMP rank test statistics is (Bell and Doksum (1966), 
Ferguson (1967),  
Nabeya and Miura (1972)):
\begin{eqnarray}
\sum_{i=1}^n \frac{R_i}{n+1}(x_i-\bar{x}) \quad \mbox{ (Wilcoxon type)},
\end{eqnarray}
where $R_i=$Rank($T_i$). The log rank statistic gives more weights to large observations, that is, in survival analysis, to those that live longer, while the Wilcoxon statistics is even handed. 


In order to compare how much relative weight is given to the small, in between, 
and large observed survival times for the PH and PEHR models, we next consider the rank likelihood for d covariates. Note that if $h(\cdot)$ is decreasing, then Rank$(h(T_i))=n+1-R_i$.  For the proportional hazard model, transform $T_i$ by $U_i=1-F_0(T_i)$, then by (\ref{LEHMANN_MODEL}) we have $f_{U_i}(u)=\triangle_i u^{\triangle_i-1}, 0 < u <1$.
Hoeffding (1951) formula shows,
\begin{eqnarray}
\mbox{Rank lkhd}= \Pi_{i=1}^n \triangle_i  \underset{0 < u_1<u_2< \cdots <u_n <1}{\int \cdots \cdots \cdots \int} \Pi_{i=1}^n u_i^{\delta_i-1}du_1\cdots du_n ,
\end{eqnarray}
where $\delta_i=\triangle_{b_i}$ and $b_i$= index on the $T$ with rank $n+1-i$= reverse anti-rank.  It follows that
\begin{eqnarray}
\mbox{Rank lkhd} \propto \Pi_{i=1}^{n} \frac{\triangle_i}{\sum_{k: T_k \geq T_{(i)}}\triangle_k },
\end{eqnarray}
that is, the familiar Cox (1972, 1975) \emph{partial likelihood} formula.
Here $\{ k: T_k \geq T_{(i)} \}$= patients at risk at time $T_{(i)}$ where $T_{(i)}$ is the $i$th ordered survival time.
Kalbfleich and Prentice (1973, 2002) 
called the rank likelihood the \emph{marginal likelihood} and extended it to censored data.

For the PEHR model, transform $T_i$ by the decreasing function  \\
$U_i$ = $a^{-1}\{\exp\{-F_0(T_i)\}-b\}$, 
then, we have 
$f_i(u)=a \tau_i (au+b)^{\theta_i-1}$,
$0<u<1$,
where
$b=e^{-1}$, $a=1-b$, and 
$\tau_i \equiv \tau(\theta_i)$. Let $\gamma_i=\theta_{b_i}$, $b_i$=index on the $T$ with rank $n+1-i$. Then,
\begin{eqnarray}
\quad \quad \quad \mbox{Rank lkhd} \propto (\Pi_{i=1}^n \tau_i) \underset{0 < u_1<u_2< \cdots <u_n <1}{\int \cdots \cdots \cdots \int} \Pi_{i=1}^n (au_i+b)^{\gamma_i-1} du_1 \cdots du_n.
\end{eqnarray}
If we perform the integration, we find that the likelihood for the PEHR model is similar to the likelihood for the Cox model except that in addition to terms involving
$\{k: T_k \geq T_{(i)}\}$, $i=1,\cdots,n$,
it includes terms involving
$\{k: T_{(i)} \leq T_k \leq T_{(j)}\}$ , $i=1,\cdots,n$, $j=1,\cdots,n$, $i \not =j$.
That is, the PEHR likelihood gives more weight to the intermediate survival times than the Cox likelihood.

Computationally, the Cox rank likelihood is easier than the PEHR rank likelihood. However, we can handle the PEHR rank likelihood with available algorithms and software (e.g. MATLAB.) More generally, $F(x)=C_{\theta}(F_0(x))$ type models, originally considered by 
Lehmann (1953),
 can be handled effectively by considering the profile NP likelihood of the next 
section (e.g. Tsodikov and Gabribotti (2007)),  Zeng and Lin (2007) ).


\section{Profile NP Likelihood}

Andersen, Borgan, Gill, and Keiding (1996),
Bickel, Klaassen, Ritov, and Wellner (1993, 1996), 
van der Vaart (1998), 
Murphy and van der Vaart (2000), Tsodikov and Garibotti (2007), Zeng and Lin (2007) 
and many others considered the problem of finding the MLE of all the parameters in a semiparametric model.  It is useful to divide the procedure into two steps by grouping parameters into two groups. 
Suppose the distribution function of T is of the form
$P(T \leq t) =$ $F(\theta, \eta(t))$, where $\theta \in R^d$ and $\eta(\cdot)$ is a nondecreasing function. If we assume temporarily that $\eta(\cdot)$ has a positive derivative $\eta'(t)$ for $t \in \{t_1, \cdots, t_n\}$, then the likelihood is
\begin{eqnarray*}
\label{}
\prod_{i=1}^{n} \eta'(t_i) f( \theta, \eta(t_i)),
\end{eqnarray*}
where $f( \theta, \eta)=\partial F(\theta, \eta)/\partial \eta$. The NP likelihood we consider is of the form
\begin{eqnarray*}
\label{}
L_{NP}(\theta, \eta)= \prod_{i=1}^{n} \eta \{t_i\} f(\theta, \eta [t_i]),
\end{eqnarray*}
where $\eta [t_i]= \sum_{j \leq i} \eta \{t_i\}$ is a step function with positive jumps 
$\eta \{t_i\}$ at the data points $t_i$, $i=1,\cdots,n$.

We assume that
\begin{eqnarray*}
a f( \theta, a) \rightarrow 0 \quad as\  a \rightarrow \infty 
\end{eqnarray*}

Next we fix  $\theta$, 
and define $\hat{\eta}_{\theta}\{t_i\}$ as
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{eqnarray}
\label{eq_PROLIK1}
\hat{\eta}_{\theta}\{.\} = ARG \ MAX_{\eta\{. \}} L_{NP}(\theta, \eta).
\end{eqnarray}
Set
\begin{eqnarray}
\label{eq_PROLIK2}
PROF\ NPLIK = l(\theta) = MAX_{\eta\{. \}} L_{NP}(\theta, \eta),
\end{eqnarray}
and solve
\begin{eqnarray}
\label{eq_PROLIK3}
\hat{\theta} = ARG \ MAX \ l(\theta).
\end{eqnarray}
Next estimate $\eta\{\cdot \}$ as 
$\hat{\eta}_{\hat{\theta}}\{\cdot \}$.
In the Lehmann model (\ref{eq1.1}), the NP likelihood is
\begin{eqnarray*}
\label{eq_LM_NP}
\prod_{i=1}^{n}  F_0\{x_i\} C'_{\theta}( \sum_{j \leq i} F_0\{x_j\}),
\end{eqnarray*}
The method is similar to finding the empirical MLE, Owen (1988, 2001), 
and profile (partial) MLE's as in Andersen, et al. (1996), and Murphy and van der Vaart (2000). \\
\\
\noindent $\bm{Remark\ 4.1:}$ Note that when
$P(T \leq t) =$ $F(\theta, \eta(t))$,
(\ref{eq_PROLIK1}),
(\ref{eq_PROLIK2}), and 
(\ref{eq_PROLIK3}) do not depend on the values of 
$t_1,\cdots,t_n$. In regression experiments, they will depend on the ranks of 
$t_1,\cdots,t_n$.
For example, see (\ref{BRESLOW_EQ}) and (\ref{eq_PEHR_PR1}).
This is in contrast to the 
Hodges and Lehmann (1963) approach that uses estimating equations based on rank test statistics to obtain estimates of parameters. In this Hodges-Lehmann ``rank inversion'' approach, estimates are functions of the ``raw'' data rather than the ranks.\\

As an example that will guide the algorithm for the PEHR model, consider the Cox model.  Set $\Lambda (t) = -\log (1-F_0(t))$, then,
\begin{eqnarray*}
L_{NP}(\bm{\beta},\Lambda) = \prod_{i=1}^{n} e^{\beta^T x_i}\Lambda\{t_i\} 
e^{-(e^{\beta^T x_i}) \Lambda[t_i]},
\end{eqnarray*}
where
\begin{eqnarray*}
\Lambda[t_i]= \sum_{j:t_j \leq t_i}\Lambda\{t_j\}.
\end{eqnarray*}
Using calculus, we find
\begin{eqnarray*}
\hat{\Lambda}_{\beta}\{t_i\} = ARG\ MAX_{\Lambda\{t_i\}} \ L_{NP}(\bm{\beta}, \Lambda) = (\sum_{j:t_j \geq  t_i} e^{\beta^T x_j})^{-1},
\end{eqnarray*}
and
\begin{eqnarray}
\label{BRESLOW_EQ}
l(\bm{\beta}) = PROF\ NPLIK = \prod_{i=1}^{n} \frac{e^{\beta^T x_i}}{ \sum_{j:t_j \geq  t_i} e^{\beta^T x_j}}.
\end{eqnarray}
This is exactly the same as the rank, the partial, and the marginal likelihood.

\section{Profile NP Likelihood for the PEHR Model}

Consider model  (\ref{SINAMI_PLUS}) with $\theta >0$. 
Set $\tau(\theta) = \theta [1-e^{-\theta}]^{-1}$, then
\begin{eqnarray}
\label{eq_PEHR_PR1}
L_{NP}(\bm{\theta}, F_0) = [\prod_{i=1}^{n} \tau(\theta_i) ] \prod_{i=1}^{n} F_0\{t_i\} e^{-\theta_i {F}_0[t_i]},
\end{eqnarray}
where ${F}_0[t_i]= \sum_{j:t_j \leq t_i}F_{0}\{t_j\}$.
Now set $p_i \equiv F_0\{t_i\}$ and maximize with respect to 
$p_1, \cdots,p_n$ with $\theta$ fixed.
The maximization problem looks very similar to Cox model maximization except for the constraint $\sum p_i =1$. We handle this constraint by writing $F_0(t)=1-\exp[-\Lambda(t)]$ with $\Lambda(t)$ unconstrained except for $\Lambda(t) \geq 0$ and by using
 a new approach based on the MM algorithm.

\subsection{The MM Algorithm} 
Lang, Hunter and Yang (2000) introduced a concept called the MM algorithm.  Its idea is that instead of maximizing a complicated original objective function, use a  simpler surrogate function so that each iteration is faster and guarantees that the original objective function increases. Given the original objective function $l(\bm{h})$ for a maximization problem,  a surrogate function $g(\bm{h}|\bm{h}_{old})$ must satisfy two properties:
\begin{eqnarray}
\label{SUR1}
l(\bm{h}_{old}) = g(\bm{h}_{old}|\bm{h}_{old}),
\end{eqnarray}
\begin{eqnarray}
\label{SUR2}
l(\bm{h}) \geq  g(\bm{h}|\bm{h}_{old}).
\end{eqnarray}
The EM algorithm is a special case of the MM algorithm.  A practical implementation issue of the MM algorithm is that we have to find a nice surrogate function case by case.  



Now we construct a surrogate function based on Tsodikov (2003).  Suppose we can write $l(\cdot)$ in the form  $l(\bm{h}) = B(\bm{h})-A(\bm{h})$ for some parameter vector $\bm{h}>\bm{0}$, where A and B are differentiable concave functions.  Then by the concavity property,
\begin{eqnarray}
\label{B_A}
 g(\bm{h}|\bm{h}_{old})=B(\bm{h})-A(\bm{h}_{old})-\nabla ^T A(\bm{h}_{old})(\bm{h}-\bm{h}_{old}),
\end{eqnarray}
where $\nabla^T A(\bm{h})=\partial A/\partial \bm{h}$ is the gradient of A, satisfies (\ref{SUR1}) and (\ref{SUR2}), and $g(\bm{h}|\bm{h}_{old})$ is a surrogate function for $l(\bm{h})$. Differentiating (\ref{B_A}) gives
%%% PP7
\begin{eqnarray}
\label{UPDATE1}
 \nabla ^T B(\bm{h}_{new})=\nabla ^T A(\bm{h}_{old}).
\end{eqnarray}
Solve (\ref{UPDATE1}) for $\bm{h}_{new}$.  Iterate the procedure until 
there is a minimal change in
$\bm{h}$ . 

\subsection{The MM Algorithm for the PEHR Model with $\theta \geq 0$ (SINAMI with $\theta \geq 0$)} 
\label{SEC_MM_PLUS}

Let $\theta_i = g(\bm{\beta}^{T} \bm{x}_i)$, for some known function $g(\cdot ) \geq 0$. When $\theta_i=0$, the distribution function of $T_i$ is $F_0(t)$. 
Let $F_0(t) = 1- \exp[-\Lambda(t)]$,
$h_k=\Lambda\{t_k\}$, 
and $\Lambda[t_i] = \sum_{k=1}^{i}{h_k}$ with $\Lambda(t) \geq 0$ and $h_k \geq 0$. Then for a temporarily fixed numerical vector $\bm{\beta}$,
\begin{eqnarray}
\begin{aligned}
\label{LOGLIK1}
l(\bm{h})=log[L_{NP}(\bm{\beta}, \bm{h})] &= 
%\displaystyle\sum_{i=1}^{n}(log(a(\theta_i))
\displaystyle\sum_{i=1}^n{\log \tau( \theta_i )}
+\displaystyle\sum_{i=1}^n{ \log h_i } \\
  & -\displaystyle\sum_{i=1}^{n} [\displaystyle\sum_{k=1}^{i}{h_k} +\theta_i (1-\exp (-\displaystyle\sum_{k=1}^{i}{h_k} ))].
\end{aligned}
\end{eqnarray}
Now we can write $l(\bm{h})=B(\bm{h})-A(\bm{h})$ with
\begin{eqnarray}
\label{B_org}
B(\bm{h})=\displaystyle\sum_{i=1}^n{ \log h_i },
\end{eqnarray}

\begin{eqnarray}
\begin{aligned}
\label{A_org}
A(\bm{h})&= \displaystyle\sum_{i=1}^{n} [\displaystyle\sum_{k=1}^{i}{h_k} +\theta_i (1-\exp (-\displaystyle\sum_{k=1}^{i}{h_k} ))].
\end{aligned}
\end{eqnarray}	
Here we may ignore $\Sigma \log \tau( \theta_i )$ because we maximize (\ref{LOGLIK1}) w.r.t. $\bm{h}$.

$B(\bm{h})$ and $A(\bm{h})$ are concave, because for $0  \leq t \leq 1$,
$B(t\bm{h}_a+(1-t)\bm{h}_b) \geq tB(\bm{h}_a)+(1-t)B(\bm{h}_b)$ and by mathematical induction,
$A(t\bm{h}_a+(1-t)\bm{h}_b) \geq tA(\bm{h}_a)+(1-t)A(\bm{h}_b)$ hold. Note that
\begin{eqnarray}
\label{B_h1}
\partial B/ \partial h_j = 1/h_j, \ j=1,...,n,
\end{eqnarray}	
\begin{eqnarray}
\label{A_h1}
\partial A/ \partial h_j = \displaystyle\sum_{i=1}^{n}(1+\theta_i (1-\exp (-\displaystyle\sum_{k=1}^{i}{h_k} )))1(j\leq i).
\end{eqnarray}	
Using (\ref{UPDATE1}), (\ref{B_h1}), and (\ref{A_h1}), update $h_j,\ j=1,\cdots,n$, at the same time,
\begin{eqnarray}
\label{MM_MAIN}
h_{j,new} = [\displaystyle\sum_{i=1}^{n}(1+\theta_i (1-\exp (-\displaystyle\sum_{k=1}^{i}{h_{k,old}} )))1(j\leq i)]^{-1}.	
\end{eqnarray} 
Iterate (\ref{MM_MAIN}) until there is a minimal change in $l(\hat{\bm{h}}_{new})$; call the result $\hat{\bm{h}}_{A}$ \linebreak ($\bm{A}$pproximated profile NPMLE).   Note that we call $\hat{\bm{h}}_{A}$ approximated profile NPMLE because $\bm{\hat{h}}_{A}$ is obtained by fixing $\bm{\beta}$. This approximation is necessary because there is no closed form $\hat{\bm{h}}$ w.r.t. $\bm{\beta}$. Next set $l(\bm{\beta})=log[L_{NP}(\bm{\beta}, \hat{\bm{h}}_A)]$ and maximize w.r.t. $\bm{\beta}$.

\subsection{The MM Algorithm for the SINAMI Model with $\theta \leq 0$} \
\label{SEC_THETA_NEG}
Consider model (\ref{SINAMI_MINUS}) with $\theta_i = g(\bm{\beta}^{T} \bm{x}_i)$, for some known function $g(\cdot) \leq 0$.
In this case we can use the algorithm of Section 5.2. To see this, suppose T satisfies model (\ref{SINAMI_MINUS}) with parameter $\theta_2 < 0$.  Set $V=1-F_0(T)$, then $V$ satisfies model (\ref{SINAMI_PLUS}) with parameter $\theta_1 = - \theta_2$. Moreover, the rank of $1-F_0(T_i)$ is $n+1-R_i$.

\subsection{The MM Algorithm for the SINAMI Model with $\theta \in R$} \
\label{SEC_THETA_BOTH}
Consider model (\ref{SINAMI_MINUS}) with $\theta_i = g(\bm{\beta}^{T} \bm{x}_i)$, for some known function $g(\cdot) \in R$.	In this case we can not use the transformation in Section \ref{SEC_THETA_NEG} because it changes the likelihood and the monotonicity of the likelihood as a function of $\theta$ does not necessarily hold. Instead, we modify the algorithm as follows:  If the value $\hat{\theta}_j$ in the jth iteration is positive, use the MM algorithm in Section  5.2. to find  $\hat{\bm{h}}_j$.  If $\hat{\theta}_j < 0$, then (\ref{LOGLIK1})
implies that finding the maximizer $\bm{h}$ is a convex optimization problem which produces $\hat{\bm{h}}_j$.
	
\subsection{Profile NPMLE Implementation} \
Successful convergence of the MM algorithm depends on a  good starting point ($\hat{\theta}, \hat{\bm{h}}$).
We consider the two sample problem:

\begin{eqnarray}
\label{}
T_{0,i} \sim  F_0(t),\ i=1,...,n_0, 
\end{eqnarray}
\begin{eqnarray}
\label{EQ6.20}
T_{1,i} \sim F(t) = \left \{
\begin{array}{l}
 \frac{1-e^{-\theta F_0(t)}}{1-e^{-\theta}}, 
\quad \quad (\theta > 0),\ i=1,...,n_1, \\  
		  F_0(t), \quad  \quad \quad \quad (\theta = 0),\ i=1,...,n_1, 
\end{array}
\right.		
\end{eqnarray}
where $F_0(\cdot)$ is an unknown distribution with density $f_0$. Note that the density of $F(t)$ is
\begin{eqnarray}
\label{SINAMI_PDF}
f(t; \theta)= \left \{
\begin{array}{l} \tau(\theta) f_0(t) e^{-\theta F_0(t)}, \quad \quad \quad (\theta > 0), \\
f_0(t), \quad \quad \quad \quad \quad \quad \quad \quad (\theta = 0).	
\end{array}
\right.
\end{eqnarray}

We use an algorithm to find ($\hat{\theta}, \hat{\bm{h}}$) where $\theta = \beta$ in this case.  For fixed $F_0$, (\ref{SINAMI_PDF}) gives an MOM estimating equation for $\theta$.  We plug in an estimate $\hat{F}_0$ for $F_0$, and use the following algorithm:
\begin{eqnarray*}
\label{SINAMI_ALG}
step\ (1):\hat{F}_0 \rightarrow \ 
step\ (2):\hat{\theta}_{} \rightarrow \ 
step\ (3):(\hat{\theta}_{A}\leftrightarrow \hat{\bm{h}}_{A}    )_{until\ \hat{\theta}_{A}\ converges.}
\end{eqnarray*}
Here $\hat{\theta}_{A}$ and $\hat{h}_{A}$ are 
approximated profile NPMLE's from
Section \ref{SEC_MM_PLUS}. The details are as follows:\\
$\bm{Step\ (1):}$ Compute the empirical distribution $\hat{F_0}(t)$ based on $T_{0,i}$ only:
\begin{eqnarray}
\label{SINAMI_ALG1}
\hat{F}_{0,[0]}(t) \equiv  \hat{F}_{0}(t) =\frac{1}{1+n_0}\displaystyle\sum_{i=1}^{n_0}{1( T_{0,i} \leq t)}.
\end{eqnarray}
Here $\hat{F}_{0}(t)\rightarrow_{a.s.} F_0(t)$ uniformly in $t$ as $n_0\ \rightarrow \infty$.
The subscript $[0]$ indicate iteration zero (starting point) for $step\ (3)$.  Note that the one-to-one relation between $\hat{F_0}$ and $\hat{\Lambda}_0$ is used to obtain $\hat{\bm{h}}$ by solving for ${\bm{h}}$ in the equations:
\begin{eqnarray}
\begin{aligned}
\label{F_H}
\hat{F}_{0}(t_i)&\equiv & \hat{F}_{0,i} &= \exp(-\hat{\Lambda}_{i}),\ where\\
\hat{\Lambda}[t_i] &\equiv & \hat{\Lambda}_{i} &= \displaystyle\sum_{k=1}^{i}\hat{h}_k,\ i=1,...,n_0.
\end{aligned}
\end{eqnarray}
\\
$\bm{Step\ (2):}$ Solve for $\hat{\theta}$ based on $\hat{F}_{0,[0]}(t)$:
\begin{eqnarray}
\begin{aligned}
	\bar{T}_1 &=& \tau(\theta) \int y e^{-\theta \hat{F}_{0,[0]}(y)} d\hat{F}_{0,[0]}(y) \\
&=& \tau(\theta) \displaystyle\sum_{i=1}^{n_0}T_{0,(i)} \exp(-\theta \frac{i}{n_0}) \frac{1}{n_0},	
\end{aligned}	
\end{eqnarray}
where $T_{0,(i)}$ is an order statistics of $T_{0,1},\cdots,T_{0,n_0}$ and $\bar{T}_1 =\displaystyle\sum_{i=1}^{n_1}T_{1,i}/n_1$. The solution  is uniquely determined because 
the distribution function (\ref{EQ6.20})
is monotone increasing in $\theta$ and hence its mean is monotone decreasing in $\theta$. If a model is $\theta\geq 0$ and $\hat{\theta}<0$, set $\hat{\theta}=0$. If a model is $\theta<0$ and $\hat{\theta}>0$, set $\hat{\theta}=0$.\\
$\bm{Step\ (3):}$ Compute $\hat{\theta}_{A}$ and $\hat{\bm{h}}_{A}$ as follows:
\begin{eqnarray}
\label{SINAMI_ALG}
(\hat{\theta}_{A}\leftrightarrow  \hat{\bm{h}}_A )_{until\ \hat{\theta}_{A}\ converges.}
\end{eqnarray}
%%%%%%%%%%%% Should move to Sec5.3   %%%%%%%%%%%%%%%%
The first iteration $\hat{\theta}_{A,[1]}$, is obtained by maximizing (\ref{LOGLIK1}) with $\hat{\theta}$ as a starting point, i.e., 
$\hat{\theta}_{A,[0]} =\hat{\theta} $ and with fixed $\hat{\bm{h}}_{[0]}$ obtained from (\ref{F_H}) and $\hat{F}_{0,[0]}$, i.e., 
\begin{eqnarray}
\hat{\theta}_{A,[1]} = \arg \max_{\theta}\ \{ l(\theta, \hat{\bm{h}}_{[0]}) : \theta \geq 0\}. 
\end{eqnarray}
Then by the MM algorithm in Section 5.2 with $\bm{\beta}_0 = \theta_0 = \hat{\theta}_{A,[1]}$ (see (\ref{MM_MAIN})), obtain
$\hat{\bm{h}}_{[1]}$ using the starting point $\hat{\bm{h}}_{[0]}$. Next obtain 
\begin{eqnarray}
\hat{\theta}_{A,[2]} = \arg \max_{\theta}\ \{ l(\theta, \hat{\bm{h}}_{A,[1]}) : \theta \geq 0\},
\end{eqnarray}
with starting point $\hat{\theta}_{A,[1]}$.  Then by the MM algorithm in Section 5.2 with $\bm{\beta}_0 = \theta_0 = \hat{\theta}_{A,[2]}$, obtain
$\hat{\bm{h}}_{A, [2]}$ with starting point $\hat{\bm{h}}_{A,[1]}$.
Repeat the procedure to get $\hat{\theta}_{A,[j]}$ and $\hat{\bm{h}}_{A,[j]}$ until convergence, i.e., $| \hat{\theta}_{A,[j]}- \hat{\theta}_{A,[j-1]}| < \epsilon$ for some small $\epsilon$.

Numerical optimizations for $\hat{\theta}$ and $\hat{\theta}_{A}$ are carried out by the MATLAB fmincon() function. \\
\\
\noindent $\bm{Remark\ 5.1:}$ For fixed $\hat{\bm{h}}$, $l(\theta, \hat{\bm{h}})$ is strictly concave and have a unique maximum.\\
\\
\noindent $\bm{Remark\ 5.2:}$ The estimate of $\bm{\beta}$ in the Cox model that we have discussed is asymptotically optimal in the semiparametric sense (
Begun, Hall, Huang, and Wellner (1983),
Bickel et al. (1993, 1998),
van der Vaart (1998)),
Murphy and van der Vaart (2000)).
These references and others give results that can be used to check the semiparametric asymptotic optimality of the profile NPMLE in the PEHR model.\\
\\
\noindent $\bm{Remark\ 5.3:}$ \emph{Transformation models}. We can show that the PEHR is a special case of 
transformation models as follows: Let $F_{\lambda}$ be the exponential $(\lambda)$ distribution function and define

\begin{eqnarray}
\label{eq_R5.3}
 G_0(y|\bm{x}) = \frac{1-e^{-\theta F_{\lambda}(y)}}{1-e^{-\theta}}, \quad y > 0,
\end{eqnarray}
where $\theta=g(\bm{x},\bm{\beta})$. Let $\psi$ be an increasing function from $[0,\infty)$ to $[0,\infty)$ and define the transformation model
\begin{eqnarray}
\label{eq_R5.3}
 G(y|\bm{x}) = G_0(\psi(y)|\bm{x}).
\end{eqnarray}

This model is of the form (\ref{SINAMI_PLUS}) with
$F_0=F_{\lambda} \psi(t)$.
Klaassen (2007) gives results
for general transformation models
that can be used to check semiparametric asymptotic efficiency of estimates of $\bm{\beta}$ in the model
(\ref{SINAMI_PLUS}).

\subsection{Estimation of the Variance of the Profile NPMLE}  
Hypothesis tests and confidence intervals require standard errors (estimates of the standard deviation) of $\hat{\theta}_{A}$.  An algorithm developed by Tsodikov and Garibotti (2007) 
combined with the preceding algorithm
allows us to compute the profile information matrix which is the observed information matrix derived from the profile likelihood.  This provides standard errors  $SE(\hat{\theta}_A)$ of $\hat{\theta}_{A}$.

\section{Simulation Results} 
\subsection{PEHR Model Estimates} 
Monte Carlo (MC) simulation for model
(\ref{EQ6.20}) with $F_0$ equal to the exponential distribution EXP(1)
is based on 1000 Monte Carlo samples with n=100, 200, and 300. $T_{0,i} \sim EXP(1)$, $i=1,\cdots,n_0$, $T_{1,i} \sim PEHR( \theta=2,\ or\ 3),\ i=1,\cdots,n_1$, $n_0=n_1=n/2$, iid. We compute 
Monte Carlo estimates of 
the expected values, standard deviations (SD's), and MSE's of 
$\hat{\theta}_{MOM}$,
$\hat{\theta}_{A}$, and
$\hat{\theta}_{PAR}$ where
$\hat{\theta}_{PAR}$ is
the parametric model MLE, obtained by assuming that $F_0$ is known and equal to the EXP(1) distribution.

We also compute the Monte Carlo estimates of the expected values 
$E[SE$ $(\hat{\theta}_{A})]$ of the standard errors computed as described in Section 5.5.
Table \ref{SINAMI_TB1} summarizes the result. 

\begin{table}[tbp]
\begin{center}


\begin{tabular}{|l||r|r|r|r|r|r|}
\hline
$\theta$ &
\multicolumn{3}{|c|}{2}& 
\multicolumn{3}{|c|}{3} \\ 
\hline
n & 100 & 200 & 300 & 100 & 200 & 300 \\
\hline
$E[\hat{\theta}_{MOM}]$ & 1.97 & 2.00 & 2.02 & 3.05 & 3.00 & 3.01 \\
$E[\hat{\theta}_{A}]$ & 1.99 & 2.00 & 2.01 & 3.07 & 3.00 & 3.01 \\
$E[\hat{\theta}_{PAR}]$ & 2.02 & 2.00 & 2.01 & 3.03 & 3.01 & 3.00 \\
\hline 
$SD(\hat{\theta}_{MOM})$ & 0.819 & 0.551 & 0.479 & 0.972 & 0.635 & 0.506 \\
$SD(\hat{\theta}_{A})$ & 0.788 & 0.541 & 0.466 & 0.960 & 0.618 & 0.497 \\
$SD(\hat{\theta}_{PAR})$ & 0.541 & 0.387 & 0.317 & 0.605 & 0.431 & 0.335 \\
$E[SE(\hat{\theta}_{A})]$ & 0.752 & 0.527 & 0.430 & 0.837 & 0.578 & 0.471 \\
\hline
$MSE[\hat{\theta}_{MOM}]$ & 0.671 & 0.303 & 0.230 & 0.949 & 0.403 & 0.256 \\
$MSE[\hat{\theta}_{A}]$ & 0.621 & 0.292 & 0.217 & 0.926 & 0.383 & 0.248 \\
$MSE[\hat{\theta}_{PAR}]$ & 0.289 & 0.150 & 0.101 & 0.367 & 0.186 & 0.112 \\
\hline
\end{tabular}
\caption[mysmallcaption]{PEHR model simulation estimates (MC=1000, $\theta=2,3$)}
\label{SINAMI_TB1}
\end{center}
\end{table}

Overall $\hat{\theta}_{MOM}$, $\hat{\theta}_{A}$, and  $\hat{\theta}_{PAR}$ have almost no bias in the estimation of $\theta=$2 or 3. 
As expected, the parametric model estimate $\hat{\theta}_{PAR}$ has the smallest MSE.  The approximated profile NPMLE  $\hat{\theta}_{A}$ has a smaller MSE than $\hat{\theta}_{MOM}$, but the difference is small.
The 
%\citeauthor{r42}%
approximation to SD($\hat{\theta}_{A}$) is very good and improves as the sample size increases. 


\subsection{Model Fit for Misspecified Model} 

Next consider a model that is neither a Cox PH model nor a PEHR model: i.e., $T_{0,i} \sim EXP(1)$ and the
true model for $T_{1,i}$ is 
: Case 1, Gamma(shape=0.5, \linebreak scale=0.5), and the
: Case 2, Weibull(shape=0.5, scale=0.2). 
Here the target values $\theta$ and $\bm{h}$ are those that minimize the Kullback-Leibler divergence between the true distribution and the model class of distributions  (Doksum, Ozeki, Kim and Neto (2007)) .

Fig. \ref{FIG_2outside} shows that PEHR gives better fit than the Cox model.
\begin{figure}[h!]
\vspace{-0pc}
\begin{tabular}{ll}
\begin{minipage}{0.6\hsize}
%\begin{center}
\includegraphics[width=60mm, height=60mm]{doksum_Gp0_5_Gm0_5_HazardRatio.pdf}
\begin{center}\hspace{-2pc} $(a)\ True \ model: \mbox{Gamma}\ (shape=0.5, $\vspace{.35pc}

\hspace{6pc} $scale=0.5)$.
\end{center}
\end{minipage}\hspace{-3pc}
\begin{minipage}{0.6\hsize}
%\begin{center}
\includegraphics[width=60mm,height=60mm]{doksum_Wm5_Wr0_5_HazardRatio.pdf}
\begin{center}\hspace{-4pc}\hspace{2pc} $(b)\ True \ model: \mbox{Weibull}\ (shape=0.5, $\vspace{.35pc}

\hspace{6pc} $scale=0.2)$.
\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_2outside}\noindent Cox and PEHR estimated hazard ratio when $F_0 \sim EXP(1)$ and $F \sim$ Gamma or Weibull.}
\end{figure}


\section{Estimation in the Normal Copula Model} 
\subsection{The One Covariate Case}
Assume that the pair (X,Y) has a joint density f(x,y) with respect to Lebesgue measure on $R^2$ and a joint distribution function F(x,y). Let $F_1$ and $F_2$ be the marginal distribution functions of X and Y, respectively, and let $\Phi$ denote the standard normal distribution function. Consider the transformations $X \rightarrow  Z$ $= \Phi^{-1}(F_1(X))$, $Y \rightarrow  W$ $= \Phi^{-1}(F_2(Y))$. Then the marginal distributions of Z and W are standard normal.  The bivariate normal copula model $\mathcal{F}$ is defined by the assumption that the joint distribution of (Z,W) is bivariate normal with zero mean, unit variance, and correlation coefficient $\rho$. That is,
\begin{eqnarray*}
\label{COPULA01}
\mathcal{F}=\{ F\ : (\Phi^{-1}(F_1(X)), \Phi^{-1}(F_2(Y))) \sim N(0,0,1,1,\rho) \},
\end{eqnarray*}
where $F_1$ and $F_2$ are the marginals of F.
Let $\{(X_i,Y_i), i=1,2, \cdots,n \}$ be independent and identically distributed with distribution function 
$F  \in \mathcal{F}$, 
and set 
$Z_i = \Phi^{-1}(F_1(X_i))$,
$W_i =$ $\Phi^{-1}$ $(F_2(Y_i))$, $i=1,2, \cdots,n $. If we (temporarily) assume that $F_1$ and $F_2$ are known, then
because E(ZW) = $\rho$, a method of moments "estimate" of $\rho$, is ${r}_{MOM} = n^{-1} \sum_{i=1}^n{Z_iW_i}$.  
The asymptotic distribution of 
$\sqrt{n} ( {r}_{MOM}-  \rho)$ 
is $N(0,1+ \rho^2)$ when $F \in \mathcal{F}$.  Assuming $F_1$ and $F_2$ known, the asymptotic distribution of
$\sqrt{n} ( {r}_{MLE}-  \rho)$, where ${r}_{MLE}$ is the maximum likelihood "estimate" of $\rho$ is $N(0, (1-\rho^2)^2/(1+ \rho^2))$.  The asymptotic variance $(1-\rho^2)^2/(1+ \rho^2)$ of
$ {r}_{MLE}$ is smaller than the asymptotic variance $(1-\rho^2)^2$ of the usual Pearson correlation coefficient  ${r}_P$ and much smaller than the asymptotic variance $1+ \rho^2$ of ${r}_{MOM}$. 

Note that $  \mathcal{F}$ is invariant under coordinate-wise increasing transformations.  That is, if $(X,Y) \sim F \in  \mathcal{F}$ and $U=h_1(X), V=h_2(Y)$ with $h_1$ and $h_2$ increasing, then the distribution G of (U,V) is in $  \mathcal{F}$. If we want methods that are invariant under such transformations, we must use statistics based on the ranks defined in Section 1.

Suppose next that $F_1$ and $F_2$ are unknown. It may then make sense to replace the ordered Z's and W's by their expected values.  This leads to the Fisher and Yates (1938) or normal scores $E(Z_{(i)})$, $i=1, \cdots,n$ where $Z_{(1)}, \cdots, Z_{(n)}$ are N(0,1) order statistics. We write $a(i)=E(Z_{(i)})$. An accurate approximation to  
$E(Z_{(i)})$ is $\Phi^{-1}[ (i-3/8)/(n+1/4)]$, e.g. Cox (2006).

Let $Z_i^{'} = a(R_i)$, $W_i^{'} = a(S_i)$ where $R_i$ and $S_i$ are the ranks of $X_i$ and $Y_i$ when the X's and Y's are ranked separately.  Then we obtain estimates $ \hat{\rho}_{MOM}$,
$ \hat{\rho}_{MLE}$, and $ \hat{\rho}_{P}$ of $\rho$ 
when $F_1$ and $F_2$ are unknown by replacing $Z_i$ and $W_i$ by 
 $Z_i^{'}$ and  $W_i^{'}$ in ${r}_{MOM}$, ${r}_{MLE}$ and ${r}_{P}$. 
In this case $ \hat{\rho}_{MOM}$, $ \hat{\rho}_{P}$ are nearly identical and asymptotically equivalent, but they are different from $ \hat{\rho}_{MLE}$.
We will use $ \hat{\rho}_{P}$ because it is slightly less biased, and denote it by $\hat{\rho}_{NS}$ where NS signifies normal scores.  Thus
\begin{eqnarray}
\hat{\rho}_{NS}= \sum Z'_i W'_i /\sum a^2(i).
\end{eqnarray}
It follows from Bhuchongkul (1964)
that based on the rank likelihood, $\hat{\rho}_{NS}$ is, uniformly in $F_1$ and $F_2$, a locally most powerful test statistics in the bivariate normal copula model.
 Zou and Hall (2002) gave an asymptotic extension of this result. They also
computed the rank likelihood estimate of $\rho$ in the 
bivariate normal copula model using an improved version of the likelihood  sampler in Doksum(1987).

Klaassen and Wellner (1997) found
$\sqrt{n} ( \hat{\rho}_{NS}-  \rho) \rightarrow _{d} N(0, (1-\rho^2)^2)$ in the copula model $\mathcal{F}  $  with $F_1$ and $F_2$ unknown;
the same as for the Pearson correlation in the bivariate normal model.  In fact, in a bivariate 
normal$(\mu_1, \mu_2, \sigma_1^2, \sigma_2^2, \rho)$ 
model, ${r}_{P}$ is the MLE, and ${r}_{P}$ and $\hat{\rho}_{NS}$ are asymptotically optimal in the parametric sense.

A fourth possible estimate is the profile NP estimate obtained by fixing $\rho$ and replacing $F_1$ and $F_2$ by step functions with jumps $\{p_i\}$ and $\{q_i\}$ at 
$(X_i,Y_i)$ in the log likelihood for the normal copula model.  That is, ignoring constants ($\rho$ is fixed), we maximize
\begin{eqnarray}
\begin{aligned}
\label{COPU_L1}
l(\bm{p}, \bm{q})&= \sum_i \{ \log{p_i} + \log{q_i}   \\
&+ \frac{1}{2}((\Phi^{-1}(\sum_{k: X_k \leq X_i} p_k) )^2 
+ \frac{1}{2}(\Phi^{-1}(\sum_{k: Y_k \leq Y_i} q_k))^2  \\
&+ \frac{1}{2} (1-\rho^2 )^{-1} 
[
(\Phi^{-1}(\sum_{k: X_k \leq X_i} p_k))^2 \\
&- 2\rho \Phi^{-1}(\sum_{k: X_k \leq X_i} p_k) \Phi^{-1}(\sum_{k: Y_k \leq Y_i} q_k) +
(\Phi^{-1}(\sum_{k: Y_k \leq Y_i} q_k))^2
]\},
\end{aligned}
\end{eqnarray}
w.r.t. $(\bm{p}, \bm{q})$ where $\sum p_i =1$ and  $\sum q_i =1$. Then given $(\hat{\bm{p}}, \hat{\bm{q}})$, maximize the log likelihood w.r.t. $\rho$, which gives $\hat{\rho}_{PROF}$, a profile NPMLE.\\
$ $\\
\noindent $\bm{Remark\ 7.1:}$ An estimate $\hat{\theta}$ of a parameter $\theta \in R$ in a semiparametric model is regular if 
$\sqrt{n}(\hat{\theta} - \theta) \rightarrow_d$
$ N(0, V_{\hat{\theta}}(\theta, \eta))$
for some asymptotic variance $V_{\hat{\theta}}(\theta, \eta)$ 
and if $\hat{\theta}$ satisfies additional regularity conditions given in Bickel et al. (1993, 1998). 
For $F \in \mathcal{F}$, $V_{{r}_P}(\rho, F)$ depends on F (e.g. Bickel and Doksum (2007), Example 5.3.6), while 
$V_{\hat{\rho}_{NS}}(\rho, F)$ does not, as shown by 
Klaassen and Wellner (1997).
Klaassen and Wellner (1997) go on to argue that $(1-\rho^2)^2$ is a semiparametric asymptotic variance lower bound for the class $\mathcal{S}$ of all regular estimates of $\rho$.  Thus
$\hat{\rho}_{NS}$ is semiparametrically optimal in the minimax sense:
\begin{eqnarray}
\label{remark7.1eq}
\sup\{V_{\hat{\rho}_{NS}}(\rho, F): F \in \mathcal{F}\} =
\inf_{\hat{\rho} \in \mathcal{S}}
\sup\{V_{\hat{\rho}}(\rho, F): F \in \mathcal{F}\} 
\end{eqnarray}

\noindent $\bm{Remark\ 7.2:}$ 
Recall that $\hat{\rho}_{NS}$ was obtained by inserting normal scores 
$Z'_i$ and $W'_i$ in the MOM estimate for the model with $F_1$ and $F_2$ known, and that the MLE ${r}_{MLE}$ for this model has variance  
 $(1-\rho^2)^2/(1+ \rho^2)$.
Klaassen and Wellner (1997) have shown that the
approximate MLE $\hat{\rho}_{MLE}$ obtained from 
${r}_{MLE}$ by replacing $(Z_i, W_i)$ with
$(Z_i', W_i')$ is semiparametrically optimal in the same sense as 
$\hat{\rho}_{NS}$.  
Because the distribution of the ranks do not depend on $F_1$ and $F_2$, this implies that 
$\hat{\rho}_{NS}$ and
$\hat{\rho}_{MLE}$ are asymptotically equivalent for every 
$F \in \mathcal{F}$.
We conjecture that 
$\hat{\rho}_{PROF}$ is also asymptotically optimal and equivalent to
$\hat{\rho}_{NS}$.\\
\\
\noindent $\bm{Remark\ 7.3:}$ 
The asymptotic distribution of 
$\hat{\rho}_{NS}$ when the distribution of (X,Y) is not in $\mathcal{F}$ can be obtained from  
Ruymgaart, Shorack, and Van Zwet (1972) and
Ruymgaart (1974) . \\

\subsection{The Multivariate Covariate Case}
\label{SEC_MCOV}
The normal copula model in the multivariate case is defined as follows: Let $Y \sim G$, $X_j \sim F_j$ ,
$\bm{h}(\bm{X})=(h_1 (X_1), \cdots, h_d (X_d))$, where $h_j$, $j=0,\cdots,d$ are increasing functions
defined by
\begin{eqnarray}
\label{EQ_Y_COPULA}
h_0(Y) = \Phi^{-1}(G(Y)),
\end{eqnarray}
\begin{eqnarray}
\label{EQ_X_COPULA}
h_j(X_j) =  \Phi^{-1}(F_j(X_j)).
\end{eqnarray}
The distribution of the untransformed variables 
$(\bm{X}, Y)$
is a \emph{copula model} if we assume that 
($\bm{h}(\bm{X})$, $h_0(Y)$) is multivariate normal with 0 means and unit variances.

\section{Transformation and NP Models}

Consider a regression experiment with response Y and a random covariate vector $\bm{X}=(X_1,\cdots,X_d)^T$. 
We will extend the normal scores estimate 
$\hat{\rho}_{NS}$ of Section 7 
to the  d dimensioned case and compare it
with estimates appropriate for parametric and nonparametric models.
In the  \emph{copula regression} model of Section \ref{SEC_MCOV}, we can write
\begin{eqnarray}
\label{EQ8.1}
h_0(Y)= \beta^T \bm{h}(\bm{X})+ \epsilon, \ \epsilon \sim N(0, \sigma^2),
\end{eqnarray}
where $\bm{\beta}$ is the set of regression coefficients when regressing $ h_0(Y) $ on $\bm{h}(\bm{X})$.
The transform both sides Box-Cox model is based on (\ref{EQ8.1}) with  $h_0(Y)=Y^{(\lambda_{d+1})}$, 
$h_j(X_j)=X_j^{(\lambda_{j})}$, $j=1, \cdots, d$, where 
$t^{(\lambda)}=(t^{\lambda}-1)/\lambda$. Thus for this case, we can write

\begin{eqnarray}
\label{EQ8.2}
Y^{(\lambda_{d+1})}=\alpha + \bm{\beta}^T \bm{X}^{(\bm{\lambda})} + \epsilon, 
\ \epsilon \sim N(0, \sigma^2).
\end{eqnarray}
We first consider a procedure for estimating the parameters in model (\ref{EQ8.2}):\\

$\bm{I}$ Profile Likelihood for a multivariate model.

Hernandez and Johnson (1980) considered the one sample multivariate Box-Cox transformation model. This was adopted to regression by
Doksum, Ozeki, Kim and Neto (2007).
We regard 
$(Y^{(\lambda_{d+1})},\bm{X}^{(\bm{\lambda})})$ as a d+1 multivariate normal $(\bm{\mu}, \Sigma)$ vector. 
Regressing 
$Y^{(\lambda_{d+1})}$ on
$\bm{X}^{(\bm{\lambda})}$
leads to (\ref{EQ8.2}).
We fix $\xi \equiv (\bm{\lambda},\lambda_{d+1})$ and estimate the parameters in the normal model by maximizing the likelihood thereby obtaining the familiar normal theory estimates
$(\hat{\bm{\mu}(}\xi), \hat{\Sigma}(\xi))$.  We plug these into the likelihood and obtain the profile likelihood $l(\xi)$, which we maximize to get $\hat{\xi}$ and the final estimates $(\hat{\bm{\mu}(}\hat{\xi}), \hat{\Sigma}\hat{(\xi)})$.  These are the usual linear model estimates with $Y_i$, $X_{ij}$ replaced by 
$Y_i^{(\hat{\lambda}_{d+1})}$,
$X_{ij}^{(\hat{\lambda}_j)}$.
Similarly, the estimate of $\bm{\beta}$ in (\ref{EQ8.2}) is the usual linear model estimate with 
$Y_i$,
$X_{ij}$ replaced by
$Y_i^{(\hat{\lambda}_{d+1})}$,
$X_{ij}^{(\hat{\lambda}_j)}$.\\
\\
\noindent $\bm{Remark\ 8.1:}$  
We also considered the maximum likelihood estimates of the parameters $\bm{\beta}, \sigma^2$, and $(\bm{\lambda},\lambda_{d+1})$ in model (\ref{EQ8.2}).
This approach has the problem that if we want to test $H_0:\beta_j=0$, then $\lambda_j$ is not identifiable under $H_0$. Approach I does not have this problem. This is one case where likelihood and profile likelihood are very different. The algorithm for this MLE often failed to converge. When it did converge, it produced results close to those of method I. We omit the details.\\
\\
\noindent $\bm{Remark\ 8.2:}$  
As pointed out by Zou and Hall (2002), when d=1, the MLE of $\rho$ in the 
Box-Cox transformation model with unknown transformation parameters and standardized transformations have 
the same efficiency as the MLE for the model with known transformation 
parameters because this Box-Cox model is between the bivariate normal model 
and the bivariate normal copula model and the MLE's in these models have 
the same asymptotic variance $(1-\rho^2)^2$.
The result that the efficiency is the same whether or not the $\lambda$'s  are known in this Box-Cox model was also obtained by Wong (1981).
This result is very different from the 
results of Bickel and Doksum (1981) regarding the estimation of regression coefficients.\\
\\
\noindent $\bm{Remark\ 8.3:}$    
Consider the transformation model
\begin{eqnarray}
\label{EQ8.3NEW}
h_0(Y)=\bm{\beta}^T \bm{X} + \epsilon, \ \epsilon \sim N(0, \sigma),
\end{eqnarray}
where $\bm{X}$ is a vector of random covariates and $h(\cdot)$ is increasing. In this case we can consider the rank estimate $\hat{\bm{\beta}}_R$ obtained by maximizing the rank likelihood $l_{\bm{r}}(\bm{\beta}) = P(\bm{R}=\bm{r})$ defined in Section \ref{SEC_RANK}. The results of Bickel and Ritov (1997) imply that in a certain sense  
$\hat{\bm{\beta}}_R$ is semiparametrically optimal for model (\ref{EQ8.3NEW}). However the normal scores estimate of 
$\hat{\bm{\beta}}= (\bm{x}^T \bm{x})^{-1} \bm{x}^T \bm{a}$, where
$\bm{a}= (a(S_1),\cdots,a(S_n))^T  $
and $\bm{x}$ is a vector of nonrandom covariates,
is not asymptotically optimal unless $|\bm{\beta}|/\sigma$ tends to zero at a certain rate as $n \rightarrow \infty$
(Doksum (1987)). MC methods for $\hat{\bm{\beta}}_R$ is introduced in Bickel and Doksum (2009), 
Section 10.5.\\


We next introduce a semiparametric approach for the copula regression model and a nonparametric regression approach.\\
\begin{figure}[t]
\vspace{-0pc}
\begin{minipage}[t]{125mm}
\centering
\includegraphics[ angle=0, width=125mm, ]{doksum_rho_boxplot_nSim1000_nSample128_L1n0_5_L2n0_5.pdf}
\includegraphics[ angle=0, width=125mm, ]{doksum_rho_boxplot_nSim1000_nSample128_L1n1_L2n1.pdf}
\end{minipage}
 \caption{ \label{FIG_rho}\noindent $\hat{\rho}$ boxplots for correctly specified model (\ref{eq_sim1}) with n=128. Top 5: $\lambda_1 = \lambda_2=0.5$. Bottom 5: $\lambda_1 = \lambda_2=1$. The horizontal line gives the true value of $\rho$.}
\end{figure}

$\bm{II}$ Normal score substitution.

The model (\ref{EQ8.1}) with $h_j$, $j=0,\cdots,d$, satisfying 
(\ref{EQ_Y_COPULA}) and (\ref{EQ_X_COPULA}) is invariant under increasing transformations.  As in the d=1 case, this leads to using the ranks $\{S_i\}$ of the Y's and the ranks
$\{R_{ij}: i=1, \cdots, n\}$ of $X_{ij}$ among 
$\{X_{ij}: i=1, \cdots, n, j=1, \cdots, d\}$.
Because the distribution of the ranks is invariant under increasing transformations, for rank methods, model (\ref{EQ8.1})
is equivalent to
\begin{equation}
\label{EQ8.3}
Y' = \bm{\alpha}^T \bm{X}' + \epsilon',
\end{equation}
where 
$X_{j}' \sim N(0,1)$ and 
$\epsilon '$ 
are independent, 
$Y' \sim N(0,1)$ and 
$\bm{\alpha}^T = \Sigma^{-1} \bm{\rho}$ with
$  \bm{\rho}= ( Corr(X_1',Y'), \cdots, Corr(X_d',Y') )^T$ and 
$\Sigma $ the correlation matrix of 
%%$\bm{X}=(X_{ij} ')$.
$\bm{X}'=(X_{j} ')_{d \times 1}$. 
Here $\Sigma$ is assumed to be nonsingular.
Based on the distribution of the ranks, $\alpha_1, \cdots, \alpha_d$ are identifiable parameters in model (\ref{EQ8.3}).
These parameters represents the relative importance of the $X_j$'s.

The normal scores $Z_{ij}' = a(R_{ij})$ and $W'_i = a(S_i)$
have approximately the same distribution as the unobservable 
$X'_{ij}$ and $Y'_i$ in model (\ref{EQ8.3}). Because
$E(Y'|\bm{x}')=\bm{\alpha}^T \bm{x}'$, if we replace $X_{ij}'$ and $Y'_i$ with $Z'_{ij}$ and $W'_i$, we find that an approximate method of moments estimate of $\alpha$ is
\begin{equation}
\label{EQ_ALPHA_LSE}
\hat{\bm{\alpha}} = (Z_D^T Z_D)^{-1} Z_D^T \bm{W'},
\end{equation}
\noindent where $Z_D$ is the no intercept design matrix $(Z'_{ij})_{n \times  d'}$,
$d'$ is the rank of the matrix $(Z'_{ij})$,
and 
$\bm{W'}=$  $(W'_1, \cdots,$ $W'_n)^T$.
%%%%%%%%%%%%%%%%%%%%%%%%%%
Any subset of variables ${X_j: j \in J}$ with the same ranks, say $R_{1J},\cdots,$ $R_{nJ}$, is collapsed into one variable denoted as $X_J$ with ranks $R_{1J},\cdots, R_{nJ}$ to avoid singularity. Based on
Klaassen and Wellner (1997), we conjecture that 
$\hat{\bm{\alpha}}$
is semiparametrically efficient for the multivariate normal copula model.\\

$\bm{III}$ Nonparametric estimation. 

We next introduce a nonparametric approach. We consider the model 
\begin{equation}
\label{NONPARA1}
 Y=m(\bm{X})+ \epsilon,
\end{equation}
where $m()$ is unknown and $\epsilon$ has median zero. To estimate $m()$, we use a cubic B-spline and the R function smooth.spline().
The number of knots are automatically selected (less than the number of observations n.)  The smoothing parameter is chosen by generalized cross validation (GCV.)

\subsection{Simulation Results} 
We consider the d=1 case and consider the properties of estimates of $\rho= Corr(h_1(X),$ \linebreak $h_0(Y))$. In this case, the method II estimate is $\hat{\rho}_{NS}$.

\subsubsection{Correctly Specified Model} 
The true model satisfies
\begin{equation}
\label{eq_sim1}
 Y^{(\lambda_2)}=\alpha_0 + \alpha_1 X^{(\lambda_1)}+ \epsilon,
\end{equation}
where   
$\epsilon \sim N(0, \sigma^2)$ and $X^{(\lambda_1)} \sim N(\mu_1, \sigma_0^2)$ are independent.
This model is a subset of the normal copula model $\mathcal{F}$ with
\begin{equation}
\label{eq_F1F2}
	F_1(x) = \Phi(\frac{x^{(\lambda_1)}-\mu_1}{\sigma_0}),\ 
	F_2(y) = \Phi(\frac{y^{(\lambda_2)}-\mu_2}{\sigma_2}),
\end{equation}
where $\mu_2=E Y^{(\lambda_2)}$, and 
$\sigma^2_2=Var Y^{(\lambda_2)}$. We use 
1000 MC trials and
take
$\sigma^2 =1$,
$\sigma_0^2=1$,
$(\lambda_1,\ \lambda_2) \in \{(0.5,0.5), (1,1) \}$,
$\alpha_0=6$,
$\alpha_1 = \in \{0, 0.1, 0.5, 1, 2\}$, and
$\mu_1 = 5$.

 

Fig. \ref{FIG_rho} shows that methods I and II have similar properties for $\alpha_1 \leq 0.5$. For larger $\alpha_1$, the normal scores estimate has a downward bias which is
negligible for $n \geq  500$ (not shown here.) Method I converges all the time 
with the constraint $-\bm{4} \leq  \bm{\lambda} \leq  \bm{4}$. Method II does not involve any optimization and hence converges all the time.

\subsubsection{Misspecified Model} 

We simulate the data from
\begin{equation}
\label{eq_sim2}
 Y^{(\lambda_2)}=(1-\gamma)(\alpha_0 + \alpha_1 X^{(\lambda_1)}) +
 \gamma[L(X)] + \epsilon,
\end{equation}
where $L()$ is a nonlinear function.
Thus the model is a Box-Cox model when $\gamma=0$, but when $\gamma>0$, we are checking the performance of the methods when the model generating the methods are misspecified.

For comparisons of methods we need a parameter that makes sense for all three methods.  One such parameter is
\begin{equation}
m(x)=Median(Y|X=x).
\end{equation}
We consider the 25th, 50th, and 75th population quantiles of X, i.e., our parameters of interest are 
$m(x_{0.25})$,
$m(x_{0.50})$, and 
$m(x_{0.75})$.

Methods I and II are based on models of the form
\begin{equation}
h_0(Y)=g(\bm{X},\bm{\beta}) + \epsilon,
\end{equation}
where $h_0(\cdot)$ is an increasing function. If $X$ and $\epsilon$ are independent and $median(\epsilon)$ $=0$, then 
\begin{equation}
\label{EQ8.4}
m(x)=h_0^{-1}(g(\bm{X},\bm{\beta}) ).
\end{equation}
For method I, the MLE of m(x)  in model (\ref{eq_sim1}) is, 
\begin{equation}
\hat{m}(x) =(\hat{\lambda}_2 [\hat{\beta}_0 +\hat{\beta}_1 \frac{x^{\hat{\lambda}_1}-1}{\hat{\lambda}_1} ] +1)^{1/\hat{\lambda}_2}.
\end{equation}
For method II, write model (\ref{EQ8.1}) as 
\begin{equation}
h_0(Y)= \rho h_1(X)+ \epsilon, \ \epsilon \sim N(0, \sigma^2),
\end{equation}
where $\rho \equiv \rho(h_1(X), h_0(Y))$ is the correlation coefficient. Then by 
(\ref{EQ8.4}), 
\begin{equation}
m(x) = h_0^{-1}(\rho h_1(x)),
\end{equation}
where
\begin{equation}
h_0^{-1}(t) = F_2^{-1}(\Phi(t)).
\end{equation}
It follows that
\begin{equation}
m(x) = {F_2}^{-1}(\Phi(\rho \Phi ^{-1}(F_1(x)) ),
\end{equation}
and by replacing $F_1$ and $F_2$ by their empiricals, a natural estimate of m(x) is
\begin{equation}
\hat{m}(x) = y_{([n\Phi(\hat{g}(x))])},
\end{equation}
where $\hat{g(}x)=\hat{\rho}_{NS} \Phi^{-1}(\hat{F_1}(x))$ and $[\ ]$ is the greatest integer function. 


For method III, we use the smoothing spline estimate of
$E(Y|X=x)$ described earlier.  In our models with normal errors, $E(Y|X=x)$ coincide with the conditional median m(x).

In the simulation we use model (\ref{eq_sim2}) with
$\epsilon \sim N(0, \sigma^2)$ and $X^{(\lambda_1)} \sim N(\mu_1, \sigma_0^2)$ independent,
\begin{equation}
L(t)=\alpha_0+ \alpha_1
 \mu_1-1.25+2.5[1+\exp(-10(t-\mu_1))]^{-1},
\end{equation}
$\sigma^2 \in \{0.01, 0.1, 0.5, 1\}$,
$\sigma_0^2=1$,
$(\lambda_1,\ \lambda_2) \in \{(0.5,0.5), (1,1) \}$,
$\alpha_0=6.25$,
$\alpha_1 = \in \{0, 0.5, 1, 2\}$, 
$\mu_1 = 5$, 
and
$\gamma \in \{0, 0.25, 0.5, 0.75, 1 \}$.
The sample size is n=512.  There are 1000 MC trials.


Fig. 
\ref{FIG_L1_a1n0.5_var0.5},   
\ref{FIG_L1_a1n0.5_var1},
\ref{FIG_L0.5_a1n0.5_var0.5}, and
\ref{FIG_L0.5_a1n0.5_var1}
 are boxplots of $\hat{m}(x_{0.25})$, $\hat{m}(x_{0.50})$, and $\hat{m}(x_{0.75})$ with the setting
$(\lambda_1,\lambda_2,$ $\alpha_1, \sigma^2)$ = $(1,1,0.5, 0.5)$,
$(\lambda_1,\lambda_2,$ $\alpha_1, \sigma^2)$ = $(1,1,0.5, 1)$,
$(\lambda_1,\lambda_2,$ $\alpha_1, \sigma^2)$ = $(0.5, 0.5,0.5, 0.5)$, and
$(\lambda_1,\lambda_2,$ $\alpha_1, \sigma^2)$ = $(0.5, 0.5,0.5, 1)$ respectively. 
Fig.8-11 give MSE's for the estimates 
$\hat{m}(x_{0.25})$,
$\hat{m}(x_{0.50})$, and
$\hat{m}(x_{0.75})$.

We see that method I is the best when model (\ref{EQ8.2}) is correct, that is, $\gamma = 0$.  However when the model is increasingly misspecified, i.e., as $\gamma$ increases, its absolute bias increases which leads to low MSE performance.


Method II is overall best in terms of MSE when $\lambda_1=\lambda_2=0.5$ (Fig 6 and 7). When $\gamma=0$, it is unbiased and its variance is between Method I and Method III (Fig 4,5,6, and 7).

Method III is overall best in terms of MSE when $\lambda_1=\lambda_2=1$  and the model is badly misspecified. It's smaller bias makes up for its large variance in this case. But its MSE suffers at and near model (\ref{EQ8.2}) (Fig 6, 7, $\gamma=0$).

In summary, the normal score procedure performs very well at and close to a copula model.  For n large, this is to be expected from the results of Klaassen and Wellner (1997).  The normal scores estimate is competitive with the Box-Cox estimate in the transform both sides Box-Cox model.



\begin{figure}[h!]
%\vspace{-0pc}
%\hspace{-1.5pc}
\begin{tabular}{l}
\begin{minipage}{1\hsize}
%\begin{center}
\includegraphics[width=125mm, height=165mm]
{doksum_boxplot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr0_5.pdf}
\noindent %%(a)$\ \theta = -3$
%\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_L1_a1n0.5_var0.5}
\noindent Boxplots of the three estimates of median regression m(x) for model (\ref{eq_sim2}) with ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(1, 1, 0.5, 0.5). I: profile MLE, II: normal scores, and III: NP, spline. The true value of m(x) is the solid line.}
\end{figure}


\begin{figure}[h!]
\vspace{-0pc}
%\hspace{-1.5pc}
\begin{tabular}{c}
\begin{minipage}{1\hsize}
\begin{center}
\includegraphics[width=125mm, height=165mm]
{doksum_boxplot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr1.pdf}
\noindent %%(a)$\ \theta = -3$
\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_L1_a1n0.5_var1}\noindent Boxplots of the three estimates of median regression m(x) for model (\ref{eq_sim2}) with ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(1, 1, 0.5, 1). I: profile MLE, II: normal scores, and III: NP, spline. The true value of m(x) is the solid line.}
\end{figure}

\begin{figure}[h!]
\vspace{-0pc}
%\hspace{-1.5pc}
\begin{tabular}{c}
\begin{minipage}{1\hsize}
\begin{center}
\includegraphics[ width=132mm, height=165mm ]
{doksum_boxplot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr0_5.pdf}
\noindent %%(a)$\ \theta = -3$
\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_L0.5_a1n0.5_var0.5}\noindent Boxplots of the three estimates of median regression m(x) for model (\ref{eq_sim2}) with ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(0.5, 0.5, 0.5, 0.5). I: profile MLE, II: normal scores, and III: NP, spline. The true value of m(x) is the solid line.}
\end{figure}

%%%%%%%%%%%%%% lambda=0.5 case
\begin{figure}[h!]
\vspace{-0pc}
%\hspace{-1.5pc}
\begin{tabular}{c}
\begin{minipage}{1\hsize}
\begin{center}
\includegraphics[ width=125mm, height=165mm ]
{doksum_boxplot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr1.pdf}
\noindent %%(a)$\ \theta = -3$
\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_L0.5_a1n0.5_var1}\noindent Boxplots of the three estimates of median regression m(x) for model (\ref{eq_sim2}) with ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(0.5, 0.5, 0.5, 1). I: profile MLE, II: normal scores, and III: NP, spline. The true value of m(x) is the solid line.}
\end{figure}
\clearpage
\newpage

\begin{figure}[h!]
%\vspace{-2pc}
\begin{tabular}{ccc}
\begin{minipage}{0.38\hsize}
%\begin{center}
%\resizebox{1.0}{1.0}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr0_5_25.pdf}
%\noindent $(a)MSE25$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr0_5_50.pdf}
%\noindent $(b)MSE50$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm  ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr0_5_75.pdf}
%\noindent $(c)MSE75$
%\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_3MSE}\noindent MSE of the three estimates of m(x) as a function of the misspecification parameter $\gamma$ for ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(1, 1, 0.5, 0.5). $\bigcirc = I$, $\bigtriangleup = II$, $+ = III$.}
\end{figure}

\begin{figure}[h!]
%\vspace{-2pc}
\begin{tabular}{ccc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm  ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr1_25.pdf}
%\noindent $(a)MSE25$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr1_50.pdf}
%\noindent $(b)MSE50$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm  ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n1_L2n1_b0n6_25_b1n0_5_varErr1_75.pdf}
%\noindent $(c)MSE75$
%\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_4MSE}\noindent MSE of the three estimates of m(x) as a function of the misspecification parameter $\gamma$ for ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(1, 1, 0.5, 1). $\bigcirc = I$, $\bigtriangleup = II$, $+ = III$.}
\end{figure}

\begin{figure}[h!]
%\vspace{1pc}
\begin{tabular}{ccc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm  ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr0_5_25.pdf}
%\noindent $(a)MSE25$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr0_5_50.pdf}
%\noindent $(b)MSE50$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm  ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr0_5_75.pdf}
%\noindent $(c)MSE75$
%\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_3MSE}\noindent MSE of the three estimates of m(x) as a function of the misspecification parameter $\gamma$ for ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(0.5, 0.5, 0.5, 0.5). $\bigcirc = I$, $\bigtriangleup = II$, $+ = III$.}
\end{figure}

\clearpage
\newpage
\begin{figure}[h!]
%\vspace{-4pc}
\begin{tabular}{ccc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr1_25.pdf}
%\noindent $(a)MSE25$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr1_50.pdf}
%\noindent $(b)MSE50$
%\end{center}
\end{minipage}\hspace{-2pc}
\begin{minipage}{0.38\hsize}
%\begin{center}
\includegraphics[width=38mm,height=52mm ]
{doksum_MSE_plot5_nSim1000_nSample512_L1n0_5_L2n0_5_b0n6_25_b1n0_5_varErr1_75.pdf}
%\noindent $(c)MSE75$
%\end{center}
\end{minipage}
\end{tabular}
 \caption{ \label{FIG_4MSE}\noindent MSE of the three estimates of m(x) as a function of the misspecification parameter $\gamma$ for ($\lambda_1, \lambda_2,\alpha_1, \sigma^2$)=(0.5, 0.5, 0.5, 1). $\bigcirc = I$, $\bigtriangleup = II$, $+ = III$.}
\end{figure}

\section*{Acknowledgements} 
We thank Peter Bickel, Aad van der Vaart and Jon Wellner for helpful comments.

%\bibliographystyle{plain}
\begin{thebibliography}{0}
 \bibitem{r1}
 \textsc{Andersen, P. K., Borgan, O., Gill, R. D.} and \textsc{Keiding N.} (1996).
	\textit{Statistical Models Based on Counting Processes}.
 Springer-Verlag, New York.
 
\bibitem{r2}
  \textsc{Begun, J. M., Hall, W. J., Huang, Wei-Min} and \textsc{Wellner, J. A.} (1983).
  Information and Asymptotic Efficiency in Parametric-Nonparametric Models.
 \textit{The Annals of Statistics},
 \textbf{11}, 2, 432--452. 
 

\bibitem{r3}
 \textsc{Bell, C. B.} and \textsc{Doksum, K. A.} (1966).
  ``Optimal'' One-Sample Distribution-Free Tests and Their Two-Sample Extensions.
 \textit{The Annals of Mathematical Statistics},
 \textbf{36}, 1, 120--132.


\bibitem{r4}
 \textsc{Bhuchongkul, S.} (1964).
A Class of Nonparametric Tests for Independence in Bivariate Populations.
 \textit{Ann. Math. Statist.},
 \textbf{35},1, 138--149.

 
\bibitem{r5}
 \textsc{Bickel, P. J.} and \textsc{Doksum, K. A.} (1981).
  An analysis of transformations revisited.
 \textit{Journal of the American Statistical Association},
 \textbf{76}, 296--311.

\bibitem{r6}
  \textsc{Bickel, P. J.} and \textsc{Doksum, K. A.} (2007).
	\textit{Mathematical Statistics: Basic Ideas and Selected Topics, 2nd ed. Vol I}, Updated Printing.  
   {Pearson Prentice Hall, Upper Saddle River, NJ.}

\bibitem{r7}
  \textsc{Bickel, P. J.} and \textsc{Doksum, K. A.} To appear (2009).
	\textit{Mathematical Statistics: Basic Ideas and Selected Topics Vol II.}  
   {Pearson Prentice Hall, Upper Saddle River, NJ.}


\bibitem{r8}
 \textsc{Bickel, P. J., Klaassen, C. A., Ritov, Y.} and \textsc{Wellner, J. A.} (1993, 1998).
  \textit{Efficient and Adaptive Estimation for Semiparametric Models}.
 {Springer-Verlag, New York}.

\bibitem{r9}
 \textsc{Bickel, P. J.} and \textsc{Ritov, Y.} (1997).
  Local asymptotic normality of ranks and covariates in transformation models. In \textit{Festschrift for Lucien Le Cam (D. Pollard and G. L. Yang, eds.)}
 {Springer, New York}.


\bibitem{r10}
 \textsc{Cox, D. R.} (1964).
  Some applications of exponential ordered scores.
 \textit{J. R. Stat. Soc.},
 \textbf{B.26}, 103--110.


\bibitem{r11}
 \textsc{Cox, D. R.} (1972).
  Regression models and life tables (with discussion).
 \textit{J. Roy. Statist. Soc.},
 \textbf{B.34}, 187--220.
 
\bibitem{r12}
 \textsc{Cox, D. R.} (1975).
  Partial likelihood.
 \textit{Biometrika.},
 \textbf{62.2}, 269--276. 
 
\bibitem{r13}
 \textsc{Cox, D. R.} (2006).
  \textit{Principles of Statistical Inference}.
 {Cambridge University Press, Cambridge},
 \textbf{} 40. 
 
\bibitem{r14}
 \textsc{Doksum, K. A.} (1987).
   An extension of partial likelihood methods from 
proportional hazard models to general transformation models.
 \textit{Annals of Statistics}, 
 \textbf{15}, 325--345. 
 
 
 
 
\bibitem{r15}
 \textsc{Doksum, K. A., Ozeki, A., Kim, J.} and \textsc{Neto, E. C.} (2007).
 Thinking outside the box: Statistical inference based on Kullback-Leibler Empirical Projections.  
 \textit{Statistics and Probability Letters}.
 \textbf{77}, 1201--1213. 

\bibitem{r16}
 \textsc{Fisher, R. A.} and \textsc{Yates, F.} (1938).
\textit{Statistical Tables for Biological, Agricultural and Medical Research}. (5th edition, 1957)
 {Oliver and Boyd, Edinburgh}.


\bibitem{r17}
 \textsc{Ferguson, T. S.} (1967).
  \textit{Mathematical Statistics. A Decision Theoretic Approach}.
 {New York and London, Academic Press}.


\bibitem{r18}
 \textsc{Hajek, J.} and \textsc{Sidak, Z.} (1967).
  \textit{Theory of Rank Tests}.
 {Academic Press, New York}.

\bibitem{r19}
 \textsc{Hodges, J. L.} and \textsc{Lehmann, E. L.} (1963).
 Estimates of location based on rank tests. 
 \textit{Ann. Math. Stat},
 \textbf{34}, 598–-611. 


 \bibitem{r20}
 \textsc{Hoeffding, W.} (1951).
  `Optimum' nonparametric tests.
 \textit{Proc. 2nd Berkeley Symposium on Mathematical Statistics and Probability}, Berkeley, Univ. Calif. Press,
 \textbf{} 83--92.
 
\bibitem{r21}
 \textsc{Hernandez, F.} and \textsc{Johnson, R. A.} (1980).
 The large-sample behavior of transformations to normality.
 \textit{Journal of the American Statistical Association},
 \textbf{75}, 855--861.
%%%
\bibitem{r22}
 \textsc{Kalbfleich, J. D.} and \textsc{Prentice, R. L.} (1973).
Marginal likelihoods based on Cox's regression and life model.
 \textit{Biometrika},
 \textbf{60}, 267--278.

\bibitem{r23}
 \textsc{Kalbfleich, J. D.} and \textsc{Prentice, R. L.} (2002).
  \textit{The Statistical Analysis of failure Time Data}, 2nd edition.
 {John Wiley and Sons, Inc., Hoboken, New Jersey}.

\bibitem{r24}
 \textsc{Klaassen, C. A. J.} (2007).
  A Sturm-Liouville Problem in Semiparametric Transformation Models.  In \textit{Advances In Statistical Modeling And Inference.  Essays in Honor of Kjell A Doksum}. Editor, Vijay Nair.
 {World Scientific Pub Co Inc., New Jersey}.


\bibitem{r25}
 \textsc{Kosorok, M. R., Lee B. L.} and \textsc{Fine J. P.} (2004).
Robust inference for univariate proportional hazards frailty regression models.
 \textit{The Annals of Statistics},
 \textbf{32}, No. 4, 1448--1491.



\bibitem{r26}
 \textsc{Klaassen, C. A. J.} and \textsc{Wellner, J. A.} (1997).
  Efficient Estimation in the Bivariate Normal Copula Model: Normal Margins Are Least Favourable.
 \textit{Bernoulli},
 \textbf{3}, No. 1, 55--77.

\bibitem{r27}
 \textsc{Lange, K., Hunter, D. R.} and \textsc{Yang, I.} (2000).
  Optimization Transfer Using Surrogate Objective Functions.
 \textit{Journal of Computational and Graphical Statistics},
 \textbf{9},  No. 1, 1--20.
 

\bibitem{r28}
 \textsc{Lehmann, E. L.} (1953).
  The power of rank tests.
 \textit{Ann. Math. Statist.},
 \textbf{}24, 23--43.

\bibitem{r29}
 \textsc{Murphy, S. A.} and \textsc{van der Vaart, A. W.} (2000).
  On Profile Likelihood.
 \textit{JASA},
 \textbf{9}, No. 450, 449--465.


\bibitem{r30}
 \textsc{Nabeya, S.} and \textsc{Miura, R.} (1972).
  Locally Most Powerful Rank Tests for Lehmann mixture alternatives.
 \textit{Technical Report}. University of California, Berkeley.
 
\bibitem{r31}
 \textsc{Oakes, D.} (1992).
Bivariate survival models induced by frailties.
 \textit{JASA},
 \textbf{84}, 487--493.

\bibitem{r32}
 \textsc{Oakes, D.} and \textsc{Jong-Hyeon Jeong} (1998).
Frailty Models and Rank Tests.
 \textit{Lifetime Data Analysis},
 \textbf{4}, 3, 209--228.


\bibitem{r33}
 \textsc{Owen, A. B.} (1988).
 Empirical likelihood ratio confidence intervals for a single functional.
 \textit{Biometrika},
 \textbf{75}, No.2,  237--249.

\bibitem{r34}
 \textsc{Owen, A. B.} (2001).
 \textit{Empirical Likelihood}.
 {Chapman and Hall.}

%%%%%%%%%
\bibitem{r35}
 \textsc{Ruymgaart, F. H.} (1974).
  Asymptotic Normality of Nonparametric Tests for Independence. 
 \textit{Ann. Statist.},
 \textbf{2}, 892--910.
 
\bibitem{r36}
 \textsc{Ruymgaart, F. H., Shorack, G. R.} and \textsc{Van Zwet, W. R.} (1972).
  Asymptotic Normality of Nonparametric Tests for Independence.
 \textit{Ann. Math. Statist.},
 \textbf{43}, 1122--1135.
 
 



\bibitem{r37}
 \textsc{Savage, I. R.} (1956).
  Contributions to the theory of rank orders statistics: the two-sample case.
 \textit{Ann. Math. Stat.},
 \textbf{27}, 590--615.
 
\bibitem{r38}
 \textsc{Savage, I. R.} (1957).
  Contributions to the Theory of Rank Order Statistics--The "Trend" Case.
 \textit{Ann. Math. Stat.},
 \textbf{28}, 4, 968--977.
 
\bibitem{r39}
 \textsc{Sibuya, M.} (1968).
  Generating Doubly Exponential Random Numbers.
 \textit{Annals of the Institute of Statistical Mathematics, Tokyo, Supplement V},
 \textbf{} 1--7.
 
%\bibitem[Shao (2003)]{SHAO2003}
% \textsc{Shao, J.} (2003).
% Mathematical Statistics 2nd ed.
% \textit{Springer-Verlag, New York, NY.}
% \textbf{}, 321.   

\bibitem{r40}
 \textsc{Sklar, A.} (1959).
 Fonctions de repartition a n dimensions et leurs marges.
 \textit{L'Institut de Statistique de L'Universite de Paris.}
 \textbf{8}, 229--231.   


\bibitem{r41}
 \textsc{Tsodikov, A.} (2003).
  Semiparametric models: a generalized self-consistency approach.
 \textit{J. R. Statist. Soc.},
 \textbf{B.65}, 3, 759--774.


\bibitem{r42}
 \textsc{Tsodikov, A.} and \textsc{Garibotti, G.} (2007).
  Profile information matrix for nonlinear transformation models.
 \textit{Lifetime Data Analysis},
 \textbf{13}, 139--159.

\bibitem{r43}
 \textsc{van der Vaart, A. W.} (1998).
  \textit{Asymptotic statistics}.
 {Cambridge University Press, Cambridge, UK}.


\bibitem{r44}
 \textsc{Zeng, D.} and \textsc{Lin, D. Y.} (2007).
  Maximum likelihood estimation in semiparametric regression models with censored data.
 \textit{J. R. Statist.},
 \textbf{69}, 4, 1--30. 
 
\bibitem{r45}
 \textsc{Zou, K. H.} and \textsc{Hall, W. J.} (2002).
  On estimating a transformation 
correlation coefficient. 
 \textit{Journal of Applied Statistics},
  \textbf{29}, 745--760. 

 

\bibitem{r46}
 \textsc{Wong, C. W.} (1981).
  Transformation of Independent Variables in Regression Models.
 \textit{Ph.D. Thesis}. Department of Statistics, University of Berkeley, CA.
\end{thebibliography}




\end{document}





