Bayesian Nonparametric Model for Weighted Data Using Mixture of Burr XII Distributions

TL;DR

This paper extends Bayesian nonparametric modeling to weighted data using Dirichlet process mixtures with Burr XII distributions, demonstrating its effectiveness on lifetime data with censored observations.

Contribution

It introduces a novel application of DPMMs with Burr XII kernels to weighted data and develops an MCMC-based method to estimate unweighted distributions from weighted data.

Findings

Effective modeling of lifetime data with censored observations

Accurate estimation of density and survival functions

Successful application to real and simulated datasets

Abstract

Dirichlet process mixture model (DPMM) is a popular Bayesian nonparametric model. In this paper, we apply this model to weighted data and then estimate the un-weighted distribution from the corresponding weighted distribution using the metropolis-Hastings algorithm. We then apply the DPMM with different kernels to simulated and real data sets. In particular, we work with lifetime data in the presence of censored data and then calculate estimated density and survival values.

Equations

$$g(x)=\frac{w(x)f(x)}{E[w(X)]},E[w(X)]<\infty ,x\ge 0,$$

$${\int }_0^{\infty }w(x{)}^{-1}g(x)dx<\infty ,$$

$$Algorithm1:Metropolis-Hastingsalgorithm$$

$$Initializewith{x}^{(0)}\sim q(x)$$

$$\text{for}i=1,2,...\text{do}$$

$$Propose{x}^{cand}\sim q(x^{(i)}|x^{(i-1)})$$

$$Calculatetheacceptanceprobability:$$

$$\alpha (x^{cand}|x^{(i-1)})=min\{1,\frac{q(x^{(i-1)}|x^{cand})\pi (x^{cand})}{q(x^{cand}|x^{(i-1)})\pi (x^{(i-1)})}\}$$

$$Generateu\sim Uniform(u;0,1)$$

$$\text{if}u<\alpha \text{then}$$

$$x^{(i)}\leftarrow x^{cand}$$

$$x^{(i)}\leftarrow x^{(i-1)}$$

$$x_{j+1}=y_{j+1}withprobabilitymin\{1,\frac{w^{-1}(y_{j+1})}{w^{-1}(x_j)}\},$$

$$x_{j+1}=x_{j}otherwise.$$

$$p(x_{j+1}|x_j)=min\{1,\frac{w^{-1}(y_{j+1})}{w^{-1}(x_j)}\}g(x_{j+1})+\{1-r(x_j)\}1(x_{j+1}=x_j),$$

$$r(x)=\int min\{1,\frac{w^{-1}(x^{\ast })}{w^{-1}(x)}\}g(x^{\ast })d{x}^{\ast }.$$

$$g_{\mathcal{P}}(y)=\int \kappa (y;\theta )P(d\theta ),$$

$${\int }_0^{\infty }{w}^{-1}(y)\kappa (y;\theta )dy<\infty .$$

$$g_{c,k,P}(y)={\int }_{\mathcal{R}}Bur{r}_{XII}(y|c,k)P(dc,dk),$$

$$y|c,k$$

$$(c,k)|P$$

$$P|\upsilon ,P_0$$

$$P_0(c,k|\varphi ,\gamma )=Unif(c|0,\varphi )\times Exp(k|\gamma ).$$

$$t_i|c_i,k_i$$

$$(c_i,k_i)|\mathcal{P}$$

$$\mathcal{P}$$

$$P_0|\gamma ,\varphi$$

$$\nu ,\gamma ,\varphi$$

$$(1)[({\theta }_i)|({\theta }_{-i},z_{-i}),\upsilon ,...,t],fori=1,...,n$$

$$(2)[({\theta }_j^{\ast })|z,n^{\ast },\upsilon ,...,t],forj=1,...,n^{\ast }$$

$$(3)[\upsilon |\{({\theta }_j^{\ast }),j=1,...,n^{\ast }\},n^{\ast },t],[...|\{({\theta }_j^{\ast }),j=1,...,n^{\ast }\},n^{\ast },t].$$

$$Algorithm2:Gibbssampler$$

$$1.Initializewith{\theta }^{(0)}\sim f(\theta )$$

$$2.Fori=1,2,...do$$

$${\theta }_1^{(i)}\sim f({\theta }_1|{\theta }_2^{(i-1)},{\theta }_3^{(i-1)},...,{\theta }_d^{(i-1)},D),$$

$$⋮$$

