The Asymmetric Power-Student-t Model for Censored and Truncated Data

Abstract

In this paper, we propose the power Student-t regression model for censored (limited) observations which extends the Student-t censored regression model. This extension is based on the asymmetric and heavy-tailed power Student-t distribution. The score functions and expected information matrix are given as well as the process for estimating the parameters in the model is discussed by using the likelihood approach. Two simulation studies are conducted to evaluate parameter recovery and properties of the model and finally, two applications to a real data set are reported to demonstrate the usefulness of this new methodology.

Key words
Censored regression model; Fisher information matrix; maximum likelihood estimation; power Student-$t$ distribution

INTRODUCTION

Regression models where the response variable is censored or limited are common in different fields: clinical essays, econometric analysis, social phenomena, engineering studies, among others. In clinical essays for example, in the first phases of development of the new vaccines, the determination of antibody concentration values often are left-censored due to detection limit by lack of sensitivity of the essay when the concentrations are near zero, see Moulton & Halsey 1995. In social phenomena, the study on extramarital behavior where the variable of interest is the number of extramarital affairs in the previous year, for example, it can result in a left-censored variable (Fair 1978). In econometrics analysis, the ordinary Tobit model (Tobin 1958) is commonly used to conduct studies of the labor force participation of married women. In this case, the observed response is the wage rate, which is typically considered as censored below zero, i.e., for working women, positive values for the wage rates are registered, whereas for the non-working women the observed wage rates are zero; see Mroz 1987.

In situations such as previously discussed, where censored regression (CR) models are proposed, it is common to assume a normal distribution for the error term, however, this assumption can not be suitable and it can be unrealistic due to the presence of atypical observations or high (or low) degree of skewness and kurtosis of the response variable, which the normal model is unable to capture, so considerable interest has centered on relaxing the assumption of normality of the errors in CR models. In this context, some authors have proposed a wide range of alternatives to the normal censored regression (NCR) model which is widely known in the literature as the Tobit model. Arellano-Valle et al. 2012 for example, extend the classical Tobit model by introducing the Student-$t$ censored regression (TCR) model that can be suitable when the response variable has heavy-tails and the kurtosis is greater than the usual normal model. Another extension of the Tobit model was proposed by Martínez-Flórez et al. 2013 by considering that random errors follow a power-normal (PN) distribution (Gupta & Gupta 2008). The novelty of this proposal is the incorporation of a shape parameter which gives flexibility to the assumption of the symmetric errors (normality assumption) and it allows to accommodating skewed forms to the left and the right for the error term in CR models. Recently Garay et al. 2017 proposed a family of censored regression models based on the family of symmetric distributions commonly known as the scale mixture of normal (SMN) distributions, which includes the TCR model proposed by Arellano-Valle et al. 2012. This family also includes Pearson type VII (PVII), slash (SL), power exponential (PE), contaminated normal (CN) and normal (N) distributions. In addition to being robust, these models have shown to be useful in detecting atypical observations in CR models.

Although some proposals that take into account the problem of atypical observations in censored regression models, most of them are based on the assumption of symmetry of the error and few studies that capture departure from symmetry in the distribution of errors as in Martínez-Flórez et al. 2013, for example, who support their work in the great virtues of the alpha-power models to fit data where distribution presents high or low asymmetry and/or kurtosis.

Within this class of alpha-power models Zhao & Kim 2016 proposed an extension of the Student-$t$ model by defining the power-Student-$t$ (PT) distribution as an alternative to the skew-$t$ model by Azzalini & Capitanio 2003 for fitting skewed and heavy-tailed data. The PT model, which extends the power-normal model by Gupta & Gupta 2008 seems to be useful in situations where the data present higher degree of skewness and kurtosis than PN model in presence of atypical observations.

In this paper, we propose a censored regression model under the assumption that errors follow a PT distribution (hereafter we will call it the PTCR model). The assumption of PT distribution gives flexibility for accommodating skew forms to the left and the right, and kurtosis greater or smaller than the Student $t$-distribution can be also accommodated, hence, PTCR model extends the TCR model. The process of inference in the model is conducted by using the maximum likelihood (ML) approach and its large sample properties. Application is implemented to real data set where it is demonstrated that the proposed model can be very useful in fitting real data sets.

The rest of this paper is organized as follows: Section “The Power Student-t Distribution“ presents a brief review of the main properties of the PT distribution. In Section “Power-Student-t Model for Censored and Truncated Data“, we introduce the censored and truncated PT models. Section “Censored Power-Student-t Regression Model“ introduces the PTCR model. Here, ML equation and the observed and expected information matrices are given. Section “Simulation Study“ presents the results of a simulation study which reveals the good performance of the estimation approach. The PTCR model is fitted to a data set of housewives wages in Section “Real Data Application“, revealing that the data set in question can be fitted by PTCR as well as by a CR model where the observational errors have a SMN distribution (SMNCR model).

THE POWER STUDENT-T DISTRIBUTION

In this section, we present the PT distribution and review some of its main characteristics and properties. The PT model was introduced by Zhao & Kim 2016 and it is an alternative to skew-t model for fitting data with high indices of asymmetry and kurtosis in addition to heavy tails.

Definition 1. The random variable $X$ is said to have a PT distribution with parameter $\alpha$ and degree of freedom $\nu$, if $X$ has probability density function (PDF) given by

$$f_{PT}(x;\alpha ,\nu )=\alpha f_T(x;\nu )[F_T(x;\nu ){]}^{\alpha -1},$$ (1)

for $x\in ℝ$ and $\alpha >0$. Functions $f_T(\cdot ;\nu )$ and $F_T(\cdot ;\nu )$ are the PDF and cumulative distribution function (CDF) of the standard Student-$t$ distribution.

Random variable having $f_{PT}(x;\alpha ,\nu )$ distribution is denoted shortly by $$X\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$. Figure 1 displays some forms of the PDF of the PN distribution for selected values of $\alpha$. Note from figure that parameter $\alpha$ controls the skewness and kurtosis of the distribution. The CDF of the PT model is given by

$$F_{PT}(x;\alpha ,\nu )=[F_T(x;\nu ){]}^{\alpha },$$ (2)

for $x\in ℝ$. Some properties of the PT distribution can be proven as result of Definition 1.

Figure 1
Density function of $$\mathrm{P}\mathrm{T}(\alpha ,10)$$ for $\alpha =5$ (solid line), $\alpha =2$ (dashed line), $\alpha =1$ (dotted line), $\alpha =0.5$ (dotted-dashed line).

Proposition 1. Let $$X\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$, then

• (i) if $\alpha =1$, $X$ follows Student-$t$ distribution and we write $$X\sim \mathrm{T}(\nu )$$,

• (ii) if $\nu \to \infty$, $X$ converges to power-normal (PN) model with parameter $\alpha$. The PDF is given by

$$f_{PN}(x;\alpha )=\alpha \varphi \left(x\right)[\Phi \left(x\right){]}^{\alpha -1},x\in ℝ.$$ (3)

• More details of PN distribution can be found in Gupta & Gupta 2008 and Pewsey et al. 2012.

• (iii) if $\alpha =1$ and $\nu \to \infty$, $X$ converges to standard normal distribution.

Proof. Proof of (i)-(iii) are directly obtained from definition of PT distribution ◻

Proposition 2. Let $$X\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$, then for $k\in \mathbb{N}$

$$E[X^k]=E[(F_T^{-1}(Y;\nu ){)}^k],$$ (4)

where $Y$ has a beta distribution, and $F_T^{-1}(\cdot ;\nu )$ denotes the inverse of the function $F_T(\cdot ;\nu )$ .

Proof. We have by definition that

$$E[X^k]={\int }_{ℝ}{x}^k\alpha f_T(x;\nu )[F_T(x;\nu ){]}^{\alpha -1}dx$$

thus, letting $y=F_T(x;\nu )$, then $x=F_T^{-1}(y;\nu )$, it follows that

$$E[X^k]={\int }_0^1(F_T^{-1}(y;\nu ){)}^k\alpha y^{\alpha -1}dy$$

which is the expected value of the function $(F_T^{-1}(Y;\nu ){)}^k$, where $Y$ follows a beta distribution with parameter $\alpha$ and 1. ◻

The expected value, variance, indices of asymmetry and kurtosis of the PT model can be found by using the expressions

• (i) $E\left[X\right]={\mu }_1$

• (ii) $V\left[X\right]={\mu }_2-{\mu }_1^2$

• (iii) ${\gamma }_1={\displaystyle \frac{{\mu }_3-3{\mu }_1{\mu }_2+2{\mu }_1^3}{{({\mu }_2-{\mu }_1^2)}^{3/2}}}$

• (iv) ${\gamma }_2={\displaystyle \frac{{\mu }_4\prime }{{\left({\mu }_2\prime \right)}^2}}$

where ${\mu }_r=E\left[X^r\right]$ and ${\mu }_r\prime =E{[X-E\left(X\right)]}^r$. Table I presents the values of the asymmetry coefficient of the PT model for some values of the $\nu$ parameter and for values of $\alpha$ in the range of 0.1 to 100000.

Table I
Skewness of the $$\mathrm{P}\mathrm{T}(\alpha ,\nu )$$ model. Values of $\alpha$ ranging of 0.1 to 100000.

$\nu$	$min$	$max$
4	-9.3423	5.5913
5	-7.5147	3.5101
6	-7.3825	2.7668
7	-4.6073	2.3721
8	-4.1468	2.1228
9	-3.8410	1.9490
10	-3.4826	1.8202
12	-2.2654	1.6409
15	-2.0131	1.4752
20	-1.8802	1.3208
25	-1.8245	1.2329
30	-1.6372	1.1760
50	-1.1158	1.0664
100	-1.0401	0.9864
500	-0.9390	0.9264

Definition 2. Let $$X\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$. The PT density of location and scale is defined as the distribution of $Z=\xi +\eta X$, for $\xi \in ℝ$ and $\eta >0$. The corresponding PDF is given by

$$f_{PT}(z;\xi ,\eta ,\alpha ,\nu )=\alpha f_T(\frac{z-\xi }{\eta };\nu )[F_T(\frac{z-\xi }{\eta };\nu ){]}^{\alpha -1},$$ (5)

for $z\in ℝ$. We denote this extension as $$Z\sim \mathrm{P}\mathrm{T}(\xi ,\eta ,\alpha ,\nu )$$, and we have that $$\mathrm{P}\mathrm{T}(0,1,\alpha ,\nu )\equiv \mathrm{P}\mathrm{T}(\alpha ,\nu )$$.

The $k$th moment of the random variable $Z$ is given by

$$E[Z^k]=\sum _{j=0}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right){\xi }^j{\eta }^{k-j}{\mu }_{k-j},$$ (6)

where ${\mu }_j$ is the $j$th moment of a random variable $$X\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$. Zhao & Kim 2016 derived the information matrix for the location-scale version and showed that it is non-singular when $\alpha =1$ for small values of the parameter $\nu$ (i.e., $\nu \ll 30$). When $\nu$ tends to $+\infty$, then PT distribution converges to PN model and here, we recall Pewsey et al. 2012 showed that PN model has non-singular information matrix. This result guarantees that regularity conditions are satisfied for the likelihood approach, hence, with PT model, symmetry can be tested by using ordinary large sample properties of the likelihood ratio statistics.

POWER-STUDENT-T MODEL FOR CENSORED AND TRUNCATED DATA

