Keywords

Stochastic process, Longitudinal, Integrated ornstein-uhlenbeck, Bayesian, Spatial, Convolution

Mathematics Subject Classification

62M05, 62M30, 62P10

Introduction

A longitudinal data consists of measurements of a single variable taken repeatedly over time from an individual. Any approach used to analyze such data must properly consider the correlations among the observations, see Li, N., et al. [1]. The typical structure of longitudinal data is numerous measurements of a possible multivariate response variable on each subject. There could also be covariates, possibly time varying, that influence the response variable [2]. The aim in the analysis of such data is to understand the changes in the mean structure of the response variable with time, to understand the effect of the covariates on the response variable, and to understand the within subject correlation structure.

Let $Y_i\left(t_{ij}\right)$ be the longitudinal process of subject $i$ at time $t_{ij}\ge \text{ 0}$. Values of $Y_i\left(t_{ij}\right)$ are measured intermittently at some set of times $t_{ij}$. The observed longitudinal data on subject $i$ may be subject to "error", thus we observed only $X_i\left(t_{ij}\right)$ whose elements may not exactly equal the corresponding $Y_i\left(t_{ij}\right)$. Given a specification for $Y_i\left(t_{ij}\right)$, the observed longitudinal data are taken to follow

$$Y_i\left(t_{ij}\right)\text{ = }{X}_i\left(t_{ij}\right)+{\epsilon }_i\left(t_{ij}\right)(1.1)$$

Where ${\epsilon }_i\left(t_{ij}\right)$ is an intra subject error.

In this paper, following Taylor, et al., Xu and Zeger, Wang and Taylor, and Abu Bakar, et al. [3-6] we consider a model of the form

$$Y_i\left(t_{ij}\right)\text{ = }{X}_i\left(t_{ij}\right)+{\epsilon }_i\left(t_{ij}\right)$$

$$X_i\left(t_{ij}\right)\text{ = }{a}_i+{\text{ bt}}_{ij}+\beta {\text{Z}}_i\left(t_{ij}\right)+{\text{ W}}_i\left(t_{ij}\right)(1.2)$$

With the notations (N) denotes for Normal Distribution, and $\left(N_{n_i}\right)$ denotes for Multivariate Normal with dimension $n_i$. $Y_i\left(t_{ij}\right)$ denotes the observed value of a continuous time-dependent covariates (disease marker or longitudinal measurements) for subject $i,\left(i=\text{ 1, }\mathrm{...}\text{, M}\right)$ at $t_{ij}\left(j=\text{ 1, }\mathrm{...}{\text{, n}}_i\right);\text{ M}$ number of the subjects in the study; $n_i$ number of repeated measurements for subject i and may be different for each subject, $X_i\left(t_{ij}\right)$ is the true value of the marker at time $t_{ij}$, ${\text{a}}_i$ are independent random intercept of subject $i$ follow the normal distribution with mean ${\mu }_a$ and variance ${\sigma }_a^2,\left\{a_i\sim N\left({\mu }_a,{\sigma }_a^2\right)\right\},\text{ b}$ is fixed slope, ${\text{Z}}_i\left(t_{ij}\right)$ is the covariates for subject $i$ at time $t_{ij}$ with corresponding unknown regression parameter $\beta ,{\epsilon }_i\left(t_{ij}\right)\sim N\left(0,{\sigma }_{\epsilon }^2\right)$ represents deviations due to measurement error and "local" biological variation that is on a sufficiently short time scale that ${\epsilon }_i\left(t_{ij}\right)$ may be taken independent across j, and $W_i\left(t_{ij}\right)\sim N_{n_i}\left(0,\wedge \right)$ is a vector of independent IOU stochastic process, the covariance matrix $\wedge$ with parameters ${\alpha }_w$ and ${\sigma }_w^2$; only depends on $i$ through the number $n_i$ of observations and through the time points $t_{ij}$ at which measurements are taken, see Henderson, Diggle and Dobson [7].

