# Estimation of Treatment Effects for Heterogeneous Matched Pairs Data with
Probit Models

Jun Wanga Wei Gao${}^{a\ *}$ Man-Lai Tangb
a Key Laboratory for Applied Statistics of MOE, School of Mathematics and
Statistics,
Northeast Normal University, Changchun, Jilin 130024, China
b Department of Mathematics and Statistics, Hang Seng Management College, Hong
Kong, China

###### Abstract

Estimating the effect of medical treatments on subject responses is one of the
crucial problems in medical research. Matched-pairs designs are commonly
implemented in the field of medical research to eliminate confounding and
improve efficiency. In this article, new estimators of treatment effects for
heterogeneous matched pairs data are proposed. Asymptotic properties of the
proposed estimators are derived. Simulation studies show that the proposed
estimators have some advantages over the famous Heckman’s estimator and
inverse probability weighted (IPW) estimator. We apply the proposed
methodologies to a blood lead level data set and an acute leukaemia data set.

Key Words: Confounding; Grouped effects; Heterogeneity; Matched pairs data;
Probit models; Treatment effects

##  1 Introduction

In many medical studies, researchers are often interested in evaluating the
treatment effects for some improved surgery or new medicine, and one popular
measurement is average treatment effects (ATE)

$ATE=E[Y(1)-Y(0)]$

where $Y(0)$ and $Y(1)$ denote the two potential outcomes given the control
treatment ($D=0$) and the active treatment($D=1$). Neyman (1923) originally
proposed potential outcomes models, Rubin (1974) generalized it, which is
usually called Rubin causal models(RCM). In order to effectively estimate
treatment effects, one challenge problem is to control confounding, and Abadie
and Imbens (2006, 2011), Díaz, et al. (2015) proposed covariates balanced
methods to estimate average treatment effects (ATE) under strong ignorability
conditions.

Another method to control confounding is the matched pairs studies (e.g.,
twins study), in which one of the subjects is randomly selected to expose to
treatments and the other is selected to receive the control treatment. Let $A$
and $B$ be the pairs with covariates ${\bm{x}}^{A}$ and ${\bm{x}}^{B}$
respectively, $D$ denote the treatment status with $D=0$ for control and $D=1$
for new treatment, and $Y^{A}$ be the response value for subject $A$ with
treatment status $D$ and $Y^{B}$ is the response value for subject $B$ with
treatment status $1-D$. When response value is linearly relationship with
treatment status, treatment effects are usually estimated by difference-in-
differences (DID) approaches (Ashenfelter and Card, 1985; Imbens and Angrist,
1994; Card and Krueger, 1994; Abadie and Gardeazabal, 2003). Lately, Sakaguchi
(2016) considered the time-varying treatment effects model for panel data and
estimated the time-varying average treatment effects by difference-in-
differences (DID) approaches.

However, in practice, in terms of matched pairs data, the relationship between
response value and treatment status is usually nonlinear, namely, the response
value is assumed to be binary with $0$ representing negative effect and $1$
representing positive effect. The following model is assumed:

$Y^{A}=I({{\bm{\beta}}}^{{}^{\prime}}{{\bm{x}}}^{A}+\lambda
D+\eta^{A}>0),\mbox{~{}and}\;\;\;Y^{B}=I({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}+\lambda(1-D)+\eta^{B}>0)$
(1)

where ${\bm{\beta}}$ and $\lambda$ are $k\times 1$ regression coefficients and
treatment parameters respectively, and $\eta^{A}$ and $\eta^{B}$ are errors.

If $\eta^{A}$ and $\eta^{B}$ in Model (1) are independent of ${\bm{x}}^{A}$,
${\bm{x}}^{B}$ and $D$, various authors have considered the problem of
estimating the treatment effect. For example, Butler and Moffitt (1982)
proposed a computationally efficient quadrature procedure to obtain the ML
(maximum likelihood) estimators for the one-factor multinomial probit model
while Bertschek and Lechner (1998) proposed a set of convenient GMM
(generalized methods of moments) estimators for the binomial probit model
based on panel data. Rubin (1974) presented a discussion of matching,
randomization, random sampling, and other methods of controlling extraneous
variation. Robins, et al. (1994), Robins and Rotnizky (1995), Hirano, et al.
(2003), Hernan and Robins (2006) proposed the inverse probability weighted
estimator for treatment effects.

In practice, matched pairs data are collected from different backgrounds and
it is not uncommon that some variables may not be observed for some pairs.
Besides, although paired observations may share common features (e.g., same
genes or environments), there may exist a great deal of heterogeneity among
paired observations. For non-linear models, estimation of treatment effects in
the presence of heterogeneity is a challenging problem. Heckman et al. (1997)
gave detailed discussions on this issue and proposed the following models

$Y^{A}=I({{\bm{\beta}}}^{{}^{\prime}}{{\bm{x}}}^{A}+\lambda
D+\tau+\epsilon^{A}>0),\mbox{~{}and}\;\;\;Y^{B}=I({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}+\lambda(1-D)+\tau+\epsilon^{B}>0)$
(2)

where ${\bm{\beta}}$ and $\lambda$ are $k\times 1$ regression coefficients and
treatment parameters respectively, $\tau$ is the unobserved grouped effects
(e.g., twins share common genes) and is independent of $\epsilon^{A}$ and
$\epsilon^{B}$, and $\epsilon^{A}$ and $\epsilon^{B}$ are errors with means
$0$ and variances $1$.

It is noteworthy that $\tau$ given in Models $(\ref{MB})$ can be considered a
constant (i.e., fixed effect model) or random variable (i.e., random effect
models). If we assume that $\tau$ is a random variable, there’s usually have a
prior information about group effects. Hence, assuming that $\epsilon^{A}$ and
$\epsilon^{B}$ are normally distributed and $\tau$ has prior normal
distributions, Heckman (1978) developed the maximum likelihood and other
estimates based on the multivariate probit model with structural shift to
estimate ${\bm{\beta}}$ and $\lambda$. However, Heckman’s estimates can be
very sensitive to the prior distribution for $\tau$. Estimation problems of
random effect models is discussed in Guilkey and Murphy(1993), Arellano and
Bonhomme (2009), Walker (1996), Dey et al. (1997), Horowitz (1992) and
Arulampalam and Stewart (2009). Recently Gao et al. (2017) proposed a new
method to estimate the parameters for dynamic probit models.

When $\tau$ is fixed effect, different estimation methods have been considered
by various authors for panel data. For instance, Arellano (2003) applied the
orthogonal technique to analyze discrete models for panel data. Carro (2007)
demonstrated that the orthogonal technique can reduce the order bias of the
maximum likelihood estimator from $O(T^{-1})$ to $O(T^{-2})$ where $T$ is the
length of individual observations. Bartolucci et al. (2016) proposed the
modified profile likelihood estimator for fixed-effects panel data models.

In this article, we will consider models (2) when $\epsilon^{A}$ and
$\epsilon^{B}$ are normally distributed with means $0$ and variances $1$
(i.e., probit model). Estimators for ${\bm{\beta}}$ and $\lambda$ are
considered and their theoretical properties are presented in Section 2.
Simulation studies will be conducted and the results will be reported in
Section 3. In Section 4, we demonstrate our methodologies with two real data
sets. Conclusions will be discussed in Section 5. Proofs of lemmas and
theorems will be presented in Appendix.

##  2 Probit Model and Estimator

By assuming that $\epsilon^{A}$ and $\epsilon^{B}$ are normally distributed
with means being $0$ and variances being $1$, Models (2) becomes the famous
Probit Models. We aim to estimate $\lambda$ and ${\bm{\beta}}$ of model (2)
when $D$ is randomly selected between $0$ and $1$.

Let
$(D_{1},Y^{A}_{1},Y^{B}_{1},{{\bm{x}}}^{A}_{1},{{\bm{x}}}^{B}_{1}),\cdots,(D_{n},Y^{A}_{n},Y^{B}_{n},{{\bm{x}}}^{A}_{n},{{\bm{x}}}^{B}_{n})$
be $n$ observations with

