Accessible Published by De Gruyter July 7, 2018

A Spatio-Temporal Model and Inference Tools for Longitudinal Count Data on Multicolor Cell Growth

PuXue Qiao, Christina Mølck, Davide Ferrari and Frédéric Hollande


Multicolor cell spatio-temporal image data have become important to investigate organ development and regeneration, malignant growth or immune responses by tracking different cell types both in vivo and in vitro. Statistical modeling of image data from common longitudinal cell experiments poses significant challenges due to the presence of complex spatio-temporal interactions between different cell types and difficulties related to measurement of single cell trajectories. Current analysis methods focus mainly on univariate cases, often not considering the spatio-temporal effects affecting cell growth between different cell populations. In this paper, we propose a conditional spatial autoregressive model to describe multivariate count cell data on the lattice, and develop inference tools. The proposed methodology is computationally tractable and enables researchers to estimate a complete statistical model of multicolor cell growth. Our methodology is applied on real experimental data where we investigate how interactions between cancer cells and fibroblasts affect their growth, which are normally present in the tumor microenvironment. We also compare the performance of our methodology to the multivariate conditional autoregressive (MCAR) model in both simulations and real data applications.

1 Introduction

Longitudinal image data based on fluorescent proteins play a crucial role for both in vivo and in vitro analysis of various biological processes such as gene expression and cell lineage fate. Assessing the growth patterns of different cell types within a heterogeneous population and monitoring their interactions enables biomedical researchers to determine the role of different cell types in important biological processes such as organ development and regeneration, malignant growth or immune responses under various experimental conditions. For example, tumor progression has been shown to be affected by bidirectional interactions among cancer cells or between cancer cells and cells from the microenvironment, including tumor-infiltrating immune cells [1]. Being able to study these interactions in a laboratory setting is therefore highly relevant, but is complicated by the difficulty of dissecting the effect of the different cell types as soon as the number of cell types exceeds two. In the present study we used longitudinal image data collected from multicolor live-cell imaging growth experiments of co-cultures of cancer cells and fibroblasts (a key cell type in the tumor microenvironment) as well as behaviourally distinct (cloned) cancer cells. Using a high-content imaging system, we were able to acquire characteristics for each individual cell at subsequent times, including fluorescent properties, spatial coordinates, and morphological features. The motivation of this work was to design a model allowing the determination of spatio-temporal growth interactions between these multiple cell populations.

In longitudinal growth experiments, the two important goals are to determine growth rates for different cell populations and to assess how interactions between cell types may affect their growth. Whilst a wide range of descriptive data analysis approaches have been used in applications, inference based on a comprehensive model of multicolor cell data is an open research area. The main challenges are related to the presence of complicated spatio-temporal interactions amongst cells and difficulties related to tracking individual cells across time from image data. Typical longitudinal experiments consist of a relatively small number of measurements (e.g. 5 to 20 images taken every few hours), which is adequate for monitoring cell growth. Tracking individual cells would typically require more frequent measurements, complicating the practicality of the experiments in terms of the storage cost of very large image files and the cytotoxicity induced by the imaging process.

Although tracking individual cell trajectories is difficult due to cell migration, overlapping cells, changes in cell morphology, image artifacts, cell death and division, obtaining cell counts by cell type (represented by a certain color) is straightforward and can be easily automated. To describe the spatial distribution for different cell types, we propose to divide an image into a number of contiguous regions (tiles) to form a regular lattice structure as shown in Figure 1(a). We then record the frequency of cells of different colors in each tile at subsequent time points, and based on which we model the spatial and temporal dependencies of the cell growth.

Figure 1: (a) Microscope images for the cancer cell growth data obtained from a high-content imager (Operetta, Perkin Elmer) at the initial and final time points of the experiment. In each image, colors for non-fluorescent fibroblasts, as well as red and green fluorescent cancer cells are merged. (b) Illustration of the local structure for the model in (1). The two planes correspond to 3×3$3 \times 3$ tiles at times t$t$ and t+1$t+1$. The average number of cells of color c$c$ in a given tile at time t+1$t+1$ is assumed to depend on the number of cells of other colors in contiguous neighboring tiles at time t$t$.

Figure 1:

(a) Microscope images for the cancer cell growth data obtained from a high-content imager (Operetta, Perkin Elmer) at the initial and final time points of the experiment. In each image, colors for non-fluorescent fibroblasts, as well as red and green fluorescent cancer cells are merged. (b) Illustration of the local structure for the model in (1). The two planes correspond to 3×3 tiles at times t and t+1. The average number of cells of color c in a given tile at time t+1 is assumed to depend on the number of cells of other colors in contiguous neighboring tiles at time t.

To model spatio-temporal data, one could choose to approximate the spatio-temporal process by a spatial process of time series, that is, to view the process as a multivariate spatial process where the multivariate dependencies are inherited from temporal dependencies. In other words, it can be seen as a temporal extension of spatial processes.

The most popular way of developing a spatial process is through the conditionally auto-regressive (CAR) model proposed by Besag [2]. Waller, Carlin, Xia, and Gelfand [3] extend the CAR model into a spatio-temporal setting by allowing spatial effects to vary across time. However, the model lacks a specification of temporal dependency, as also noted by Knorr-Held [4]. More recently, Quick, Waller, and Casper [5] proposed a multivariate space-time CAR (MSTCAR) model, which is essentially a multivariate CAR model, where both temporal and between group dependencies are modelled as multivariate dependencies. Other works related to spatial process of time series include Sans, Schmidt, Nobre, et al. [6] and Quick, Waller, and Casper [7].

Alternatively, one also think of the process as a time series of spatial process, or a spatial extension of time series. This is the approach we take in our spatio-temporal modelling. The underlying notion is that “the temporal dependence is more natural to model than the spatial dependence” [8].

Following Cox et al. [9], it is useful to distinguish two modelling approaches for the analysis of time series data commonly seen in spatial-temporal modelling literature: the parameter-driven and observation-driven model. In a parameter-driven model, the dependence between subsequent observations is modelled by a latent stochastic process, which evolves independently of the past history of the observation process. In contrast, in an observation-driven model, time dependence arises because the conditional expectation of the outcome given the past depends explicitly on the past values.

For multivariate count data, the advantage of parameter-driven models is that one can easily assume that the conditional expectation of the observed process (on log-scale), as a latent process, is (multivariate) normal. There are extensive works related to latent spatio-temporal models under the Bayesian framework, including models with Gaussian data modelled by (multivariate) Gaussian process with an additive error [10, 11, 12, 13], Poisson data with conditional expectation modelled by Gaussian latent process ([14, 15] and Chapter 7 of [8]) and Poisson data with multivariate log-gamma latent process [16]. However, estimation of parameters in parameter-driven models requires considerable computational effort, as does prediction of the latent process.

On the other hand, in observation-driven models, inference is possible in a (penalized) maximum likelihood framework and therefore can be easily fitted even for quite complex regression models [17]. Schrödle, Held, and Rue [18] proposed a parameter-driven spatio-temporal model and compared it with a similar observation-driven model proposed by Paul, Held, and Toschke [19]. They conclude that the parameter-driven models perform slightly better in terms of prediction in some cases, however, while the computation time for the observation-driven model is mostly less than a second, fitting a parameter-driven model takes several hours if it ever converges, because of the complexity with the latent autoregressive process. Besides, their model contains only five parameters, while in our application, the number of parameters of interest grows quadratically with the number of cell populations, which will make the parameter-driven models intractable even with a moderate number of cell populations.