In this paper, we introduce an approximation of the IOU process which still gives an efficient inference of model parameters and reducing the dimensionality and complexity of the model. The details of this simple approach are given in the following sections, including an approximate formulation, the likelihood and parameter estimations. The usefulness of this modeling approach is then demonstrated by simulations.

The Stochastic Process and Approximation

The main disadvantage of the IOU process is that it is not stationary; hence it is necessary to have a natural time zero for each individual. In some applications, it may be that there is no natural time zero, or that time zero is not exactly known. Large Longitudinal datasets are often defined on naturally heterogeneous fields or have other inherently spatially varying conditions. Therefore, it is unreasonable to expect a response variable to be well-modeled by a stationary process over a large domain space. However, using non-stationary models is difficult in practice due to the conceptual challenges in specifying the model and the computational challenges of fitting the model when the data is so large that memory constraint prevent formation of the covariance matrix.

We propose to approximate the IOU stochastic process $W_i$ into flexible spatial model that can be constructed by convolving a very simple and independent, process with a kernel function. This approximation for constructing a spatial process introduces a number of advantages over specification through a spatial covariogram. In particular, this process convolution specification leads to computational simplifications and easily extends beyond simple stationary models. Our modeling approach is similar to that in Higdon [8], provide simple representations of such model by convolving continuous white noise with a kernel, whose shape determines the covariance structure of the resulting process. This approach is an alternative to traditional geostatistical techniques, where a covariance function is specified directly, but allows for increased flexibility, since the choice of the kernel also allows for features such as non-stationary, anisotropy, and edge effects. Moreover, model (1.2) is temporal longitudinal model, by applying the proposed approximation, the dimensionality of the complex temporal process significantly reduced.

The model we propose for $X_i\left(t_{ij}\right)$ at time $t_{ij}$ is

$$X_i\left(t_{ij}\right)\text{ = }{a}_i+b{t}_{ij}+\beta Z_i\left(t_{ij}\right)+U_i\left(t_{ij}\right),$$

Where U is a mean zero Gaussian process, and is an approximation of the IOU stochastic process W. Hence, model in (1.2) becomes

$$Y_i\left(t_{ij}\right)\text{ = }{X}_i\left(t_{ij}\right)+{\epsilon }_i\left(t_{ij}\right)$$

$$X_i\left(t_{ij}\right)\text{ = }{a}_i+b{t}_{ij}+\beta Z_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)(2.1)$$

But rather than specify $U\left(t\right)$ through its covariance function, it is determined by the latent process $U^{\ast }\left(t\right)$ and the smoothing kernel $K\left(t\right)$.The restriction on the latent process $U^{\ast }\left(t\right)$ to be nonzero at the spatial sites $t^{\ast }=\left\{t_1^{\ast }=t_{i1},t_2^{\ast }=t_{i2},\mathrm{...}\text{,}{t}_{m_i}^{\ast }=t_{i{m}_i}\right\}$ in a space time t and define $U^{\ast }\left(t^{\ast }\right)=\left\{u^{\ast }\left(t_1^{\ast }\right),\mathrm{...}\text{, }{u}^{\ast }\left(t_{m_i}^{\ast }\right)\right\}$. Each $u^{\ast }$ is then modeled as independent draws from a mean zero Gaussian distribution with variance ${\sigma }_u^2$.

Hence, $U^{\ast }$ will follow a multivariate normal with mean zero and variance covariance matrix ${\displaystyle \sum =}{\text{ I}}_{m_i}{\sigma }_u^2$, here $m_i$ is the number of space times over spatial sites $t^{\ast }$.

The new Gaussian process $U_i\left(t\right)$ is then

$U_i\left(t_{ij}\right)\text{ = }{\displaystyle {\sum }_{r=1}^{m_i}K\left(t_{ij}-t_r^{\ast }\right)}{U}_i^{\ast }\left(t_r^{\ast }\right)(2.2)$

Where $K\left(.-t_r^{\ast }\right)$ is a kernel centered at $t_r^{\ast }$.