$${\theta }_d^{(i)}\sim f({\theta }_d|{\theta }_1^{(i)},{\theta }_2^{(i)},...,{\theta }_{d-1}^{(i)},D),$$

$$f(c_i,k_i|\{(c_{i^{\prime }},k_{i^{\prime }});i\ne i^{\prime }\},\nu ,\gamma ,\varphi ,t_{io})=\frac{q_0^o{h}^o(c_i,k_i|\varphi ,\gamma ,t_{io})+{\sum }_{j=1}^{n^{\ast (i)}}{n}_j^{\ast (i)}{q}_j^o{\delta }_{c_j^{\ast },k_j^{\ast }}}{q_0^o+{\sum }_{j=1}^{n^{\ast (i)}}{n}_j^{\ast (i)}{q}_j^o},$$

$$q_0^o$$

$$h^o(c_i,k_i|\gamma ,\varphi ,t_{io})\propto k_B(t_{io}|c_i,k_i)P_0(c_i,k_i|\gamma ,\varphi )\propto [c_i|\gamma ,\varphi ,t_{io}][k_i|c_i,\gamma ,\varphi ,t_{io}]$$

$$[c_i|\gamma ,\varphi ,t_{io}]\propto c_i{t}_{io}^{c_i-1}{I}_{(0,\varphi )}(c_i)i=1,...,n$$

$$[k_i|c_i,\gamma ,\varphi ,t_{io}]\propto Gamma(.|2,\frac{1}{[\frac{1}{\gamma }+ln(1+t_{io}^{c_i})]}).$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Let $X$ be a non-negative random variable with density function $f(x)$ and $w(x)$ be a non-negative function of $x$. A new random variable $X^{w}$ with density function $g(x)$ that is defined as bellow,

[TABLE]

is called weighted random variable with respect to X, and $g(x)$ is called the weighted density function with respect to $f(x)$. Especially, if $w(x)=x$, the resulted weighted distribution is called length-biased distribution which has various applications in many areas. Zelen and Feinleib(1969) applied the length-biased distribution to detect breast cancer early, and also Patil et al. (1977) used this weighted distribution to study human families and wild-life population. Patil et al.(1978) introduced distributions of the type given in equation 1 with an arbitrary non-negative weight function $w(x)$ and gave practical examples. For more examples of weighted distributions and their applications, see [2],[21],[22],[16],[9].

In this work we consider wighted distributions nonparametric Bayesian methods by using Dirichlet process mixture model(DPMM) which is a papular Bayesian nonparametric model and then, appply these idea to survival analysis.

In [10], DPMM is used for density estimation under length-biased data. They consider a log-normal distribution as the kernel function which its shape parameter has the determined distribution. Here we consider a DPMM with Burr(XII) distribution as the kernel function with two parameters (scale parameters) whose distributions are random in the model. Since the support of the Burr(XII) distribution is $R^{+}$, it is suitable for survival study ([1],[13],[24],[25]). In [1] Burr(XII) distribution is used as the kernel in DPMM, and then, survival function and hazard rate are calculated for simulated and real data.

In the next section, preliminary concepts and the methodology are presented, and section 3 describes the nonparametric Bayesian approach that allows us to define the model which will be used in the next section. Section 4 is contained the modeling and algorithm that is used for the sampling. In the next section, we consider data illustrations for applying the model to the data. Finally, we summarized our results in the conclusion section.

2 Preliminary and Methodology

We want to estimate the density and other survival functions by considering a general case of weight function $w(x)$. The strategy for avoiding computing the intractable normalizing constant would be to model $g(x)$ directly and then make inference about $f(x)$ by considering this fact that $g(x)\propto w(x)f(x)$. If we set $f(x;\theta)$ as a parametric family, so $f(x;\theta)$ and $g(x;\theta)$ are known except the normalizing constant that may be not tractable.

Let $w(.)$ be a general weight function, an essential condition to model F(.) through G(.) (F(.) and G(.) denote the distribution functions of f(.) and g(.) respectively.) is

[TABLE]

because f is a distribution function.

Through invertibility, equation (2) enables us to reconstruct F from G.

In the Bayesian nonparametric framework, we assign an appropriate nonparametric prior distribution on g, providing relation (2). The question that now arises is how the posterior structures obtained after modeling g directly can convert to the posterior structures from f.

The first step is to construct a method to convert a weighted sample to an un-weighted one. Then it is possible to inference about the posteriors.