Therefore, we choose to work with a spatial extension of observation-driven time series. Zeger and Qaqish [20] review various observation-driven time series models with a quasi-likelihood estimation. Fokianos and Tjøstheim [21] develop and study the probabilistic properties of a log-linear autoregressive time series model for Poisson data, as an extension of the model considered by Fokianos, Rahbek, and Tjøstheim [22]. See Scott, et al. and Kedem and Fokianos [23, 24] for a complete review.

Literature about observation-driven spatio-temporal models, however, is relatively sparse. Held, Höhle, and Hofmann [25] propose a multivariate time series model where parameters are allowed to vary across space. Paul et al. [19] extended the model such that spatial dependences are captured by additional parameters that quantify the “directed influence” of neighboring areas at previous time points on the observation of interest. Paul and Held [26] further extend the model by introducing random effects. Note that these approaches model directly the conditional expectation of the count data, meaning they are using an identity link function, instead of the canonical log-link. Thus, it is required that the parameters are positive to ensure that the resulting conditional expectation is positive. Knorr-Held and Richardson [27] propose a space-time model for surveillance data, apart from separate seasonal and spatial components, they include an autoregressive term with a latent indicator.

In this paper, we develop a conditional spatial-temporal model for multivariate count data on tiled images, and provide its application on tiled images in the context of longitudinal cancer cell monitoring experiments. Our model enables us to measure the effect on the growth rate of each cell population and changes due to local cross-population interactions. Specifically, we consider a multivariate Poisson model with intensity modeled as a log-linear form similar to those in [27] and [21], and we quantify spatio-temporal impacts of different cell populations in neighboring tiles through model parameters, as illustrated in Figure 1(b). Impacts are allowed to be positive or negative, and unlike those models that describe between group dependence through a covariance matrix, influences do not have to be symmetrical in our model. Another main advantage of the proposed framework is that it enables one to accommodate spatio-temporal cell interactions for heterogeneous cell populations within a relatively parsimonious statistical model.

Since the model complexity can be potentially very large in the presence of many cell types, it is also important to address the question of how to select an appropriate model by retaining only the meaningful spatio-temporal interactions between cell populations We cary out model selection using the common model selection criteria for parametric models, the Akaike and the Bayesian information criteria (AIC and BIC).

The remainder of the paper is organized as follows. In Section 2, we introduce the conditional spatio-temporal lattice model for multivariate count data and develop maximum likelihood inference tools. In the same section, we discuss the asymptotic properties of our estimator and standard errors. In Section 3, we study the performance of our methodology using simulated data, and compare it to that of the multivariate conditional autoregressive (MCAR) model. In Section 4, we apply our method, as well as the MCAR model to analyze datasets from an in-vitro experiment, where cancer cells are co-cultured with fibroblasts. In Section 5, we conclude and give final remarks.

2 Methods

2.1 Multicolor spatial autoregressive model on the lattice

Let LN2 be a discrete lattice. In the context of our application, the lattice is obtained by tiling a microscope image into nL tiles, denoted by Ln(L). The total number of tiles nL is a monotonically increasing function of n. One can choose various forms of lattice, for example, the regular or hexagonal lattices. For simplicity, we tile the image into n×n regular rectangular tiles, which makes nL=n2. An example of a tiled image with n=10 is shown in Figure 1(a). Denote a pair of neighboring tiles {i,j} with ij, if tiles i and j share the same border or coincide (i=j). Each tile may contain cells of different colors; thus, we let C={1,,nC} be a finite set of colors and denote by nC the total number of colors. Let Y={Yt,t=1,,T} be the sample of observations where Yt={Yt(c),cC} is the collection of observations at time point t, and Yt(c)=(Y1,t(c),,YnL,t(c)) is the vector of observed frequencies for color c on the lattice Ln at time t. The joint distribution for the spatio-temporal process on the lattice is difficult to specify, due to local spatial interactions for neighboring tiles and global interactions occurring at the level of the entire image. An additional issue is that cells tend to be clustered together due to the cell division process and other biological mechanisms; thus it is not uncommon to observe low counts in a considerable portion of tiles. In typical longitudinal experiments, the number of time points seldom go beyond 50 due to experimental, storage and processing cost, while nL can be relatively large. So we work under the framework where T is assumed to be finite, while nL is allowed to grow to infinity.

We suppose that the count for the ith tile Yi,t(c) follows a marginal Poisson distribution Yi,t(c)|Yt1Pois(λi,t(c)), with intensity modeled by the canonical log-link vi,t(c)=logλi,t(c), where vi,t(c) takes the following spatial autoregressive form:


for all cC,t=1,,T, with ni={#j:ij,jLn} being the number of tiles in a neighborhood of tile i. Although we are adopting the regular grids for simplicity, the model is readily applicable to other tiling strategies. Changing the tiling strategy would only change the realisations of Si,t1(c) in (2).

Here, we assume that the conditional count for different tiles at time t is independent conditioning on information from t1, i.e.


for all c,cC,t=1,,T, and i,jLn,ij. This does not suggest that they (Yi,t(c) and Yj,t(c)) are independent, but rather that their spatio-temporal dependence is due to the structure of intensity λi,t(c) in (1). Conditional independence is a commonly used assumption for spatio-temporal models in a non-gaussian setting [3, 28], since it’s exceedingly difficult to work with multivariate non-Gaussian distribution [8].

The elements of the parameter vector α=(α(1),,α(nC)) are main effects corresponding to a baseline average count for cells of different colors. The spatio-temporal interactions are measured by the statistic Si,t1(c) in (2), which essentially counts the number of cells of color c in the neighborhood of tile i at time t1. Hence, the autoregressive parameter β(c|c) is interpreted as positive or negative change in the average number of cells with color c, due to interactions with cells of color c in neighbouring tiles. A positive (or a negative) sign of β(c|c) means that the presence of cells of color c in neighboring tiles promotes (or inhibits) the growth of cells of color c. The spatio-temporal effects βc|c,c,cC, are collected in the nC×nC weighted incidence matrix B. This may be used to generate weighted directed graphs, as shown in the example of Figure 2, where the nodes of the directed graph correspond to cell types, and the directed edges are negative or positive spatio-temporal interactions between cell types.

Equation (1) could be extended to some more specific form, for example, vi,t(c)=α(c)+cCβ1(c|c)Si,t1(c)+β0(c|c)log1+Yi,t1(c), where β1(c|c) are interpreted as the effect of cells of color c from neighbouring (but not the same) tiles have on the growth of cells with color c, while β0(c|c) as the effect of cells of color c from the same tile. However, we stick to the model in (1) because we have no evidence showing that the more complex model is advantageous from a model selection view point.

We choose to work with a log-linear form for the autoregressive equation of vi,t(c) in eq. (1), where we apply a logarithmic transform and add 1 to the counts at time t1, Yi,t1(c). It offers several advantages compared to the more commonly used linear form. First, λi,t(c) and Yi,t1(c) are transformed on the same scale. Moreover, this model can accommodate both positive and negative correlations, while it is not possible to account for positive association in a stationary model if past counts are directly included as explanatory variables. For example, with the model vi,t=α+βYi,t1 for a single color, the intensity would be λi,t=expαexpβYi,t1, which may lead to instability of the Poisson means if β>0 since λi,t is allowed to increase exponentially fast. Finally, adding 1 to Yi,t1(c) is for coping with zero data values, since log(Yi,t1(c)) is not defined when Yi,t1(c)=0, which arises often, and it maps zeros of Yi,t1(c) into zeros of log(1+Yi,t1(c)).

2.2 Likelihood inference

Let θ be the overall parameter vector θ=(α,vec(B))Rp, where α is a nC-dimensional vector defined in Section 2.1 and B is a nC×nC matrix of color interaction effects, p=nC(1+nC) is the total number of parameters. In this section, we develop a weighted maximum likelihood estimator for our model,


where λi,t(c)(θ) is the expected number of cells with color c in tile i at time t, defined in (1) and the weights wi,t(c) are given constants. The weighted maximum likelihood estimator (MLE), θˆ, is obtained by maximizing the weighted log-likelihood function


where vi,t(c)(θ)logλi,t(c)(θ). Equivalently, θˆ is formed by solving the weighted estimating equations


where γi,t(θ)=yi,t(1)expvi,t(1)(θ),,yi,t(nC)expvi,t(nC)(θ), denotes the Kronecker product, is the gradient operator with respect to θ and vi,tvi,t(c)(θ)=(1,Si,t1(1),,Si,t1(nC)).

Specific weights could be used to address the presence of outliers. Following Ferrari and Vecchia [29] and La Vecchia, Camponovo, and Ferrari [30], the influence of strong outliers could be avoided by taking weights of form wi,t(c)=exp{(1q)[Yi,t(c)v(c)(θ)exp(v(c)(θ))]} with q being a tuning constant smaller than 1. However, for the current application we use constant weights all equal to 1.

Our empirical results show that this choice performs reasonably well in terms of estimation accuracy in all our numerical examples and guarantees optimal variance for the estimator θˆ under correct model specification. The solution to eq. (5) is obtained by a standard Fisher scoring algorithm, which is found to be stable and converges fast in all our numerical examples.

Finally, in practical applications it is also important to address the question of how to select an appropriate model by retaining only the meaningful spatio-temporal interactions between cell populations, and avoid over-parametrized models. Model selection plays an important role by balancing goodness-of-fit and model complexity. Here, we select non-zero model parameters based traditional model selection approaches: the Akaike Information criterion, AIC=2(θˆ)+2p, and the Bayesian information criterion, BIC=2(θˆ)+plog(|nLT|).

2.3 Asymptotic properties and standard errors

In this section, we overview the asymptotic behavior of the estimator introduced in Section 2.2. In our setting we consider a fixed number of time points, T, whilst the lattice Ln is allowed to increase. This reflects the notion that the statistician is allowed to choose an increasingly fine tiling grid as the number of cells increases. If the regularity conditions stated in the Appendix hold, then nLHn(θ0)1/2(θˆnθ0) converges in distribution to a p-variate normal distribution with zero mean vector and identity variance, as nL, with Hn(θ) given in (6). Asymptotic normality of θˆn follows by applying the limit theorems for M-estimators for nonlinear spatial models developed by Jenish and Prucha [31]. One condition required to ensure this behaviour is that Yt has constant entries at the initial time point t=0, which is quite realistic since typically cells are seeded randomly at the beginning of the experiment. Our proofs mostly check α-mixing conditions and L2-Uniform Integrability of the score functions ui,t(θ) ensures a pointwise law of large numbers, with additional stochastic equicontinuity, a uniform version of the law of large numbers required by Jenish and Prucha [31].

The asymptotic variance of θˆ is Vn(θˆ)=Hn1(θ0), where Hn(θ) is the p×p Hessian matrix


with ui(θ)=ui,1(θ)++ui,T(θ) being the partial score function for the ith tile. Direct evaluation of H(θ) may be challenging since the expectations in (6) is intractable. Thus, we estimate Hn(θ) by the empirical counterpart




Note that the above estimators approximate the quantities in formula (6) by conditional expectations. Our numerical results suggest that the above variance approximation yields confidence intervals with coverage close to the nominal level (1α). Besides the above formulas, we also consider confidence intervals obtained by a parametric bootstrap approach. Specifically, we generate B bootstrap samples Y(1),,Y(B) by sampling at subsequent times from the conditional model specified in eqs. (1) and (2) with θ=θˆ. From such bootstrap samples, we obtain bootstrapped estimators, θˆ(1),,θˆ(B), which are used to estimate var(θˆ0) by the usual covariance estimator Vˆboot(θˆ)=b=1B(θˆ(b)θ)2/(B1), where θ=b=1Bθˆ(b)/B. Finally, a (1α)100% confidence interval for θj is obtained as θˆj±z1α/2{Vˆ}jj1/2, where zq is the q-quantile of a standard normal distribution, and Vˆ is an estimate of var(θˆ) obtained by either eq. (7) or bootstrap resampling.

3 Monte Carlo simulations

In our Monte Carlo experiments, we generate data from a Poisson model as follows. At time t=0, we populate nL tiles using equal counts for cells of different colors. For t=1,,T, observations are drawn from the multivariate Poisson model Yi,t(c)|Yt1Poisson(λi,t(c)),cC. Recall that the rate λi,t(c) defined in Section 2.1 contains autoregressive coefficients β(c|c), which are collected in the nC×nC matrix B.

We assess the performance of MLE under different settings concerning the size and sparsity of B. Consider the three models with the following choices of B:


Denote Model i as the model corresponding to Bi,i=1,2,3. In Model 1, all the effects in B have the same size; in Model 2, the effects have decreasing sizes; Model 3 is the same as Model 1, but with some interactions exactly equal to zero.

We set α(1)==α(nC)=0.1 for all three models. The above parameter choices reflect the situation where the generated process Y has a moderate growth.

In Table 1 and Table 2, we show results based on 1000 Monte Carlo runs generated from Models 1-3, for n=25,nC=3 and T=10 and 25. In Table 1, we show Monte Carlo estimates of squared bias and variance of θˆ. Both squared bias and variance of our estimator are quite small in all three models, and decrease as T gets larger. The variances of Model 2 are slightly larger than those in the other two models due to the increasing difficulty in estimating parameters close to zero.

Table 1:

Monte Carlo estimates for squared bias (×106) and variance (×104) of the MCLE for three models with time points T=10,25. Simulation standard errors are shown in parenthesis. The three models differ in terms of the coefficients β(c|c),c,cC, as described in Section 3: Non-zero equal effects (Model 1), non-zero decreasing interactions (Model 2), and sparse effects (Model 3). For all models, α(c)=0.1,c=1,2,3. Estimates are based on 1000 Monte Carlo runs.

Model 10.45(0.57)5.75(0.26)0.29(0.32)2.36(0.11)
Model 20.64(0.91)9.66(0.42)0.67(0.71)4.45(0.20)
Model 30.77(0.97)8.09(0.36)0.52(0.51)3.47(0.16)

In Table 2, we report the coverage probability for symmetric confidence intervals of the form θˆ±z1α/2sdˆ(θˆ), where zq is the qquantile for a standard normal distribution, with α=0.01,0.05,0.10. The standard error, sdˆ(θˆ), is obtained by the squared root of diagonal elements of Vn(θˆ) and the parametric bootstrap estimate, Vˆest and Vˆboot, described in Section 2.3. The coverage probability of the confidence intervals are very close to the nominal level for both methods.

Table 2:

Monte Carlo estimates for the coverage probability of (1α)% confidence intervals θˆ±z1α/2sdˆ(θˆ), with sdˆ(θˆ) obtained using bootstrap (Vˆboot) and sandwich (Vˆest) estimators in Section 2 and 3. The three models differ in terms of the coefficients β(c|c),c,cC as described in Section 3: Non-zero equal effects (Model 1), non-zero decreasing interactions (Model 2), and sparse effects (Model 3). For all models, α(c)=0.1,c=1,2,3, estimates are based on 1000 Monte Carlo runs.

Model 198.699.098.999.0
Model 398.999.098.998.9
Model 395.495.594.995.1
Model 189.290.390.190.3
α=0.10Model 290.690.089.790.0
Model 390.690.690.290.2

In Table 3, we show results for the model selection based on 1000 Monte Carlo samples from Model 3 using the AIC and the BIC given in Section 2 for n=25 and T=10,25. We report Type A error (a term is not selected when it actually belongs to the true model ) and Type B error (a term is selected when it is not in the true model ). For both AIC and BIC model selection is more accurate for large T. As expected AIC tends to over select, and BIC outperforms AIC, with zero Type A error, and very low Type B error.

Table 3:

Monte Carlo estimates for % Type A error (a term is not selected when it actually belongs to the true model) and % Type B error (a term is selected when it is not in the true model) using AIC and BIC criteria. Results are based on 1000 Monte Carlo samples generated from Model 3 with n=25 and T=10,25.

Type AType BType AType B

Finally, we compare the performance of our model with the following Multivariate conditional autoregressive (MCAR) model proposed by Leroux, Lei, and Breslow [32]:


where Zi,iLn are random effects with conditional distribution


where ρ is a spatial autocorrelation parameter, with ρ=0 corresponding to independence, while ρ=1 corresponds to the intrinsic model, and ΣZ is a nCT×nCT between variable covariance matrix, which is assumed to have no fixed structure, and ni is the number of tiles in a neighborhood of tile i as defined in Section 2.1. Let β=(αT,vec(B)T)T be a vector of regression parameters, where B is defined in Section 2.1 and α is the intercept. Let the covariate xi,t be a nC2-dimensional vector consists of nC vectors: (Si,t1(1),,Si,t1(nC)), where Si,t1(c) carries the information from the neighbouring tiles on the previous time point, defined in eq. (2).

An independent Gaussian prior, N(0,100000), is specified for each regression parameter in β. A uniform prior on the unit interval, U(0,1), is specified for ρ. For covariance matrix ΣZ, assume an inverse Wishart distribution with identity scale matrix and nCT degree of freedom.

To evaluate the performance of MLE under our model and estimators obtained by the MCAR model, we generate 1000 set of data from Model 1. Estimation of the MCAR model is done by MCMC sampling, using R package CARBayes by Lee [33]. Table 4 show Monte Carlo estimates of squared bias, variance, the coverage probability of 95% confidence intervals and computation time for n,T{10,25} and nC=1,2,3. Two of the settings are the same as those shown for Model 1 in Table 1: n=25,nC=3,T=10 and n=25,nC=3,T=25 . In estimation of MCAR, we also show results of two MCMC settings: 1. MCAR1: 1000 MCMC samples generated and 200 discarded as the burn-in period; 2. MCAR2: 5000 samples with 100 discarded. Coverage probabilities of our model is computed as θˆ±z0.975sdˆ(θˆ), where zq is the qquantile for a standard normal distribution. The standard error, sdˆ(θˆ), is obtained by taking the squared root of diagonal elements of Vn(θˆ) described in Section 2.3.

Table 4:

Monte Carlo estimates for squared bias (×106), variance (×104), the coverage probability of 95% confidence intervals as well as computation time for n,T{10,25} and nC=1,2,3 of MLE of our model, and MCAR, where in MCAR1, 1000 MCMC samples generated and 200 discarded as the burn-in period; and in MCAR2, 5000 samples with 100 discarded. True values of regression parameters are shown as B1. Estimates are obtained from 1000 Monte Carlo runs.

n = 10Biasˆ2Varˆ95%TimeBiasˆ2Varˆ95%Time
nC= 1Model 12.89(3.40)21.22(0.98)94.51s1.29(1.62)10.18(0.46)94.63s
nC= 2Model 13.86(3.87)43.81(1.95)95.32s3.51(1.96)33.64(1.52)94.63s
nC= 3Model 14.83(5.10)34.66(1.59)94.93s4.4(2.48)14.21(0.65)94.77s
n = 25T=10T=25
nC= 1Model 10.51(0.56)3.18(0.14)94.910s0.40(0.37)1.73(0.08)94.523s
nC= 2Model 10.76(0.34)13.22(0.56)94.48s0.59(0.67)5.69(0.25)94.419s
nC= 3Model0.67(0.66)5.91(0.27)94.624s0.31(0.44)2.35(0.14)94.955s

In overall, our method performs better than MCAR at analysing the kind of data that we generate, especially when n and/or T is small, with much smaller bias and variance, as well as computation time. The performance of MCAR improves significantly as the model gets more complicated (i.e. larger nC), and when n and T increases. In the case where n=25,T=25 and nC=3, it almost performs equally well with our model, however, it takes almost an hour to obtain the estimates, while our method requires less than a minute. Besides, for the coverage probabilities to reach the nominal level, it seems that MCAR requires larger MCMC sample size as the model gets more complicated, while those of our model has been stable and close to the nominal level in all cases.

4 Analysis of the cancer cell growth data

Cancer cell behaviour is believed to be determined by several factors including genetic profile and differentiation state. However, the presence of other cancer cells and non-cancer cells has also been shown to have a great impact on overall tumor behaviour [34, 35]. It is therefore important to be able to dissect and quantify these interactions in complex culture systems. The data sets in this section represent a cancer cell-fibroblast co-culture experiment. The data sets analyzed consist of counts of cell types (different cancer cell populations expressing different fluorescent proteins, and non-fluorescent fibroblasts) from 9 subsequent images taken at an 8-hour frequency over a period of 3 days using the Operetta high-content imager (Perkin Elmer). Information regarding cell type (fluorescent profile) and spatial coordinates for each individual cell were extracted using the associated software (Harmony, Perkin Elmer).

Each image was subsequently tiled using a 25×25 regular grid.

We choose the number of tiles for a balance between the fit of the model and capturing the local impact between cell populations. More specifically, decreasing tile sizes enables one to detect local impacts between cell populations, which is one of the objectives of our analysis. However, if the tiles are too small, we will end up with mostly no cells in most tiles. In this situation the conditional Poisson model would not fit well the data. On the other hand, when the tiles are too large the model would fit the data well (the conditional Poisson would be approximately a conditional normal model), but we lose information on local impacts. We recommend 0 to 20 average cells per tile, since for such choice our diagnostic and goodness-of-fit analyses suggest that the conditional Poisson model fits well the data whilst enabling us to measure local correlation effects between populations.

4.1 Cancer cell-fibroblast co-culture experiment

In this experiment, cancer cells are co-cultured with fibroblasts, a predominant cell type in the tumor microenvironment, believed to affect tumor progression, partly due to interactions with and activation by cancer cells [34]. In this experiment, fibroblasts (F) are non-fluorescent whereas cancer cells fluoresce either in the red (R) or green (G) channels due to the experimental expression of mCherry or GFP proteins, respectively. Cells were initially seeded at a ratio of 1:1:2 (R:G:F).

Model selection and inference.   We applied our methodology to quantify the magnitude and direction of the impacts have on growth for the considered cell types. To select the relevant terms in the intensity expression (1), we carry out model selection using the BIC model selection criterion. In Table 5, we show estimated parameters for the full and the BIC models, with bootstrap 95% confidence intervals in parenthesis. Figure 2 illustrates estimated spatio-temporal impacts between cell types using a directed graph. The solid and dashed arrows represent respectively significant and not significant impacts between cell types at the 95% confidence level. Significant impacts coincide with parameters selected by BIC.

The interactions within each cell type (βˆ(c|c),c=R,G,F) are significant, which is consistent with healthy growing cells. As anticipated, the effects βˆ(c|c) for the cancer cells are larger than those for the slower growing fibroblasts. The validity of the estimated parameters is also supported by the similar sizes of the parameters for the green and red cancer cells. This is expected, since the red and green cancer cells are biologically identical except for the fluorescent protein they express. Interestingly, the size of the estimated effects within both types of cancer cells (βˆ(c|c),c=R,G) are larger than the impact they have on one another (βˆ(G|R) and βˆ(R|G)). This is not surprising, since βˆ(c|c)(c=R,G) reflects not only impacts between cells from the same cell population, but also cell proliferation. The fact that we are able to detect the impacts between the red and green cancer cells confirms that our methodology is sensitive enough to detect biologically relevant impacts even though no interactions were found between the cancer cells and the fibroblasts. This might be due to the fact that we used normal fibroblasts that had not previously been in contact with cancer cells and thus had not been activated to support tumor progression as is the case with cancer-activated fibroblasts.

Figure 2: Directed graph showing fitted spatio-temporal interactions between GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). The solid and dashed arrows represent respectively the significant and not significant interactions between cell types at the 95%$95\%$ confidence level.

