Modelling Inverse Gaussian Data with Censored Response Values: EM versus MCMC

Abstract

Low detection limits are common in measure environmental variables. Building models using data containing low or high detection limits without adjusting for the censoring produces biased models. This paper offers approaches to estimate an inverse Gaussian distribution when some of the data used are censored because of low or high detection limits. Adjustments for the censoring can be made if there is between 2% and 20% censoring using either the EM algorithm or MCMC. This paper compares these approaches.

1. Introduction

Partial missing data is common in environmental applications, where measurement systems have low or high detection limits. The measurement processes at times are incapable of measuring below certain values; for example, near zero counts of colony form units, below a certain lower limit, denoted 𝑐, are considered too unreliable when the measurement process involves diluting solution samples. An example of a high detection limit (𝑐) occurs for flows when flow is measured based on the height of a river, yet the river height measurement cannot exceed the height of the bank. All that is known in such cases is that the measurement is below or above 𝑐, respectively. In this paper, accurate model fitting is the aim. We advocate the EM algorithm [1] to establish realistic penalised maximum likelihood estimates of regression parameters in such cases. Other approaches apply MCMC methods [2]. The next section introduces the EM algorithm approach for dealing with low/high detection limits. A simulation study for high detects (the case that is more unstable because it is unbounded on the high side) is used to demonstrate the capability of the methodology. Section 3 looks at the MCMC approach. Section 4 uses an example to demonstrate the EM methodology. Section 5 concludes with a discussion.

2. The EM Algorithm

The expectation step of the EM algorithm is used to replace censored values based on the current parameter estimates and the detection limit. These replaced values are, in turn, used to establish penalised maximum likelihood estimates of regression parameters. These two steps are applied iteratively until estimates converge.

In this paper, we will use the GAMLSS model to fit the location and spread for the data using penalised maximum likelihood estimation (see [3, 4]). The library GAMLSS has the capability of fitting models with low detects, but this is experimental, while our approach has been tested at least for the applications considered in this paper.

Let the location of the response 𝑦_𝑖 given the vector of explanatory variables 𝑥_𝑖 be denoted by 𝑓(𝜇_{𝑦⋅𝑥𝑖})=𝑥^𝑡_𝑖𝛽, where 𝑓 is the link function and 𝜇_{𝑦⋅𝑥𝑖}=E(𝑦_𝑖∣𝑥_𝑖) is the mean. We are interested in the inverse Gausian distribution using the usual logarithm link function givingE𝑦_𝑖∣𝑥_𝑖=𝜇_{𝑦⋅𝑥𝑖}=𝑒^{𝑥𝑡𝑖𝛽}. (2.1) The inverse Gaussian distribution is given by ℎ𝑦∣𝜇_{𝑦⋅𝑥𝑖},𝜙=𝜙𝑦2𝜋^−3/2𝑒^{−𝜙(𝑦−𝜇𝑖𝑦⋅𝑥)2/(2𝜇2𝑖𝑦⋅𝑥𝑦)}, (2.2) where 𝜙 is the shape parameter and the variance 𝜎²_𝑖 is given by 𝜇³_{𝑦⋅𝑥𝑖}/𝜙. The GAMLSS model for 𝜙 is 𝜙=𝑒^{𝑥𝑡𝑖𝛼}, (2.3) so that 𝜎²_𝑖=𝑒^{3𝑥𝑡𝑖𝛽−𝑥𝑡𝑖𝛼}. To simplify the notation, we drop subscripts; that is, let 𝜇=𝑒^𝑥𝑡𝛽. When 𝑦 is a low detect, we want to evaluate E(𝑦∣𝑥,𝑦<𝑐) which is given by (see Appendix A) E𝑦∣𝑥,𝑦<𝑐=𝜇−2𝜇𝑒^2𝜙/𝜇Φ−𝑐/𝜇+1𝜙/𝑐𝐻𝑐,𝜇,𝜙, (2.4) where 𝐻(𝑐,𝜇,𝜙)=∫^𝑐₀ℎ(𝑦,𝜇,𝜙)𝑑𝑦 and Φ(𝑧)=∫^𝑧_−∞𝑒^{(−1/2)𝑥2}/√2𝜋𝑑𝑥.