The resulting covariance function for $U\left(t\right)$ depends only on the displacement vector d = t − s and is given by

$$C\left(d\right)=\text{ Cov}\left(U\left(t\right),U\left(s\right)\right)={\displaystyle {\int }_{s}K\left(u-t\right)}K\left(u-s\right)du(2.3)$$

Table 1 shows kernels that give standard Gaussian, exponential, and spherical covariograms for the process $U\left(t\right)$, [9]. In addition, the covariogram induced by the biwieght kernel Cleveland [10] is also shown.

The process convolution approach gives an approach to build dependent spatial processes, see Ver Hoef and Barry [11]. The basic idea is to build processes $U\left(t\right)$ that share part of a common latent process in their construction. Perhaps the biggest attraction to these process convolution models is that they give a framework for developing new classes of space and space-time models that allow for more realistic space-time dependence while maintaining some analytic tractability. Generally, one can construct a space-time process by first defining a simple, possibly discrete, process over space and time, and then smoothing it out with one or more kernels, giving a smooth process over space and time.

This constructive approach is appealing since the resulting models can be extended to allow for generalizations such as non-stationarity, non-Gaussian models, and non-separable space-time dependence structures. See Wolpert and Ickstadt [12], and Higdon, et al. [13] for some purely spatial applications, and Higdon [14] for a space-time model. In addition, models can be constructed in such a way to facilitate computation - such as restricting the underlying process to reside on a lattice so that fast Fourier transforms can be employed.

The Likelihood and Priors

Based upon our previous assumptions, the unknown parameters in model (2.1) are

$\Omega =\left\{a_i,\mathrm{....}\text{, }{a}_n,b,\beta ,{\sigma }_{\epsilon }^2,{\mu }_a,{\sigma }_a^2,{\sigma }_u^2\right\}$. The contribution of subject $i$ to the conditional likelihood function is given by