Based on the goodness of the PT distribution to fit data with high indices of asymmetry and kurtosis, in this section, we introduce the censored PT and the truncated PT models which we will be denoted by CPT and TPT, respectively.

Definition 3 (Censored PT Model). Suppose that random variable $Y$ follows a $$\mathrm{P}\mathrm{T}(\xi ,\eta ,\alpha ,\nu )$$ distribution. Let $Y_1,\dots ,Y_n$ a random sample of size $n$ of $Y$, where only those values of $Y_i$ greater than constant $k_i$ are recorded; and for values $Y_i\le k_i$ only the value $k_i$ is registered. The observed values $Y_i^o$ can be written as

$$Y_i^o=\{\begin{array}{cc}{k}_i,\hfill & \mathrm{\text{if}}{Y}_i\le k_i\mathrm{\text{,}}\hfill \\ Y_i,\hfill & \mathrm{\text{if}}{Y}_i>k_i\mathrm{\text{,}}\hfill \end{array}$$

for $i=1,\dots ,n$. The resulting sample is said to be a censored power-Student-$t$ (CPT).

From Definition 3 it follows that $P(Y_i^o=k_i)=P(Y_i\le k_i)={\left\{F_T(\frac{k_i-\xi }{\eta };\nu )\right\}}^{\alpha }$, and for the observations $Y_i^o=Y_i$ the distribution of $Y_i^o$ is the same of $Y_i$, i.e., $$Y_i^o\sim \mathrm{P}\mathrm{T}(\xi ,\eta ,\alpha ,\nu )$$. For convenience, we choose to work with the case of left-censored data, however, the followings results can be extended to other types of censorship.

Maximum Likelihood Estimation for CPT Model

Let $Y_1^o,\dots ,Y_n^o$ be a random sample of the censored $$\mathrm{P}\mathrm{T}(\xi ,\eta ,\alpha ,\nu )$$ distribution (censored in $k_i$). To perform statistical inference for parameter vector $𝛉={(\xi ,\eta ,\alpha ,\nu )}^{\top }$ by using the ML method, we use the reparameterization in Olsen 1978. Thus, let $\gamma =\sigma \xi$ and $\sigma =1/\eta$, the log-likelihood function for the new vector $𝛗={(\gamma ,\sigma ,\alpha ,\nu )}^{\top }$, given the observed sample $${\mathbf{Y}}^o$$ can be written as

$$\begin{array}{rl}\ell (\mathit{\phi };{\mathbf{Y}}^o)& amp;=\alpha \sum _{i=1}^n(1-d_i)\log \{F_T(c_i;\nu )\}\\ & amp;+\sum _{i=1}^n{d}_i\{\log \alpha +\log \sigma +\log f_T(z_i;\nu )+(\alpha -1)\log F_T(z_i;\nu )\}\end{array}$$ (7)

where $c_i=\sigma k_i-\gamma$, $z_i=\sigma y_i-\gamma$, and $d_i$ is an indicator variable defined as $d_i=1$ if $Y_i^o=Y_i$; and $d_i=0$ if $Y_i^o=k_i$. The components of the score function $𝐔\left(𝛗\right)$ are obtained by deriving partially the log-likelihood function given in (7) with regard to components $\gamma$, $\sigma$, $\alpha$ and $\nu$, we obtain the following equations

$$\begin{array}{rl}U(\gamma )& =-\alpha \sum _{i=1}^n(1-d_i)r(c_i;\nu )\\ & +\sum _{i=1}^n{d}_i\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i-(\alpha -1)r(z_i;\nu )\}\\ U(\sigma )& =\alpha \sum _{i=1}^n(1-d_i)r(c_i;\nu )k_i\\ & +\sum _{i=1}^n{d}_i\{\frac{1}{\sigma }-\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i{y}_i+(\alpha -1)r(z_i;\nu )y_i\}\\ U(\alpha )& =\sum _{i=1}^n(1-d_i)\log \{F_T(c_i;\nu )\}+\sum _{i=1}^n{d}_i\{\frac{1}{\alpha }+\log F_T(z_i;\nu )\}\\ U(\nu )& =\frac{\alpha }{2}\sum _{i=1}^n(1-d_i)\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )-\frac{c_{i}r(c_i;\nu )}{\nu }\}\\ & +\frac{1}{2}\sum _{i=1}^n{d}_i\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-\frac{1}{\nu }\\ & -\log (1+\frac{z_i^2}{\nu })+\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }\}\\ & +\frac{\alpha -1}{2}\sum _{i=1}^n{d}_i\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(z_i;\nu )-\frac{z_{i}r(z_i;\nu )}{\nu }\}\end{array}$$

where $r(x;\nu )=f_T(x;\nu )/F_T(x;\nu )$, $\psi \left(x\right)=\frac{d}{dx}log\Gamma \left(x\right)$ is the digamma function and $b_{mn}(c_m;\nu +m)$ is the truncated moment defined as

$$b_{mn}(c_m;\nu +m)={\int }_{-\infty }^{c_m}{s}^m{\{log(1+\frac{s^2}{\nu +m})\}}^n\frac{f_T(s;\nu +m)}{F_T(c_m;\nu +m)}ds,$$ (8)

with $c_m=\sqrt{\frac{\nu +m}{\nu }}c$ and $c_0=c$. The moments $b_{mn}$ in (8) are obtained by numerical integration, for example, by using integrate function of R Development Core Team 2018. The ML estimates of the parameters $\xi$, $\eta$, $\alpha$ and $\nu$ in censored PT model are obtained using iterative algorithms based on the score functions and by applying the inverse transformation $\xi =𝛄/\sigma$ and $\eta =1/\sigma$. For obtaining the standard errors of the ML estimates one should compute the information matrix $${\mathbf{I}}_{\mathit{\phi }}$$. It is well known that the elements of $${\mathbf{I}}_{\mathit{\phi }}$$ are given by

$${\mathbf{I}}_{\mathit{\phi }}(i,k)=-E\left[\frac{{\mathrm{\partial }}^2\ell (\mathit{\phi };{\mathbf{Y}}^o)}{\mathrm{\partial }{\mathit{\phi }}_i\mathrm{\partial }{\mathit{\phi }}_k}\right],i,k=1,\dots ,4$$

where $𝛗={(\gamma ,\sigma ,\alpha ,\nu )}^{\top }$. Since expectation over PT distribution and second-order derivatives are not straightforward, numerical methods should be performed to obtain the explicit form of the information matrix $${\mathbf{I}}_{\mathit{\phi }}$$. Thus, we use the observed information matrix for calculating the standard errors in the rest of the paper. To recover the information matrix $${\mathbf{I}}_{\mathit{\theta }}$$ of the original parameterization $𝛉={(\xi ,\eta ,\alpha ,\nu )}^{\top }$, we use

$${\mathbf{I}}_{\mathit{\theta }}=(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta }{)}^{\mathrm{\top }}{\mathbf{I}}_{\mathit{\phi }}(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta }),$$

where the Jacobian matrix is

$$\frac{\partial 𝛗}{\partial 𝛉}=\left(\begin{array}{cccc}\hfill \frac{1}{\eta }\hfill & \hfill -\frac{1}{{\eta }^2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -\frac{1}{{\eta }^2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)$$ (9)

Definition 4 (Truncated PT Model). Let $X$ be a random variable with distribution $$X\sim \mathrm{P}\mathrm{T}(\xi ,\eta ,\alpha ,\nu )$$. Let $a,b\in ℝ$ with $a<b$, such that $P(a<X<b)>0$. It is said that random variable $Y$ has a truncated power-Student-$t$ (TPT) distribution in the interval $(a,b)$, if $Y$ has the same distribution as $X\mid X\in (a,b)$. In this case, we write $$Y\sim {\mathrm{T}\mathrm{P}\mathrm{T}}_{(a,b)}(\xi ,\eta ,\alpha ,\nu )$$.

As a consequence of the Definition 4, the PDF of TPT distribution can be obtained as

$$\begin{array}{cc}\hfill f_{PT}(y\mid (a,b))& =\frac{\alpha }{\eta }{f}_T(\frac{y-\xi }{\eta };\nu ){\left[F_T(\frac{y-\xi }{\eta };\nu )\right]}^{\alpha -1}\hfill \\ \hfill & \times {\{{\left[F_T\left(\frac{b-\xi }{\eta }\right)\right]}^{\alpha }-{\left[F_T\left(\frac{a-\xi }{\eta }\right)\right]}^{\alpha }\}}^{-1}\hfill \end{array}$$

if $a<y<b$, and $f_{PT}(y\mid (a,b))=0$ in otherwise. Now, we consider that before the sample to be selected, the distribution of $Y$ is truncated at the value $k$, so that we can only choose observations such that $Y\le k$. Then, random variable $Y$ has PDF given by

$$\begin{array}{cc}\hfill f_{PT}(y\mid y\le k)& =\frac{\alpha }{\eta }{f}_T(\frac{y-\xi }{\eta };\nu ){\left[F_T(\frac{y-\xi }{\eta };\nu )\right]}^{\alpha -1}{\left[F_T(\frac{k-\xi }{\eta };\nu )\right]}^{-\alpha }\hfill \end{array}$$ (10)

for $-\infty <y\le k$.

Maximum Likelihood Estimation for TPT Model

Given a sample $$\mathbf{Y}=(Y_1,\dots ,Y_n{)}^{\mathrm{\top }}$$ of the TPT distribution in the value $k$, the log-likelihood function for vector $𝛗={(\gamma ,\sigma ,\alpha ,\nu )}^{\top }$, where $\gamma =\sigma \xi$ and $\sigma =1/\eta$, is given by

$$\begin{array}{rl}\ell (\mathit{\phi };\mathbf{Y})& =-n_1\alpha \log \{F_T(c;\nu )\}\\ & +\sum _{y_i\le k}^n\{\log \alpha +\log \sigma +\log f_T(z_i;\nu )+(\alpha -1)\log F_T(z_i;\nu )\}\end{array}$$ (11)

where $c=\sigma k-\gamma$, $z_i=\sigma y_i-\gamma$, and $n_1$ is the number of observations in the sample such that $-\infty <y\le k$. Deriving partially the log-likelihood function (11) with respect to the components of the vector $𝛗$ the following elements of the score function are obtained

$$\begin{array}{rl}U(\gamma )& =n_1\alpha r(c;\nu )+\sum _{y_i\le k}\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i-(\alpha -1)r(z_i;\nu )\}\\ U(\sigma )& =-n_1\alpha r(c;\nu )k+\sum _{y_i\le k}\{\frac{1}{\sigma }-\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i+(\alpha -1)r(z_i;\nu )\}\\ U(\alpha )& =-n_1\log \{F_T(c;\nu )\}+\sum _{y_i\le k}\{\frac{1}{\alpha }+\log F_T(z_i;\nu )\}\\ U(\nu )& =-\frac{n_1\alpha }{2}\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c;\nu )-\frac{cr(c;\nu )}{\nu }\}\\ & +\frac{1}{2}\sum _{y_i\le k}\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-\frac{1}{\nu }\\ & -\log (1+\frac{z_i^2}{\nu })+\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }\}\\ & +\frac{\alpha -1}{2}\sum _{y_i\le k}\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(z_i;\nu )-\frac{z_{i}r(z_i;\nu )}{\nu }\}\end{array}$$