High detect limits are less common in the application considered later but do occur in practice. In this situation, it can be shown that the expectation step is given by (see [5] for a different result that is often numerically equivalent to the result below)𝐸𝑦∣𝑥,𝑦>𝑐=𝜇+2𝜇𝑒^2𝜙/𝜇Φ−𝑐/𝜇+1𝜙/𝑐1−𝐻𝑐,𝜇,𝜙. (2.5)

A simulation study, of 1000 runs per table entry, generated data from the inverse Gaussian distribution with density 𝑓(𝑥,𝜇=2,𝜙=1) and truncated values at high detection limits 𝑐=3, 4, 5, and 6. The application of the EM algorithm is outlined in Whitemore [5]. This was used to find the maximum likelihood estimates of hidden values 𝜇 and 𝜙. The results are reported in Table 1. Similar results were found for inverse Gaussian distributions with low detects.

Table 1 

EM simulation results for fitting Inverse Gaussian models with high detects.

It is clear from the standard errors in Table 1 that the EM algorithm is unstable for 50 observations with standard errors more than 25% of the parameter estimated values. However, for this few observations, we are unlikely to diagnose the data as Inverse Gaussian, and any methodology is going to be suspicious.

If the presence of high detects is ignored so that data with actual values greater than 𝑐 are set equal to 𝑐, then the resulting parameter estimates are biased. For 𝑐=3, 4, 5, and 6, the 1000 observations example with high-truncated values in Table 1, produced average estimates of 𝜇 and 𝜙, without correction for the truncation, equal to (𝜇=) 1.374, 1.528, 1.639, and 1.716 (actual 𝜇=2) and (𝜙=) 1.257, 1.169, 1.119, and 1.092 (actual 𝜙=1), respectively. These parameter estimates are significantly biased. A simulation exercise for low detects was carried out but is not presented, because the results are less extreme (e.g., if 𝑐≪𝜇, then little is lost by replacing the low detects by 𝑐/2). This will be demonstrated in the application in Section 4.

3. Markov Chain Monte Carlo

In this section, the results of a Markov chain Monte Carlo analysis of the simulation study of high detects are presented (see previous section for details). The mean and dispersion parameters are estimated using Gibbs sampling (e.g., [6]). Defining the mean parameter as 𝜇=𝑒^𝛽, the MCMC analysis is undertaken with a uniform prior distribution for 𝛽 and a gamma prior distribution for 𝜙, with shape parameter equal to 0.001 and scale parameter equal to 1,000 (i.e., prior variance = 1000). That is, the mean of the prior distribution for 𝜙 is equal to 1. Based on these prior distributions, the conditional distribution for 𝜙 is Gamma (𝑛/2+0.001, ∑^𝑛_𝑖=1(𝑦_𝑖−𝜇)²/(2𝜇²𝑦_𝑖)+0.001), so samples of 𝜙 can be drawn from this distribution, while samples for 𝛽, and hence 𝜇, are drawn using rejection sampling (e.g., [7, Chapter  3]).

The results of the MCMC analysis of the simulations are presented in the Table 2 below. These results are generated by using MCMC to draw 7,000 parameter samples for each of the 1,000 simulated data sets of the 16 different combinations of cutoff and number of observations in Table 2. That is, a total of 16,000 MCMC analyses are run and each of these is run for 7,000 sample iterations. The first 2,000 samples of each MCMC run are discarded, and of the remaining 5,000 samples, every fifth sample is saved and used to estimate 𝜙 and 𝜇. The medians of these 1,000 realizations is taken as the estimate of 𝜙 and 𝜇 for that simulated data set. in Table 2 presents summary statistics for the 1,000 median estimates for 𝜇 and 𝜙 from the 1,000 simulated data sets for each cutoff and number of observations combinations. An examination of the results in this table and those for the EM analysis in Table 1 reveals that the MCMC estimates are generally less biased for the dispersion parameter 𝜙 than the corresponding EM estimates, especially for the smaller sample sizes. For the mean parameter 𝜇 the means of the 1,000 median estimates for 𝜇 are generally more biased than the EM estimates.

Table 2 