An indirect method to simulate samples of complex distribution is Monte Carlo Markov chain (MCMC) approach. Metropolise-Hastings algorithm [metro95] is one of the MCMC methods which simulates samples from a probability distribution by making use of the full joint density function and proposal distributions for each of the variables of interest.

In general form, the Metropolis-Hastings algorithm is as the following form,

[TABLE]

Hatjispyros and et al. used the Metropolis-Hastings algorithm to convert a length biased sample to the unbiased version. Now, we want to apply this algorithm by using general weight function satisfied in relation (2) to convert the sample from weighted distribution to the un-weighted version.

Suppose that $y_{1},y_{2},...,y_{N}$ denote a random sample of $g$. The Metropolis-Hastings algorithm converts this sample to a sample from $f(x)\propto w(x)^{-1}g(x)$. We assume $g(.)$ is replaced by $q(.)$ in the algorithm 1 with acceptance probability $min\{1,\frac{w^{-1}(y_{j+1})}{w^{-1}(x_{j})}\}$. If $x_{j}$ denotes the current sample from f(x), then

[TABLE]

The transition density is,

[TABLE]

where

[TABLE]

We can have the following outline methodology in a general form:

1- $(y_{1},...,y_{n})$ is a sample from g that we are going to assign a suitable nonparametric prior to it.

2- Using MCMC methods, the posterior values of the random measure $\Pi(dg|y_{1},...,y_{n})$ and other relevant parameters will obtain. So a sequence $\{y^{l}_{n+1}\},l=1,2,...$, from the posterior predictive density $g(y|y_{1},...,y_{n})$ will be generated.

3- $\{y^{l}_{n+1}\}$ will form a sequence of proposal values of a Metropolis-Hastings chain with the stationary density of the weighted posterior predictive i.e. $\{y^{l}_{n+1}\}\propto w(y)^{-1}g(y|y_{1},...,y_{n})$. By the equation (3) we generate the $\{x^{l}_{n+1}\}$ values at the level l.

4- Then $\{x^{l}_{n+1}\}$ values are a sample from the posterior of predictive f (un-weighted density).

3 The model and inference

Modeling $g(x)$, the weighted distribution, in the Bayesian nonparametric framework is based on infinite mixture model [15] as the following form

[TABLE]

where P is a discrete probability measure and $\kappa(y;\theta)$ is a kernel density on $(0,\infty)$ for all $\theta$’s in the parameter space, satisfying the following condition,

[TABLE]

By choosing Burr(XII) density with two parameters c and k as the kernel of the mixture model, so we have

[TABLE]

where $P$ is a discrete random probability measure. Suppose $P\sim DP(\upsilon,P_{0})$ where $DP(\upsilon,P_{0})$ denotes the Dirichlet process with precision parameter $\upsilon>0$ and base measure $P_{0}$ [7]. We named this mixture model by Dirichlet process Burr(XII) mixture model (DPBMM).

Hierarchical representation of DPBMM can be presented as the following form,

[TABLE]

Suppose that the base distribution $P_{0}$ is the prior distribution for the joint distribution of $c$ and $k$. By choosing Burr(XII) distribution as the kernel, $P_{0}$ that yields closed-form expression for $\int k_{B}(.|c,k)P_{0}(dc,dk)$ is not available. Moreover, we choose multiple distributions of uniform$(0,\phi)$ and exponential with the parameter $\gamma$ for the $P_{0}$, i.e.

[TABLE]

This choice achieves determined goals. By considering hyper-parameters, $\gamma$ and $\phi$ are random, we choose prior distributions $Pareto(a_{\phi},b_{\phi})$ and $IGamma(a_{\gamma},b_{\gamma})$ for them respectively.

Set the $a_{\phi}=a_{\gamma}=d$ and also d=2 since this value makes the variance of Pareto distribution infinite that cover all values in R. $b_{\phi}$ and $b_{\gamma}$ are determined by the data [1].

Finally, for any $t_{i},i=1,...,n$, lifetime data in a sample of n observations, by considering DPBMM and selecting priors for parameters of the model we have,

[TABLE]

After determining the model, we want to formulate how to sample from DPMMs by Gibbs sampling. According to [12], Gibbs sampling for drawing a sample from $[(\theta_{1},...,\theta_{n}),\upsilon,...|t]$ based on the following full conditional distributions ( bracket is used to show the conditional and marginal distributions):

[TABLE]

where t is the vector of failure time data. Here, $\theta_{i}$’s are parameters of the kernel in DPMMs that will be analyzed.