$$\begin{array}{l}\left[Y_i\left(t_{ij}\right)|\Omega ,{\text{ X}}_i\right]=\left({\displaystyle {\prod }_{j=1}^{n_i}\left[Y_i\left(t_{ij}\right)|\Omega ,{\text{ X}}_i\right]}\right)\\ =\left({\displaystyle {\prod }_{j=1}^{n_i}\left[Y_i\left(t_{ij}\right)|a_i,b,\beta ,{\sigma }_{\epsilon }^2,Z_i\left(t_{ij}\right),U_i\left(t_{ij}\right)\right]}\right)\\ ={\displaystyle {\prod }_{j=1}^{n_i}\frac{1}{\sqrt{2\pi {\sigma }_{\epsilon }^2}}}\exp \left(\frac{-{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta Z_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}{2{\sigma }_{\epsilon }^2}\right)(3.1)\end{array}$$

Where, $[.]$ and $\left[.|.\right]$ denote marginal and conditional densities, respectively, $\Omega$ denotes all model unknown parameters.

For the prior density of $\Omega$ we assume tha $b,\beta ,{\sigma }_{\epsilon }^2,{\mu }_a,{\sigma }_a^2,and{\sigma }_u^2$ have independent prior densities, so that $\left\{a_i,i=1,\mathrm{...},n\right\}$ are independent and normally distributed with parameters ${\mu }_a$ and ${\sigma }_a^2$, also $\left\{U_i,i=1,\mathrm{...},n\right\}$ are independent stochastic process with parameter ${\sigma }_u^2$.

From the independency assumptions, the posterior density of all unknown model parameters $\left[\Omega |data\right]$ is proportional of

$$\underset{\left[U_i\left(t\right)|,{\alpha }_w,{\sigma }_w^2\right]\left[a_i|{\mu }_a,{\sigma }_a^2\right]\}\times \left[b\right]\left[{\mu }_a\right]\left[{\sigma }_a^2\right]\left[\beta \right]\left[{\sigma }_u^2\right]\left[{\sigma }_{\epsilon }^2\right]}{\{{\displaystyle {\prod }_{i=1}^N\left({\displaystyle {\prod }_{j=1}^{n_i}\left[Y_i\left(t_{ij}\right)|a_i,b,\beta ,{\sigma }_{\epsilon }^2,X_i,U_i\left(t\right)\right]}\right)}}(3.2)$$

To fit the full model and make inference about the population parameters, Adaptive Rejection Metropolis Sampling (ARMS), and Gibbs and ARMS sampling techniques are used. These methods are a MCMC technique for drawing dependent samples from complex high dimensional distributions, see Waezizadeh, and Mehrpooya [15]. The posterior distribution converges was checked by Gelman-Rubin convergence statistic R (posterior consistency), as modified by Brooks and Gelman [16]. In order to apply one of these methods on our model, posterior for each parameter must be derived, and then a proposed prior density for each of these parameters must be chosen. Based on the likelihood functions in (3.1) and (3.2), with the notations (IG) denotes for Inverse Gamma, and (N) denotes for Normal, the conditional densities of the unknown model parameters are given as follows:

For the error parameter ${\sigma }_{\epsilon }^2$

$$\begin{array}{l}\left[{\sigma }_{\epsilon }^2|.\right]\alpha {\displaystyle {\prod }_{i=1}^n{\displaystyle {\prod }_{j=1}^{n_i}\left[Y_i\left(t_{ij}\right)|a_i,b,{\sigma }_{\epsilon }^2,Z_i\left(t_{ij}\right),U_i\right]}}\left[{\sigma }_{\epsilon }^2\right]\\ \alpha {\displaystyle {\prod }_{i=1}^n{\displaystyle {\prod }_{j=1}^{n_i}\frac{1}{\sqrt{{\sigma }_{\epsilon }^2}}}}\exp \left(\frac{-{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}{2{\sigma }_{\epsilon }^2}\right)\left[{\sigma }_{\epsilon }^2\right]\\ \alpha {\displaystyle {\prod }_{i=1}^n{\displaystyle {\prod }_{j=1}^{n_i}\frac{1}{\sqrt{{\sigma }_{\epsilon }^2}}}}\exp \left(\frac{-{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}{2{\sigma }_{\epsilon }^2}\right)\left[{\sigma }_{\epsilon }^2\right]\\ \alpha \frac{1}{{\left({\sigma }_{\epsilon }^2\right)}^{\frac{{\displaystyle {\sum }_{i=1}^n{n}_i}}{2}}}\exp \left(\frac{-{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}}}{2{\sigma }_{\epsilon }^2}\right)\left[{\sigma }_{\epsilon }^2\right]\\ \alpha {\left({\sigma }_{\epsilon }^2\right)}^{-\left(\left(\frac{{\displaystyle {\sum }_{i=1}^n{n}_i}}{2}-1\right)+1\right)}\exp \left(\frac{-{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}}/2}{{\sigma }_{\epsilon }^2}\right)\left[{\sigma }_{\epsilon }^2\right]\\ \alpha IG\left({\alpha }_0,{\beta }_0\right)\left[{\sigma }_{\epsilon }^2\right]\end{array}$$

Where

${\alpha }_0=\frac{{\displaystyle {\sum }_{i=1}^m{n}_i}}{2}-1$, and ${\beta }_0=\frac{{\displaystyle {\sum }_{i=1}^m{\displaystyle {\sum }_{j=1}^{n_i}{\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}^2}}}{2}$

In the same manner we found:

The intercept variance ${\sigma }_a^2:\left[{\sigma }_a^2|.\right]\alpha \text{IG }\left({\alpha }_0,{\beta }_0\right)\left[{\sigma }_a^2\right]$, where ${\alpha }_0=\frac{n}{2}-1$, and ${\beta }_0=\frac{{\displaystyle {\sum }_{i=1}^n{\left(a_i-{\mu }_a\right)}^2}}{2}$

The intercept mean ${\mu }_a:\left[{\mu }_a|.\right]\alpha N\left({\alpha }_0,{\beta }_0\right)\left[{\mu }_a\right]$, where ${\alpha }_0=\frac{{\displaystyle {\sum }_{i=1}^n{a}_i}}{N}$, and ${\beta }_0=\frac{{\sigma }_a^2}{n}$

The random intercept $a_i$, for $i=\text{ 1,2, }\mathrm{....}\text{n,}\left[a_i|.\right]\alpha \text{ N}\left({\alpha }_0,{\beta }_0\right)$, where

$${\alpha }_0=\left(\frac{{\displaystyle {\sum }_{j=1}^{n_i}\left(Y_{ij}-\left(b{t}_{ij}+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}}{{\sigma }_{\epsilon }^2}+\frac{{\mu }_a}{{\sigma }_a^2}\right){\left(\frac{n_i}{{\sigma }_{\epsilon }^2}+\frac{1}{{\sigma }_a^2}\right)}^{-1}and{\beta }_0={\left(\frac{n_i}{{\sigma }_{\epsilon }^2}+\frac{1}{{\sigma }_a^2}\right)}^{-1}$$

The average rate of the slope $b:\left[b|.\right]\alpha N\left(m_b,{\beta }_0\right)\left[b\right]$, where

$m_b=\frac{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{t}_{ij}\left(Y_i\left(t_{ij}\right)-\left(a_i+\beta X_i\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}}}{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{t}_{ij}^2}}}and{\beta }_0=\frac{{\sigma }_{\epsilon }^n}{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{t}_{ij}^2}}}$