MCMC simulation results for fitting inverse Gaussian models with high detects.

4. Application

The data include analytes (e.g., total nitrogen, chlorophyll-a, e. coli, turbidity, etc.) and flows routinely measured at 35 different sites in the waterways that supply Sydney with drinking water for a period going back 20 years. These measures are used to assess the quality of water in the various reservoirs, rivers, and dams that are used to supply Sydney and greater Sydney with drinking water. Some measurements are taken daily (temperature, dissolved oxygen, turbidity, and pH), while others are collected roughly either fortnightly or monthly. The fitted models were used to assess the risk that particular analytes pose to the drinking water supply to Sydney residents.

We used AIC to select the distribution and model which was most appropriate for several water quality variables. The distributions considered were normal, log-normal, inverse gaussian, and zero-adjusted inverse Gaussian. The analytes that appeared to select the inverse Gaussian as the preferred distribution in more than 75% of the sites were firstly phosphorus filterable, nitorogen ammoniacal, and e. coli. Other analytes that where selected less frequently as inverse Gaussian distributed were entrococci, toxic cyanobacterial count, chloride, nitrogen oxide, and phosphorus total. The number of measurements ranged from several thousand to a few hundred with different number of low detects. All of these were successfully fitted by the approach advocated in the paper. As an example, we model phosphorus total(mg/L) at Bendeela Pondage in the Sydney water catchment supply area. This example is illustrated below.

Table 3 presents the result of the fitted model for phosphorus total (mg/L) at Bendeela Pondage. This was measured for the past 15 years with frequency of just over one per month (giving 213 measurements). There are seven low detects with one being below 0.001 and six being below 0.01. Not adjusting for this can adversely influence the model fit, because 0.01 is close to the mean value. Applying the above EM algorithm to this example results in a nonhomogeneous mean and variance, which was fitted in R [8] by GAMLSS with family equal to the inverse Gaussian. The EM algorithm ran in a fraction of the time of the MCMC approach.

Table 3 

Fitted models replacing low detects with min(𝑦)/2, 𝑐/2, 𝑐 or the EM E-step.

In addition, qq-plots of the theoretical values versus the quantile residuals were closer to each other after convergence using EM algorithm. In Table 3, setting the low detects equal to 𝑐/2 produces a fitted model that is closest to the EM fit. However, even with the propotionally few low detects in the sample, there are significant differences in some of the time-of-the-day harmonic coefficients for the trend and spread (i.e., sin/cos(2𝜋𝑡) coefficients).

5. Discussion

The paper offers a way of improving models with the modelled response having an inverse Gaussian distribution when the response variable has low or high detects. Although only one example is reported in the paper, the methodology has been successfully applied to many more examples for the Sydney Catchment Authority. Problems with convergence were encountered in a few applications. The MCMC approach was computationally too intensive to run on all applications, but little appears to be lost by applying the EM algorithm. We used the AIC criterion to select the family from normal, log-normal, inverse Gaussian, or zero-adjusted inverse Gaussian (for which low detects were assumed zeros). This latter model is only realistic when the low detects are below all other measured values, hence, it is inappropriate in the first example in Section 4, where some low detects (i.e., 0.01) are close to the average response value. The high detection limit is seldom an issue in applications, but its theory was included, because it was encountered on a few occasions in the Sydney Catchment Authority work. The EM algorithm in this paper makes no attempt to achieve a maximum likelihood estimate for the variance; however, this limitation is of little consequence, since the variance estimate is reasonably comparable to the MCMC approach, which does properly account for the variance.

Numerical instabilities in the approaches were experienced in simulations, but results can be achieved by small adjustments to the estimated variances.

Appendix A

Our aim in this appendix is to find 𝐄(𝐱∣𝐱<𝐜) and 𝐄(𝐱∣𝐱>𝐜) using the moment generating function (mgf) of 𝐱 given 𝐱<𝐜, denoted 𝐌_{𝐱∣𝐱<𝐜}(𝐭) or 𝐌_{𝐱∣𝐱>𝐜}(𝐭), respectively; for example, 𝐄(𝐱∣𝐱<𝐜)=(𝜕𝐭𝐌_{𝐱∣𝐱<𝐜}(𝐭))/(𝜕𝐭)|_𝐭=𝟎.