Figure 2:

Directed graph showing fitted spatio-temporal interactions between GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). The solid and dashed arrows represent respectively the significant and not significant interactions between cell types at the 95% confidence level.

Goodness-of-fit and one-step ahead prediction   To illustrate the goodness-of-fit of the estimated model, we generate cell counts for each type in each tile, yˆi,t(c), from the Pois(λˆi,t(c)) distribution for t1, where λˆi,t(c) is computed using observations at time t1, with parameters estimated from the entire dataset. In Figure 4, we compare the actually observed and generated cell counts for GFP cancer cells (G) and mCherry cancer cells (R) and fibroblasts (F) across the entire image. The solid and dashed curves for all cell types are close, suggesting that the model fits the data reasonably well. As anticipated, the overall growth rate for the red and green cancer cells are similar, and sensibly larger than the growth rate for fibroblasts.

To assess the prediction performance of our method, we consider one-step-ahead forecasting using parameters estimated from a moving window of five time points. In Figure 3, we show quantiles of observed cell counts against predicted counts for each tile. The upper and lower 95% confidence bounds are computed non-parametrically by taking Fˆ11(Fˆ0(yt(c))0.95) and Fˆ11(Fˆ0(yt(c))+0.95), where Fˆ0 and Fˆ1 are the empirical distributions of the observations and predictions at time t respectively [36]. The identity line falls within the confidence bands in each plot, indicating a satisfactory prediction performance.

