Date Geändert

Introduction

Ever since I have read the refresher of Tjalling Jager on differential equations and likelihood functions, I have been irritated about the common way to express the likelihood of data following a normal distribution using the density function of the normal distribution. The statement that made me wonder was "note that we have seamlessly gone from probabilities to probability densities; this should not concern us here" (p. 23).

It seemed to me that the likelihood should rather be expressed as the integral of the probability density over the interval that is covered by a certain representation of a numeric value. For example, if an observation was reported to be "1.2", the associated interval would range from 1.15 to 1.25.

Only recently, when reading in the documentation of the mlt package by Hothorn (2018), I noticed that a similar argument about the approximation of the likelihood of continuous random variables has been mad earlier (Lindsey, 1999).

The case of a single observation ...

In order to put the interval idea to work, we start with a single observation. For a continuous random variable \(Y\) with a known normal distribution, the probability of a single observation to be within the interval \(] \underline{y}, \overline{y} ]\) is

$$ \mathbb{P}(\underline{y} < Y \le \overline{y} \mid \mu, \sigma) = \Phi\left(\frac{\overline{y} - \mu}{\sigma}\right) - \Phi\left(\frac{\underline{y} - \mu}{\sigma}\right) $$

Using our example of an observed value of 1.2, assuming that it implicitly means that the true value is in the interval ]1.15,1.25], we can define a function to calculate this probability in R as:

prob_obs <- function(lower, upper, mean = 0, sd = 1) {
  pnorm(upper, mean, sd) - pnorm(lower, mean, sd)
}
(p_1.2 <- prob_obs(1.15, 1.25))
## [1] 0.01942216
log(p_1.2)
## [1] -3.94134

The probability of this observation, assuming a standard normal distribution for the underlying statistical population is obtained to be about 0.0194, which is very close to the probability density at the value of 1.2 times the width of the interval. If we simply approximate the probability by the probability density, the factor by which this is off (and thus the shift in the log-likelihood) will depend on the units chosen for the observation.

Note that we have not distinguished between inclusive and exclusive bounds of the interval, as the numeric representation of the bounds in R will generally be of sufficient precision so that the difference between the closest match to a bound (e.g. 1.15) and the next possible value can be considered negligible for this analysis.

In order to obtain the parameters of the normal distribution that is most consistent with our observation, we can then proceed to interpret the probability of the observation given above as the likelihood of the model parameters

$$ L(\mu, \sigma | y) = \mathbb{P}(\underline{y} < Y \le \overline{y} \mid \mu, \sigma) $$

and maximise it to obtain maximum likelihood estimates (MLE) \(\hat{\mu}\) and \(\hat{\sigma}\) for the distribution parameters by minimizing the negative log-likelihood.

log_lik <- function(mean, sd, data) {
  sum(log(prob_obs(data$lower, data$upper, mean, sd)))
}
objective <- function(theta, data) {
  - log_lik(theta["mean"], theta["sd"], data)
}
min_1.2 <- nlminb(start = c(mean = 0, sd = 1),
  objective, data = data.frame(lower = 1.15, upper = 1.25),
  lower = c(mean = -Inf, sd = 0),
  control = list(trace = 1))
##   0:     3.9413405:  0.00000  1.00000
##   1:     3.5358156: 0.939180  1.34343
##   2:     3.0483807:  1.41780 0.810734
##   3:     2.9329842:  1.39233 0.722755
##   4:   0.014035849:  1.19062 0.0183273
##   5:    -0.0000000:  1.22768  0.00000
##   6:    -0.0000000:  1.22768  0.00000
print(min_1.2$par)
##     mean       sd 
## 1.227683 0.000000

We see that the optimization converges to a mean somewhere in the specified interval and a standard deviation of zero.

In order to better understand this result, we can visualize the likelihood profile:

mu_plot = seq(1, 1.4, by = 0.01)
sigma_plot = seq(0, 0.4, by = 0.01)
z <- outer(mu_plot, sigma_plot,
  function(mean, sd) prob_obs(1.15, 1.25, mean, sd))
persp(mu_plot, sigma_plot, z,
  theta = 35, phi = 20, ticktype = "detailed",
  xlab = "Mean", ylab = "Sigma", zlab = "Likelihood")

plot of chunk profile

Thus our model suggests that the most likely normal distribution given our single observation has a standard deviation close to zero and a mean somewhere in the interval between 1.15 and 1.25.

Several, potentially censored observations

In order to conveniently deal with several such observations, we define a function that converts certain character representations of our observations into a dataframe containing the lower and upper bounds of the corresponding observational intervals. Also, we want to be able to specify left and right censored observations using representations like '< 0.2' or '> 5.0'.

as.obs <- function(x) {
  if (any(sapply(x, function(a) grepl("e", a)))) stop ("Scientific notation is not supported")

  left <- grepl("^<", x)
  right <- grepl("^>", x)
  x_num <- as.numeric(gsub("^[<>]", "", x))
  decimals <- gsub(".*\\.", "", x)
  n_decimals <- nchar(decimals)
  interval_range <- 10^-n_decimals
  lower <- ifelse(left, -Inf, ifelse(right, x_num, x_num - 0.5 * interval_range))
  upper <- ifelse(right, Inf, ifelse(left, x_num, x_num + 0.5 * interval_range))
  ret <- data.frame(lower, upper)
  return(ret)
}
x <- c("1.2", "3.273", "< 0.2", "> 2.7")
x_obs <- as.obs(x)
print(x_obs)
##    lower  upper
## 1 1.1500 1.2500
## 2 3.2725 3.2735
## 3   -Inf 0.2000
## 4 2.7000    Inf

