Poisson-Gamma and Spatial-Temporal Models: with Application to Cervical Cancer in Kenya’s Counties

In Africa, Cancer is an emerging health problem where in 2012 new cancer cases were about 847,00 and around 519,00 deaths, three quarters of those deaths occurred in sub-Saharan region. In 2018, cancer was ranked as the third leading cause of deaths in Kenya after infectious and cardiovascular diseases. In 2018 cancer incidences were estimated to be 47,887 new cancer cases and 32,987 deaths. According to data from World Health Organization in 2020, cervical cancer is the second most prevalent cancer among women while breast cancer is the first. In this study, data collected by the Nairobi Cancer Registry (NCR) was used to produce spatial-temporal distribution of the cervical cancer in counties in Kenya. The results showed that counties where data was available among them Embu, Meru, Machakos, Mombasa, Nyeri, Kiambu, Kakamega, Nairobi and Bomet respectively had high risk of cervical cancer. Availability of county-based estimates and spatial-temporal distribution of cervical cancer cases will aide development of targeted county strategies, enhance early detection, promote awareness and implementation of universal coverage of major control interventions which will be crucial in reducing and halting the rising burden of the cancer cases in Kenya. In counties where data was not available the model showed relative risks for cervical cancer disease was minute but it was present, therefore spatial temporal models are very appropriate to estimate relative risks of diseases even when there is a small sample (and possibly without a sample) in a given area by borrowing information from other neighboring regions.


Introduction
Cancer arises when normal cells transforms into tumour cell in various stages from pre-cancerous lesion to a malignant tumour. According to the International Agency for Research on Cancer (IARC), in 2018 the global new cancer cases was estimated to be 18.1 million and approximately 9.6 million deaths. In Africa, Cancer is an emerging health problem where in 2012 new cancer cases were about 847,00 and 519,00 deaths, three quarters of those deaths occurred in sub-Saharan region. In 2018, cancer was ranked as the third leading cause of deaths in Kenya after infectious and cardiovascular diseases.
The annual incidence of cancer in Kenya was estimated to be 47,887 new cancer cases, with an annual mortality 32,987 in 2018 [6].
Among women cervical cancer is the fourth prevalent cancer worldwide [27]. According to data from World Health Organization in 2020, in Kenya breast cancer in the most prevalent followed by cervical cancer among women [1]. Human papillomavirus (HPV) infection which is transmitted through direct contact is the cause of almost all cervical cancers [9]. According to Schiffman et al. [25] sexually transmitted HPV genotypes, notably HPV16 cause virtually all cervical cancers world-wide if not controlled immunologically or by screening. The control strategies for Cervical Cancer in Kenya's Counties cervical cancer includes early screening, vaccination against HPV, treatment of pre-cancerous lesions, diagnosis and treatment of invasive cervical cancer and palliative care [28].
Cervical cancer screening aides in detection of abnormalities which can be treated and pre cancers which may progress into actual cancer thus reducing cervical cancer incidences, deaths and morbidity related to treatment [25].
Spatial temporal components of events associated with time and location can be utilized to reveal aspects related to where and when the events occurred. Cervical cancer is an event associated with time and space (location) therefore analysis can be conducted to determine the trends, spread and patterns to develop ways to halt its spread [20]. This research project presented a spatial temporal approach to analyze spatial dynamics of the cervical cancer cases in Kenya over two years (2015 and 2016) in ten counties namely Bomet, Embu, Kakamega, Kiambu, Machakos, Meru, Mombasa, Nairobi, Nakuru, and Nyeri County.
Generating spatial or spatial-temporal maps by mapping cancer rates may help to identify distribution of cervical cancer in different geographical locations. Identification of areas with high disease burden helps in prioritization of control efforts and interventions leading to change of risk behaviours [12].
The study main aim was to determine cervical cancer hotspots by density analysis, and evaluate the trend of spatial correlation for the distribution of cases in Kenya from counties with available data.