Figure 3: QQ-plots for cell growth, comparing observed (horizontal axis) and one-time ahead predicted (vertical axis) cell counts per tile on the entire image at times t=6,7,8$t=6,7,8$ for GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). One-time ahead predictions are based on the model fitted using a moving window of five time points.

Figure 3:

QQ-plots for cell growth, comparing observed (horizontal axis) and one-time ahead predicted (vertical axis) cell counts per tile on the entire image at times t=6,7,8 for GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). One-time ahead predictions are based on the model fitted using a moving window of five time points.

Comparison with MCAR model   Next, we compare the estimates as well as the goodness-of-fit on the real data with the MCAR model. Parameter estimates are shown in Table 5, with 95% confidence intervals given in parenthesis. Results from both models are mostly consistent with each other, specifically, both models show that impacts within each cell type (βˆ(c|c),c=R,G,F) are significant, the effects βˆ(c|c) for cancer cells are larger than those for the slower growing fibroblasts, the green and red cancer cells have positive impact on each other, and cancer cells have no impact on fibroblasts. The only difference is, the MCAR model shows a negative impact of fibroblasts on the green cancer cells only, while our model detect no significant impact on either cancer cells. Since the red and green cancer cells are biologically identical except for the fluorescent protein they express, we expect a symmetrical result with both cancer cells.

Table 5:

Estimated parameters for the full, the BIC models and the MCAR model based on the cancer cell growth data described in Section 4. Bootstrap 95% confidence intervals based on 50 bootstrap samples are given in parenthesis.

Full model
αˆ(c)−0.99 (−1.19, −0.79)−0.50 (−0.70, −0.30)−0.26 (−0.45, −0.06)
βˆ(G|c)1.23 (1.10, 1.35)0.34 (0.21, 0.48)0.12 (−0.03, 0.27)
βˆ(R|c)0.28 (0.17, 0.38)1.09 (0.96, 1.21)0.02 (−0.09, 0.13)
βˆ(F|c)0.10 (−0.01, 0.21)0.02 (−0.07, 0.12)0.92 (0.81, 1.03)
BIC model
αˆ(c)−0.88 (−1.04, −0.72)−0.49 (−0.66, −0.31)−0.19 (−0.36, −0.02)
βˆ(G|c)1.24 (1.11, 1.37)0.35 (0.21, 0.48)/
βˆ(R|c)0.28 (0.17, 0.39)1.09 (0.96, 1.21)/
βˆ(F|c)//0.93 (0.82, 1.04)
αˆ(c)−0.45 (−0.54, −0.38)−0.45 (−0.54, −0.38)−0.45 (−0.54, −0.38)
βˆ(G|c)1.06 (0.93, 1.16)0.16 (0.09, 1.16)−0.15 (−0.22, −0.09)
βˆ(R|c)0.25 (0.15,0.31)1.01 (0.92,1.08)0.05 (−0.02,0.10)
βˆ(F|c)0.03 (−0.06,0.20)0.03 (−0.07,0.19)0.96 (0.83,1.08)

In Figure 4, apart from the observed (solid curve) and generated (dashed curve) cell counts from our model, we also show the generated cell counts from the MCAR model (dotted curve) for the green cancer cells (G), red cancer cells (R) and fibroblasts (F) across the entire image. Compared to the dotted curves, the dashed curves are slightly closer to the solid ones, which means our model seems more appropriate for analysing this type of data than the MCAR model.

Figure 4: Goodness-of-fit of the estimated models. Observed (solid) and predicted (dashed for our model and dotted for the MCAR model) number of GFP cancer cells (G), mCherry cancer cells (R) cancer cells and fibroblasts (F) for the entire image. Predicted cell counts for each cell type in each tile yˆi,t(c)$\hat{y}_{i,t}^{{(c)}}$ is generated from the conditional Poisson model with intensity λˆi,t(c)$\hat{\lambda}_{i,t}^{{(c)}}$ defined in eqs. (1) and (2), where the coefficients βˆ(c|c′)$\hat{\beta}^{(c|c\prime)}$ are estimated from the entire dataset.

Figure 4:

Goodness-of-fit of the estimated models. Observed (solid) and predicted (dashed for our model and dotted for the MCAR model) number of GFP cancer cells (G), mCherry cancer cells (R) cancer cells and fibroblasts (F) for the entire image. Predicted cell counts for each cell type in each tile yˆi,t(c) is generated from the conditional Poisson model with intensity λˆi,t(c) defined in eqs. (1) and (2), where the coefficients βˆ(c|c) are estimated from the entire dataset.

5 Conclusion and final remarks

