8.2 Regression with censoring

In defining whether an observation is censored or not, we introduced indicator variables to identify whether a particular measurement arises due to censoring. This then raises the question: can the censoring indicator variables be used as an explanatory variable in a linear regression model?

To examine this consider the prisoners data set that contains height measurements of 1000 prisoners according to the wall-chart illustrated in Figure 8.3 (Link). The data set also contains two categorical variables: the prisoner’s gender and whether the prisoner’s height was ‘O’bserved, ‘L’eft censored if they are shorter than 60 inches or ‘R’ight censored of they are taller than 72 inches.

Let $Y_i^{*}$ denote the true height of prisoner $i$, which we wish to model using simple linear regression model

$$Y_i^{*}\sim \mathrm{N}({\mu }_i,{\sigma }^2)\mathit{\hspace{1em}}\mathrm{with}\mathit{\hspace{1em}}{\mu }_i={𝜷}^{\prime }{𝐱}_i\mathit{\hspace{1em}}\mathrm{for}i=1,\mathrm{\dots },n,$$

where ${𝐱}_i$ is a vector of explanatory variables for prisoner $i$. However, the recorded height $y_i$ is restricted to the interval $[a,b]$ such that:

Here, $a=60$ inches is the left censoring point and $b=72$ inches is the right censoring point.

In the following extract, the data set prisoners is read into R and a few summaries are evaluated.

prisoners <- read.table("prisoners.dat")
range(prisoners$height)
[1] 60 72
table( prisoners$gender, prisoners$censor )
      O   L   R
  M 564   0 127
  F 304   5   0

Here we see that there are 868 observed height measurements with 5 left censored events at 60 inches, all occurring for female prisoners, and 127 right censored events at 72 inches, all occurring for male prisoners.

If we then wish to use gender and the censor categorical variable as an explanatory variable in a regression model, then we can estimate the co-efficients as follows:

Model <- glm(height ~ 1 + gender + censor, data = prisoners,
  family=Gaussian)
summary(Model)

Deviance Residuals:
   Min      1Q  Median      3Q     Max
-6.502  -0.940   0.000   1.166   6.882

Coefficients:
            Estimate Std. Error  t value Pr(>|t|)
(Intercept) 69.40877    0.06824 1017.185  < 2e-16 ***
genderF     -5.21103    0.11530  -45.195  < 2e-16 ***
censorL     -4.19773    0.73065   -5.745 1.22e-08 ***
censorR      2.59123    0.15917   16.280  < 2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

(Dispersion parameter for gaussian family taken to be 2.626082)

    Null deviance: 10470.5  on 999  degrees of freedom
Residual deviance:  2615.6  on 996  degrees of freedom

 
Exercise 8.65
Calculate the expected height for all gender and censoring combinations. Comment on what the expected values show.

 

Using the censor variable in a regression model results in identifying the censoring points. Furthermore, the estimates of interest (expected height of male and female prisoners) are only derived from data that is not censored. This demonstrates the difference between identifying the best fitting model verses the usefulness of the model.

To overcome these issues, we must examine the likelihood function. Let $I_i^a$ and $I_i^b$ be indicator variables identifying whether the measurement is left or right censored respectively:

If the $i$th prisoner’s height is observed ($I_i^a=0$ and $I_i^b=0$) then the contribution to the likelihood is:

$$L_i(𝜷)=\frac{1}{\sigma }\varphi \left(\frac{y_i-{\mu }_i}{\sigma }\right)=\frac{1}{\sigma }\varphi \left(\frac{y_i-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right),$$

where $\varphi (.)$ denotes the standard normal pdf. If instead the measured height of the $i$th prisoner is left-censored ($I_i^a=1$ and $I_i^b=0$), then their true height could take any value less than $a$. The corresponding contribution to the likelihood is:

$$L_i(𝜷)=\mathrm{Pr}(Y_i^{*}\le a)=\mathrm{\Phi }\left(\frac{a-{\mu }_i}{\sigma }\right)=\mathrm{\Phi }\left(\frac{a-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right),$$

where $\mathrm{\Phi }(.)$ denotes the standard normal cdf. Finally, if the measured height of the $i$th prisoner is instead right-censored ($I_i^a=0$ and $I_i^b=1$), then their true height could take any value greater than $b$. The corresponding contribution to the likelihood is:

$$L_i(𝜷)=\mathrm{Pr}(Y_i^{*}>b)=1-\mathrm{Pr}(Y_i^{*}\le b)=1-\mathrm{\Phi }\left(\frac{b-{\mu }_i}{\sigma }\right)=\mathrm{\Phi }\left(-\frac{b-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right).$$

Combining these components obtains the likelihood function:

$$L(𝜷)=\prod _{i=1}^n{\left\{\mathrm{\Phi }\left(\frac{a-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right)\right\}}^{I_i^a}\times {\left\{\mathrm{\Phi }\left(-\frac{b-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right)\right\}}^{I_i^b}\times {\left\{\frac{1}{\sigma }\varphi \left(\frac{y_i-{𝜷}^{\prime }{𝐱}_i}{\sigma }\right)\right\}}^{1-I_i^a-I_i^b}$$

It is not possible to evaluate the maximum likelihood estimates analytically, so we must adopt a numerical approach.

  censReg( formula, left = 0, right = Inf, data, start = NULL )

library("censReg")
Model <- censReg(height ~ -1 + gender, left=60, right=72, data=prisoners)
summary(Model)

Observations:
         Total  Left-censored     Uncensored Right-censored
          1000              5            868            127

Coefficients:
         Estimate Std. error t value Pr(> t)
genderM  70.09065    0.08057  869.93  <2e-16 ***
genderF  64.11728    0.11780  544.30  <2e-16 ***
logSigma  0.72694    0.02476   29.36  <2e-16 ***
---
Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

Newton-Raphson maximisation, 4 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -2003.058 on 3 Df