where $r(x;\nu )=f_T(x;\nu )/F_T(x;\nu )$, $\psi \left(x\right)=\frac{d}{dx}log\Gamma \left(x\right)$ is the digamma function and $b_{mn}(c_m;\nu +m)$ is given by (8). To obtain the ML estimates of the parameters $\xi$, $\eta$, $\alpha$ and $\nu$ in the TPT model, we proceed in a similar way to the CPT model, and iterative methods based on the Newton-Rapshon algorithm are used with the score functions. We use the observed information matrix for calculating the standard errors and to recover the information matrix $${\mathbf{I}}_{\mathit{\theta }}$$ of the original parameterization $𝛉={(\xi ,\eta ,\alpha ,\nu )}^{\top }$, we use $${\mathbf{I}}_{\mathit{\theta }}=(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta }{)}^{\mathrm{\top }}{\mathbf{I}}_{\mathit{\phi }}(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta })$$, where $\partial 𝛗/\partial 𝛉$ is given in (9).

CENSORED POWER-STUDENT-T REGRESSION MODEL

In this section, we introduce the censored power-Student-$t$ regression model, which is denoted by PTCR. This model results from the consideration of the observed random variable $Y_i^o=D_i{Y}_i$, where $D_i=I_{(0,+\infty )}\left(Y_i\right)$ and $$Y_i={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }+{\epsilon }_i$$, with $${\epsilon }_i\stackrel{\mathrm{i}\mathrm{i}\mathrm{d}}{\sim }\mathrm{P}\mathrm{T}(0,\eta ,\alpha ,\nu )$$, for $i=1,\dots ,n$; i.e.,