In this paper, we introduced a conditional spatial autoregressive model and accompanying inference tools for multivariate spatio-temporal cell count data. The new methodology enables one to measure the overall cell growth rate in longitudinal experiments and spatio-temporal interactions with either homogeneous or heterogeneous cell populations. The proposed inference approach is computationally tractable and strikes a good balance between computational feasibility and statistical accuracy. Numerical findings from simulated and real data in Sections 3 and 4 confirm the validity of the proposed approach in terms of prediction, goodness-of-fit and estimation accuracy.

The data sets described in this paper serve as a proof-of-concept that the proposed methodology works. However, the potential applications and the relevant questions that the methodology can help to answer in cancer cell biology are plentiful. To build on from the examples given in this paper, the methodology can be used to study interactions between cancer cells and a wide range of cancer-relevant cell types such as cancer-activated fibroblasts, macrophages, and other immune cells when co-cultured. Since a substantial proportion of cancer cells in tumors are in close proximity to other cell types that have been shown to affect tumor progression, using these co-cultures is more representative of the situation in a patient compared to studying cancer cells on their own. In addition to just giving the final cell number, the presented approach can dissect which cell types affect the growth of others and to what extent in complex heterogeneous populations. This could be relevant in a drug discovery setting to determine if a drug affects cancer cell growth due to internal effects (on other cancer cells) or by interfering with the interaction between the cancer cells and other cell types. Finding drugs with different targets and mechanisms of action are particularly sought after as they provide a wider target profile, increasing the chance of patients responding as well as reducing the risk of tumors becoming resistant. The impact of different genes and associated pathways in different cell types in relation to inter-cellular interactions can also be studied by genetically modifying the cell type(s) in question before mixing the cells together. This could be beneficial to identify new potential drug targets. Our approach is also applicable in other kinds of studies where local spatial cell-cell interactions are believed to affect cell growth such as studies of neurodegenerative diseases [37] and wound healing/tissue re-generation [38]. In addition to evaluating cell growth, our approach can also be used to study transitions between cellular phenotypes upon interaction with other cell types, provided that the different phenotypes studied can be distinguished from one another based on the image data. Finally, it is worth noting that issues may arise when cells become too confluent/dense, this may lead to segmentation problems of the imaging system. If they become completely confluent, they are likely to progressively stop growing. If one wants to measure for longer period of time, experiments can be performed in larger wells/plates or with smaller starting cell numbers.

Our method offers several practical advantages to researchers interested in analysing multivariate count data on heterogeneous cell populations. First, the conditional Poisson model does not require tracking individual cells across time, a process that is often difficult to automate due to cell movement, morphology changes at subsequent time points, and additional complications related to storage of large data files. Second, we are able to quantify local spatio-temporal interactions between different cell populations from a very simple experimental set-up where the different cell populations are grown together in a single experimental condition (co-culture). An alternative, solely experimentally-based strategy would require monitoring the different cell types alone and together at different cell densities (number of cells per condition) in order to make inferences in terms of potential interactions. However, such an approach would give no possibility of evaluating the spatial relations in the co-culture conditions and would still restrict the number of simultaneously tested cell types to two.

In the future, we foresee several useful extensions of the current methodology, possibly enabling the treatment of more complex experimental settings. First, complex experiments involving a large number of cell populations, nC, would imply an over-parametrized model. Clearly, this large number of parameters would be detrimental to both statistical accuracy and reliable optimization of the likelihood objective function n(θ) in (4). To address these issues, we plan to explore a penalized likelihood of form n(θ)penλ(θ), where pen(θ) is a nonnegative sparsity-inducing penalty function. For example, in a different likelihood setting, Bardic et al. [39] consider the L1-type penalty pen(θ)=λ|θ|,λ>0.

Second, for certain experiments, it would be desirable to modify the statistics in eq. (2) to include additional information on cell growth such as the distance between heterogeneous cells, and covariates describing cell morphology.

Thirdly, it would be useful to develop a more principled way to select the tile sizes/number, and consider tiling the microscope image into a hexagonal lattice, which is a more natural choice in real application, since the distance between neighboring tiles would be more even than that of a regular lattice.

Finally, although numerical results (results not reported here) show that our method are quite robust in the presence of mild outliers (with around 5% of contaminated data), for more severe situations, we expect that severe or numerous outliers will have some influence on the estimates since the Poisson score function is unbounded. To address this problem, the log-likelihood scores in eq. magenta (5) should be replaced by some other robust alternative. Following Ferrari and Vecchia [29] and La Vecchia et al. [30], robustness can be obtained by the so-called q-entropy estimation method simply obtained by replacing the usual logarithm in the log-likelihood estimating equation by the q-logarithm logarithm function logq(x)=(x1q1)/(1q)if q1, and logq(u)=log(x) if q=1, for all x>0. This ensures a bounded influence function for the implied estimator and therefore guarantees control of the bias under contamination.

Funding statement: This article was supported by the Australian National Health and Medical Research Council (1049561, 1064987 and 1069024)


The authors wish to acknowledge support from the Australian National Health and Medical Research Council grants 1049561, 1064987 and 1069024 to Frédéric Hollande. Christina Mølck is supported by the Danish Cancer Society.


In the first part of this section, we provide technical lemmas required to prove asymptotic properties of the estimator θˆn.

Denote Et[] as the expectation with respect to Yt={Yi,t,iLn}, and E[] as the expectation of Y={Yt,t=1,,T}. Let Ni,r be the set of tiles in the neighborhood of tile i, with radius r. Specifically, for two locations i and j, we say jNi,r if ijr. Thus, the neighborhood defined in Section 2 is of radius 1, i.e. {j:ji}={j:jNi,1}. Denote nr=maxiLn|Ni,r|=r2+r+1. Actually, for any tile i that is not on the boundary of the image, |Ni,r|=nr.

In the remainder of this paper we use the following assumptions:

  1. A.1: The parameter space Θ is a compact subset of Rp, and that θ0 is the unique maximiser of (θ)=limnLn(θ).

  2. A.2: The (nC+1)×nLT matrix (v1,1,v1,2,,v1,T,v2,1,,vn,T) is full rank.

Lemma 1.

Let Y1,,Yn be independent Poisson random variables with mean λ1,,λn respectively, where N is a finite positive integer. Then for any positive integer h,




Lemma 2

Denote Y˜Ni,r,t=maxjNi,r,cCYj,t(c), with corresponding observation y˜Ni,r,t and conditional mean λ˜Ni,r,t, then




the {} denotes Stirling number of the second kind, α˜=maxcCα(c),B=maxc(cCβ(c|c))n1.



Similarly, for any cC, we have λNi,r,t(c)eα˜y˜Ni,r+1,t1+1B, since {jNj,1;jNi,r,iLn,r>0}={jNi,r+1;iLn,r>0}.

Next, we proceed by induction.

For T=1, by the conditional independence assumption and Lemma 1, we have


Since T1=0 and Yt has constant entries at time point 0, ET11+Y˜Ni,r+1,T1Bk=1+y˜Ni,r+1,0Bk.

Suppose eq. (8) is true for T=t, then for T=t+1, we have


Lemma 3.

Given Assumption A.1, for any finite constant a,b0 and θΘ,E(λi,t(c)aSi,t1(c)b)<,c,cC,iLn,t=1,,T.


By the definition of ft(k) given in Lemma 2, we know that ft(k) is bounded for all bounded t under assumption A.1. Thus, Lemma 2 implies


For simplicity, define the distance between tile i and j as d(i,j)=r if r1<ijr.

Lemma 4.

For anyiLn,t1=1,,T,

Cov(Yi,t1,Yj,t2)=0,for jLn,t2=1,,T,if d(i,j)>t1+t2.