$Y^{A}_{i}=I({\bm{\beta}}^{{}^{\prime}}{{\bm{x}}}_{i}^{A}+\tau_{i}+\lambda
D_{i}+\epsilon^{A}_{i}>0),\mbox{~{}and}\;\;Y_{i}^{B}=I({\bm{\beta}}^{{}^{\prime}}{{\bm{x}}}_{i}^{B}+\tau_{i}+\lambda(1-D_{i})+\epsilon^{B}_{i}>0)$
(3)

where ${\bm{\beta}}$ and $\lambda$ are $k\times 1$ regression coefficients and
treatment parameters respectively, $\tau_{i}$ is unobserved grouped effects
and independent of $\epsilon^{A}_{i}$ and $\epsilon^{B}_{i}$, and
$\epsilon^{A}_{i}$ and $\epsilon^{B}_{i}$ are independently normally
distributed with means $0$ and variances $1$.

When there are no covariates, estimating the treatment effects of $D$ is of
the major interest and one usually considers the average treatment effect
(ATE), i.e.,

$ATE=E[Y^{A}D+Y^{B}(1-D)-(Y^{A}(1-D)+Y^{B}D)].$

The average treatment effect (ATE) can be simply estimated by the sample
$\\{Y_{i}^{A},{\bf x}^{A}_{i},D_{i}\\}_{i=1}^{n}$, $\\{Y_{i}^{B},{\bf
x}^{B}_{i},1-D_{i}\\}_{i=1}^{n}$, i.e.,

$\widehat{ATE}=\frac{1}{n}\sum\limits_{i}^{n}\Bigg{(}Y_{i}^{A}D_{i}+Y_{i}^{B}(1-D_{i})-(Y_{i}^{A}(1-D_{i})+Y_{i}^{B}D_{i})\Bigg{)}.$

When the response is a linear function of $D$ or individuals are homogeneous,
the average treatment effect is an efficient estimator for measuring the
treatment effect of $D$. Under the Probit model (3), it is noteworthy that the
treatment effect is unidentifiable if individuals demonstrate substantial
heterogeneity. For example,

$\displaystyle ATE$ $\displaystyle=$ $\displaystyle
E[Y^{A}D+Y^{B}(1-D)-(Y^{A}(1-D)+Y^{B}D)]$ $\displaystyle=$ $\displaystyle
E[Y^{A}\mid D=1]p(D=1)+E[Y^{B}\mid D=0]p(D=0)$ $\displaystyle-E[Y^{B}\mid
D=1]p(D=1)-E[Y^{A}\mid D=0]p(D=0)$ $\displaystyle=$
$\displaystyle\int\Phi(\lambda+\tau)f(\tau)d\tau-\int\Phi(\tau)f(\tau)d\tau=\Phi\left(\frac{\lambda-\mu_{\tau}}{\sqrt{1+\sigma_{\tau}^{2}}}\right)-\Phi\left(\frac{-\mu_{\tau}}{\sqrt{1+\sigma_{\tau}^{2}}}\right)$

when $\tau$ is normally distributed with mean $\mu_{\tau}$ and standard
deviation $\sigma_{\tau}$, or

$\displaystyle ATE$ $\displaystyle=$ $\displaystyle
E[Y^{A}D+Y^{B}(1-D)-(Y^{A}(1-D)+Y^{B}D)]$ $\displaystyle=$ $\displaystyle
E[Y^{A}\mid D=1]p(D=1)+E[Y^{B}\mid D=0]p(D=0)$ $\displaystyle-E[Y^{B}\mid
D=1]p(D=1)-E[Y^{A}\mid D=0]p(D=0)$ $\displaystyle=$
$\displaystyle\int\Phi(\lambda+\tau)f(\tau)d\tau-\int\Phi(\tau)f(\tau)d\tau$
$\displaystyle=$ $\displaystyle
p\left[\Phi\left(\frac{\lambda-\mu_{\tau_{1}}}{\sqrt{1+\sigma_{\tau_{1}}^{2}}}\right)-\Phi\left(\frac{-\mu_{\tau_{1}}}{\sqrt{1+\sigma_{\tau_{1}}^{2}}}\right)\right]+(1-p)\left[\Phi\left(\frac{\lambda-\mu_{\tau_{2}}}{\sqrt{1+\sigma_{\tau_{2}}^{2}}}\right)-\Phi\left(\frac{-\mu_{\tau_{2}}}{\sqrt{1+\sigma_{\tau_{2}}^{2}}}\right)\right]$

when $\tau$ is a mixture of two normal variates, e.g.,
$p*N(\mu_{1},\sigma_{\tau_{1}}^{2})+(1-p)*N(\mu_{2},\sigma_{\tau_{2}}^{2})$.
The average treatment effect is close to zero and unidentifiable if the mean
or variance of the group effect $\tau$ is large.

In practical, collected data specially for observational studies exist
heterogeneous among groups. Following Gao et al. (2017), we propose new
estimators for the treatment effect $\lambda$ and ${\bm{\beta}}$ in model (3)
which can overcome the unidentifiable problem when the variance of individual
effects (i.e., $\tau$) is large. Under the condition that $D_{i}$ is randomly
selected between $0$ and $1$ and is independent of individual effect
$\tau_{i}$, the following estimation procedure for $\lambda$ and
${\bm{\beta}}$ is proposed

$(\hat{{\bm{\beta}}},\hat{\lambda})=\arg\max\prod_{i=1}^{n}\left\\{p_{i}^{Y^{A}_{i}(1-Y_{i}^{B})}[1-p_{i}]^{(1-Y^{A}_{i})Y_{i}^{B}}\right\\},$
(4)

where

$p_{i}=p({{\bm{\beta}},\lambda};{\bm{x}}_{i}^{A},{\bm{x}}_{i}^{B},D_{i})=\frac{G(\lambda(1-2D_{i})+{\bm{\beta}}^{{}^{\prime}}({\bm{x}}^{B}_{i}-{\bm{x}}^{A}_{i}))}{G(\lambda(1-2D_{i})+{\bm{\beta}}^{{}^{\prime}}({\bm{x}}^{B}_{i}-{\bm{x}}^{A}_{i}))+G(-\lambda(1-2D_{i})-{\bm{\beta}}^{{}^{\prime}}({\bm{x}}^{B}_{i}-{\bm{x}}^{A}_{i}))}$

and

$G(x)=-\sqrt{\pi}x\Phi(-\frac{x}{\sqrt{2}})+\exp\\{-\frac{x^{2}}{4}\\}.$

Let

$L({\bm{\beta}},\lambda)=\prod_{i=1}^{n}\left\\{p_{i}^{Y^{A}_{i}(1-Y_{i}^{B})}[1-p_{i}]^{(1-Y^{A}_{i})Y_{i}^{B}}\right\\}.$
(5)

Here, (5) can be regarded as a conditional likelihood functions with
$Y^{A}_{i}+Y^{B}_{i}=1$ while (4) can be viewed as a conditional maximum
likelihood estimation problem. The identifiable problem for (4) is reported in
the following theorem and its proof can be obtained by similar method given by
Gao et al. (2017).

Theorem 1. There exists an unique solution (i.e.,
$(\hat{{\bm{\beta}}},\hat{\lambda}$)) to (4) if the following conditions hold

(a) the rank of

$\left\\{({\bm{x}}^{B}_{1}-{\bm{x}}^{A}_{1})I(Y^{A}_{1}+Y^{B}_{1}=1),\cdots,({\bm{x}}^{B}_{n}-{\bm{x}}^{A}_{n})I(Y^{A}_{n}+Y^{B}_{n}=1)\right\\}$

is $k$ (i.e., the number of covariates), and

(b) there exists some $j$ such that $D_{j}=0$ and
${\bm{x}}^{B}_{j}-{\bm{x}}^{A}_{j}=\sum\limits_{i\not=j}d_{i}({\bm{x}}^{B}_{i}-{\bm{x}}^{A}_{i})I(D_{i}=0)$
where $d_{1},\cdots,d_{j-1},d_{j+1},\cdots,d_{n}$ are non-positive numbers.

If individuals $A$ and $B$ are sampled from the same population, condition of
(b) given in Theorem 1 is satisfied with probability close to $1$ for large
$n$. Condition (b) given in Theorem 1 can also be rewritten as