$$Y_i^o=\{\begin{array}{ll}{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }+{\epsilon }_i,& \text{if }{Y}_i>0\text{, }\\ 0,& \text{if }{Y}_i\le 0\text{,}\end{array}$$ (12)

where $𝛃$ is a vector of dimension $p$ of unknown parameters, $${\mathbf{x}}_i=(x_{i1},\dots ,x_{ip}{)}^{\mathrm{\top }}$$, for $i=1,\dots ,n$, are vectors of known covariates, and ${\epsilon }_i$, for $i=1,\dots ,n$, are independent random variables with PT distribution with location parameter 0, scale parameter $\eta$, shape parameter $\alpha$, and degrees of freedom $\nu$. This assumption is equivalent to considering that unobserved random variables $Y_1,\dots ,Y_n$ are independent with $$Y_i\sim \mathrm{P}\mathrm{T}({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta },\eta ,\alpha ,\nu )$$, that is, with PDF given by

$$g(y_i;\mathit{\theta })=\frac{\alpha }{\eta }{f}_T(\frac{y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu ){\left[F_T(\frac{y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu )\right]}^{\alpha -1}$$

for $i=1,\dots ,n$, where $𝛉={({𝛃}^{\top },\eta ,\alpha ,\nu )}^{\top }$. The contribution to likelihood function for observations $Y_i^o=0$ is given by

$$P(Y_i^o=0)={\left\{F_T(-\frac{{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu )\right\}}^{\alpha },$$

and for observations $Y_i^o\ge 0$, we have that $${}_i^o\sim \mathrm{P}\mathrm{T}({\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta },\eta ,\alpha ,\nu )$$. Therefore, the likelihood function of the PTCR model based on the observed sample $${\mathbf{Y}}^o=(Y_1^o,\dots ,Y_n^o{)}^{\mathrm{\top }}$$ is given by

$$\begin{array}{rl}L(\mathit{\theta };{\mathbf{Y}}^o)=& \prod _{i=1}^n{\left\{F_T(-\frac{{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu )\right\}}^{\alpha (1-d_i)}\{\frac{\alpha }{\eta }{f}_T(\frac{y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu )\\ & \times {\left[F_T(\frac{y_i-{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }}{\eta };\nu )\right]}^{\alpha -1}{\}}^{d_i}\end{array}$$ (13)

where

$$d_i=\{\begin{array}{cc}1,\hfill & \mathrm{\text{if}}{Y}_i^o>0\mathrm{\text{,}}\hfill \\ 0,\hfill & \mathrm{\text{if}}{Y}_i^o=0\mathrm{\text{,}}\hfill \end{array}$$

Model (13) can be extended to the situation where the value of the censorship associated with the observation $i$ is replaced by the value $k_i$ (a known value), i.e.,

$$Y_i^o=\{\begin{array}{ll}{\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }+{\epsilon }_i,& \text{if }{Y}_i>k_i\text{, }\\ k_i,& \text{if }{Y}_i\le k_i\text{,}\end{array}$$ (14)

for $i=1,\dots ,n$. Note that, by making $Y_i^{o*}=Y_i^o-k_i$, $${\mathbf{x}}_i^{\ast }=({\mathbf{x}}_i^{\mathrm{\top }},k_i{)}^{\mathrm{\top }}$$ and ${𝛃}^{*}={({𝛃}^{\top },-1)}^{\top }$, we have the previous model in (12), hence, the results of the inference based on the ML method can be used to fit the more general model in (14).

Proposition 3. Consider the model (12) with assumption $${\epsilon }_i\stackrel{\mathrm{i}\mathrm{i}\mathrm{d}}{\sim }\mathrm{P}\mathrm{T}(0,\eta ,\alpha ,\nu )$$, for $i=1,\dots ,n$, then

1. (i) if $\alpha =1$, model (12) is reduced to $t$-Student censored regression (TCR) model

2. (ii) if $\nu \to \infty$, model (12) converges to power-normal censored regression (PNCR) model see Martínez-Flórez et al. 2013

3. (iii) if $\alpha =1$ and $\nu \to \infty$, model (12) converges to usual normal censored regression (NCR) model, that is, Tobit model.

Proof. Proof of (i)-(iii) are directly obtained from definition of PTCR model. ◻

Moments

Proposition 4. The mean and variance for the $i$th observed response in PTCR model are given by

$$\begin{array}{rl}E[Y_i^o]& =\eta (c_i+m_1(c_i;\nu ))(1-[F_T(-c_i;\nu ){]}^{\alpha })\end{array}$$ (15)

and

$$\begin{array}{rl}Var[Y_i^o]& ={\eta }^2(c_i^2+2c_i{m}_1(c_i;\nu ))([F_T(-c_i;\nu ){]}^{\alpha })(1-[F_T(-c_i;\nu ){]}^{\alpha })\\ & +{\eta }^2(m_2(c_i;\nu )-{[m_1(c_i;\nu )]}^2(1-[F_T(-c_i;\nu ){]}^{\alpha }))(1-[F_T(-c_i;\nu ){]}^{\alpha })\end{array}$$ (16)

respectively, where $c_i={\mu }_i/\eta$ with $${\mu }_i={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }$$ and

$$m_r(c;\nu )=\frac{1}{1-\{F_T(-c;\nu ){\}}^{\alpha }}{\int }_{F_T(-c;\nu )}^1{\left[F_T^{-1}(u;\nu )\right]}^r\alpha u^{\alpha -1}du,r=1,2.$$

Note that $m_r(c;\nu )$ is the moment $E\left[{\left(F_T^{-1}(U;\nu )\right)}^r\right]$ of a random variable $U$ with distribution $$\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(\alpha ,1)$$, truncated in the interval $(F_T(-c;\nu ),1)$.

Proof. The mean and variance for the $i$th observed response in PTCR model can be obtained by noting that $Y_i^o=D_i{Y}_i$, where $D_i=I_{(0,+\infty )}\left(Y_i\right)$ and $Y_i={\mu }_i+\eta Z_i$, with $${\mu }_i={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\beta }$$ and $$Z_i\stackrel{\text{iid}}{\sim }\mathrm{P}\mathrm{T}(\alpha ,\nu )$$, $i=1,\dots ,n$. We have

$$\begin{array}{rl}E[Y_i^o]& =E[Y_i\mid Y_i>0]P(Y_i>0)\\ & =E[{\mu }_i+\eta Z_i\mid {\mu }_i+\eta Z_i>0]P({\mu }_i+\eta Z_i>0)\\ & =({\mu }_i+\eta E[Z_i\mid Z_i>-\frac{{\mu }_i}{\eta }])P(Z_i>-\frac{{\mu }_i}{\eta })\\ & =\eta (\frac{{\mu }_i}{\eta }+E[Z_i\mid Z_i+\frac{{\mu }_i}{\eta }>0])(1-P(Z_i\le -\frac{{\mu }_i}{\eta }))\\ & =\eta (c_i+E[Z_i\mid Z_i+c_i>0])(1-[F_T(-c_i;\nu ){]}^{\alpha })\end{array}$$

By Lemma 1 in Appendix A, we have

$$E[Y_i^o]=\eta (c_i+m_1(c_i;\nu ))(1-[F_T(-c_i;\nu ){]}^{\alpha })$$

To obtain the variance of $Y_i^o$, note that

$$\begin{array}{rl}E[(Y_i^o{)}^2]& =E[Y_i^2\mid Y_i>0]P(Y_i>0)\\ & =E[({\mu }_i+\eta Z_i{)}^2\mid {\mu }_i+\eta Z_i>0]P({\mu }_i+\eta Z_i>0)\\ & =({\mu }_i^2+2\eta {\mu }_{i}E[Z_i\mid Z_i>-\frac{{\mu }_i}{\eta }]+{\eta }^{2}E[Z_i^2\mid Z_i>-\frac{{\mu }_i}{\eta }])P(Z_i>-\frac{{\mu }_i}{\eta })\\ & ={\eta }^2(\frac{{\mu }_i^2}{{\eta }^2}+2\frac{{\mu }_i}{\eta }E[Z_i\mid Z_i>-\frac{{\mu }_i}{\eta }]+E[Z_i^2\mid Z_i>-\frac{{\mu }_i}{\eta }])(1-P(Z_i\le -\frac{{\mu }_i}{\eta }))\\ & ={\eta }^2(c_i^2+2c_{i}E[Z_i\mid Z_i+c_i>0]+E[Z_i^2\mid Z_i+c_i>0])(1-[F_T(-c_i;\nu ){]}^{\alpha })\end{array}$$

Using Lemma 1 in Appendix, it follows that

$$E[(Y_i^o{)}^2]={\eta }^2(c_i^2+2c_i{m}_1(c_i;\nu )+m_2(c_i;\nu ))(1-[F_T(-c_i;\nu ){]}^{\alpha })$$

thus, by calculating $Var\left[Y_i^o\right]=E[{\left(Y_i^o\right)}^2]-(E\left[Y_i^o\right]{)}^2$ and after some algebraic manipulations, we obtain

$$\begin{array}{rl}Var[Y_i^o]& ={\eta }^2(c_i^2+2c_i{m}_1(c_i;\nu ))([F_T(-c_i;\nu ){]}^{\alpha })(1-[F_T(-c_i;\nu ){]}^{\alpha })\\ & +{\eta }^2(m_2(c_i;\nu )-{[m_1(c_i;\nu )]}^2(1-[F_T(-c_i;\nu ){]}^{\alpha }))(1-[F_T(-c_i;\nu ){]}^{\alpha })\end{array}$$

◻

It is important to note that, if $\nu$ tends to infinity, then (15) and (16) converge to the mean and variance of the PNCR model (Martínez-Flórez et al. 2013), i.e., when $\nu \to \infty$

$$\begin{array}{rl}E[Y_i^o]& =\eta (c_i+m_1(c_i))(1-[\mathrm{\Phi }(-c_i){]}^{\alpha })\\ Var[Y_i^o]& ={\eta }^2(c_i^2+2c_i{m}_1(c_i))([\mathrm{\Phi }(-c_i){]}^{\alpha })(1-[\mathrm{\Phi }(-c_i){]}^{\alpha })\\ & +{\eta }^2(m_2(c_i)-{[m_1(c_i)]}^2(1-\{\mathrm{\Phi }(-c_i){\}}^{\alpha }))(1-[\mathrm{\Phi }(-c_i){]}^{\alpha })\end{array}$$

where

$$m_r\left(c\right)=\frac{1}{1-\{\Phi (-c){\}}^{\alpha }}{\int }_{\Phi (-c)}^1{\left[{\Phi }^{-1}\left(u\right)\right]}^r\alpha u^{\alpha -1}du.$$

with $\Phi (\cdot )$ the CDF of the standard normal distribution, and ${\Phi }^{-1}(\cdot )$ the inverse function of $\Phi (\cdot )$. Also, worth noting that, if $\alpha =1$ and $\nu \to \infty$, we have $m_1(c_i;\nu )\to \varphi \left(c_i\right)/\Phi \left(c_i\right)$ and $m_2(c_i;\nu )\to (\Phi \left(c_i\right)-c_i\varphi \left(c_i\right))/\Phi \left(c_i\right)$, thus

$$\begin{array}{rl}E[Y_i^o]& =\eta (c_i\mathrm{\Phi }(c_i)+\varphi (c_i))\\ Var[Y_i^o]& ={\eta }^2(c_i^2\mathrm{\Phi }(c_i)+2c_i\varphi (c_i))(1-\mathrm{\Phi }(c_i))\\ & +{\eta }^2(\mathrm{\Phi }(c_i)-\varphi (c_i)[c_i+\varphi (c_i)])\end{array}$$

which are the mean and variance, respectively, of the Tobit model (Tobin 1958).

Maximum Likelihood Estimation

The ML method is considered by using the reparameterization of Olsen 1978. Let $𝛄=\sigma 𝛃$, $\sigma =1/\eta$ and

$$d_i=\{\begin{array}{ll}1,& \text{if }{Y}_i^o>k_i\text{,}\\ 0,& \text{if }{Y}_i^o=k_i\text{,}\end{array}$$

the log-likelihood function for $𝛗={({𝛄}^{\top },\sigma ,\alpha ,\nu )}^{\top }$ obtained from (13) under the new parameterization is given by

$$\begin{array}{rl}\ell (\mathit{\phi };{\mathbf{Y}}^o)=& \alpha \sum _{i=1}^n(1-d_i)\log \{1-F_T(c_i;\nu )\}\\ & +\sum _{i=1}^n{d}_i\{\log \alpha +\log \sigma +\log f_T(z_i;\nu )+(\alpha -1)\log F_T(z_i;\nu )\}\end{array}$$ (17)

where $z_i=\sigma y_i-c_i$, with $$c_i={\mathbf{x}}_i^{\mathrm{\top }}\mathit{\gamma }$$. The components of the score function $𝐔\left(𝛗\right)$ are obtained by deriving $$\ell (\mathit{\phi };{\mathbf{Y}}^o)$$ partially in relation to the components $𝛄$, $\eta$, $\alpha$ and $\nu$. After some algebraic manipulations the following components of the score function are obtained

$$\begin{array}{rl}U(\gamma )& =-\alpha \sum _{i=1}^n(1-d_i)r(-c_i;\nu ){\mathbf{x}}_i\\ & +\sum _{i=1}^n{d}_i\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i-(\alpha -1)r(z_i;\nu )\}{\mathbf{x}}_i\end{array}$$ (18)

$$\begin{array}{rl}U(\sigma )& =\sum _{i=1}^n{d}_i\{\frac{1}{\sigma }-\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i{y}_i+(\alpha -1)r(z_i;\nu )y_i\}\end{array}$$ (19)

$$\begin{array}{rl}U(\alpha )& =\sum _{i=1}^n(1-d_i)\log \{1-F_T(c_i;\nu )\}+\sum _{i=1}^n{d}_i\{\frac{1}{\alpha }+\log F_T(z_i;\nu )\}\end{array}$$ (20)

$$\begin{array}{rl}U(\nu )& =-\frac{\alpha }{2}\sum _{i=1}^n(1-d_i)\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )-\frac{c_{i}r(c_i;\nu )}{\nu }\}R(c_i;\nu )\\ & +\frac{1}{2}\sum _{i=1}^n{d}_i\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-\frac{1}{\nu }-\log (1+\frac{z_i^2}{\nu })\}\\ & +\frac{1}{2}\sum _{i=1}^n{d}_i\left\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }\right\}\\ & +\frac{\alpha -1}{2}\sum _{i=1}^n{d}_i\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(z_i;\nu )-\frac{z_{i}r(z_i;\nu )}{\nu }\}\end{array}$$ (21)

where $r(x;\nu )=f_T(x;\nu )/F_T(x;\nu )$, $R(x;\nu )=F_T(x;\nu )/(1-F_T(x;\nu ))$, $\psi \left(x\right)=\frac{d}{dx}log\Gamma \left(x\right)$ is the digamma function and $b_{mn}(c_m;\nu +m)$ is the truncated moment defined by (8). Note that, if $\alpha =1$ the equations (18)-(21) are reduced to the functions of the TCR model (Arellano-Valle et al. 2012), while, if $\alpha =1$ and $\nu \to \infty$, then $r(c;\nu )\to r\left(c\right)=\varphi \left(c\right)/\Phi \left(c\right)$, $R(c;\nu )\to R\left(c\right)=\Phi \left(c\right)/(1-\Phi \left(c\right))$, and $U\left(\nu \right)\to 0$, therefore, the equations (18) and (19) are reduced to score functions of the Tobit model.

The elements of the observed information matrix $${\mathbf{J}}_{\mathit{\phi }}$$ for PTCR model, which are denoted by $j_{{𝛗}_i{𝛗}_j}$, can be obtained by calculating the second partial derivative of the log-likelihood function (17), i.e., $$j_{{\mathit{\phi }}_i{\mathit{\phi }}_j}=-{\mathrm{\partial }}^2\ell (\mathit{\phi };{\mathbf{Y}}^o)/\mathrm{\partial }{\mathit{\phi }}_i\mathrm{\partial }{\mathit{\phi }}_j$$, while the expected information matrix is obtained as $${\mathbf{I}}_{\mathit{\phi }}=E[{\mathbf{J}}_{\mathit{\phi }}]$$, which involves the calculation of truncated expected values that have no closed form and must be obtained numerically. The Appendix B presents the expressions for the elements of the matrices $${\mathbf{I}}_{\mathit{\phi }}$$ and $${\mathbf{J}}_{\mathit{\phi }}$$. The expected information matrix $${\mathbf{I}}_{\mathit{\theta }}$$ of the original parameterization $𝛉={(𝛃,\eta ,\alpha ,\nu )}^{\top }$ can be recovered by using $${\mathbf{I}}_{\mathit{\theta }}=(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta }{)}^{\mathrm{\top }}{\mathbf{I}}_{\mathit{\phi }}(\mathrm{\partial }\mathit{\phi }/\mathrm{\partial }\mathit{\theta })$$, where

$$\frac{\mathrm{\partial }\mathit{\phi }}{\mathrm{\partial }\mathit{\theta }}=\left(\begin{array}{cccc}\frac{1}{\eta }{\mathbf{I}}_p& -\frac{1}{{\eta }^2}\mathit{\beta }& 0& 0\\ 0& -\frac{1}{{\eta }^2}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$$

Finally, the ML estimates for $𝛉={(𝛃,\eta ,\alpha ,\nu )}^{\top }$ can be obtained using iterative methods based on the Newton-Rapshon algorithm from the score function (18) - (21) and applying the inverse transformation $𝛃=𝛄/\sigma$ and $\eta =1/\sigma$. Estimates of the variances of the estimator can be obtained by evaluating the inverse of the observed information matrix $${\mathbf{J}}_{\mathit{\phi }}^{-1}$$ at the ML estimators $\widehat{𝛗}={(\widehat{𝛄},\hat{\sigma },\hat{\alpha },\hat{\nu })}^{\top }$ and by using the previous result.

Model Selection and Residual Analysis

In this section some criteria for the selection of the best-fitted model and a methodology for residual analysis are proposed.

Model Selection

Many model selection tools are generally used, such as the Akaike information criteria (AIC), (Akaike 1974), Bayesian information criterion (BIC) (Schwarz 1978), and the AIC corrected (AICc) (Sugiura 1978), which are defined by

$$-2\ell (\mathit{\theta }\mid \mathbf{y})+c(k,n)$$

where the term $c(k,n)$ is a quantity that depends on the number of free parameters that are estimated in the model $k$, and the number of observations in the sample $n$. For the AIC one has $c(k,n)=2k$, for BIC $c(k,n)=klog\left(n\right)$ and for AICc, $c(k,n)=2k(k+1)/(n-(k+1))$. To choose the best-fitted model, the criteria AIC, BIC and AICc are used.

Residual Analysis

The residual analysis has the purpose of detecting the presence of atypical observations and to evaluate the assumptions of the model, being able to include formal tests to detect departures from the assumptions of the considered model, as well as informal graphs to present general characteristics of the residuals.

Following Garay et al. 2017 and Arellano-Valle et al. 2012, in this work we considered the transformed martingal residuals $r_{MTi}$ proposed by Barros et al. 2010 as diagnostic tool to evaluate deviations from the postulated model for the response variable, as well as to detect the presence of atypical observations. The residuals are defined as

$$r_{MTi}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(r_{Mi})\sqrt{-2[r_{Mi}+{\delta }_i\log ({\delta }_i-r_{M_i})]},i=1,\dots ,n$$

where $r_{Mi}={\delta }_i+logS(y_i;\widehat{𝛉})$ is the martingal residual proposed by Ortega et al. 2003, where ${\delta }_i=0,1$ indicates whether the $i$th observation is censored or not, respectively, $$\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(r_{Mi})$$ denotes the sign of $r_{Mi}$ and $S(y_i;\widehat{𝛉})=P_{\widehat{𝛉}}(Y_i>y_i)$ represents the survival function evaluated at $y_i$, where $\widehat{𝛉}$ are the MLE for $𝛉$.

As suggested by Garay et al. 2017, this type of standardized residuals is used due to the fact that they are symmetrically distributed around zero, which facilitates the construction of the simulated envelopes with little computational effort and will be useful to detect an incorrect specification of the model, as well as the presence of observations atypical.

SIMULATIONS STUDIES

Simulation Study 1: Robustness of the Maximum Likelihood Estimates

In this section, we compare the performance of the estimates for PTCR model in the presence of outliers on the response variable. Following Garay et al. 2017 and Mattos et al. 2018 we performed a simulation study based on the NCR model. Specifically, we considered (12) with $${\mathbf{x}}_i^{\mathrm{\top }}=(1,x_i)$$ and $${\epsilon }_i\sim \mathrm{N}(0,{\eta }^2)$$ for $i=1,\dots ,n$. As in Garay et al. 2017 we generated 1000 artificial samples of size $n=300$, considering $𝛃={({\beta }_1,{\beta }_2)}^{\top }={(2,3.5)}^{\top }$, $\eta =2$ and fixing the left censoring level at $p=10,20$ and 30% (that is, 10, 20 and 30% of the observations in each data set were left censored, respectively). We generated independently the values $x_i$, for $i=1,\dots ,n$, from a uniform distribution on the interval (2, 20). These values were fixed throughout the simulations.

To assess how much the ML estimates are influenced by the presence of outliers, we replaced the observation $y_{150}$ by $y_{150}\left(\delta \right)=y_{150}+\delta$, with $\delta =1,2,\dots ,10$. Let ${\hat{\beta }}_i\left(\delta \right)$ and ${\hat{\beta }}_i$ be the ML estimates of ${\beta }_i$ with and without contamination, respectively, $i=1,2$. We are particularly interested in the relative changes

$$RC({\hat{\beta }}_i\left(\delta \right))=\left|({\hat{\beta }}_i\left(\delta \right)-{\hat{\beta }}_i)/{\hat{\beta }}_i\right|$$

We define the relative changes for $\eta$ analogously. For each replication we obtained the parameter estimates with and without outliers, under the PTCR model. Table II and Figure 2 depict the average values of the relative changes across all samples and different censoring levels.

Table II
Average relative changes on estimates for different contaminations d and censoring level p.

Cens.	$RC(\cdot )$	$\delta$
Level		1	2	3	4	5	6	7	8	9	10
	$RC\left({\beta }_1\right)$	0.00194	0.02610	0.01427	0.05336	0.08081	0.09474	0.12469	0.15654	0.15260	0.15643
10%	$RC\left({\beta }_2\right)$	0.00010	0.00004	0.00009	0.00022	0.00024	0.00032	0.00045	0.00060	0.00064	0.00080
	$RC\left(\eta \right)$	0.00046	0.01137	0.00780	0.02271	0.03341	0.03780	0.04594	0.05558	0.05101	0.04596
	$RC\left({\beta }_1\right)$	0.01337	0.01963	0.00945	0.06445	0.12935	0.10080	0.09735	0.07929	0.04515	0.03117
20%	$RC\left({\beta }_2\right)$	0.00010	0.00008	0.00030	0.00031	0.00002	0.00018	0.00040	0.00075	0.00124	0.00157
	$RC\left(\eta \right)$	0.00219	0.00020	0.00160	0.00849	0.00984	0.00615	0.00417	0.00024	0.00910	0.01333
	$RC\left({\beta }_1\right)$	0.02466	0.00670	0.10623	0.06087	0.09179	0.15043	0.05653	0.19465	0.03994	0.16353
30%	$RC\left({\beta }_2\right)$	0.00062	0.00190	0.00016	0.00109	0.00006	0.00066	0.00019	0.00015	0.00012	0.00123
	$RC\left(\eta \right)$	0.00176	0.00860	0.01234	0.00391	0.00881	0.01215	0.00324	0.01675	0.00467	0.01425

Figure 2
Simulation study 1. Average relative changes on estimates for different contaminations d and censoring level.

We observe that influence increases dramatically when $\delta$ increases for $p=10\%$, specially for the parameter. However, for the $p=20$ and 30%, these measures vary little, which indicates that PTCR model is more robust in these cases in the presence of discrepant observations.

Simulation Study 2: Asymptotic properties

To study the performance of the ML estimator $\widehat{𝛉}={({\widehat{𝛃}}^{\top },\hat{\eta },\hat{\alpha },\hat{\nu })}^{\top }$, a Monte Carlo simulation study with sample sizes $n=$150, 300, 750, and 1000 is presented. We considered the PTCR model defined in Section “Censored Power-Student-t Regression Model” with $${\mathbf{x}}_i^{\mathrm{\top }}=(1,x_i)$$ for $i=1,\dots ,n$. The true values of the parameters were taken as $𝛃=$(2,1.5)${}^{\top }$, $\eta =$1.5, $\alpha =$ 0.4, 2.5 and $\nu =$ 3.0. We also consider levels of censorship equal to $p=$ 0, 10, 25 and 45$\%$. The covariate $x_i$ was generated from a uniform distribution $$\mathrm{U}$$(0.1,20) as considered in Garay et al. 2017. For each combination of parameters, sample sizes and censorship levels, 2000 samples of the PTCR model were generated with errors $${\epsilon }_i\sim \mathrm{P}\mathrm{T}(0,\eta ,\alpha ,\nu )$$. To evaluate the performance of the estimators, the absolute value of the relative bias (RB) and the mean square error (MSE) were considered, they are given by

$$\mathrm{\text{RB}}\left({\hat{\theta }}_i\right)=\frac{1}{2000}\sum _{j=1}^{2000}(\frac{{\hat{\theta }}_i^{\left(j\right)}}{{\theta }_i}-1),\mathrm{\text{MSE}}\left({\hat{\theta }}_i\right)=\frac{1}{2000}\sum _{j=1}^{2000}{({\hat{\theta }}_i^{\left(j\right)}-{\theta }_i)}^2,$$

respectively, where ${\hat{\theta }}_i^{\left(j\right)}$ is the estimator of ${\theta }_i$ for the $j$th sample, for ${\theta }_i\in 𝛉={({𝛃}^{\top },\eta ,\alpha ,\nu )}^{\top }$. The ML estimates of the parameters were calculated by using the optim function of R Development Core Team 2018. The optimization of the likelihood function was done using iterative methods based on the Newton-Rapshon algorithm by using the score functions.

It can be seen from the Table III that RB and MSE tend to decrease when the value of $n$ increases, indicating that estimates based on the ML method have good asymptotic properties. That pattern is the same for the different levels of censorship of $p$ under consideration. Note that, when the sample sizes is $n=150$, the estimates for the ${\beta }_0$ parameter are unstable (in terms of MSE) because it is affected by the bias of the asymmetry parameter $\alpha$, however, when the sample size increases, the estimates become more stable. In general, this problem is very common in these types of models, see for example, Martínez-Flórez et al. 2013, so we recommend moderate and large sample sizes in these types of models

Table III
Performance of the ML estimators of ${\beta }_0$, ${\beta }_1$, η, α and ν for the PTCR model.

		$\alpha =0.40$								$\alpha =2.5$
		$p=0\%$		$p=10\%$		$p=25\%$		$p=45\%$		$p=0\%$		$p=10\%$		$p=25\%$		$p=45\%$
$n$	$\widehat{𝛉}$	RB	MSE	RB	MSE	RB	MSE	RB	MSE	RB	MSE	RB	MSE	RB	MSE	RB	MSE
150	${\hat{\beta }}_1$	0.010	0.982	0.019	1.127	0.022	0.746	0.050	0.657	0.092	0.700	0.132	0.813	0.138	0.913	0.155	1.200
	${\hat{\beta }}_2$	0.002	0.083	0.002	0.070	0.002	0.047	0.001	0.041	0.002	0.014	0.002	0.026	0.003	0.034	0.005	0.055
	$\hat{\eta }$	0.030	0.278	0.035	0.260	0.028	0.225	0.033	0.215	0.022	0.207	0.018	0.231	0.015	0.248	0.004	0.275
	$\hat{\alpha }$	0.039	0.119	0.040	0.112	0.041	0.102	0.044	0.098	0.309	1.778	0.301	1.878	0.297	1.913	0.315	2.020
	$\hat{\nu }$	0.309	2.673	0.296	2.428	0.255	2.125	0.223	1.980	0.188	1.735	0.196	1.840	0.211	2.057	0.216	2.228
300	${\hat{\beta }}_1$	0.010	0.839	0.025	0.787	0.014	0.504	0.025	0.441	0.061	0.552	0.086	0.601	0.094	0.689	0.122	0.922
	${\hat{\beta }}_2$	0.001	0.062	0.001	0.051	0.001	0.033	0.001	0.029	0.001	0.011	0.001	0.018	0.001	0.024	0.002	0.038
	$\hat{\eta }$	0.020	0.203	0.020	0.184	0.014	0.155	0.016	0.147	0.018	0.150	0.018	0.168	0.016	0.178	0.016	0.205
	$\hat{\alpha }$	0.035	0.091	0.031	0.080	0.026	0.072	0.026	0.068	0.195	1.223	0.184	1.343	0.206	1.472	0.262	1.715
	$\hat{\nu }$	0.186	1.739	0.160	1.610	0.123	1.263	0.102	1.127	0.093	0.973	0.107	1.025	0.109	1.026	0.128	1.175
750	${\hat{\beta }}_1$	0.002	0.588	0.003	0.481	0.006	0.303	0.009	0.263	0.024	0.219	0.037	0.346	0.044	0.415	0.064	0.574
	${\hat{\beta }}_2$	0.001	0.046	0.001	0.031	0.001	0.021	0.001	0.018	0.001	0.009	0.001	0.012	0.001	0.014	0.001	0.024
	$\hat{\eta }$	0.008	0.136	0.007	0.111	0.005	0.097	0.006	0.092	0.005	0.077	0.006	0.099	0.007	0.108	0.008	0.128
	$\hat{\alpha }$	0.009	0.042	0.009	0.045	0.008	0.042	0.008	0.040	0.084	0.599	0.070	0.640	0.089	0.781	0.128	1.024
	$\hat{\nu }$	0.052	0.715	0.041	0.602	0.034	0.534	0.029	0.502	0.020	0.512	0.035	0.498	0.040	0.531	0.048	0.603
1000	${\hat{\beta }}_1$	0.001	0.417	0.001	0.394	0.004	0.262	0.008	0.227	0.010	0.167	0.026	0.289	0.031	0.347	0.044	0.474
	${\hat{\beta }}_2$	0.001	0.031	0.001	0.025	0.001	0.017	0.001	0.015	0.001	0.009	0.001	0.010	0.001	0.012	0.001	0.020
	$\hat{\eta }$	0.001	0.091	0.002	0.095	0.003	0.082	0.004	0.078	0.003	0.069	0.005	0.083	0.005	0.093	0.006	0.107
	$\hat{\alpha }$	0.008	0.035	0.001	0.040	0.006	0.036	0.006	0.035	0.068	0.489	0.051	0.519	0.063	0.609	0.089	0.778
	$\hat{\nu }$	0.035	0.605	0.036	0.513	0.027	0.444	0.004	0.036	0.013	0.293	0.028	0.401	0.031	0.435	0.038	0.487

REAL DATA APPLICATIONS

Application 1: Wage rate

To illustrate the proposed model, we consider a data set described by Mroz 1987. The data set consists of 753 married white women with ages between 30 and 60 years old in 1975, with 428 women that worked at some point during that year. The response variable used in this application is the wage rate, which represents a measure of the wage of the housewife known as the average hourly earnings. In data set, we have that if the wage rates are set equal to zero, these wives did not work in 1975. Therefore, these observations are considered left-censored at zero.

The considerated covariates were the wife’s age $\left(X_1\right)$, years of schooling $\left(X_2\right)$, the number of children younger than six years old in the household $\left(X_3\right)$, and the number of children between six and nineteen years old $\left(X_4\right)$. These data were analyzed previously by Arellano-Valle et al. 2012 using a TCR model and later by Garay et al. 2017 using the Scale Mixture of Normal Censored Regression (SMNCR) models. We analyzed the data set by fitting a PTCR model and we compare our proposal with SMNCR models by Garay et al. 2017: Student-$t$ censored regression model (TCR) (Arellano-Valle et al. 2012), Slash censored regression model (SLCR), and normal censored regression model (NCR), that is, the usual tobit model. Table IV shows skewness and kurtosis index for complete data and also for uncensored observations. Notice that values for the skewness and kurtosis indexes justify using the PTCR.

Table IV
Statistical summary for wage rate data.

	$n$	Mean	SE	Skewness	Kurtosis
Complete Data	753	2.3746	3.2418	2.7778	15.7967
Non Censored Data	428	4.1777	3.3103	3.0801	16.8212

Table V presents parameter estimates together with their corresponding standard errors (SE) for the PTCR, TCR, SLCR and NCR models. To fit model of the PTCM, we use the R code in the Appendix C. Table VI presents some model selection criteria, together with the values of the log-likelihood. According to the AIC, BIC and AICc criteria, the PTCR model seems to yield a better fit to the Mroz’s data than the SMNCR models (TCR and SLCR models) and the usual Tobit model (NCR model), supporting the contention of a departure from symmetry of the errors. Also, the SE of the PTCR model are smaller than SE of the SMNCR and NCR models.

Table V
Parameters and standard errors (SE) of the PTCR, TCR, SLCR and NCR models fitted to Wage rate data.

	PTCR		TCR		SLCR		NCR
$𝛉$	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE
${\beta }_0$	-8.8412	0.3785	-1.0474	1.1549	-1.1445	1.4377	-2.7510	1.7334
${\beta }_1$	-0.1331	0.0004	-0.1108	0.0204	-0.1084	0.0232	-0.1046	0.0276
${\beta }_2$	0.4887	0.0050	0.6475	0.0658	0.6434	0.0635	0.7281	0.0831
${\beta }_3$	-3.4010	0.1869	-3.1636	0.3878	-3.0958	0.3785	-3.0264	0.4406
${\beta }_4$	-0.3053	0.0149	-0.2964	0.1247	-0.2946	0.1210	-0.2143	0.1527
$\eta$	4.1538	0.1320	3.2616	0.1731	2.5958	0.3071	4.5760	0.1697
$\nu$	6.5899	0.9014	4.2000	0.8068	1.4061	0.4361	–	–
$\alpha$	12.8760	0.9871	–	–	–	–	–	–

Table VI
Wage rate data. Model selection criteria.

Criterion	PTCR	TCR	SLCR	NCR
$\ell \left(𝛉\right)$	$-1413.942$	-1440.145	-1436.286	-1481.655
AIC	2843.883	2894.291	2886.573	2975.311
AICc	2828.077	2880.441	2872.723	2963.424
BIC	2880.876	2926.659	2918.941	3003.055

A more emphatic indication that an asymmetric model should be considered comes from testing the hypothesis a TCR model against an asymmetric (PTCR model), that is,

$$H_0:\alpha =1\mathrm{\text{versus}}{H}_1:\alpha \ne 1,$$

by using the likelihood ratio (LR) statistics, $-2log\left(\Lambda \right)=-2({\ell }_{TCR}\left(\widehat{𝛉}\right)-{\ell }_{PTCR}\left(\widehat{𝛉}\right))$, which for the data set under study, leads to $-2log\left(\Lambda \right)=52.41$, which is greater than the critical 5% value with one degree of freedom which is given by ${\chi }_{1,95\%}^2=3.84$. This is an indication that the PTCR model fits Mroz’s data better than the ordinary TCR model.

Finally, in order to verify if there is any incorrect specification in the assumptions of the fitted model, the simulated envelope graphs for the transformed martingal residuals are shown in Figure 3. This figure indicates that the PTCR model is, apparently, more suitable for the adjustment of this data than the SMNCR models. It can also be observed that the SMNCR models with heavy tails fit the data better than the NCR model, since there are few observations that are outside the envelopes.

Figure 3
Wage rate data. Envelopes of transformed martingale residuals for PTCR, TCR, SLCR and NCR models.

Application 2: Stellar Abundances Data

The second censored dataset is described in Santos et al. 2002 and are available in the R package astrodatR (Feigelson 2014) under the name Stellar abundances. These data were analyzed Mattos et al. 2018 by using the Scale Mixture of Skew Normal Censored Regression (SMSNCR) models. We analyzed the data set by fitting a PTCR model and again, we compare our proposal with SMNCR models by Garay et al. 2017.

The dataset consists of measurements for 68 solar-type stars and for our analysis we followed Mattos et al. 2018 and consider:

• log N(Be) as the response variable, which represents the log of the abundance of beryllium scaled to Sun’s abundance (i.e. the Sun has $logN\left(Be\right)=0.0$)

• $T_{\mathrm{\text{eff}}}/1000$ as the explanatory variable, which represents the effective stellar surface temperature (in kelvin).

In astronomical research, a previously identified sample of objects (stars, galaxies, quasars, X-ray sources, etc.) is observed at some new wavebands. According to Feigelson 2014, due to limited sensitivities, some objects may be undetected, leading to upper limits in their derived luminosities. For this dataset we have 12 left-censored data points, i.e. 12 undetected beryllium measurement, that represents 19.35% of observations.

Table VII presents the ML estimates for the parameters of the four models, i.e. PTCR, TCR, SLCR and NCR models, together with their corresponding standard errors. Table VIII compares the fit of the four models using the model selection criteria (AIC, AIC${}_c$ and BIC). Note that again the PTCR model with heavy tails have better fit than the TCR, SLCR and NCR models. The QQ-plots and envelopes for the martingale residuals are shown in Figure 4. This figure clearly indicates that the PTCR, TCR and SLCR models are more suitable for modeling the current data than the NCR model, since there are not observations falling outside the envelope.

Table VII
Parameters and standard errors (SE) of the PTCR, TCR, SLCR and NCR models fitted to stellar abundances data.

	PTCR		TCR		SLCR		NCR
$𝛉$	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE
${\beta }_0$	-1.9692	0.0658	-1.8923	0.1482	-0.5441	0.8733	-1.6478	0.5941
${\beta }_1$	0.5353	0.0021	0.5091	0.0046	0.2719	0.0262	0.4372	0.0180
$\eta$	0.0608	0.0002	0.0873	0.0004	0.1987	0.0254	0.3843	0.0014
$\nu$	2.0891	0.4530	0.8335	0.0419	1.0567	0.1652	–	–
$\alpha$	0.3265	0.0084	–	–	–	–	–	–

Table VIII
Stellar abundances data. Model selection criteria.

Criterion	PTCR	TCR	SLCR	NCR
$\ell \left(𝛉\right)$	3.688	-6.635	-25.096	-38.281
AIC	2.624	21.269	58.192	82.562
AICc	-6.408	13.905	50.827	76.937
BIC	13.722	30.148	67.069	89.220

Figure 4
Stellar abundances data. Envelopes of transformed martingale residuals for PTCR, TCR, SLCR and NCR models.

CONCLUSIONS

In this paper, an asymmetric alternative for the Student-$t$ censored regression model by Arellano-Valle et al. 2012 and SMNCR by Garay et al. 2017 has been developed. It is based on the new family of asymmetric and heavy-tailed power-$t$ distribution (Zhao & Kim 2016). Moreover, it follows that the ordinary tobit model and the Student-$t$ censored regression models are special cases. The observed and expected information matrix is analytically obtained, allowing for the direct implementation of the inference on this type of models. The problem of estimating the parameters in the model is dealt by using the maximum likelihood approach which is also used for developing large sample properties for the estimators. The likelihood ratio statistics can be used for testing the PTCR null hypothesis since the TCR model is special case of the model entertained. Applications to Wage rate data and Stellar Abundances Data indicate that the PTCR model can be a useful alternative to the TCR and SMNCR models.

The proposed PT distribution can be considered in the statistical models based on the scale mixtures of normal family to improve the fit the models such as Maleki & Nematollahi 2017. Also, the methodology of constructing the asymmetric distribution on the symmetric version of the Skew-Reflected-Gompertz distribution which recently introduced by Hosseinzadeh et al. 2019, can be considered as a future work for researchers.

ACKNOWLEDGMENTS

• AKAIKE H. 1974.A new look at statistical model identification. IEEE Trans Automat Contr AU-19(4): 716–722.
• ARELLANO-VALLE RB, CASTRO LM, GONZÁLEZ-FARÍAS G & MUNÕZ-GAJARDO KA. 2012. Student-$t$ Censored Regression Model: Properties and Inference. Stat Methods Appt 21(4): 453–473. doi:10.1007/s10260-012-0199-y
• AZZALINI A & CAPITANIO A. 2003. Distributions generated by perturbation of symmetry with emphasis on a multivariate skew $t$-distribution. J R Stat Soc Series B Stat Methodol 65(2): 367–389.
• BARROS M, GALEA M, GONZÁLEZ M & LEIVA V. 2010. Influence Diagnostics in the Tobit Censored Response Model. Stat Methods Appt 19(3): 379–397. doi:10.1007/s10260-010-0135-y.
• FAIR RC. 1978.A theory of extramarital affairs. J Polit Econ 86(1): 45–61.
• FEIGELSON ED. 2014. astrodatR: Astronomical Data. Available at. URL https://cran.r-project.org/web/packages/astrodatR/ R package v. 0.1.
  » https://cran.r-project.org/web/packages/astrodatR/
• GARAY AM, LACHOS VH, BOLFARINE H & CABRAL CRB. 2017. Linear Censored Regression Models with Scale Mixtures of Normal Distributions. Stat Pap 58(1): 247–278. doi:10.1007/s00362-015-0696-9.
• GUPTA RD & GUPTA RC. 2008. Analyzing skewed data by power–normal model. Test 17: 197–210.
• HOSSEINZADEH A, MALEKI M, KHODADADI Z & CONTRERAS-REYES JE. 2019. The Skew-Reflected-Gompertz distribution for analyzing the symmetric and asymmetric data. J Comput Appl Math 349: 132–141.
• MALEKI M & NEMATOLLAHI AR. 2017. Autoregressive Models with Mixture of Scale Mixtures of Gaussian innovations. Iran J Sci Technol Trans A Sci 41(4): 1099–1107.
• MARTÍNEZ-FLÓREZ G, BOLFARINE H & GÓMEZ HW. 2013. The alpha–power tobit model. Commun Stat Theory Methods 42(4): 633–643.
• MATTOS T, GARAY AM & LACHOS VH. 2018. Likelihood-based inference for censored linear regression models with scale mixtures of skew-normal distributions. J Appl Stat 45(11): 2039–2066.
• MOULTON LH & HALSEY NA. 1995. A Mixture Model With Detection Limits for Regression Analyses of Antibody Response to Vaccine. Biometrics 51: 1570–1578.
• MROZ TA. 1987. The Sensitivity of an Empirical Model of Married Women’s Hours of Work to Economic and Statistical Assumptions. Econometrica 55(4): 765–799. doi:10.2307/1911029.
• OLSEN RJ. 1978. Note on the Uniqueness of the Maximum Likelihood Estimator for the Tobit Model. Econometrica 46(5): 1211–1215.
• ORTEGA EM, BOLFARINE H & PAULA GA. 2003. Influence Diagnostics in Generalized Log-Gamma Regression Models. Comput Stat Data Anal 42: 165–186.
• PEWSEY A, GÓMEZ HW & BOLFARINE H. 2012. Likelihood–based inference for power distributions. Test 21(4): 775–789.
• R DEVELOPMENT CORE TEAM. 2018. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. URL http://www.R-project.org ISBN 3-900051-07-0.
  » http://www.R-project.org
• SANTOS N, LÓPEZ RG, ISRAELIAN G, MAYOR M, REBOLO R, GARCÍA-GIL A, TAORO MP DE & RANDICH S. 2002. Beryllium Abundances in Stars Hosting Giant Planets. Astron Astrophys 386: 1028–1038.
• SCHWARZ G. 1978. Estimating the dimension of a model. Ann Stat 6(2): 461–464.
• SUGIURA N. 1978. Further Analysis of the Data by Akaike’s Information Criterion and the Finite Corrections. Commun Stat Theory Methods 7(1): 13–26. doi:10.1080/03610927808827599.
• TOBIN J. 1958.Estimation of relationship for limited dependent variables. Econometrica 26(1): 24–36.
• ZHAO J & KIM HM. 2016.Power $t$ distribution. Ommun Stat Appl Methods 23(4): 321–334.

APPENDIX A LEMMAS

Lemma 1. Let $$Z\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$, then $E[Z^k\mid Z+c>0]=m_k(c;\nu )$, with

$$m_k(c;\nu )=\frac{1}{1-\{F_T(-c;\nu ){\}}^{\alpha }}{\int }_{F_T(-c;\nu )}^1{\left[F_T^{-1}(u;\nu )\right]}^k\alpha u^{\alpha -1}du.$$

where $F_T^{-1}(\cdot ;\nu )$ is the inverse of $F_T(\cdot ;\nu )$.

Lemma 2. Let $$Z\sim \mathrm{P}\mathrm{T}(\alpha ,\nu )$$, and define $r(Z;\nu )=f_T(Z;\nu )/F_T(Z;\nu )$. Then

• (i)$E\{{\left[r(Z;\nu )\right]}^k{(1+\frac{Z^2}{\nu })}^{-m/2}{Z}^n\mid Z+c>0\}={\displaystyle \frac{{\left[f_T(0;\nu )\right]}^k}{1-{\left\{F_T(-c;\nu )\right\}}^{\alpha }}}{a}_{kmn}(-c;\nu )$, where

$$\begin{array}{cc}\hfill a_{kmn}(-c;\nu )=& {\int }_{F_T(-c;\nu )}^1{\left[F_T^{-1}(u;\nu )\right]}^n{\{1+\frac{{\left[F_T^{-1}(u;\nu )\right]}^2}{\nu }\}}^{-\frac{k(\nu +1)+m}{2}}\alpha u^{\alpha -k-1}du\hfill \end{array}$$

• (ii)$E\{{\left[r(Z;\nu )\right]}^{k}log(1+\frac{Z^2}{\nu })Z^n\mid Z+c>0\}={\displaystyle \frac{{\left[f_T(0;\nu )\right]}^k}{1-{\left\{F_T(-c;\nu )\right\}}^{\alpha }}}{a}_{kn}(-c;\nu )$, where

  $$a_{kn}(-c;\nu )={\int }_{F_T(-c;\nu )}^{1}log\{1+\frac{{\left[F_T^{-1}(u;\nu )\right]}^2}{\nu }\}{\left[F_T^{-1}(u;\nu )\right]}^n\alpha u^{\alpha -k-1}du$$

• (iii) $E\{{\left[r(Z;\nu )\right]}^k{\left[b_{01}(Z;\nu )\right]}^m{Z}^n\mid Z+c>0\}={\displaystyle \frac{{\left[f_T(0;\nu )\right]}^k}{1-{\left\{F_T(-c;\nu )\right\}}^{\alpha }}}{b}_{kmn}(-c;\nu )$, where

$$b_{kmn}(-c;\nu )={\int }_{F_T(-c;\nu )}^1{\left[b_{01}\left(F_T^{-1}(u;\nu )\right)\right]}^m{\left[F_T^{-1}(u;\nu )\right]}^n\alpha u^{\alpha -k-1}du$$

• con $b_{01}(c;\nu )={\int }_{-\infty }^c\{log(1+\frac{s^2}{\nu })\}\frac{f_T(s;\nu )}{F_T(c;\nu )}ds$.

The proof of the Lemmas 1 and 2 are straightforward and they follow directly for the definition of expected value.

APPENDIX B INFORMATION MATRICES FOR PTCR MODEL

Observed Information Matrix

The elements of the observed information matrix $𝐉\left(𝛗\right)$ for the PTCR model are given by

$$\begin{array}{rl}\mathit{j}\mathit{\gamma }{\mathit{\gamma }}^{\mathrm{\top }}& =\alpha \sum _{i=1}^n(1-d_i)\{r(-c_i;\nu )-\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i^2}{\nu })}^{-1}{c}_i\}r(-c_i;\nu ){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ & +\left(\frac{\nu +1}{\nu }\right)\sum _{i=1}^n{d}_i\{{(1+\frac{z_i^2}{\nu })}^{-1}-\frac{2}{\nu }{(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^2\}{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ & +(\alpha -1)\sum _{i=1}^n{d}_i\{r(z_i;\nu )+\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i\}r(z_i;\nu ){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ \mathit{j}\mathit{\gamma }\sigma & =-\left(\frac{\nu +1}{\nu }\right)\sum _{i=1}^n{d}_i\{{(1+\frac{z_i^2}{\nu })}^{-1}-\frac{2}{\nu }{(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^2\}{y}_i{\mathbf{x}}_i\\ & -(\alpha -1)\sum _{i=1}^n{d}_i\{r(z_i;\nu )+\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i\}\\ & \times r(z_i;\nu )y_i{\mathbf{x}}_i\\ \mathit{j}\mathit{\gamma }\nu & =\frac{\alpha }{2}\sum _{i=1}^n(1-d_i)\{(\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right))(1+R(c_i;\nu ))\\ & -\frac{1}{\nu }(1+c_{i}r(c_i;\nu )R(c_i;\nu ))-R(c_i;\nu )b_{01}(c_i;\nu )-\log (1+\frac{c_i^2}{\nu })\\ & +\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i^2}{\nu })}^{-1}\frac{c_i^2}{\nu }\}r(-c_i;\nu ){\mathbf{x}}_i\\ & +\frac{1}{{\nu }^2}\sum _{i=1}^n{d}_i\{{(1+\frac{z_i^2}{\nu })}^{-1}{z}_i-\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^3\}{\mathbf{x}}_i\\ & +\frac{\alpha -1}{2}\sum _{i=1}^n{d}_i\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }+b_{01}(z_i;\nu )\\ & -\log (1+\frac{z_i^2}{\nu })-\frac{1}{\nu }(1-z_{i}r(z_i;\nu ))\}r(z_i;\nu ){\mathbf{x}}_i\end{array}$$

$$\begin{array}{rl}\mathit{j}\mathit{\gamma }\alpha & =\sum _{i=1}^n(1-d_i)r(-c_i;\nu ){\mathbf{x}}_i+\sum _{i=1}^n{d}_{i}r(z_i;\nu ){\mathbf{x}}_i\\ j\sigma \sigma & =\sum _{i=1}^n{d}_i\{\frac{1}{{\sigma }^2}+\left(\frac{\nu +1}{\nu }\right)[{(1+\frac{z_i^2}{\nu })}^{-1}-\frac{2}{\nu }{(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^2]y_i^2\}\\ & +(\alpha -1)\sum _{i=1}^n{I}_i\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}{z}_i+r(z_i;\nu )\}r(z_i;\nu )y_i^2\\ j\sigma \nu & =-\frac{1}{{\nu }^2}\sum _{i=1}^n{d}_i\{{(1+\frac{z_i^2}{\nu })}^{-1}{z}_i-\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^3\}{y}_i\\ & -\frac{\alpha -1}{2}\sum _{i=1}^n{d}_i\{\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }+b_{01}(z_i\nu )\\ & -\log (1+\frac{z_i^2}{\nu })-\frac{1}{\nu }(1-z_{i}r(z_i;\nu ))\}r(z_i;\nu )y_i\\ j\sigma \alpha & =-\sum _{i=1}^n{d}_{i}r(z_i;\nu )y_i\\ j\nu \alpha & =-\frac{1}{2}\sum _{i=1}^n(1-d_i)\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )-\frac{c_i}{\nu }r(c_i;\nu )\}\\ & \times R(c_i;\nu )+\frac{1}{2}\sum _{i=1}^n{d}_i\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(z_i;\nu )\\ & -\frac{z_i}{\nu }r(z_i;\nu )\}\\ j\alpha \alpha & =-\sum _{i=1}^n\frac{d_i}{{\alpha }^2}\end{array}$$