Let Ni,t={j:Cov(Yj,0,Yi,t)0;jL} be the collection of counts in tiles at time 0 that are correlated with the count in tile i at time t (Yi,t). Due to the neighborhood structure in the autoregressive term described in Section 2, one can easily tell that Ni,t is a neighbourhood around tile i, with the radius equal to t.

Due to the condition that Yt has constant entries at time 0, we have Cov(Yi,t1,Yj,t2)=0 if Ni,t1Nj,t2=, which is true when d(i,j)>t1+t2.

For any (i,t1)Dn, {(j,t2):Ni,t1Nj,t2} is a neighborhood around tile i, with a radius t1+t2.

Since nr=2r2+2r+1, we have


In the second part of this section, we study the asymptotic properties of the estimator θˆn.

Proposition 1

(Existence and uniqueness) If assumption A.3 holds, then there exist unique maximizer of n(θ), denoted by θˆn.


First, since Θ is compact and n(θ) is continuous, at least one maximiser of n(θ) exist. Next, we wish to prove that the maximiser is unique. The p×p Hessian matrix of n(θ) can be written as a block matrix


where Hn(c)(θ)=iLnt=1Texpvi,t(c)(θ)vi,tvi,t is a (nC+1)×(nC+1) matrix. Matrix vi,tvi,t is positive semidefinite with rank 1. By Assumption A.2, iLnt=1Tvi,tvi,t is full rank, which means Hn(c)(θ) is positive definite for all cC and θΘ, since expvi,t(c)(θ)>0. This shows that n(θ) is strictly convex, which implies θˆn is unique.

Proposition 2

[Consistency] If the regularity assumption A.1 holds, then θˆnpθ0 with probability tending 1, as nL.


We proceed by verifying the conditions of Theorem 2 in [31]. First we show that the score functions are Lp-Uniform Integrable for p<3, i.e.

(10)limnsupiLnt=1,,TsupθΘE[ui,tp(θ)I(ui,t(θ)>k)]0,as k.

The general form of each entry of ui,t(θ) is (λi,t(c)yi,t(c))Si,t1c, take p=3, we have


which is finite by lemma 3. This gives us the L3boundedness of ui,t(θ), i.e.


which implies Lp-Uniform Integrability, for p<3.

Second, we show the stochastic equicontinuity of ui,t(y;θ), i.e.


The ui,t(θ) is a p×p matrix, with each column being either γi,t(θ)β(c|c)vi,t or γi,t(θ)α(c)vi,t, and


Thus, the non-zero entries of EsupθΘui,t(θ) have the general form: EsupθΘλi,t(c)Si,t(c)Si,t(c), which are bounded by an equivalent analogous to Lemma 3.

Thirdly, we check αmixing conditions. Let U and V be two subsets of Dn, and let σ(U)=σ{Yi,t;(i,t)U} be the σalgebra generated by random variables Yi,t,(i,t)U.



Then the αmixing coefficient for the random field {Yi,t,iLn,t=1,,T} is defined as


Following Bai et al. [40], in an adimensional space, we need (a) δ>0s.t.m=1ma1α(1,1,m)δ/(2+δ)<, (b) For k+l4,m=1ma1α(k,l,m)<, (c) ϵ>0s.t.α(1,,m)=O(maϵ), where k,l,mN and d(U,V)=min{ij:iU,jV} is the distance between sets U and V.

For any fixed i1,,ikLn,k< and t1=0,,T,

consider U={Yi,t1=yi,t1,,Yik,t1=yik,t1} and V={Yj,t2=yj,t2;jLn,t2=0,,T}, then |U|=k and |V| as n. By Lemma 4, we have P(Yi,t1=yi,t1,Yj,t2=yi,t1)P(Yi,t1=yi,t1)P(Yj,t2=yj,t2)=0, if d(i,j)>t1+t2. Thus, α(U,V)=0 for any |U|=k, provided that d(U,V)>2T, that is, α(k,,m)=0 if m>2T.

This implies all three mixing conditions.

Finally, by Theorem 3 in Jenish and Prucha [31], Uniform Integrability in eq. (10) and mixing condition (a) ensure that the score functions ui,t(y;θ) satisfy a point wise law of large numbers in the sense that

1nLiLnt=1TsupθΘ(ui,t(y,θ)Eui,t(y;θ))p0,as nL,

for all θΘ.

Proposition 3.

If the regularity assumptions A.1 and A.2 hold, we havenLVn(θ)1/2(θˆnθ0)converges in distribution to apvariate Normal with zero mean vector and identity variance, asnL.


First, we show the uniform law of large numbers for un(θ):

(11)supθun(θ)Eun(θ)p0,as nL,

where un(θ)=n(θ)/nL as defined in Section 2. Note that


The first term in eq. (12) is O(nL1), since Varui,t(θ)Eui,t(θ)2, which is shown to be finite in the proof of Proposition 2.

For the second term in eq. (12), by Lemma 2 we have


where Covui,t1(θ),uj,t2(θ)Eui,t1(θ),uj,t2(θ)Eui,t1(θ)2+Eui,t2(θ)2 is finite by Lemma 2. Thus, the second term in eq. (12) is also of order O(nL1) element wise, which means Varun(θ)0 as n. Therefore, eq. (11) follows by Chebyshev’s inequality.

Second, Vn(θ)=1/nLVariLnt=1Tuit(θ)=1/nLEiLnt=1Tuit(θ)=1/nLHn(θ), which is shown to be positive definite under Assumption A.2 in Proposition 1. Thus, together with uniform Integrability in eq. (10) and the mixing conditions, by Theorem 1 in [31], we have


Finally, by Taylor’s expansion,


where θ˜n is a vector with elements between θˆn and θ0. Since θˆn=θ0+op(1) by Proposition 2, we have (θ˜nθ0)2=(θˆnθ0)op(1). The second derivative 2un(θ) is a p×p×p matrix, with entries being either 0 or λit(c)Si,t1(c1)Si,t1(c2)Si,t1(c3), where i=1,,n, and t=1,,T, and c,c1,c2,c3C. Due to the structure of λit(c) and Si,t1(c) in Section 2, all non-zero elements in 2un(θ) are monotone with respect to θ. Thus, there exists θsΘ such that 2un(θs)2un(θ) for all θΘ. Therefore, we have EsupθΘ2un(θ)=supθΘE2un(θ),

which can be shown to be finite by an equivalent analogous to Lemma 3.

Thus, eq. (13) can be written as


By eq. (11), un(θ0)pEun(θ0)=Vn(θ0), since n(θ) is the full likelihood. Therefore, by eqs. (12) and (13), we have



[1] Medema JP, Vermeulen L. Microenvironmental regulation of stem cells in intestinal homeostasis and cancer. Nature. 2011;474:318–326.2167774810.1038/nature10212 Search in Google Scholar

[2] Besag J. Spatial interaction and the statistical analysis of lattice systems. J Royal Stat Soci Series B Methodol. 1974;192–236. Search in Google Scholar

[3] Waller LA, Carlin BP, Xia H, Gelfand AE. Hierarchical spatio-temporal mapping of disease rates. J Am Stat Assoc. 1997;92:607–617.10.1080/01621459.1997.10474012 Search in Google Scholar

[4] Knorr-Held L. Bayesian modelling of inseparable space-time variation in disease risk. 1999. Search in Google Scholar

[5] Quick H, Waller LA, Casper M. A multivariate space–time model for analysing county level heart disease death rates by race and sex. J R Stat Soc: Ser C. Appl Stat. 2017. Search in Google Scholar