The effect of the regression parameter on the marker ${\beta }_l,l=1,\mathrm{...},p:\left[{\beta }_l|.\right]\alpha N\left(m_{{\beta }_l},{\beta }_0\right)\left[{\beta }_l\right]$, where

-l

$$m_{{\beta }_l}=\frac{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{X}_{il}\left(t_{ij}\right)\left(Y_i\left(t_{ij}\right)-\left(a_i+b{t}_{ij}+\beta \left(-l\right)X_i\left(-l\right)\left(t_{ij}\right)+U_i\left(t_{ij}\right)\right)\right)}}}{{\displaystyle {\sum }_{i=1}^n{\displaystyle {\sum }_{j=1}^{n_i}{X}_{il}^2\left(t_{ij}\right)}}}and$$

$${\text{β}}_{\text{0}}\text{= }\frac{{\text{σ}}_{\text{ε}}^2}{{\displaystyle {\sum }_{\text{i=1}}^{\text{n}}{\displaystyle {\sum }_{\text{j=1}}^{{\text{n}}_{\text{i}}}{X}_{il}^2\left(t_{ij}\right)}}}$$

$\beta \left(-l\right)=\left({\beta }_1,\mathrm{...},{\beta }_{l-1},{\beta }_{l+1},\mathrm{...},{\text{β}}_{\text{0}}\right)$ and $X_i\left(-l\right)$ are the remaining covariates after $l^{th}$ covariate $X_{il}$ is excluded from $X_i$.

For the stochastic process $U\left(t\right)$ parameter ${\sigma }_u^2:\left[{\sigma }_u^2|.\right]\alpha \text{IG }\left({\alpha }_0,{\beta }_0\right)\left[{\sigma }_u^2\right]$, where

${\alpha }_0=\frac{{\displaystyle {\sum }_i^m{t}_i^{\ast }}}{m}-1$, and ${\beta }_0=\frac{{\displaystyle {\sum }_{i=1}^m{U}_i^T{A}^{-1}{U}_i}}{m}$, ${\text{U}}_i=\left\{u_i\left(t_i^{\ast }\right),\mathrm{....},u_i\left(t_{mi}^{\ast }\right)\right\}$, $A^{-1}=\frac{{\sum }^{-1}}{{\sigma }_u^2}$, and $\sum$ is given in (2.2).

Since all posterior densities are in standard form, then it is easy to choose conjugate priors for all model (2.1) parameters, drawing random variates using Gibbs sampler from their full conditional distributions is straightforward since their full conditional densities are standard distributions. Therefore, we use the full conditional density as proposal density. At each updating step for these parameters, a new draw from the full conditional density is always accepted.