Model (4) and discreteness property of Dirichlet process, exhibit a clustering in $\theta$’s. We present $n^{*}$ as the number of the clusters between $\theta_{i}$’s that denote by $\theta_{j}^{*}$’s. The vector of indicators $z=(z_{1},...,z_{n})$ indicates the clustering configuration such that, $z_{i}=j$ when $\theta_{i}=\theta_{j}^{*}$. Also, the $\theta_{-i}$ that used in (8), will be defined by $\theta_{-i}=(\theta_{1},\theta_{2},...,\theta_{i-1},\theta_{i+1},...,\theta_{n})$.

4 Modeling

For modeling the un-weighted density $f(x)$ from the weighted density $g(x)$, we apply the following algorithm. At first, to generate a sample from $g(x)$, it needs to estimate the parameters of the model. To this aim we draw a sample from $(c_{i},k_{i})$ and update $z_{i}$ for each $t_{i}$.

In simulation-based parameter estimation, we use the Gibbs sampler that it includes two steps to reach the goal.

[TABLE]

where $\theta_{1},...,\theta_{d}$ are model parameters and D is the vector of observations. The values of iteration i would be sampled from the distribution with last version of the other parameter values.

Now, the model will be applied for lifetime data with the presence of right censored data, that is very common in the survival study. To calculate the related distributions we divide data into censored and uncensored observations.

1- Uncensored data:

For uncensored data($t_{io}$), the conditional posterior density of $(c_{i},k_{i})$ is a mixed distribution [20],

[TABLE]

where $q_{j}^{o}=k_{B}(t_{io}|c_{j}^{*},k_{j}^{*})$ and

[TABLE]

in which the last integration can be computed numerically and

[TABLE]

where

[TABLE]

and

[TABLE]

2- right censored data:

For right censored data $(t_{ic})$, the conditional posterior density of $(c_{i},k_{i})$ is

[TABLE]

where $q_{j}^{c}=1-K_{B}(t_{ic}|c_{j}^{*},k_{j}^{*})$, and

[TABLE]

Also, the last integration in above equation is computed numerically. By using the property of the censored data we have,

[TABLE]

To sample from the first part of the last equation, we use the slice sampling method. Therefore, by using this MCMC method, we can have a sample from $h^{c}(c_{i},k_{i}|\phi,\gamma,t_{ic})$. Now for observed and censored data, $(c_{i},k_{i})$ for $i=1,...,n$ can be updated and improved.

In a general form, $(c^{*}_{j},k_{j}^{*})^{\prime}s$ can be updated on $\phi,\gamma$ and $t$ as the following

[TABLE]

where $B^{*}=\sum_{\{io:s_{io}=j\}}(\frac{1}{\gamma}+ln(1+t_{io}^{c_{j}^{*}}))+\sum_{\{ic:s_{ic}=j\}}ln(1+t_{ic}^{c_{j}^{*}})$ and $n_{j}^{o}$ is the number of observed data which located in cluster j. The important task to generate a sample from equation (9), is drawing from the first part of the equation. Sampling from the gamma distribution is simple.

To sample from

[TABLE]

it needs to consider auxiliary variables W=\bigl{\{}w_{io};\{io:s_{io}=j\}\bigr{\}} such that

[TABLE]

By marginalization over the auxiliary variables, then $[c_{j}^{*}|\phi,t_{io}]$ for $j=1,...,n^{*}$ will be obtained. Moreover, $w_{io}$’s are uniform variables on $(0,\frac{t_{io}^{c_{j}^{*}-1}}{1+t_{io}^{c_{j}^{*}}})$. Therefore we have

[TABLE]

where $B=max\{0,\frac{ln(w_{io})}{1+t_{io}}\}$. Now drawing from $[c_{j}^{*}|\phi,t]$ is straightforward.

Afterward, using the method that is applied in [6], $\phi,\gamma$ and $\nu$ will be updated. If we take $u$ as a latent variable such that

[TABLE]

then

[TABLE]

where $p=\frac{a_{\nu}+n^{*}-1}{n(b_{\nu}-log(u))+a_{\nu}+n^{*}-1}$.

And finally to update $\phi$ we have

[TABLE]

where $b^{*}=max\{b_{\phi},max_{1\leq j\leq n^{*}}c_{j}^{*}\}$. So

[TABLE]

Repeating this technique can update $\gamma$

[TABLE]

Now, the all conditional distributions on the equation (8) will be computed.

5 Data illustrations

In this section, we consider two kinds of data set, simulated data and real data. For a given sample $(x_{1},...,x_{n})$, the density function is estimated and compared with the following two density estimators which are used in [10]:

i) The classical kernel density estimation,

[TABLE]

ii) The kernel density estimation for indirect data,

[TABLE]

where $\hat{\mu}$ is harmonic mean of $(x_{1},...,x_{n})$.

As we will see, these estimators have good fitness for these type of data and values of $\tilde{g}_{h}$ and $\hat{f}_{J,h}$ are close to the exact values of the true density. To simulate the required samples, the Gibbs sampler iterates 60,000 times with a burn-in period of 10,000 times.

5.1 Simulated data

5.1.1 Length biased distribution of log-normal

The first data set simulated from the log-normal distribution with parameter $(\mu,\sigma^{2})=(0.5,0.5)$. We use this fact that Length biased distribution of a log-normal with parameters $\mu+\sigma^{2}$ and $\sigma^{2}$ is again a log-normal with parameters $\mu$ and $\sigma^{2}$ [22].

By choosing the log-normal distribution as the kernel, we can show the preference of the model and algorithm. This model tested in [10] for length-biased data using simulated data from the gamma distribution with DPMM when the gamma distribution considered as the kernel.

Figure 1 shows this simulated data. In part (a) histogram of simulated data is drawn and also the curve of data is depicted with the solid line and $\tilde{g}_{h}$ is shown with the dashed line.

In part (b) we see the histogram of predictive distribution and also real curve and $\tilde{g}_{h}$ with solid and dashed line respectively.

Part (c) depict the histogram of convert data to the un-weighted version using the Metropolis-Hastings algorithm. Real data here is the log-normal distribution with parameters (0,0.5), that is drawn with the solid line. Kernel density estimator for indirect data is shown with the dash-dot line.

As we see in figure 1, the predictive distribution for simulated data have good fitness and its curve is very similar to the classical kernel density estimator . Distribution of un-weighted data is also close to the real distribution, log-normal distribution with parameters (0,0.5).

5.1.2 Weighted distribution of Gamma

Here we consider a $Gamma(\alpha,\beta)$ distribution with the weight function $w(x|a,b)=x^{a}exp(\frac{-x}{b})$. So by computing the weighted distribution, we see resulted distribution is again a gamma distribution with parameters $(\alpha+a,\frac{\beta+b}{b\beta})$.

This data set simulated from the $gamma(1,2)$ as the weighted distribution with $w(x)=exp(-x)$(a=0,b=1), and then its un-weighted version is a gamma distribution with (1,1) as the parameters.

In figure 2 part (a), we depict the histogram of the data with its real curve and $\tilde{g}_{h}$.

In part (b), the histogram of predictive values of data, $\tilde{g}_{h}$ and real curve are shown.

In part (c), we can see the histogram of un-weighted distribution that is obtained from the Metropolis-Hastings algorithm and $\hat{f}_{J,h}$ for this data.

5.2 Real data

5.2.1 Widths of shrubs data

For real data, we consider the data that can find in [19]. This data consists of 46 measurements of widths of shrubs that are sampled by line-transect. In this method of sampling, the probability of inclusion in the sample is proportional to the width of the shrub that it makes it a case of length-biased sampling.

The un-weighted version of values of data and $\hat{f}_{J,h}$ depicted in part (a) of figure 3, also we can see the predictive values of the DPBM model with histogram and $\tilde{g}_{h}$ with the dashed line in part (b).

5.2.2 Bladder cancer data

The next real data set is a survival data which included censored values. This data is taken from [14] which corresponding to remission times (in the month) of a random sample of 138 bladder cancer patients. For fitting the model, we divided data to 10. At first, we apply the DPBM model to the data and in figure 4 draw the estimated density function base on the model and histogram of data.

We consider data that comes from weighted distribution (as in [11] a weighted Lomax distribution is fitted to this data set.), then by considering $w(x)=e^{-x}$, the histogram of un-weighted values with the Metropolis-Hastings method and curve is obtained as the figure 5.

6 Conclusion

In this article, we apply the Bayesian nonparametric approach to model weighted data. We use the Dirichlet process mixture model (DPMM) with Burr(XII) distribution as the kernel function in mixing models. We assumed weighted distribution with arbitrary weight function that satisfies in equation(2). By using the Metropolis-Hastings algorithm, the weighted distribution converted to the un-weighted one. As an application, we fit the DPMM with the different kernels and weight functions for real and simulated data sets. As application in the survival study, a real lifetime data set which contained censored observations are used and density and survival functions are calculated.