(b${}^{{}^{\prime}}$) there exists some $j$ such that $D_{j}=1$ and
${\bm{x}}^{B}_{j}-{\bm{x}}^{A}_{j}=\sum\limits_{i\not=j}d_{i}({\bm{x}}^{B}_{i}-{\bm{x}}^{A}_{i})I(D_{i}=1)$
where $d_{1},\cdots,d_{j-1},d_{j+1},\cdots,d_{n}$ are non-positive numbers.

Gao et al (2017) developed estimators for the dynamic parameters in dynamic
models and they are consistent and asymptotically normal distributed. We will
show that $\hat{{\bm{\beta}}}$ and $\hat{\lambda}$ are also consistent and
asymptotically normal distributed under some regular conditions. For the sake
of simplicity, denote
$\hat{{\bm{\theta}}}=(\hat{{\bm{\beta}}}^{{}^{\prime}},\hat{\lambda})^{{}^{\prime}}$
and ${\bm{\theta}}=({\bm{\beta}}^{{}^{\prime}},\lambda)^{{}^{\prime}}$.

Assumption: The individual effect is continuously distributed with the density
function $f(x)$ which satisfies

$f(x)=\frac{1}{\sigma_{\tau}}h\left(\dfrac{x-\mu_{\tau}}{\sigma_{\tau}}\right).$
(6)

Under Assumption (6), $\hat{{\bm{\theta}}}$ is consistent and asymptotically
normally distributed and its proof is given in Appendix A.

Theorem 2. (Consistency) Under Assumption (6) with $\sigma_{\tau}=an^{\alpha}$
$(0<\alpha<1)$ for some constant $a>0$, $\hat{{\bm{\theta}}}$ converges to
${\bm{\theta}}$ in probability $1$ as $n\rightarrow\infty$.

Let

$\displaystyle K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)$ $\displaystyle=$
$\displaystyle\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}+\lambda
D+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}-\lambda(1-D)-u)du$
$\displaystyle+\int\Phi(-{\bm{\beta}}^{{}^{\prime}}x^{A}-\lambda
D-u)\Phi({\bm{\beta}}^{{}^{\prime}}x^{B}+\lambda(1-D)+u)du.$

Noting that

$p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})=\frac{\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)du}{K({\bm{\theta}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}$

and

$1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})=\frac{\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)du}{K({\bm{\theta}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})},$

we have the following result

Theorem 3. (Asymptotic Normality) Let $h(x)$ be continuously derivative and
$\sigma_{\tau}=a\sqrt{n}$($a>0$). Suppose that the covariance matrix of
${\bm{x}}^{A}-{\bm{x}}^{B}$ is positive definite, the conditional distribution
of ${\bm{x}}^{A}-{\bm{x}}^{B}$ given $D=0$ is symmetrical at the origin and
$0<P[D=1\mid x]<1$. Under Assumption (6) , we have
$k_{n}(\hat{{\bm{\theta}}}-{{\bm{\theta}}})$ is asymptotically distributed as
$N(0,c\Sigma^{-1})$ as $n\rightarrow\infty$, where

$k_{n}=\left\\{\sum\limits_{i=1}^{n}I{(Y^{A}_{i}+Y^{B}_{i}=1)}\right\\}^{1/2}$

and

$\Sigma=E\left\\{\frac{K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}}\Big{)}\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}\partial}\Big{)}^{{}^{\prime}}}{p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)\times(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))}\right\\},\mbox{~{}and}\;\;c=E[K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)].$

##  3 Simulations

In this section, we evaluate the finite sample performances of the treatment
effects estimator $\hat{\lambda}$ which is proposed in Section 2. Several
scenarios will be considered.

In the first simulation study, we mainly consider models $(\ref{ME})$ without
covariates. In this case, the error terms, i.e., $\epsilon_{i}^{A}$ and
$\epsilon_{i}^{B}$, are assumed to be i.i.d. normally distributed with means
being 0 and variances being 1 and $D_{i}$ follows binomial distribution with
success probability being $p=2/3$. We assume that the individual effects come
from different distributions (e.g., uniform, normal, Student t and Cauchy
distributions). We obtain the mean bias (BIAS) and root mean squared errors
(RMSEs) of the estimates of different $\lambda$ for different distributions of
$\tau$ with sample sizes being 1000 or 5000. The results based on 100
repetitions are reported in Table 1. From Table 1, we observe that our
proposed estimator is robust to different prior distributions for $\tau$. In
particular, it is noteworthy that our proposed estimator works well for long-
tailed distribution (i.e., Student-t) and distribution with no moments (i.e.,
t(1)(Cauchy)).

Table 1: Simulated RMSEs and BIAS of the new estimator for treatment effect $\lambda$ in model $(\ref{ME})$ without covariates based on 100 replications. sample size $\tau_{i}$ | parameter | U(-4,4) |  | N(0,4) |  | t(1) |
---|---|---|---|---|---|---|---
| $\lambda$ | BIAS | RMSE | BIAS | RMSE | BIAS | RMSE
n=1000 |  | 2
---
1.5
1
0.5
0
-0.5
-1
-1.5
-2
| 0.022
---
0.011
0.018
0.002
-0.018
-0.006
-0.053
-0.030
-0.026
| 0.184
---
0.150
0.126
0.103
0.099
0.091
0.122
0.158
0.242
| 0.019
---
0.000
-0.016
-0.010
-0.011
0.004
-0.019
-0.003
-0.053
| 0.256
---
0.178
0.124
0.126
0.114
0.107
0.144
0.171
0.254
| 0.026
---
-0.013
-0.046
-0.012
-0.008
0.022
0.038
0.025
-0.025
| 0.236
---
0.143
0.109
0.093
0.092
0.082
0.109
0.128
0.206
sample size $\tau_{i}$ | $\lambda$ | U(-10,10) |  | N(0,25) |  | t(1) |
n=5000 |  | 2
---
1.5
1
0.5
0
-0.5
-1
-1.5
-2
| 0.038
---
-0.005
0.006
0.001
0.003
0.013
-0.002
-0.006
-0.030
| 0.180
---
0.114
0.079
0.062
0.062
0.071
0.071
0.097
0.151
| 0.039
---
0.055
0.028
0.024
-0.011
0.014
-0.015
-0.011
-0.056
| 0.283
---
0.202
0.160
0.140
0.123
0.115
0.136
0.228
0.271
| -0.008
---
-0.008
-0.035
-0.005
0.003
0.018
0.029
0.014
0.011
| 0.075
---
0.066
0.054
0.037
0.033
0.037
0.054
0.038
0.073

Next, we consider models $(\ref{ME})$ without covariates and with the
individual effect $\tau_{i}\sim N(0,\pi^{2}/3)$ where $\pi$ = 3.14. The
results for standard error (SE), BIAS and RMSE of our proposed estimate for
various $\lambda$-values based on sample sizes being 200, 500, 1000, 2000,
3000 and 5000 are reported in Table 2. We find that as the sample size
increases, all BIAS, RMSE and SE for $\hat{\lambda}$ decrease, which is
consistent with Theorem 2 that $\hat{\lambda}$ is a consistent estimator of
$\lambda$. Besides, we observe that both SE and RMSE are almost identical. In
addition, we consider model $(\ref{ME})$ with covariates in which $x_{i}^{A}$
being sampled from N(0, 1) and $x_{i}^{B}=x_{i}^{A}+N(0,1)$. Mean BIASs and
RMSEs of the estimates of $\lambda$ and $\beta$ based on 100 repetitions for
different values of $\lambda$ and $\beta$ with sample size being 1000 are
reported in Table 3. We find that the biases for both $\hat{\lambda}$ and
$\hat{\beta}$ are close to zero, which again confirms the result of Theorem 2
that $\hat{{\bm{\theta}}}$ is a consistent estimator of ${{\bm{\theta}}}$.