Simulation

To illustrate our proposed model, we setup our simulation study represents a randomize clinical trial, in which M = 500 subjects are randomized. Each longitudinal marker in model (2.1), $Y_i\left(t_{ij}\right)$, $i=1,\mathrm{...},M,j=1,\mathrm{...},n_i$ was simulated as the sum of the trajectory function $X_{il}\left(t_{ij}\right)$ and the error terms ${\epsilon }_i\left(t_{ij}\right)$, each subject has its observed longitudinal measured n i =10 at time points $t=\left\{t_1=0.1,\mathrm{...},t_{10}=1\right\}$. For the simulations considered in this research the smoothing kernel $K(.)$ will be a radially symmetric kernel, such that $K\left(t\right)\alpha \exp \left(-\frac{1}{2}{\Vert t\Vert }^2\right)$, with variance covariance function $C\left(d\right)\alpha \exp \left(-\frac{1}{2}{\left(\frac{t-s}{\sqrt{2}}\right)}^2\right)$, and defining the latent process support so that the $t_j^{\ast }$ are m = 5 equally spaced points ranging from -0.1 to 1.2. Note that this combination of kernel width and spacing of the $t_j^{\ast }$ yields a spatial process $U\left(t\right)$ via (2.3) that is nearly stationary. If the spacing become much larger, or if the kernel width is reduced, the covariance structure $C\left(d\right)$ for $U\left(t\right)$ becomes unduly influenced by sparseness artifacts.

Since the calculations for the simulation study were highly computationally intensive, we have used the cluster with about 20 nodes with AMD Quad-Core Opteron 835X, 4 × 2G Hz, and 16 GB RAM per node. We obtained 3000 iterations after a burn-in of 1000 iterations. Convergence was checked by monitoring histories of sampled quantities with several different starting points. The histogram, the time series plots of one sequence of Gibbs samples for different number of iterations and the average number of these iterations for the parameter ${\mu }_a$ are presented in Figure 1.

We also used the Gelman-Rubin convergence statistic, as modified by Brooks and Gelman [15]. They emphasize that one should be concerned both with convergence of R is the ratio pooled-width within-width (the ratio of width of the central 80% interval based on pooled runs and the average width of 80% intervals within the individual runs), to one, and with convergence of both the pooled and within interval widths to stability. For our analyses, all R's converged to 1 within 3000 iterations and hence the burn-in of 1000. The analysis for the 500 simulated data sets for a single scenario took approximately 2:30 hours to run the model under the approximate model. While it took almost 6 hours to run the analysis without approximations.

With the initial values of the parameters for which the data are generated considered as the truth values of the parameters, estimate Monte Carlo Summary statistics, Monte Carlo Standard Deviation (MCSD), Mean Squared Error (MSE), 95% Confidence Converge Rate (CCR), and Bias in Percentage Terms (BPT) are presented in Table 2, where, MCE stand for Monte Carlo Error and it can be evaluated as follows : In our simulate study we used 500 data replications, thus the resulting estimates are subject to sampling variation (Monte Carlo Error), this variation for the point estimate can be calculated as $\widehat{p}=MCSD/\sqrt{500}$, the MCE then can be found by $MCE=\sqrt{\frac{\widehat{p}\left(1-\widehat{p}\right)}{500}}$

Results in Table 2 assert the convergence of the Markov Chain and the samplers reached the convergence after 2,000 iterations after 1,000 iterations are burn-in. Posterior means, posterior standard deviations, Bias as percent of true parameter and 95% highest posterior density intervals for each parameter in the proposed and IOU models, are represented in Table 3. The estimates of all parameters from the approximate modelling analysis are quite accurate and efficient.

Discussion and Conclusion