$$\begin{array}{rl}j\nu \nu & =+\frac{\alpha }{4}\sum _{i=1}^n(1-d_i)\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )\\ & {-\frac{c_i}{\nu }r(c_i\mid \nu )\}}^{2}R(c_i;\nu )(1+R(c_i;\nu ))\\ & -\frac{\alpha }{4}\sum _{i=1}^n(1-d_i)\{{\psi }_1\left(\frac{\nu }{2}\right)-{\psi }_1\left(\frac{\nu +1}{2}\right)-\frac{2}{\nu (\nu +1)}+b_{01}^2(c_i;\nu )\\ & -\frac{1}{\nu }{b}_{01}(c_i;\nu )+\frac{1}{\nu }(2b_{01}(c_i;\nu )-\frac{\nu +3}{\nu (\nu +1)})c_{i}r(c_i;\nu )\\ & +\frac{1}{\nu +2}\frac{F_T(c_{2i};\nu +2)}{F_T(c_i;\nu )}{b}_{21}(c_{2i};\nu +2)-b_{02}(c_i;\nu )\\ & +\frac{1}{\nu }[\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i^2}{\nu })}^{-1}\frac{c_i^2}{\nu }-\log (1+\frac{c_i^2}{\nu })+\frac{c_i}{\nu }r(c_i;\nu )]\\ & \times c_{i}r(c_i;\nu )\}R(c_i;\nu )\\ & -\frac{1}{4}\sum _{i=1}^n{d}_i\{{\psi }_1\left(\frac{\nu +1}{2}\right)-{\psi }_1\left(\frac{\nu }{2}\right)+\frac{2}{{\nu }^2}-\frac{4}{{\nu }^3}{(1+\frac{z_i^2}{\nu })}^{-1}{z}_i^2\\ & +\frac{2}{{\nu }^3}\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-2}{z}_i^4\}\\ & +\frac{\alpha -1}{4}\sum _{i=1}^n{d}_i\{{\psi }_1\left(\frac{\nu }{2}\right)-{\psi }_1\left(\frac{\nu +1}{2}\right)-\frac{2}{\nu (\nu +1)}+b_{01}^2(z_i;\nu )\\ & -\frac{1}{\nu }{b}_{01}(z_i;\nu )+\frac{1}{\nu }(2b_{01}(z_i;\nu )-\frac{\nu +3}{\nu (\nu +1)})z_{i}r(z_i;\nu )\\ & +\frac{1}{\nu +2}\frac{F_T(z_{2i};\nu +2)}{F_T(z_i;\nu )}{b}_{21}(z_{2i};\nu +2)-b_{02}(z_i;\nu )\\ & +\frac{1}{\nu }[\left(\frac{\nu +1}{\nu }\right){(1+\frac{z_i^2}{\nu })}^{-1}\frac{z_i^2}{\nu }\\ & -\log (1+\frac{z_i^2}{\nu })+\frac{z_i}{\nu }r(z_i;\nu )]z_{i}r(z_i;\nu )\}\end{array}$$