Materials and Methods
The study focused on cervical cancer data collected from different counties for a period of two years (2015 and 2016) by Kenya National Cancer Registry, which is a national Population-Based Cancer Registry (PBCR) that provides high quality cancer surveillance data which includes incidences, mortality, and survival information of cancer patients.
Kenya has area of 582,650 km square and is divided into 47 counties (See Figure 1) and has population of 47.5 Million people as per the census conducted in 2019 by Kenya National Bureau of Statistics (KNBS).

Methodology
According to Chandra [7] Small Area Estimation is a twofold problem. First is how to obtain reliable estimates based on data obtained from small areas where some have empty samples and others with very small samples. The second question is how estimation error is to be assessed. Survey data can be used to estimate indicators of small areas with small samples sizes or empty samples. Direct estimators of indicators in small areas might have large sampling errors, their estimation is improved by introducing regression models which provides a relation between independent and dependent variables of interest. Small area estimation is applied in modelling data from this form of non-planned domains (small areas) [2]. Small Area Estimation methods are divided into area and unit level small area models [22].
Gómez-Rubio et al. [11] noted that, fitting Bayesians models using available in-samples data (samples where data is available) can provide reliable estimates for off-samples areas (where data was not available), area level covariates can be included in calculation of the estimates when available. The estimates maybe less accurate as compared when the survey data is available for each area, but they are reasonable since they have lower bias. To model large-scale spatial patterns and with very sparse data, spatial random effects are included at regional level.
In this study the Bayesian Hierarchical Generalized Linear Mixed Models (BHGLMMs) which are used in small area estimation because of their ability to incorporate multiple levels of model dependencies [8] were considered. Posterior distributions of model parameters are estimated using Bayes method through combination of observed data, prior distributions and the available covariates [13]. Strong spatial autocorrelation is generally exhibited in small area level data According to Lawson [15], introduction of spatially structured random effects in the model residual spatial autocorrelation is accounted. Conditional autoregressive priors first proposed by Besag et al. [4] are applied in modelling of spatially structured random effects. The equations in the paper were written using Math Type software [19].

Generalized Linear Mixed Model
In Generalized Linear Mixed Models (GLMMs), the major assumption is that the response variable comes from an exponential family of the form The exponential family is defined as, ε are error terms.

Poisson-Gamma Model
A typical Poisson model cannot model extra variance, Poisson-gamma (PG) model is also known as Negative Binomial model which incorporates gamma distributed random-effects is used as an alternative. Poisson-gamma (PG) model, has the following two-level formulation: In first level the assumption is that the random variable i y (count data) has Poisson distribution or can be written as y Poisson e θ µ with probability density function: is a vector of covariates and based on equation above the joint probability density function is obtained as follows: The marginal probability density function can be obtained as follows The distribution of equation is negative binomial with mean and variance for i y given as follows: The posterior distribution for i θ is estimated as follows The posterior distribution for i θ is Gamma based on equation 7 or can be written as: The Bayes estimate for i θ , generates the posterior mean and posterior variance a, as follows: Spatial correlation of the data is not taken in to account in this model since it only introduces a spatially-unstructured over-dispersion [21]. Due to the stated disadvantage and inability to include covariates Poisson-Gamma model has been criticized and has shown to be inferior to Conditional Autoregressive (CAR) convolution models which are more complex [16].