Table 2: Simulated RMSEs, BIAS and SE of the new estimator for treatment effect $\lambda$ in models $(\ref{ME})$ without covariates and $\tau_{i}\ \sim\ N(0,\pi^{2}/3)$ based on 100 replications. sample size | $\lambda$ | BIAS | SE | RSME | sample size | $\lambda$ | BIAS | SE | RMSE
---|---|---|---|---|---|---|---|---|---
| 200
---
| 1
---
0.5
0
-0.5
-1
| 0.015
---
0.003
0.003
-0.029
-0.029
| 0.295
---
0.259
0.236
0.248
0.293
| 0.294
---
0.258
0.235
0.248
0.293
2000 |  | 1
---
0.5
0
-0.5
-1
| -0.007
---
0.004
-0.006
0.009
0.011
| 0.088
---
0.068
0.069
0.076
0.087
| 0.088
---
0.068
0.069
0.076
0.087
500 |  | 1
---
0.5
0
-0.5
1
| 0.016
---
0.001
-0.018
0.001
-0.029
| 0.164
---
0.133
0.127
0.143
0.177
| 0.164
---
0.132
0.128
0.143
0.178
3000 |  | 1
---
0.5
0
-0.5
-1
| -0.008
---
0.002
-0.004
-0.001
0.009
| 0.070
---
0.059
-0.004
0.062
0.077
| 0.070
---
0.059
0.059
0.062
0.077
1000 |  | 1
---
0.5
0
-0.5
-1
| -0.005
---
-0.016
-0.010
0.003
0.001
| 0.116
---
0.101
0.093
0.105
0.126
| 0.115
---
0.102
0.093
0.105
0.125
5000 |  | 1
---
0.5
0
-0.5
-1
| -0.006
---
-0.001
0.004
-0.001
0.004
| 0.054
---
0.041
0.042
0.047
0.053
| 0.054
---
0.041
0.042
0.047
0.053
Table 3: Simulated RMSEs and BIAS of the new estimator of the parameter in models $(\ref{ME})$ with covariates, sample size n=1000, $\tau_{i}\ \sim\ N(0,\pi^{2}/3)$ based on 100 replications. $\lambda$ | ${\bm{\beta}}$ | BIAS$(\hat{\lambda})$ | RMSE$(\hat{\lambda})$ | BIAS$(\hat{{\bm{\beta}}})$ | RMSE$(\hat{{\bm{\beta}}})$ | $\lambda$ | ${\bm{\beta}}$ | BIAS$(\hat{\lambda})$ | RMSE$(\hat{\lambda})$ | BIAS$(\hat{{\bm{\beta}}})$ | RMSE$(\hat{{\bm{\beta}}})$
---|---|---|---|---|---|---|---|---|---|---|---
| 1
---
0.5
0
-0.5
-1
| 1
---
0.5
0
-0.5
-1
| 0.029
---
0.011
-0.001
-0.016
-0.015
| 0.029
---
0.011
-0.001
-0.016
-0.015
| 0.025
---
-0.011
0.041
-0.025
-0.027
| 0.170
---
0.108
0.105
0.109
0.179
| 1
---
0.5
-0.5
-1
| -1
---
\- 0.5
0.5
1
| 0.000
---
-0.017
0.010
-0.015
| 0.177
---
0.118
0.117
0.193
| -0.019
---
0.001
0.011
0.042
| 0.172
---
0.095
0.093
0.186
| 0
---
0
0
0
| 1
---
0.5
-0.5
-1
| 0.005
---
0.088
0.015
-0.024
| 0.133
---
0.103
0.102
0.138
| 0.016
---
0.020
0.001
-0.014
| 0.165
---
0.088
0.094
0.145
| 1
---
0.5
-0.5
-1
| 0
---
0
0
0
| 0.024
---
0.016
0.001
-0.005
| 0.145
---
0.109
0.105
0.126
| 0.016
---
0.008
-0.007
-0.014
| 0.097
---
0.081
0.078
0.094

For normally distributed group effects with mean being 0 and variance being
$\sigma^{2}$ in models $(\ref{ME})$, Heckman (1978) proposed the maximum
likelihood estimators for the parameter $\lambda$ and $\sigma^{2}$ (denoted as
$\hat{\lambda}_{H}$ and $\hat{\sigma}_{H}$, respectively). Here, we compare
our proposed estimators with Heckman estimators for the random effect in
models $(\ref{ME})$ without covariates. For this purpose, we consider the
group effects coming from N(0,1), N(0,4), and symmetrical mixed normal
distribution $0.5*N(-6,9)+0.5*N(6,9)$. BIASs and RMSEs of the estimates (based
on 100 repetitions) for different $\lambda$ values with sample size being 1000
are reported in Table 4. We also consider group effect coming from asymmetric
mixed normal distribution $0.5*N(-6,9)+0.5*N(6,3)$, $0.5*N(-6,3)+0.5*N(6,9)$
and the results are reported in Table 4. From Table 4, we find that our
proposed estimator and Heckmans estimator are comparable. However, for
asymmetrically distributed group effects, our proposed estimator is
significantly better that Heckmans estimator in terms of smaller BIASs and
RMSEs.

Table 4: Comparison of RMSEs and BIAS of our proposed estimator $\hat{\lambda}_{new}$ and Heckmans estimator $\hat{\lambda}_{H}$ for various group effect distribution with sample size being 1000 based on 100 repetitions. Distribution of the $\tau_{i}$ | $\lambda$ | BIAS$(\hat{\lambda}_{new})$ | RMSE$(\hat{\lambda}_{new})$ | BIAS$(\hat{\lambda}_{H})$ | RMSE$(\hat{\lambda}_{H})$ | BIAS$(\hat{\sigma}_{H})$ | RMSE$(\hat{\sigma}_{H})$
---|---|---|---|---|---|---|---
N(0,1) |  | 1
---
0.5
0
-0.5
-1
| -0.027
---
-0.025
0.000
0.021
0.031
| 0.092
---
0.061
0.059
0.083
0.088
| 0.018
---
0.008
0.003
0.004
0.002
| 0.077
---
0.064
0.051
0.058
0.070
| 0.011
---
0.009
0.005
0.002
0.002
| 0.084
---
0.084
0.082
0.082
0.095
N(0,4) |  | 1
---
0.5
0
-0.5
-1
| -0.016
---
0.010
-0.011
0.004
-0.019
| 0.124
---
0.126
0.114
0.107
0.144
| -0.005
---
0.010
-0.008
0.018
-0.011
| 0.138
---
0.121
0.121
0.105
0.138
| 0.057
---
0.072
0.062
0.078
0.089
| 0.470
---
0.472
0.337
0.431
0.471
| 0.5*N(-6,9)
---
+0.5*N(6,9)
| 1
---
0.5
0
-0.5
-1
| 0.043
---
0.014
0.002
-0.014
-0.014
| 0.274
---
0.174
0.186
0.200
0.251
| 0.036
---
-0.029
0.010
-0.015
-0.093
| 0.251
---
0.206
0.185
0.202
0.285
| 1.033
---
1.083
0.756
0.943
1.173
| 2.462
---
2.426
2.110
2.034
2.850
| 0.5*N(-6,9)
---
+0.5*N(6,3)
| 1
---
0.5
0
-0.5
-1
| 0.052
---
0.011
-0.001
-0.010
-0.051
| 0.303
---
0.221
0.206
0.211
0.314
| 0.327
---
0.303
0.273
0.232
0.201
| 0.427
---
0.387
0.340
0.251
0.323
| 5.688
---
5.511
5.530
5.122
5.325
| 6.262
---
6.162
6.234
5.776
6.167
| 0.5*N(-6,3)
---
+0.5*N(6,9)
| 1
---
0.5
0
-0.5
-1
| 0.021
---
0.018
0.024
-0.010
-0.023
| 0.253
---
0.190
0.175
0.239
0.283
| -0.206
---
-0.232
-0.271
-0.263
-0.380
| 0.350
---
0.302
0.344
0.346
0.481
| 4.930
---
5.882
5.379
5.022
5.757
| 5.862
---
6.584
5.592
5.590
6.645
Table 5: Compared the proposed method with the CML when the error term is standard normal distribution. With the distribution of individual effect $\tau_{i}$ distribution is N(0,1),N(0,$\pi^{2}/3$), $N(0,6)$, $0.5*N(-4,6)+0.5*N(4,6)$ for sample n=500 based on 100 repetitions. |  | New method | CML
---|---|---|---
Distribution of the $\tau_{i}$ | $\lambda$ | BIAS$(\hat{\lambda}_{new})$ | RMSE$(\hat{\lambda}_{new})$ | BIAS$(\hat{\lambda}_{CML})$ | RMSE$(\hat{\lambda}_{CML})$
$N(0,1)$ |  | 1
---
0.5
0
-0.5
-1
| -0.011
---
-0.006
0.004
-0.002
-0.001
| 0.119
---
0.102
0.083
0.110
0.118
| 0.194
---
0.077
-0.002
-0.088
-0.205
| 0.227
---
0.131
0.122
0.134
0.236
$N(0,\pi^{2}/3)$ |  | 1
---
0.5
0
-0.5
-1
| 0.011
---
0.006
0.001
0.002
-0.029
| 0.184
---
0.140
0.133
0.139
0.207
| -0.504
---
-0.239
-0.015
0.238
0.500
| 0.510
---
0.248
0.070
0.247
0.505
$N(0,6)$ |  | 1
---
0.5
0
-0.5
-1
| 0.056
---
0.010
-0.007
0.008
-0.018
| 0.248
---
0.199
0.191
0.199
0.232
| -0.712
---
-0.350
0.008
0.359
0.716
| 0.714
---
0.354
0.046
0.363
0.718
$0.5*N(-4,6)+0.5*N(4,6)$ |  | 1
---
0.5
0
-0.5
-1
| 0.019
---
0.002
-0.005
-0.010
-0.004
| 0.269
---
0.219
0.202
0.222
0.249
| -0.763
---
-0.385
0.008
0.183
0.761
| 0.765
---
0.385
0.050
0.386
0.762