Expected Information Matrix

The expected information matrix $${\mathbf{I}}_{\mathit{\phi }}$$ is obtained by taking $${\mathbf{I}}_{\mathit{\phi }}=E[{\mathbf{J}}_{\mathit{\phi }}]$$, which involves the calculation of truncated moments that have no closed form and must be obtained numerically. The elements $I_{{𝛗}_i{𝛗}_j}$ of the expected information matrix are given by

$$\begin{array}{rl}{\mathit{I}}_{\mathit{\gamma }\mathit{\gamma }}=& \sum _{i=1}^n\{r(-c_i;\nu )-\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i}{\nu })}^{-1}{c}_i\}{f}_{PT}(-c_i;\alpha ,\nu ){\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ & +\left(\frac{\nu +1}{\nu }\right)\sum _{i=1}^n\{a_{020}(-c_i;\nu )-\frac{2}{\nu }{a}_{042}(-c_i;\nu )\}{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ & +(\alpha -1)f_T(0;\nu )\sum _{i=1}^n\{\left(\frac{\nu +1}{\nu }\right)a_{121}(-c_i;\nu )+f_T(0;\nu )a_{200}(-c_i;\nu )\}{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}\\ {\mathit{I}}_{\mathit{\gamma }\sigma }=& -\frac{1}{\sigma }\left(\frac{\nu +1}{\nu }\right)\sum _{i=1}^n\{a_{021}(-c_i;\nu )+c_i{a}_{020}(-c_i;\nu )\}{\mathbf{x}}_i\\ & +\frac{2}{\sigma }\left(\frac{\nu +1}{{\nu }^2}\right)\sum _{i=1}^n\{a_{043}(-c_i;\nu )+c_i{a}_{042}(-c_i;\nu )\}{\mathbf{x}}_i\\ & -(f_T(0;\nu ){)}^2\left(\frac{\alpha -1}{\sigma }\right)\sum _{i=1}^n\{a_{201}(-c_i;\nu )+c_i{a}_{200}(-c_i;\nu )\}{\mathbf{x}}_i\\ & -f_T(0;\nu )\left(\frac{\alpha -1}{\sigma }\right)\left(\frac{\nu +1}{\nu }\right)\sum _{i=1}^n\{a_{122}(-c_i;d\nu )+c_i{a}_{121}(-c_i;\nu )\}{\mathbf{x}}_i\\ {\mathit{I}}_{\mathit{\gamma }\nu }=& \frac{\alpha }{2}\sum _{i=1}^n\{(\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right))(1+R(c_i;\nu ))\\ & -\frac{1}{\nu }(1+c_{i}r(c_i;\nu )R(c_i;\nu ))-R(c_i;\nu )b_{01}(c_i;\nu )-\log (1+\frac{c_i^2}{\nu })\\ & +\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i^2}{\nu })}^{-1}\frac{c_i^2}{\nu }\}{f}_{PT}(-c_i;\alpha ,\nu ){\mathbf{x}}_i\\ & +\frac{1}{{\nu }^2}\sum _{i=1}^n\{a_{021}(-c_i;\nu )-\left(\frac{\nu +1}{\nu }\right)a_{043}(-c_i;\nu )\}{\mathbf{x}}_i\\ & +f_T(0;\nu )\left(\frac{\alpha -1}{2}\right)\sum _{i=1}^n\{\left(\frac{\nu +1}{\nu }\right)a_{122}(-c_i;\nu )-a_{10}(-c_i;\nu )\\ & +b_{110}(-c_i;\nu )-\frac{1}{\nu }(a_{100}(-c_i;\nu )-f_T(0;\nu )a_{201}(-c_i;\nu ))\}{\mathbf{x}}_i\end{array}$$