We have proposed a new model for repeated measures that incorporated random effects, correlated stochastic process and measurements error, we proposed to approximate the IOU stochastic process $W_i$ into a spatial model that can be constructed by convolving a very simple and independent, process with a kernel function. This approach offered some advantages over specification through a spatial process of computing covariance, variogram, and extremal coefficient functions, also added to the extremal coefficient plots the empirical estimates. This approximation is attractive because it facilitates calculations, especially that contain a huge amount of data also it easily extends beyond simple stationary flexible models. Moreover, this structure can be used to significantly reduce the dimensionality of a complex temporal process.

The limitation in the proposed model is that the latent process $U^{\ast }\left(t\right)$ must be nonzero at the spatial sites $t^{\ast }$ in a space time t. To overcome this problem with the covariance function of $U\left(t\right)$ and $U\left(s\right)$ the differences through the displacement vector d = t – s was performed.

A Bayesian approach was taken to fit the proposed model through a simulation study. The numerical results demonstrate that the propose modelling method results in efficient estimates and good coverage for the population parameters. The estimates are close to the true values of the parameters and have good coverage rates. The small biases of the estimates are due to Monte Carlo simulation error. Compared to the IOU model, the proposed model results in improved estimates almost for all parameters. The proposed model demonstrated significant reductions in execution time. This approach effectively eliminates the deficiency of non-spatial huge data access by replacing such patterns in hotspots of applications with spatial sites data at a space time t at runtime. The execution time reduced by 60% illustrate the efficiency of proposed model.

References

1. Li N, Li Y, Liu Y (2015) Empirical bayes inference for the parameter of pareto distribution based on longitudinal data. International Journal of Statistics & Economics 16.
2. Wulfsohn MS, Tsiatis AA (1997) A joint model for survival and longitudinal data measured with error. Biometrics 53: 330-339.
3. Taylor JMG, Cumberland WG, Sy JP (1994) A Stochastic Model for Analysis of Longitudinal data. Journal of the American Statistical Association 89: 727-736.
4. Xu J, Zeger SL (2001) Joint analysis of longitudinal data comprising repeated measures and times to events. Applied Statistics 50: 375-387.
5. Wang Y, Taylor JMG (2001) Jointly modeling longitudinal and event time data: Application in AIDS studies. JASA 96: 895-905.
6. Abu Bakar MR, Salah KA, Ibrahim NA, Haron K (2009) A stochastic joint model for longitudinal and survival data with cure patients. International Journal of Tomography & Statistics 11: 48-67.
7. Henderson R, Diggle P, Dobson A (2000) Joint modeling of longitudinal measurements and event time data. Biostatistics 4: 465-480.
8. Higdon D (2002) Space and Space-time Modeling Using Process Convolutions. In: Quantitative methods for current environmental issues, Springer Verlag, 37-56.
9. Silverman BW (1986) Density Estimation for Statistics and Data Analysis. (1^st edn), London: Chapman and Hall.
10. Cleveland WS (1979) Robust Locally Weighted Regression and Smoothing Scatterplots. Journal of the American Statistical Association 74: 829-836.
11. Ver Hoef J, Barry RP (1998) Constructing and fitting models for cokriging and multivariable spatial prediction. Journal of Statistical Planning and Inference 69: 275-294.
12. Wolpert RL, Ickstadt K (1998) Poisson/gamma random field models for spatial statistics. Biometrika 85: 251-267.
13. Higdon DM, Swall J, Kern J (1999) Non-stationary Spatial Modeling. In: Bayesian Statistics 6, Proceedings of the Sixth Valencia International Meeting, Oxford University Press, 761-768.
14. Higdon DM (1998) A process-convolution approach to modeling temperatures in the north atlantic ocean. Journal of Environmental and Ecological Statistics 5: 173-190.
15. Waezizadeh T, Mehrpooya A (2016) A stochastic model for dynamics of two populations and its stability. Proceedings of the 47^th Annual Iranian Mathematics Conference (AIMC47) (Karaj), 10681071.
16. Brooks S, Gelman A (1998) General Methods for Monitoring Convergence of Iterative Simulations. Journal of Computational and Graphical Statistics 7: 434 -455.

---