When the error term is standard logistic distribution for models $(\ref{ME})$,
Chamberlain (1980) proposed conditional maximum likelihood (CML) estimator for
the parameter $\lambda$ and Bartolucci and Pigini (2017) developed an R
package which is named as cquad to implement the CML approach. Here, we
compare the proposed method with the CML when the individual effects
$\tau_{i}$ is distributed by N(0,1), $N(0,\pi^{2}/3)$, $N(0,6)$,
$0.5*N(-4,6)+0.5*N(4,6)$, respectively, $\epsilon_{i}^{A}$ and
$\epsilon_{i}^{B}$ are independently standard normal distribution and $D_{i}$
follows binomial distribution with success probability being $p=2/3$.
Simulation results are listed in Table 5 with sample n=500 and 100
repetitions. Simulation studies show that the biases of the CML are large,
furthermore the biases are increasing with the increase of variance of
individual effects.

Finally, we compare the proposed estimator with the inverse probability
weighted estimator developed by Robins et al. (1994). Here simulated data are
generated by the following models:

$Y^{A}_{i}=I({\beta^{{}^{\prime}}{\bm{x}}}_{i}^{A}+\tau_{i}+\lambda
D_{i}+\epsilon^{A}_{i}>0),\mbox{~{}and}\;\;Y_{i}^{B}=I({\beta^{{}^{\prime}}{\bm{x}}}_{i}^{B}+\tau_{i}+\lambda(1-D_{i})+\epsilon^{B}_{i}>0).$
$D_{i}=I\\{0.75{\bm{x}}_{i}^{A}+0.25{\bm{x}}_{i}^{B}+\varepsilon_{i}>0\\}$

where $\beta=1$, ${\bm{x}}_{i}^{A}$ being sampled from the standard normal
distribution, ${\bm{x}}_{i}^{B}$ being sampled from the uniform distribution
$U(-1,1)$, and $\varepsilon_{i}$ following standard logistic distributions and
the group effects being normally distributed with mean 0 and variance
$\pi^{2}/3$. BIASs and RMSEs are reported in Table 6 with the sample size
$n=1000$ and $100$ repetitions. It is shown that both the proposed estimators
and Heckman’s estimators perform better than IPW estimators.

Table 6:  Compared the proposed estimator $\hat{\lambda}_{new}$ with IPW estimator $\widehat{TE}_{IPW}$ and Heckman estimator $\hat{\lambda}_{H}$ for group effects being distributed as $N(0,\pi^{2}/3)$ and sample size being 1000 based on 100 repetitions. $\lambda$ | BIAS$(\hat{\lambda}_{new})$ | RMSE$(\hat{\lambda}_{new})$ | BIAS$(\widehat{TE}_{IPW})$ | RMSE$(\widehat{TE}_{IPW})$ | BIAS$(\hat{\lambda}_{H})$ | RMSE$(\hat{\lambda}_{H})$ | RMSE$(\hat{\sigma}_{H})$
---|---|---|---|---|---|---|---
| 1
---
0.5
0
-0.5
-1
| 0.015
---
-0.007
-0.006
0.003
-0.022
| 0.159
---
0.136
0.119
0.147
0.163
| -0.883
---
-0.441
0.004
0.449
0.894
| 0.883
---
0.441
0.014
0.448
0.894
| 0.030
---
0.009
0.001
0.013
-0.047
| 0.152
---
0.110
0.101
0.120
0.143
| 0.410
---
0.354
0.354
0.330
0.385

##  4 Applications

In this section, we will demonstrate our methodologies with two real datasets
- one is about lead levels in children’s blood and the other is about acute
leukaemia patients from a clinical trial.

###  4.1  Lead Levels in Children’s Blood

The data set shown in the Table 7 are extracted from the observational study
being studied in Morton et al. (1982), Pruzek & Helmreich (2009). Briefly,
children of parents who had worked in a lead-unrelated industry where lead was
used in making batteries were matched by age, exposure to traffic, and
neighborhood with children whose parents did not work in lead-related
industries. We apply our proposed method to estimate the influence of
environmental sources on children of lead absorption. here, we define that $Y$
= 1 if lead level (ug/dl) exceeds 16, 0 otherwise. Our proposed method yields
$\hat{\lambda}=0.9992$ with $\hat{\sigma}^{2}_{\lambda}=0.1733$ while
Heckman’s method produces $\hat{\lambda}_{H}=1.2278$ with
$\hat{\sigma}^{2}_{H}=0.1257$. According to our proposed method, the odd of
treatment group and control group is given by
$P(Y^{A}=1\mid\tau)/P(Y^{B}=1\mid\tau)=\int\Phi(\tau+0.9992)f(\tau)d\tau/\int\Phi(\tau)f(\tau)d\tau=1.653$
. Hence, children of employees in a battery manufacturing plant have higher
prevalence of high level of blood lead than those children whose parents are
not employed in a lead-related industry.

Table 7: Blood lead levels of case and control children and the difference between blood lead levels of the two groups by matched pairs. Matched | Lead levels(ug/dl)
---|---
pair no | case | control | Difference
| 1
---
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
| 38
---
23
41
18
37
36
23
62
31
34
24
14
21
17
16
20
15
10
45
39
22
35
49
48
44
35
43
39
34
13
73
25
27
| 16
---
18
18
24
19
11
10
15
16
18
18
13
19
10
16
16
24
13
9
14
21
19
7
18
19
12
11
22
25
16
13
11
13
| 22
---
5
23
-6
18
25
13
47
15
16
6
1
2
7
0
4
-9
-3
36
25
1
16
42
30
25
23
32
17
9
-3
60
14
14
mean | 31.84848 | 15.87879 |

###  4.2  Acute leukaemia patients of clinical trial

The data reported in Table 8 are coming from a clinical trial about acute
leukaemia patients, see Freireich et al. (1963), Gehan (1965). Briefly,
treatment group which uses 6-mercaptopurine (6-MP) to therapy is compared to
control group which employs placebo to treatment in the maintenance of
remissions in acute leukemia. One year after the start of the study, times of
remission in weeks are recorded. We define $Y$ = 1 if the time of remission is
more than 12 weeks, 0 otherwise. Our proposed method yields
$\hat{\lambda}=0.617$ with $\hat{\sigma}^{2}_{\lambda}=0.1377$ while Heckman’s
method produces $\hat{\lambda}_{H}=0.180$ with $\hat{\sigma}^{2}_{H}=0.1390$.
Obviously, simple Wald tests will suggest significant treatment effect from
our proposed method and insignificant treatment effect from Heckman’s method.
According to our method, it can be concluded that the odd of treatment group
and control group is
$P(Y^{A}=1\mid\tau)/P(Y^{B}=1\mid\tau)=\int\Phi(\tau+0.617)f(\tau)d\tau/\int\Phi(\tau)f(\tau)d\tau=1.337$.
Consequently, there is strong evidence that patients receiving 6-MP have
longer remissions than those receiving placebo.