A. Inverse Gaussian Distribution

The inverse Gaussian distribution for random variable 𝑥 is given by density ℎ𝑥,𝜇,𝜙=𝜙𝑥2𝜋^−3/2𝑒^{−𝜙(𝑥−𝜇)2/(2𝜇2𝑥)} (A.1) and has distribution function 𝐻(𝑐,𝜇,𝜙)=∫^𝑐₀ℎ(𝑥,𝜇,𝜙)𝑑𝑥=Φ(((𝑐/𝜇)−1)𝜙/𝑐)+𝑒^2𝜙/𝜇Φ(−((𝑐/𝜇)+1)𝜙/𝑐), where Φ(𝑐)=∫^𝑐_−∞(1/2𝜋)𝑒^−𝑧2/2𝑑𝑧.

The mgf for the Inverse Gaussian distribution is given by 𝑀_𝑥(𝑡)=^∞₀𝑒^𝑥𝑡ℎ𝑥,𝜇,𝜙𝑑𝑥=𝑒^{(𝜙/𝜇)(1−1−2𝜇2𝑡/𝜙)}. (A.2) By completing the square for the exponent terms in (A.1) and setting 𝜇_∗=𝜇/1−2𝜇²𝑡/𝜙, we have−𝜙𝑥−𝜇²2𝜇²𝑥+𝑥𝑡=−𝜙𝑥−𝜇_∗²2𝜇²_∗𝑥+𝜙𝜇1−1−2𝜇²𝑡𝜙=−𝜙𝑥−𝜇_∗²2𝜇²_∗𝑥+𝜙𝜇−𝜙𝜇_∗, (A.3) so that for low detects, where we know that 𝑥<𝑐, the mgf becomes 𝑀_{𝑥∣𝑥<𝑐}(𝑡)=^𝑐₀𝑒^𝑥𝑡ℎ𝑥,𝜇,𝜙𝑑𝑥𝐻𝑐,𝜇,𝜙=𝑒^{𝜙/𝜇−𝜙/𝜇∗𝑐}₀ℎ𝑥,𝜇_∗,𝜙𝑑𝑥𝐻𝑐,𝜇,𝜙=𝑒^{𝜙/𝜇−𝜙/𝜇∗}𝐻𝑥,𝜇_∗,𝜙𝐻.𝑐,𝜇,𝜙 (A.4) Since 𝜇_∗=𝜇 when 𝑡=0, we have𝐸(𝑥∣𝑥<𝑐)=𝜕𝜇_∗/𝜕𝑡𝜕𝑒^{𝜙/𝜇−𝜙/𝜇∗}𝐻𝑥,𝜇_∗,𝜙/𝜕𝜇_∗𝐻𝑐,𝜇,𝜙_𝑡=0=𝜇³_∗𝜙𝑒^{𝜙/𝜇−𝜙/𝜇∗}𝜙𝐻𝑥,𝜇_∗,𝜙/𝜇²_∗+𝜕𝐻𝑥,𝜇_∗,𝜙/𝜕𝜇_∗𝐻𝑐,𝜇,𝜙_𝑡=0.𝜕𝐻𝑐,𝜇,𝜙=𝜕Φ𝜕𝜇𝑐/𝜇−1𝜙/𝑐+𝑒^2𝜙/𝜇Φ−𝑐/𝜇+1𝜙/𝜇𝜕𝜇=−2𝜇𝑒^2𝜙/𝜇Φ−𝑐𝜇+1𝜙𝑐. (A.5) Therefore𝐸(𝑥∣𝑥<𝑐)=𝜇−2𝜇𝑒^2𝜙/𝜇Φ−𝑐/𝜇+1𝜙/𝑐𝐻𝑐,𝜇,𝜙. (A.6) Similarly, for high detects with 𝑥>𝑐,𝐸(𝑥∣𝑥>𝑐)=𝜇+2𝜇𝑒^2𝜙/𝜇Φ−𝑐/𝜇+1𝜙/𝑐1−𝐻𝑐,𝜇,𝜙. (A.7) See also Whitemore [5] which seems to often agree numerically with the result.