Conditional Autoregressive (CAR) Models
CAR models have been extensively used to model spatial data in various areas such as in epidemiology, agriculture, demography, economy and image analysis as model for both latent and observed variables To obtain latent Gaussian models, Gaussian priors are assigned to 0 , (. where Θ is unobserved multivariate Gaussian random variable, whose density π (Θ / φ) is controlled by a vector of hyper parameters Φ [17], may not follow Gaussian distribution [18].
The key components of the Latent Gaussian model are, the likelihood of the data π ( y/ Θ ), the Gaussian density of the random vector Θ, π ( Θ / Φ ) and the prior distribution of the parameter vector π ( Φ ).
The posterior is therefore defined as The posterior marginal's for i x and posterior marginal's for Φ or some j Φ can be obtained by applying integrated Nested Laplace approximation (INLA) [17].

Integrated Nested Laplace Application (INLA) Methodology
This is an appropriate inference based method for approximating the posterior marginal's of the latent Gaussian field ( / ), 1,..., i x y i n π = in three steps. The posterior marginal's of the latent effects Θ are written as The posterior marginal's can be approximated using the Laplace approximation. The first approximation to ( / ) y π Φ using Gaussian distributions is constructed as follows; The marginals of the interest can be computed using numerical integration over a multidimensional grid of values of Φ where the sum is over the values of Φ with area weights ∆ [24]. The first step in INLA computation involves approximating the posterior marginal of Φ by using Laplace approximation in equation (13).
The second step involves computing the Laplace approximation of ( / , ) for selected values of Φ which improves the Gaussian approximation in equation (10) *( , ) An improved version of ( / , ) LA i x y π Φ ɶ known as Simplified Laplace approximation was developed by Rue et al. [23]. It involves a series of expansion of ( / ,

corrects for skewness and location and it is
also less computationally expensive [23]. The third step involves combining steps 1 and 2 using numerical integration in equation 13 [23].
Cervical Cancer Distribution Generalized Linear Mixed Model with spatial and temporal random effects and assuming that the response variable was generated by a Poisson process (count data) was used to model cervical cancer cases.
Poisson-gamma (PG) model, has the following two-level formulation: Where represents observed cervical cancers cases for each county and are the expected cervical cancer cases.
In convolution model, the Poisson regression model is In this study cervical cancer cases i y 's were modelled as shown in equation (17) below.
Where j β 's are coefficients, ij X is vector of the covariates, iii. j β 's are coefficients iv. Correlated random time effects, trend f , to account for time dependence, were modelled via first order random walk [14,15]. The assumption is values in the latest year depend on values of the previous year in a specific county. The correlated temporal random effect, trend f , which has a random walk prior distribution, with precision, 1 ϕ τ was assigned Gamma (1, 0.001) prior. v. The spatial effects ( ) str i f S , estimated at county level where households were located. The spatial effects by county to account for strong spatial autocorrelation, and were modelled via normal conditionally autoregressive priors (CAR) [4]. vi. Where 1, , , i m = … and 1, , , j T = … are counties and years respectively;

Deviance Information Criteria
The Deviance Information Criterion (DIC) [26] similar to Akaike Information Criterion (AIC) and is applied in hierarchical Bayesian models. It is defined for improper priors and provides the effective number of parameters in the Bayesian models.
The deviance is The model with smaller DIC value [26] was used to calculate the relative risks. The calculated relative risks are subsequently mapped to geographical areas to produce spatial Cervical Cancer in Kenya's Counties temporal maps.

Results
Data in this study was analyzed using spatial temporal model R-packages. The packages contain functions for Generalized Linear Mixed Model (GLMM) and INLA methodology.

Spatial Temporal Maps
This section displays the distribution of notified cases, Standard Incidence Ratios (SIR) map, for all counties with notified cervical cancer cases over two years (2015-2016).  Clearly, from Figure 2 in most counties there was greater risk of cervical cancer cases than expected from the standard population since all counties where data was available had a SIR value greater than 1. Bomet=1.59, Embu =7.13, Kakamega =2.02, Kiambu =2.42, Machakos =3.44, Meru=4.82, Mombasa =5.51, Nairobi=1.66, Nakuru =2.26, Nyeri =3.07.
The deep purple areas indicated higher risks (SIR>1) while the light shaded areas are low risk (SIR<1). The highest burden of cervical cancer cases was in Embu, Mombasa, Meru, Machakos and Nyeri counties respectively.

Assessing the Presence of Spatial Autocorrelation
Spatial-autocorrelation in spatial-temporal models is modelled via random effects. Therefore, assessment of presence of spatial autocorrelation was conducted by computing the residuals of a fitted simple Poisson log-linear model. The over dispersion parameter is equal to 31.202 indicating a relatively high overdispersion, the Poisson model has equal mean and variance assumption.
Moran's I statistic was computed and permutation test for each year of the data to check for spatial-autocorrelation in the model residuals.
The null hypothesis is; there is no spatial autocorrelation while the alternative is; there is positive spatial autocorrelation. The estimated Moran's I statistic was 0.0399 and the p-value was 0.2104>0.05, suggesting there was no unexplained spatial autocorrelation in the residuals.

Poisson-gamma Model
A Poisson-gamma (Negative-Binomial) model that takes care of over dispersion model and zero inflated variables was fitted. The dispersion parameter for random effect was 1.8692 which was close to 1 and the Akaike Information Criterion (AIC) was 262.9605. Relative risk estimates for counties where data was available were presented in Table 5. In counties where data was not available the relative risks ranged from 0.01 to 0.06 (Lamu). The relative risks in Table 5 revealed that Embu county had the highest risk of cervical cancer, followed by Mombasa, Meru, Machakos, Nyeri, Kiambu, Kakamega, Nairobi, Nakuru and Bomet county respectively.

Spatio-temporal Modelling
Standardized Incidence Ratios (SIRs) sometimes can be useful, but in areas with rare diseases or with small (possibly empty) samples SIRS may be misleading and not very reliable as measure of risk since expected counts may be low as compared to actual observations. Therefore, estimating disease risks using models which borrow information from neighbouring areas and incorporate covariates information in shrinking or smoothing of extreme values based on small samples is appropriate [10].
The parametric formulation for spatial-temporal models was introduced by Bernardinelli et al. [3], with assumption that the linear predictor can be written as: % is the spatial structured and & unstructured components with a main linear trend β , which represents the overall time effect.
This The assumption of linearity in the i can be realized using a dynamic non parametric formulation for the linear predictor it i i t α υ ν = + + + γ η (21) Here, α , i υ and i ν have the same parameterisation as in (20); however, the term t γ denotes the temporally structured effect, where it is modelled using a first order random walk. In This second model was specified in R-INLA as: formula2 <-y~1+f (Area, model="bym", graph=Kenya.adj)+ f(year, model= "rw1")+ f(Area1, model="iid" )…Model 2 A third model expanding Model 1 to explain the time trend of cervical cancer cases and to include space-time interaction was specified as follows: α υ ν η = + + + γ + δ (22) In The third model R-INLA code was formulated as follows: formula3 <-y~1+f (Area, model="bym", graph=Kenya.adj) + f (year, model= "rw1")+ f(Area. Year, model=iid" )… (Model 3) The DIC values for the three models are presented in Table  6. Despite the added complexity the third model with dynamic time trend and space-time interaction had a smaller DIC value indicating a that the model fitted well to the data and was the most appropriate and the relative risks were obtained from the model. The disease risks are usually presented as relative risks in Poisson models. Posterior relative risk distributions greater than 1 indicates an elevated risk of the disease.  The relative risks in Table 7 indicates that Embu county had the highest risk of cervical cancer, followed by Mombasa, Meru, Machakos, Nyeri, Kiambu, Kakamega, Nakuru, Nairobi and Bomet county respectively.
In Figures 3-4 elevated risks are manifested by values greater than 1 in some parts of the country, and a posterior probabilities greater than 0.8.

Discussion and Conclusion
The results show that counties where data was available among them Embu, Mombasa, Meru, Machakos and Nyeri counties had very high risk of cervical cancer. The national and county institutions can utilize spatial temporal tools to identify various cancer hot spots and improve screening and treatment facilities based on specific cancer case.
In counties where data was not available the model showed relative risks of cervical cancers was not high but the risk was present, therefore spatial temporal models are very appropriate to estimate relative risks of diseases even when there is a small sample (and possibly an empty sample) in a given area by borrowing information from other neighboring regions. Based on the study findings we recommend the counties with high relatives to create awareness, provide screening services and provide vaccines to the groups which are at higher risk of cervical cancer.
Despite success of this study, the biggest impediment in spatial temporal study is non-availability of adequate county data which will provide more insight on the distribution of cervical cancer cases in Kenya. Therefore the National Cancer Registry in collaboration with counties health departments should work closely to enhance cancer data collection. This will facilitate research and inform the appropriate measures to be implemented in mitigation of the increase of cancer cases.