$$\begin{array}{rl}{\mathit{I}}_{\mathit{\gamma }\alpha }=& \frac{1}{\alpha }\sum _{i=1}^n{f}_{PT}(-c_i;\alpha ,\nu ){\mathbf{x}}_i+f_T(0;\nu )\sum _{i=1}^n{a}_{100}(-c_i;\nu ){\mathbf{x}}_i\\ I_{\sigma \sigma }=& \frac{1}{{\sigma }^2}\sum _{i=1}^n\{1+\left(\frac{\nu +1}{\nu }\right)[a_{022}(-c_i;\nu )+2c_i{a}_{021}(-c_i;\nu )\\ & +c_i^2{a}_{020}(-c_i;\nu )-\frac{2}{\nu }(a_{044}(-c_i;\nu )+2c_i{a}_{043}(-c_i;\nu )\\ & +c_i^2{a}_{042}(-c_i;\nu )\left)\right]\}\\ & +f_T(0;\nu )\left(\frac{\alpha -1}{{\sigma }^2}\right)\sum _{i=1}^n\{\left(\frac{\nu +1}{\nu }\right)[a_{123}(-c_i;\nu )+2c_i{a}_{122}(-c_i;\nu )\\ & +c_i^2{a}_{121}(-c_i;\nu )]+f_T(0;d\nu )[a_{202}(-c_i;\nu )+2c_i{a}_{201}(-c_i;\nu )\\ & +c_i^2{a}_{200}(-c_i;\nu )\left]\right\}\\ I_{\sigma \alpha }=& -\frac{f_T(0;\nu )}{\sigma }\sum _{i=1}^n\{a_{101}(-c_i\nu )+c_i{a}_{100}(-c_i;\nu )\}\\ I_{\alpha \alpha }=& \frac{1}{{\alpha }^2}\sum _{i=1}^n(1-{\{F_T(-c_i;\nu )\}}^{\alpha })\\ I_{\sigma \nu }=& -\frac{1}{\sigma {\nu }^2}\sum _{i=1}^n\{a_{022}(-c_i;\nu )+c_i{a}_{021}(-c_i;\nu )\\ & -\left(\frac{\nu +1}{\nu }\right)[a_{044}(-c_i;\nu )+c_i{a}_{043}(-c_i;\nu )]\}{\{F_T(-c_i;\nu )\}}^{\alpha }\\ & -f_T(0;\nu )\left(\frac{\alpha -1}{2\sigma }\right)\sum _{i=1}^n\{\left(\frac{\nu +1}{{\nu }^2}\right)[a_{123}(-c_i;\nu )+c_i{a}_{122}(-c_i;\nu )]\\ & -a_{11}(-c_i;\nu )-c_i{a}_{10}(-c_i;\nu )+b_{111}(-c_i;\nu )+c_i{b}_{110}(-c_i;\nu )\\ & -\frac{1}{\nu }[a_{101}(-c_i;\nu )+c_i{a}_{100}(-c_i;\nu )]\\ & +\frac{1}{\nu }[f_T(0;\nu )a_{202}(-c_i;\nu )+c_i{f}_T(0;\nu )a_{201}(-c_i;\nu )]\}\\ I_{\nu \alpha }=& \frac{1}{2\alpha }\sum _{i=1}^n\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )-\frac{c_i}{\nu }r(c_i;\nu )\}\\ & \times {(r(-c_i;\nu ))}^{-1}{f}_{PT}(-c_i;\alpha ,\nu )\\ & -\frac{1}{2}\sum _{i=1}^n\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-\frac{b_{010}(-c_i;\nu )}{1-{\{F_T(-c_i;\nu )\}}^{\alpha }}\\ & -\frac{f_T(0;\nu )a_{101}(-c_i;\nu )}{\nu (1-{\{F_T(-c_i;\nu )\}}^{\alpha })}\}(1-{\{F_T(-c_i;\nu )\}}^{\alpha })\end{array}$$