Table 8: Length of remission of 42 acute leukaemia patients with 21 patients being treated with the drug 6-mercaptopurine and 21 being assigned to control. Pair | Remission time(weeks) | Censoring | Group
---|---|---|---
| 1
---
1
2
2
3
3
4
4
5
5
6
6
7
7
8
8
9
9
10
10
11
11
12
12
13
13
14
14
15
15
16
16
17
17
18
18
19
19
20
20
21
21
| 1
---
10
22
7
3
32
12
23
8
22
17
6
2
16
11
34
8
32
12
25
2
11
5
20
4
19
15
6
8
17
23
35
5
6
11
13
4
9
1
6
8
10
| 1
---
1
1
1
1
0
1
1
1
1
1
1
1
1
1
0
1
0
1
0
1
0
1
0
1
0
1
1
1
0
1
0
1
1
1
1
1
0
1
0
1
0
| control
---
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP
control
6-MP

##  5 Conclusion

In this paper, we consider a matched pairs data design which allows for
unobserved individual heterogeneity (groups effects) with treatment effects on
the dependent variable. A new method for estimating treatment effects based on
the probit model is proposed. Here, group effects with large variances are
allowed. We demonstrate that our proposed parameter estimators are consistent
and asymptotically normally distributed. In general, our proposed estimators
are superior to existing Heckman’s estimators and inverse probability weighted
(IPW) estimators.

## Appendix. Proofs

We assume that $({\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(i=1,\cdots,n)$ are
mutually independent with a common distribution function
$F({\bm{x}}^{A},{\bm{x}}^{B},D)$. Let $\sigma_{\tau}=an^{\alpha}$ with $a>0$
and $0<\alpha<1$.

Lemma A1. Under Assumption (6) with $h(x)$ having continuous derivatives, we
have

$\displaystyle\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\xrightarrow{\
p\ }C_{1},$

where $\displaystyle
C_{1}=\frac{h(0)}{a}E[K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)]$.

Proof:

Letting $W_{i}=I(Y_{i}^{A}+Y_{i}^{B}=1)$ yields

$\displaystyle E(W_{i}|{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})$
$\displaystyle=$
$\displaystyle\bigg{[}\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)\frac{1}{\sigma_{u}}h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du$
$\displaystyle+\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)\frac{1}{\sigma_{\tau}}h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{1}{an^{\alpha}}\times\bigg{[}\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du$
$\displaystyle+\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{h(0)}{an^{\alpha}}\times\bigg{[}\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)du$
$\displaystyle+\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)du+O(\sigma_{\tau}^{-1})\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{h(0)}{an^{\alpha}}[K({\bm{\theta}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})+O(n^{-\alpha})]$

and

$E(W_{i})=\frac{C_{1}}{n^{\alpha}}+O(n^{2\alpha}),\;\;\frac{1}{n^{1-\alpha}}\sum_{i=1}^{n}E(W_{i})=C_{1}+O(n^{-\alpha})=C_{1}+o(1).$

Therefore,

$E(W_{i})^{2}=EW_{i}=\frac{C_{1}}{n^{\alpha}}+O(n^{-2\alpha}),\;\;Var(W_{i})=E(W_{i})^{2}-(E(W_{i}))^{2}=O(n^{-\alpha})$

and

$Var(\frac{1}{n^{1-\alpha}}\sum_{i=1}^{n}W_{i})^{2}=\frac{1}{n^{2-2\alpha}}nVar(W_{i})^{2}=\frac{1}{n^{2-2\alpha}}n(\frac{C_{1}}{n^{\alpha}}+O(n^{-2\alpha}))=O(n^{1-\alpha}).$

By Chebyshev’s inequality, we have

$\displaystyle
P\bigg{(}\bigg{|}\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}W_{i}-E(\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}W_{i})\bigg{|}\geq\varepsilon\bigg{)}$
$\displaystyle=$ $\displaystyle
P\bigg{(}\bigg{|}\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}W_{i}-(C_{1}+o(1))\bigg{|}\geq\varepsilon\bigg{)}$
$\displaystyle\leq$
$\displaystyle\frac{Var\left(\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}W_{i}\right)^{2}}{\varepsilon^{2}}\longrightarrow
0,$

which implies Lemma A1.

Lemma A2. Under Assumption (6), we have

$\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\left\\{{Y_{i}^{A}\log
p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}+{Y_{i}^{B}\log(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\right\\}\xrightarrow{p}C_{2},$

where

$\displaystyle C_{2}({{\bm{\theta}}}^{*})$ $\displaystyle=$
$\displaystyle\frac{h(0)}{a}E\bigg{\\{}K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)\bigg{[}p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)\log
p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D)$
$\displaystyle+(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))\log(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D))\bigg{]}\bigg{\\}}.$

Proof: Let

$U_{i}=I(Y_{i}^{A}+Y_{i}^{B}=1)\left\\{{Y_{i}^{A}\log
p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}+{Y_{i}^{B}\log(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\right\\}.$

We have

$\displaystyle E(U_{i})$ $\displaystyle=$ $\displaystyle E\left[E(U_{i}\mid
x^{A}_{i},x_{i}^{B},D_{i})\right]$ $\displaystyle=$
$\displaystyle\dfrac{h(0)}{an^{\alpha}}\Big{\\{}E\Big{[}\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}+\lambda
D+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}-\lambda(1-D)-\tau)du\log
p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D)$
$\displaystyle+{\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}-\lambda(1-D)-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}+\lambda
D+u)du}\log(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D))\Big{]}$
$\displaystyle\;\;\;+o(1)\Big{\\}}$ $\displaystyle=$
$\displaystyle\frac{1}{n^{\alpha}}\\{C_{2}({{\bm{\theta}}}^{*})+o(1)\\}$

and

$\frac{1}{n^{1-\alpha}}E\left(\sum\limits_{i=1}^{n}U_{i}\right)=\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}\frac{1}{n^{\alpha}}\\{C_{2}({{\bm{\theta}}}^{*})+o(1)\\}=C_{2}({{\bm{\theta}}}^{*})+o(1).$

Similarly, one can obtain

$\displaystyle E(U^{2}_{i})$ $\displaystyle=$ $\displaystyle
E\left[E(U^{2}_{i}|{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})\right]$
$\displaystyle=$
$\displaystyle\dfrac{h(0)}{an^{\alpha}}\Big{\\{}E\Big{[}{\log^{2}p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)\int\Phi({{\bm{\beta}}}^{{}^{\prime}}{\bm{x}}^{A}+\lambda
D+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}-\lambda(1-D)-u)du}$
$\displaystyle+{\log^{2}(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}-\lambda
D-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}+\lambda(1-D)+u)du}\Big{]}+o(1)\Big{\\}}$
$\displaystyle=$ $\displaystyle O(n^{-\alpha}),$
$Var(U_{i})=E(U^{2}_{i})-[E(U_{i})]^{2}=O(n^{-\alpha}),\mbox{~{}and}$

$Var\left(\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i}\right)=\frac{1}{n^{2(1-\alpha)}}\sum\limits_{i=1}^{n}Var(U_{i})=O(n^{1-\alpha}).$

For $\varepsilon>0$, as $n\rightarrow\infty$ we have

$\displaystyle
P\bigg{(}\bigg{|}\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i}-E(\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i})\bigg{|}\geq\varepsilon\bigg{)}$
$\displaystyle=$ $\displaystyle
P\bigg{(}\bigg{|}\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i}-(C_{2}+o(1))\bigg{|}\geq\varepsilon\bigg{)}$
$\displaystyle\leq$
$\displaystyle\frac{Var\left(\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i}\right)}{\varepsilon^{2}}\longrightarrow
0$

which implies that Lemma A2 holds.

Proof of Theorem 2: Let

$L({{\bm{\theta}}}^{*})=\prod\limits_{i=1}^{n}p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})^{Y_{i}^{A}(1-Y_{i}^{B})}(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))^{Y_{i}^{B}{(1-Y_{i}^{A})}}$