We can now estimate the parameters of the supposed underlying normal distribution and calculate the log-likelihood of the solution.

min_x_obs <- nlminb(start = c(mean = 0, sd = 1),
  objective, data = x_obs,
  lower = c(mean = -Inf, sd = 0),
  control = list(trace = 1))
##   0:     23.334776:  0.00000  1.00000
##   1:     16.546212: 0.351102  1.93634
##   2:     15.687217:  1.31191  1.73492
##   3:     15.636877:  2.14422  1.56714
##   4:     15.517024:  1.74168  1.70304
##   5:     15.494396:  1.77305  1.73891
##   6:     15.422966:  1.88078  1.96280
##   7:     15.412944:  1.90232  2.06494
##   8:     15.411206:  1.90644  2.12210
##   9:     15.411132:  1.90409  2.13473
##  10:     15.411128:  1.90225  2.13624
##  11:     15.411127:  1.90111  2.13621
##  12:     15.411127:  1.90094  2.13602
##  13:     15.411127:  1.90094  2.13598
print(min_x_obs$par)
##     mean       sd 
## 1.900940 2.135979
log_lik(min_x_obs$par["mean"], min_x_obs$par["sd"], x_obs)
## [1] -15.41113

In order to check what we have done, we can do the same using the survreg function in the survival package. At this point I would like to thank Torsten Hothorn for pointing out that this is possible and how it is done.

library(survival)
x_Surv <- Surv(x_obs$lower, x_obs$upper, type = "interval2")
x_survreg <- survreg(x ~ 1, data = data.frame(x = x_Surv), dist = "gaussian")
summary(x_survreg)
## 
## Call:
## survreg(formula = x ~ 1, data = data.frame(x = x_Surv), dist = "gaussian")
##             Value Std. Error    z     p
## (Intercept) 1.901      1.144 1.66 0.097
## Log(scale)  0.759      0.567 1.34 0.181
## 
## Scale= 2.14 
## 
## Gaussian distribution
## Loglik(model)= -15.4   Loglik(intercept only)= -15.4
## Number of Newton-Raphson Iterations: 4 
## n= 4

Voilá, we obtain the same results for the distribution parameters as well as for the log-likelihood!

Comparison with fitdistrplus

The package fitdistrplus also provides a way to fit distributions to censored data, namely via the fitdistcens function. However, the package vignette gives a likelihood function that mixes probabilities and probability densities. We can assess the effect of this in our example case. At first, we use the same interpretation of the observed data that we have used above.

library(fitdistrplus)
x_censdata <- data.frame(
  left = ifelse(x_obs$lower == -Inf, NA, x_obs$lower),
  right = ifelse(x_obs$upper == Inf, NA, x_obs$upper))

x_fitdistcens <- fitdistcens(x_censdata, "norm")
summary(x_fitdistcens)
## Fitting of the distribution ' norm ' By maximum likelihood on censored data 
## Parameters
##      estimate Std. Error
## mean 1.901758   1.144381
## sd   2.136302   1.211426
## Loglikelihood:  -15.41113   AIC:  34.82226   BIC:  33.59484 
## Correlation matrix:
##            mean         sd
## mean 1.00000000 0.01217819
## sd   0.01217819 1.00000000

As expected, we get pretty much the same results here as well. However, if we interpret the simple observations as exact (or non-censored) observations, we obtain a slightly different picture.

x_partcens <- data.frame(
  left = c(1.2, 3.273, NA, 2.7),
  right = c(1.2, 3.273, 0.2, NA))
x_fitpartcens <- fitdistcens(x_partcens, "norm")
summary(x_fitpartcens)
## Fitting of the distribution ' norm ' By maximum likelihood on censored data 
## Parameters
##      estimate Std. Error
## mean 1.901134   1.144231
## sd   2.136102   1.211229
## Loglikelihood:  -6.200705   AIC:  16.40141   BIC:  15.174 
## Correlation matrix:
##            mean         sd
## mean 1.00000000 0.01162075
## sd   0.01162075 1.00000000

While the parameter estimates are, as expected, approximately the same, we do obtain a different value for the log-likelihood. Ronald A. Fisher, when he introduced the term likelihood, stated that there is no absolute measure of likelihood (Fisher and Russell, 1922). In fact, the term was introduced in order to avoid the necessity to be able to interpret it as a probability or "inverse probability". However, in view of potential further interpretations of the likelihood it seems prudent to use a consistent definition for the likelihood of each observation in a set of observations, which must be based on an integral of the underlying probability density in the presence of censored data.


Bibliography

R. A. Fisher and Edward John Russell. On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 222(594-604):309–368, 1922. URL: https://royalsocietypublishing.org/doi/abs/10.1098/rsta.1922.0009, doi:10.1098/rsta.1922.0009.

Torsten Hothorn. Most likely transformations: the mlt package. Journal of Statistical Software, 2018. Accepted for publication 2018-03-05. URL: https://cran.r-project.org/package=mlt.

J. K. Lindsey. Some statistical heresies. The Statistician, 48:1–40, 1999.