$$\begin{array}{rl}{I}_{\nu \nu }=& \frac{\alpha }{4}\sum _{i=1}^n{\{\psi \left(\frac{\nu +1}{2}\right)-\psi \left(\frac{\nu }{2}\right)-b_{01}(c_i;\nu )-\frac{c_i}{\nu }r(c_i;\nu )\}}^2\\ & \times R(c_i;\nu )(1+R(c_i;\nu )){\{F_T(-c_i;\nu )\}}^{\alpha }\\ & -\frac{\alpha }{4}\sum _{i=1}^n\{{\psi }_1\left(\frac{\nu }{2}\right)-{\psi }_1\left(\frac{\nu +1}{2}\right)-\frac{2}{\nu (\nu +1)}+b_{01}^2(c_i;\nu )\\ & -\frac{1}{\nu }{b}_{01}(c_i;\nu )+\frac{1}{\nu }(2b_{01}(c_i;\nu )-\frac{\nu +3}{\nu (\nu +1)})c_{i}r(c_i;\nu )\\ & -b_{02}(c_i;\nu )+\frac{F_T(c_{2i};\nu +2)b_{21}(c_{2i};\nu +2)}{(\nu +2)F_T(c_{2i};\nu )}\\ & +\frac{1}{\nu }[\left(\frac{\nu +1}{\nu }\right){(1+\frac{c_i^2}{\nu })}^{-1}\frac{c_i^2}{\nu }\\ & -\log (1+\frac{c_i^2}{\nu })+\frac{c_i}{\nu }r(c_i;\nu )]c_{i}r(c_i;\nu )\}R(c_i;\nu ){\{F_T(-c_i;\nu )\}}^{\alpha }\\ & -\frac{1}{4}\sum _{i=1}^n\{{\psi }_1\left(\frac{\nu +1}{2}\right)-{\psi }_1\left(\frac{\nu }{2}\right)+\frac{2}{{\nu }^2}-\frac{4a_{022}(-c_i;\nu )}{{\nu }^3(1-{\{F_T(-c_i;\nu )\}}^{\alpha })}\\ & +\frac{2}{{\nu }^3}\left(\frac{\nu +1}{\nu }\right)\frac{a_{044}(-c_i;\nu )}{(1-{\{F_T(-c_i;\nu )\}}^{\alpha })}\}(1-{\{F_T(-c_i;\nu )\}}^{\alpha })\\ & +\frac{\alpha -1}{4}\sum _{i=1}^n\{{\psi }_1\left(\frac{\nu }{2}\right)-{\psi }_1\left(\frac{\nu +1}{2}\right)-\frac{2}{\nu (\nu +1)}\\ & +\frac{b_{020}(-c_i;\nu )}{1-{\{F_T(-c_i;\nu )\}}^{\alpha }}-\frac{b_{010}(-c_i;\nu )}{\nu (1-{\{F_T(-c_i;\nu )\}}^{\alpha })}\\ & +\frac{1}{\nu }[\frac{2f_T(0;\nu )b_{111}(-c_i;\nu )}{1-{\{F_T(-c_i\mid \nu )\}}^{\alpha }}-\frac{\nu +3}{\nu (\nu +1)}\frac{f_T(0;\nu )a_{101}(-c_i;\nu )}{(1-{\{F_T(-c_i;\nu )\}}^{\alpha })}]\\ & -\mathrm{E}[b_{02}(z_i;\nu )]+\left(\frac{1}{\nu +2}\right)\mathrm{E}[\frac{F_T(z_{2i};\nu +2)}{F_T(z_i;\nu )}{b}_{21}(z_{2i};\nu +2)]\\ & +\frac{f_T(0;\nu )}{\nu (1-{\{F_T(-c_i;\nu )\}}^{\alpha })}[\left(\frac{\nu +1}{{\nu }^2}\right)a_{123}(-c_i;\nu )-a_{11}(-c_i;\nu )\\ & +\frac{1}{\nu }{f}_T(0;\nu )a_{202}(-c_i;\nu )\left]\right\}(1-{\{F_T(-c_i;\nu )\}}^{\alpha })\end{array}$$

where $b_{01}$, $b_{02}$ and $b_{21}$ are given by (8) and must be obtained numerically. The components $a_{kmn}$, $a_{kn}$ and $b_{kmn}$, are given in the Lemman 2.

APPENDIX C R CODE TO FIT THE PTCR MODEL

Publication Dates

• Publication in this collection
  11 Oct 2021
• Date of issue
  2021

History

• Received
  11 Aug 2019
• Accepted
  22 Oct 2019