and its logarithm be

$l({{\bm{\theta}}}^{*})=\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\left\\{{Y_{i}^{A}\log
p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}+{Y_{i}^{B}\log(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\right\\}.$

By Lemmas A1 and A2, we have

$\frac{1}{n^{1-\alpha}}l({{\bm{\theta}}}^{*})=\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}U_{i}\xrightarrow{\
p\ }C_{2}({{\bm{\theta}}}^{*})$

and

$\frac{1}{n^{1-\alpha}}\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\xrightarrow{\
p\ }C_{1}.$

Thus,

$\frac{1}{\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)}l({{\bm{\theta}}}^{*})\xrightarrow{\
p\ }\frac{C_{2}({{\bm{\theta}}}^{*})}{C_{1}}.$

Since

$p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)\log\frac{p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D)}+(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))\log\frac{(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))}{(1-p({{\bm{\theta}}}^{*};{\bm{x}}^{A},{\bm{x}}^{B},D))}\geq
0,$

we have

$P\left\\{\frac{1}{\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)}l({{\bm{\theta}}}^{*})\leq\frac{1}{\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)}l({{\bm{\theta}}})\right\\}\xrightarrow{\
p\ }1,$

which implies that $\hat{{\bm{\theta}}}$ is a consistent estimator of
${{\bm{\theta}}}$.

Lemma A3. Under Assumption (6) with $h(x)$ having continuous derivatives and
$\sigma_{\tau}=a\sqrt{n}(a>0)$,

$\sum\limits_{i=1}^{n}\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)}{\sqrt{\sum\limits_{j=1}^{n}I(Y_{j}^{A}+Y_{j}^{B}=1)}}\left\\{\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\right\\}$

has the asymptotic distribution $N(0,\Sigma/c)$, where

$\Sigma=\displaystyle
E\left\\{\frac{K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}}\Big{)}\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}}\Big{)}^{{}^{\prime}}}{p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)(1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D))}\right\\},\;\;c=E[K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)].$

Proof:  Let

$\displaystyle Q_{i}$ $\displaystyle=$ $\displaystyle
I(Y_{i}^{A}+Y_{i}^{B}=1)\Big{\\{}\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D)}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\Big{\\}}$
$\displaystyle=$
$\displaystyle\Big{\\{}\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)Y_{i}^{A}}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}-\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)Y_{i}^{B}}{1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}\Big{\\}}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}.$

Since the first derivative of $h(x)$ is continuous, we have

$\displaystyle
E\left(\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)Y_{i}^{A}}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}\mid{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}\right)$
$\displaystyle=\frac{\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du}{\sigma_{\tau}\times
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}$
$\displaystyle=\frac{\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)[h(0)+h^{{}^{\prime}}(0)\frac{u-\mu_{\tau}}{\sigma_{\tau}}+o(\sigma_{\tau}^{-1})]du}{\sigma_{\tau}\times
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}$ $\displaystyle=$
$\displaystyle\frac{h(0)\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)du}{\sigma_{\tau}\times
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}+O(\sigma^{-2}_{\tau})+o(\sigma^{-2}_{\tau})$

and

$\displaystyle
E\left(\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)Y_{i}^{B}}{1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}\mid{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}\right)$
$\displaystyle=\frac{\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)h(\frac{u-\mu_{\tau}}{\sigma_{\tau}})du}{\sigma_{\tau}\times[1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})]}$
$\displaystyle=\frac{\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)[h(0)+h^{{}^{\prime}}(0)\frac{u-\mu_{\tau}}{\sigma_{\tau}}+o(\sigma_{\tau}^{-1})]du}{\sigma_{\tau}\times[1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}]$
$\displaystyle=$
$\displaystyle\frac{h(0)\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)du}{\sigma_{\tau}\times[1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})]}+O(\sigma^{-2}_{\tau})+o(\sigma^{-2}_{\tau}).$

As

$\frac{\int\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}+\lambda
D_{i}+u)\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}-\lambda(1-D_{i})-u)du}{\int\Phi(-{\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{A}_{i}-\lambda
D_{i}-u)\Phi({\bm{\beta}}^{{}^{\prime}}{\bm{x}}^{B}_{i}+\lambda(1-D_{i})+u)du}=\frac{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})},$

we have

$E(Q_{i})=E(E(Q_{i}\mid{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))=O(\sigma^{-2}_{\tau})+o(\sigma^{-2}_{\tau})=O(\sigma^{-2}_{\tau})=O(n^{-1}).$

Similarly,

$\displaystyle E[Q_{i}Q_{i}^{{}^{\prime}}]$ $\displaystyle=$
$\displaystyle\frac{h(0)}{a\sqrt{n}}\times
E\left\\{\frac{K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}}\Big{)}\Big{(}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}{\partial{{\bm{\theta}}}}\Big{)}^{{}^{\prime}}}{[1-p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)]p({{\bm{\theta}}};{\bm{x}}^{A},{\bm{x}}^{B},D)}\right\\}$
$\displaystyle+o(n^{-1/2})$ $\displaystyle=$
$\displaystyle\frac{h(0)}{a\sqrt{n}}\times\Sigma+o(n^{-1/2}).$

Let ${\bf{t}}=(t_{1},\cdots,t_{k+1})^{{}^{\prime}}$ be a k+1 dimensional row
vector. We have

$\displaystyle E[\exp\\{i{\bf{t}}^{{}^{\prime}}Q_{j}/n^{1/4}\\}]$
$\displaystyle=$ $\displaystyle
1+\frac{1}{n^{1/4}}E(i{\bf{t}}^{{}^{\prime}}Q_{j})-\frac{1}{2\sqrt{n}}E(-{\bf{t}}^{{}^{\prime}}Q_{j}Q^{{}^{\prime}}_{j}{\bf{t}})+E\Big{[}o\Big{(}\frac{{\bf{t}}^{{}^{\prime}}Q_{j}Q^{{}^{\prime}}_{j}{\bf{t}}}{n^{1/2}}\Big{)}\Big{]}$
$\displaystyle=$ $\displaystyle
1-\frac{1}{2\sqrt{n}}E({\bf{t}}^{{}^{\prime}}Q_{j}Q^{{}^{\prime}}_{j}{\bf{t}})+o(n^{-1})$
$\displaystyle=$ $\displaystyle
1-\frac{h(0){\bf{t}}^{{}^{\prime}}\Sigma{\bf{t}}}{2an}+o(n^{-1}).$

The characteristic function of $\sum\limits_{j=1}^{n}Q_{j}/n^{1/4}$ is given
by

$\displaystyle\varphi_{n}({\bf{t}})$ $\displaystyle=$ $\displaystyle
E\Big{[}\exp\\{i\sum\limits_{j=1}^{n}{\bf{t}}^{{}^{\prime}}Q_{j}/n^{1/4}\\}\Big{]}=\Big{[}E(\exp\\{i{\bf{t}}^{{}^{\prime}}Q_{j}/n^{1/4}\\})\Big{]}^{n}$
$\displaystyle=$
$\displaystyle\Big{[}1-\frac{h(0){\bf{t}}^{{}^{\prime}}\Sigma{\bf{t}}}{2an}+o(n^{-1})\Big{]}^{n}\longrightarrow\exp\Big{[}-\frac{h(0){\bf{t}}^{{}^{\prime}}\Sigma{\bf{t}}}{2a}\Big{]}.$

Hence,

$\frac{1}{n^{1/4}}\sum\limits_{j=1}^{n}Q_{j}\xrightarrow{\ d\
}N(0,\frac{h(0)}{a}\Sigma).$

By Lemma A1, we have

$\frac{1}{\sqrt{n}}\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\xrightarrow{\
p\ }\frac{h(0)}{a}E[K({\bm{\theta}};{\bm{x}}^{A},{\bm{x}}^{B},D)].$

As a result,

$\displaystyle\sum\limits_{i=1}^{n}\frac{I{(Y_{i}^{A}+Y_{i}^{B}=1)}}{\sqrt{\sum\limits_{j=1}^{n}I{(Y_{j}^{A}+Y_{j}^{B}=1)}}}\big{[}\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\big{]}\xrightarrow{\
d\ }N(0,\Sigma/c)$