[6] Sans Ó, Schmidt AM, Nobre AA, et al. Bayesian spatio-temporal models based on discrete convolutions. Can J Stat. 2008;36:239–258.10.1002/cjs.5550360205 Search in Google Scholar

[7] Quick H, Waller LA, Casper M. Hierarchical multivariate space-time methods for modeling counts with an application to stroke mortality data. arXiv preprint arXiv:1602.04528. 2016. Search in Google Scholar

[8] Cressie N, Wikle CK. Statistics for spatio-temporal data. John Wiley & Sons, 2011. Search in Google Scholar

[9] Cox DR, Gudmundsson G, Lindgren G, Bondesson L, Harsaae E, Laake P, Juselius K, Lauritzen SL. Statistical analysis of time series: Some recent developments [with discussion and reply]. Scand J Stat. 1981;93–115. Search in Google Scholar

[10] Bradley JR, Holan SH, Wikle CK. Multivariate spatio-temporal models for high-dimensional areal data with application to longitudinal employer-household dynamics. Ann Appl Stat. 2015;9:1761–1791.10.1214/15-AOAS862 Search in Google Scholar

[11] Bradley JR, Holan SH, Wikle CK. Multivariate spatio-temporal survey fusion with application to the american community survey and local area unemployment statistics. Stat. 2016;5:224–233.10.1002/sta4.120 Search in Google Scholar

[12] Shaddick G, Wakefield J. Modelling daily multivariate pollutant data at multiple sites. J R Stat Soc: Ser C. Appl Stat. 2002;51:351–372.10.1111/1467-9876.00273 Search in Google Scholar

[13] Wikle CK, Berliner LM, Cressie N. Hierarchical bayesian space-time models. Environ Ecol Stat. 1998;5:117–154.10.1023/A:1009662704779 Search in Google Scholar

[14] Holan S, Wikle C. Hierarchical dynamic generalized linear mixed models for discrete-valued spatio-temporal data. Handbook of Discrete–Valued Time Series, 2015. Search in Google Scholar

[15] Mugglin AS, Cressie N, Gemmell I. Hierarchical statistical modelling of influenza epidemic dynamics in space and time. Stat Med. 2002;21:2703–2721.1222888610.1002/sim.1217 Search in Google Scholar

[16] Bradley JR, Holan SH, Wikle CK. Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data. Bayesian Anal. 2017. Search in Google Scholar

[17] Davis RA, Dunsmuir WT, Streett SB. Observation-driven models for poisson counts. Biometrika. 2003;90:777–790.10.1093/biomet/90.4.777 Search in Google Scholar

[18] Schrödle B, Held L, Rue H. Assessing the impact of a movement network on the spatiotemporal spread of infectious diseases. Biometrics. 2012;68:736–744.2217162610.1111/j.1541-0420.2011.01717.x Search in Google Scholar

[19] Paul M, Held L, Toschke AM. Multivariate modelling of infectious disease surveillance data. Stat Med. 2008;27:6250–6267.10.1002/sim.344018800337 Search in Google Scholar

[20] Zeger SL, Qaqish B. Markov regression models for time series: a quasi-likelihood approach. Biometrics. 1988;1019–1031.3148334 Search in Google Scholar

[21] Fokianos K, Tjøstheim D. Log-linear poisson autoregression. J Multivariate Anal. 2011;102:563–578.10.1016/j.jmva.2010.11.002 Search in Google Scholar

[22] Fokianos K, Rahbek A, Tjøstheim D. Poisson autoregression. J Am Stat Assoc. 2009;104:1430–1439.10.1198/jasa.2009.tm08270 Search in Google Scholar

[23] Dunsmuir WT, Scott DJ, et al. The glarma package for observation driven time series regression of counts. J Stat Softw. 2015;67:1–36. Search in Google Scholar

[24] Kedem B, Fokianos K. Regression models for time series analysis, vol. 488. John Wiley & Sons, 2005. Search in Google Scholar

[25] Held L, Höhle M, Hofmann M. A statistical framework for the analysis of multivariate infectious disease surveillance counts. Stat Modell. 2005;5:187–199.10.1191/1471082X05st098oa Search in Google Scholar

[26] Paul M, Held L. Predictive assessment of a non-linear random effects model for multivariate time series of infectious disease counts. Stat Med. 2011;30:1118–1136.21484849 Search in Google Scholar

[27] Knorr-Held L, Richardson S. A hierarchical model for space–time surveillance data on meningococcal disease incidence. J R Stat Soc: Ser C. Appl Stat. 2003;52:169–183.10.1111/1467-9876.00396 Search in Google Scholar

[28] Wikle CK, Anderson CJ. Climatological analysis of tornado report counts using a hierarchical bayesian spatiotemporal model. J Geophys Res Atmos. 2003;108. Search in Google Scholar

[29] Ferrari D, Vecchia. On robust estimation via pseudo-additive information. Biometrika. 2011;99:238–244. Search in Google Scholar

[30] La Vecchia D, Camponovo L, Ferrari D. Robust heart rate variability analysis by generalized entropy minimization. Comput Stat Data Anal. 2015;82:137–151.10.1016/j.csda.2014.09.001 Search in Google Scholar

[31] Jenish N, Prucha IR. Central limit theorems and uniform laws of large numbers for arrays of random fields. J Econom. 2009;150:86–98.10.1016/j.jeconom.2009.02.00920161289 Search in Google Scholar

[32] Leroux BG, Lei X, Breslow N. Estimation of disease rates in small areas: a new mixed model for spatial dependence. In: Statistical models in epidemiology, the environment clinical trials, 179–191. Springer, 2000. Search in Google Scholar

[33] Lee D. Carbayes: An r package for bayesian spatial modeling with conditional autoregressive priors. J Stat Softw. 2013;55:1–24. Search in Google Scholar

[34] Kalluri R, Zeisberg M. Fibroblasts in cancer. Nat Rev Cancer. 2006;6:392–401.1657218810.1038/nrc1877 Search in Google Scholar

[35] Tabassum DP, Polyak K. Tumorigenesis: it takes a village. Nat Rev Cancer. 2015;15:473–483.10.1038/nrc397126156638 Search in Google Scholar

[36] Koenker R. Quantile regression. No. 38, Cambridge university press, 2005. Search in Google Scholar

[37] Garden GA, La Spada AR. Intercellular (mis) communication in neurodegenerative disease. Neuron. 2012;73:886–901.10.1016/j.neuron.2012.02.01722405200 Search in Google Scholar

[38] Leoni G, Neumann P, Sumagin R, et al. Wound repair: role of immune–epithelial interactions. Mucosal Immunol. 2015;8:959–968.10.1038/mi.2015.6326174765 Search in Google Scholar

[39] Bradic J, Fan J, Wang W. Penalized composite quasi-likelihood for ultrahigh dimensional variable selection. J R Stat Soc Series B Stat Methodol. 2011;73:325–349.10.1111/j.1467-9868.2010.00764.x21589849 Search in Google Scholar

[40] Bai Y, Song PX, Raghunathan T. Joint composite estimating functions in spatiotemporal models. J R Stat Soc Series B Stat Methodol. 2012;74:799–824.10.1111/j.1467-9868.2012.01035.x Search in Google Scholar

Received: 2017-07-20
Revised: 2018-05-08
Accepted: 2018-06-19
Published Online: 2018-07-07

© 2018 Walter de Gruyter GmbH, Berlin/Boston