which implies the Lemma A3 holds.

Proof of Theorem 3: Let

$l({\bm{\theta}})=\sum\limits_{i=1}^{n}Y_{i}^{A}(1-Y_{i}^{B})\log
p({\bm{\theta}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})+(1-Y_{i}^{A})Y_{i}^{B}\log(1-p({\bm{\theta}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))$

$\hat{\bm{\theta}}$ is its maximum point which satisfies

$\frac{\partial
l(\hat{{\bm{\theta}}})}{\partial{{\bm{\theta}}}}=\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\left\\{\frac{Y_{i}^{A}-p(\hat{{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p(\hat{{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p(\hat{{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p(\hat{{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\right\\}=0.$

Expanding the above function around the true value, i.e., ${{\bm{\theta}}}$,
yields

$\displaystyle\dfrac{\partial
l(\hat{{\bm{\theta}}})}{\partial{{\bm{\theta}}}}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\left\\{\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\right\\}$
$\displaystyle+\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\frac{\partial\left[\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\right]}{\partial{{\bm{\theta}}}}_{|{{\bm{\theta}}}={{\bm{\theta}}}}(\hat{{\bm{\theta}}}-{{\bm{\theta}}})$
$\displaystyle+[\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)]o(\parallel\hat{{\bm{\theta}}}-{{\bm{\theta}}}\parallel).$

By Lemma A3, it can be proved that

$\frac{1}{\sum\limits_{j=1}^{n}I(Y_{j}^{A}+Y_{j}^{B}=1)}\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)\frac{\partial\left[\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\cdot\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\right]}{\partial{{\bm{\theta}}}}$

converges in probability to $-\Sigma/c$ and

$\sum\limits_{i=1}^{n}\frac{I(Y_{i}^{A}+Y_{i}^{B}=1)}{\sqrt{\sum\limits_{j=1}^{n}I(Y_{j}^{A}+Y_{j}^{B}=1)}}\Big{\\{}\frac{Y_{i}^{A}-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})(1-p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i}))}\frac{\partial
p({{\bm{\theta}}};{\bm{x}}^{A}_{i},{\bm{x}}^{B}_{i},D_{i})}{\partial{{\bm{\theta}}}}\Big{\\}}\xrightarrow{d}N(0,\Sigma/c)$

Thus, one can obtain

$\displaystyle\sqrt{\sum\limits_{i=1}^{n}I(Y_{i}^{A}+Y_{i}^{B}=1)}(\hat{{{\bm{\theta}}}}-{{\bm{\theta}}})\xrightarrow{\
d\ }N(0,c\Sigma^{-1})$

and Theorem 3 is proved.

## References

Abadie, A., Gardeazabal, J., 2003. The economic costs of conflict: a case
study of the basque country. The American Economic Review. 93(1), 113–132.

Abadie, A., Imbens, G.W., 2006. Large sample properties of matching estimators
for average treatment effects. Econometrica.  74, 235–267.

Abadie, A., Imbens, G.W., 2011. Bias-corrected matching estimators for average
treatment effects. Journal of Business and Economic Statistics.  29, 1–11.

Arellano, M., 2003. Discrete choices with panel data. Investigaciones
Economicas. 27(3), 423–458.

Arellano, M., Bonhomme, S., 2009. Robust priors in nonlinear panel data
models. The Econometric Society. 77(2), 489–536.

Arulampalam, W., Stewart, M.B., 2009. Simplified implementation of the Heckman
estimator of the dynamic probit model anda comparison with alternative
estimators. Oxford Bulletin of Economics and Statistics. 71, 659–681.

Ashenfelter, O., Card, D., 1985. Using the longitudinal structure of earnings
to estimate the effect of training programs. The Review of Economics and
Statistics.  67(4), 648–660.

Bartolucci, F., Bellio, R., Salvan, A., Sartori, N., 2016. Modified profile
likelihood for fixed-effects panel data models. Econometric Reviews. 35(7),
1271–1289.

Bartolucci, F., Pigini, C., 2017. cquad: An R and Stata package for
conditional maximum likelihood estimation of dynamic binary panel data models.
Journal of Statistical Software. 78, 1–26.

Bertschek, I., Lechner, M., 1998. Convenient estimators for the panel probit
model. Journal of Econometrics. 87(2), 329–371.

Butler, I., Moffitt, R., 1982. A computationally efficient quadrature
procedure for the one-factor multinomial probit model. Econometrica. 50(3),
761–764.

Card, D., Krueger, A.B., 1994. Minimum wages and employment: a case study of
the fast-food industry in New Jersey and Pennsylvania. American Economic
Association. 84(4), 772–793.

Carro, J.M., 2007. Estimating dynamic panel data discrete choice models with
fixed effects. Journal of Econometrics. 140(2), 503–528.

Chamberlain, G., 1980. Analysis of covariance with qualitative data. The
Review of Economic Studies. 47(1), 225–238.

Dey, K.D., Chen, M,H., Chang, H., 1997. Bayesian approach for nonlinear random
effects models. Biometrics. 53, 1239–1252.

Díaz, J., Rau, T., Rivera, J., 2015. A matching estimator based on a bi-level
optimization problem. Review of Economics and Statistics. 97(4), 803–812.

Freireich, E.J., Gehan, E.A., et al., 1963. The effect of 6-mercaptopurine on
the duration of steroid-induced remissions in acute leukemia. Blood. 21(6),
699–716.

Gao, W., Bergsma, W., Yao, Q., 2017. Estimation for dynamic and static panel
probit models with large individual effects. Jouranl of Time Series Analysis.
38, 266-284.

Gehan, E., 1965. A generalized Wilcoxon test for comparing arbitrarily single-
censored samples. Biometrika. 52, 203–233.

Guilkey, D.K., Murphy, J.L., 1993. Estimation and testing in the random
effects probit model. Journal of Econometrics. 59(3), 301–317.

Heckman, J., 1978. Dummy endogenous variables in a simultaneous equation
system. Econometrica. 46(4), 931–958.

Heckman, J., Smith, J., Clements, N., 1997. Making the most out of programme
evaluations and social experiments: accounting for heterogeneity in programme
impacts. The Review of Economic Studies. 64(4), 487-535.

Hernan, M.A., Robins, J.M., 2006. Estimating Causal Effects From
Epidemiological Data. Journal of Epidemiology and Community Health. 60, 578
C-596.

Hirano, K., Imbens, G.M., Ridder, G., 2003. Efficient estimation of average
treatment effects using the estimated propensity score. Econometrica. 71,
1161–1189.

Horowitz, J.L., 1992. A smoothed maximum score estimator for binary response
model. Econometrica. 60, 505–531.

Imbens, G.M., Angrist, D.J., 1994. Identification and estimation of local
average treatment effects. Econometrica. 62(2), 467–475.

Morton, D.E., Saah, A.J., et al., 1982. Lead absorption in children of
employees in a lead related industry. American Journal of Epidemiology.
115(4), 549–555.

Neyman, J, 1923 [1990]. On the Application of Probability Theory to
Agricultural Experiments. Essay on Principles. Section 9. Statistical Science.
5(4), 465–472. Trans. Dorota M. Dabrowska and Terence P. Speed.

Pruzek, R.M., Helmreich, J.E., 2009. Enhancing dependent sample analysis with
graphics. Journal of Statistics Education. 17 (1).

Robins, J.M., Rotnitzky, A., Zhao, L.P., 1994. Estimation of regression
coefficients when some regressors are not always observed. Journal of the
American Statistical Association. 89(427), 846-866.

Robins, J.M., Rotnitzky, A., 1995. Semiparametric efficiency in multivariate
regression models with missing data. Journal of the American Statistical
Association. 90, 122–129.

Rubin, D.B., 1974. Estimating causal effects of treatments in randomized and
nonrandomized studies. Journal of Educational Psychology. 66(1), 688-701.

Sakaguchi, S., 2016. Estimation of time-varying average treatment effects
using panel data when unobserved fixed effects affect potential outcomes
differently. Economics Letters. 146, 82–84.

Walker, S., 1996. An EM algorithm for nonlinear random effects models.
Biometrics. 52, 934–944.