Multiple testing in spatial epidemiology: a Bayesian approach

In this work we aim to propose a new approach for preliminary epidemiological studies on Standardized Mortality Ratios (SMR) collected in many spatial regions. A preliminary study on SMRs aims to formulate hypotheses to be investigated via individual epidemiological studies that avoid bias carried on by aggregated analyses. Starting from collecting disease counts and calculating expected disease counts by means of reference population disease rates, in each area an SMR is derived as the MLE under the Poisson assumption on each observation. Such estimators have high standard errors in small areas, i.e. where the expected count is low either because of the low population underlying the area or the rarity of the disease under study. Disease mapping models and other techniques for screening disease rates among the map aiming to detect anomalies and possible high-risk areas have been proposed in literature according to the classic and the Bayesian paradigm. Our proposal is approaching this issue by a decision-oriented method, which focus on multiple testing control, without however leaving the preliminary study perspective that an analysis on SMR indicators is asked to. We implement the control of the FDR, a quantity largely used to address multiple comparisons problems in the eld of microarray data analysis but which is not usually employed in disease mapping. Controlling the FDR means providing an estimate of the FDR for a set of rejected null hypotheses. The small areas issue arises diculties in applying traditional methods for FDR estimation, that are usually based only on the p-values knowledge (Benjamini and Hochberg, 1995; Storey, 2003). Tests evaluated by a traditional p-value provide weak power in small areas, where the expected number of disease cases is small. Moreover tests cannot be assumed as independent when spatial correlation between SMRs is expected, neither they are identical distributed when population underlying the map is heterogeneous. The Bayesian paradigm oers a way to overcome the inappropriateness of p-values based methods. Another peculiarity of the present work is to propose a hierarchical full Bayesian model for FDR estimation in testing many null hypothesis of absence of risk.We will use concepts of Bayesian models for disease mapping, referring in particular to the Besag York and Mollié model (1991) often used in practice for its exible prior assumption on the risks distribution across regions. The borrowing of strength between prior and likelihood typical of a hierarchical Bayesian model takes the advantage of evaluating a singular test (i.e. a test in a singular area) by means of all observations in the map under study, rather than just by means of the singular observation. This allows to improve the power test in small areas and addressing more appropriately the spatial correlation issue that suggests that relative risks are closer in spatially contiguous regions. The proposed model aims to estimate the FDR by means of the MCMC estimated posterior probabilities b i's of the null hypothesis (absence of risk) for each area. An estimate of the expected FDR conditional on data (\FDR) can be calculated in any set of b i's relative to areas declared at high-risk (where thenull hypothesis is rejected) by averaging the b i's themselves. The\FDR can be used to provide an easy decision rule for selecting high-risk areas, i.e. selecting as many as possible areas such that the\FDR is non-lower than a prexed value; we call them\FDR based decision (or selection) rules. The sensitivity and specicity of such rule depend on the accuracy of the FDR estimate, the over-estimation of FDR causing a loss of power and the under-estimation of FDR producing a loss of specicity. Moreover, our model has the interesting feature of still being able to provide an estimate of relative risk values as in the Besag York and Mollié model (1991). A simulation study to evaluate the model performance in FDR estimation accuracy, sensitivity and specificity of the decision rule, and goodness of estimation of relative risks, was set up. We chose a real map from which we generated several spatial scenarios whose counts of disease vary according to the spatial correlation degree, the size areas, the number of areas where the null hypothesis is true and the risk level in the latter areas. In summarizing simulation results we will always consider the FDR estimation in sets constituted by all b i's selected lower than a threshold t. We will show graphs of the\FDR and the true FDR (known by simulation) plotted against a threshold t to assess the FDR estimation. Varying the threshold we can learn which FDR values can be accurately estimated by the practitioner willing to apply the model (by the closeness between\FDR and true FDR). By plotting the calculated sensitivity and specicity (both known by simulation) vs the\FDR we can check the sensitivity and specicity of the corresponding\FDR based decision rules. For investigating the over-smoothing level of relative risk estimates we will compare box-plots of such estimates in high-risk areas (known by simulation), obtained by both our model and the classic Besag York Mollié model. All the summary tools are worked out for all simulated scenarios (in total 54 scenarios). Results show that FDR is well estimated (in the worst case we get an overestimation, hence a conservative FDR control) in small areas, low risk levels and spatially correlated risks scenarios, that are our primary aims. In such scenarios we have good estimates of the FDR for all values less or equal than 0.10. The sensitivity of\FDR based decision rules is generally low but specicity is high. In such scenario the use of\FDR = 0:05 or\FDR = 0:10 based selection rule can be suggested. In cases where the number of true alternative hypotheses (number of true high-risk areas) is small, also FDR = 0:15 values are well estimated, and \FDR = 0:15 based decision rules gains power maintaining an high specicity. On the other hand, in non-small areas and non-small risk level scenarios the FDR is under-estimated unless for very small values of it (much lower than 0.05); this resulting in a loss of specicity of a\FDR = 0:05 based decision rule. In such scenario\FDR = 0:05 or, even worse,\FDR = 0:1 based decision rules cannot be suggested because the true FDR is actually much higher. As regards the relative risk estimation, our model achieves almost the same results of the classic Besag York Molliè model. For this reason, our model is interesting for its ability to perform both the estimation of relative risk values and the FDR control, except for non-small areas and large risk level scenarios. A case of study is nally presented to show how the method can be used in epidemiology.

Rainfall spatial predictions: a two-part model and its assessment

Spatial prediction of hourly rainfall via radar calibration is addressed. The change of support problem (COSP), arising when the spatial supports of different data sources do not coincide, is faced in a non-Gaussian setting; in fact, hourly rainfall in Emilia-Romagna region, in Italy, is characterized by abundance of zero values and right-skeweness of the distribution of positive amounts. Rain gauge direct measurements on sparsely distributed locations and hourly cumulated radar grids are provided by the ARPA-SIMC Emilia-Romagna. We propose a three-stage Bayesian hierarchical model for radar calibration, exploiting rain gauges as reference measure. Rain probability and amounts are modeled via linear relationships with radar in the log scale; spatial correlated Gaussian effects capture the residual information. We employ a probit link for rainfall probability and Gamma distribution for rainfall positive amounts; the two steps are joined via a two-part semicontinuous model. Three model specifications differently addressing COSP are presented; in particular, a stochastic weighting of all radar pixels, driven by a latent Gaussian process defined on the grid, is employed. Estimation is performed via MCMC procedures implemented in C, linked to R software. Communication and evaluation of probabilistic, point and interval predictions is investigated. A non-randomized PIT histogram is proposed for correctly assessing calibration and coverage of two-part semicontinuous models. Predictions obtained with the different model specifications are evaluated via graphical tools (Reliability Plot, Sharpness Histogram, PIT Histogram, Brier Score Plot and Quantile Decomposition Plot), proper scoring rules (Brier Score, Continuous Rank Probability Score) and consistent scoring functions (Root Mean Square Error and Mean Absolute Error addressing the predictive mean and median, respectively). Calibration is reached and the inclusion of neighbouring information slightly improves predictions. All specifications outperform a benchmark model with incorrelated effects, confirming the relevance of spatial correlation for modeling rainfall probability and accumulation.

Spatial variability patterns of some Vertisol properties at a field scale using standardized data

Spatial variability of Vertisol properties is relevant for identifying those zones with physical degradation. In this sense, one has to face the problem of identifying the origin and distribution of spatial variability patterns. The objectives of the present work were (i) to quantify the spatial structure of different physical properties collected from a Vertisol, (ii) to search for potential correlations between different spatial patterns and (iii) to identify relevant components through multivariate spatial analysis. The study was conducted on a Vertisol (Typic Hapludert) dedicated to sugarcane (Saccharum officinarum L.) production during the last sixty years. We used six soil properties collected from a squared grid (225 points) (penetrometer resistance (PR), total porosity, fragmentation dimension (Df), vertical electrical conductivity (ECv), horizontal electrical conductivity (ECh) and soil water content (WC)). All the original data sets were z-transformed before geostatistical analysis. Three different types of semivariogram models were necessary for fitting individual experimental semivariograms. This suggests the different natures of spatial variability patterns. Soil water content rendered the largest nugget effect (C0 = 0.933) while soil total porosity showed the largest range of spatial correlation (A = 43.92 m). The bivariate geostatistical analysis also rendered significant cross-semivariance between different paired soil properties. However, four different semivariogram models were required in that case. This indicates an underlying co-regionalization between different soil properties, which is of interest for delineating management zones within sugarcane fields. Cross-semivariograms showed larger correlation ranges than individual, univariate, semivariograms (A ≥ 29 m). All the findings were supported by multivariate spatial analysis, which showed the influence of soil tillage operations, harvesting machinery and irrigation water distribution on the status of the investigated area.

Ross River virus disease clusters and spatial relationship with mosquito biting exposure in Redland Shire, southern Queensland, Australia

The spatial heterogeneity in the risk of Ross River virus (family Togaviridae, genus Alphavirus, RRV) disease, the most common mosquito-borne disease in Australia, was examined in Redland Shire in southern Queensland, Australia. Disease cases, complaints from residents of intense mosquito biting exposure, and human population data were mapped using a geographic information system. Surface maps of RRV disease age-sex standardized morbidity ratios and mosquito biting complaint morbidity ratios were created. To determine whether there was significant spatial variation in disease and complaint patterns, a spatial scan analysis method was used to test whether the number of cases and complaints was distributed according to underlying population at risk. Several noncontiguous areas in proximity to productive saline water habitats of Aedes vigilax (Skuse), a recognized vector of RRV, had higher than expected numbers of RRV disease cases and complaints. Disease rates in human populations in areas which had high numbers of adult Ae. vigilax in carbon dioxide- and octenol-baited light traps were up to 2.9 times those in areas that rarely had high numbers of mosquitoes. It was estimated that targeted control of adult Ae. vigilax in these high-risk areas could potentially reduce the RRV disease incidence by an average of 13.6%. Spatial correlation was found between RRV disease risk and complaints from residents of mosquito biting. Based on historical patterns of RRV transmission throughout Redland Shire and estimated future human population growth in areas with higher than average RRV disease incidence, it was estimated that RRV incidence rates will increase by 8% between 2001 and 2021. The use of arbitrary administrative areas that ranged in size from 4.6 to 318.3 km2, has the potential to mask any small scale heterogeneity in disease patterns. With the availability of georeferenced data sets and high-resolution imagery, it is becoming more feasible to undertake spatial analyses at relatively small scales.

Do predators aggregate in response to pest density in agroecosystems? Assessing within-field spatial patterns

1. The spatial heterogeneity of predator populations is an important component of ecological theories pertaining to predator-prey dynamics. Most studies within agricultural fields show spatial correlation (positive or negative) between mean predator numbers and prey abundance across a whole field over time but generally ignore the within-field spatial dimension. We used explicit spatial mapping to determine if generalist predators aggregated within a soybean field, the size of these aggregations and if predator aggregation was associated with pest aggregation, plant damage and predation rate. 2. The study was conducted at Gatton in the Lockyer Valley, 90 km west of Brisbane, Australia. Intensive sampling grids were used to investigate within-field spatial patterns. The first row of each grid was located in a lucerne field (10 m from interface) and the remaining rows were in an adjacent soybean field. At each point on the grid the abundance of foliage-dwelling and ground-dwelling pests and predators was measured, predation rates [using sentinel Helicoverpa armigera (Hubner) egg cards] and plant damage were estimated. Eight grids were sampled across two summer cropping seasons (2000/01, 2001/02). 3. Predators exhibited strong spatial patterning with regions of high and low abundance and activity within what are considered to be uniform soybean fields. Ground-dwelling and foliage-dwelling predators were often aggregated in patches approximately 40 m across. 4. Lycosidae (wolf spiders) displayed aggregation and were consistently more abundant within the lucerne, with a decreasing trap catch with distance from the lucrene/soybean interface. This trend was consistent between subsequent grids in a single field and between fields. 5. The large amount of spatial variability in within-field arthropod abundance (pests and predators) and activity (egg predation and plant damage) indicates that whole field averages were misleading. This result has serious implications for sampling of arthropod abundance and pest management decision-making based on scouting data. 6. There was a great deal of temporal change in the significant spatial patterns observed within a field at each sampling time point during a single season. Predator and pest aggregations observed in these fields were generally not stable for the entire season. 7. Predator aggregation did not correlate consistently with pest aggregation, plant damage or predation rate. Spatial patterns in predator abundance were not associated consistently with any single parameter measured. The most consistent positive association was between foliage-dwelling predators and pests (significant in four of seven grids). Inferring associations between predators and prey based on an intensive one-off sampling grid is difficult, due to the temporal variability in the abundance of each group. 8. Synthesis and applications. This study demonstrated that generalist predator populations are rarely distributed randomly and field edges and adjacent crops can have an influence on within-field predator abundance. This must be considered when estimating arthropod (pest and predator) abundance from a set of samples taken at random locations within a field.

Spatial correlations between beta-amyloid (Abeta) deposits and blood vessels in early-onset Alzheimer's disease

In cases of late-onset Alzheimer’s disease (AD), there is a spatial correlation between the classsic ‘cored’ type of Beta-amyloid (Abeta) deposit and the large vertically penetrating arterioles in the cerebral cortex suggesting that blood vessels are involved in the pathogenesis of the classic deposits. In this chapter, the spatial correlations between the diffuse, primitive, and classic Abeta deposits and blood vessels were studied in 10 cases of early-onset AD in the age range 40 – 65 years. Sections of frontal cortex were immunostained with antibodies against Abeta?and with collagen IV to reveal the Abeta deposits and blood vessel profiles. In the early-onset cases as a whole, all types of Abeta? deposit and blood vessel profiles were distributed in clusters. There was a positive spatial correlation between the clusters of the diffuse Abeta deposits and the larger (>10µm) and smaller diameter (<10?m) blood vessel profiles in one and three cases respectively. The primitive and classic Abeta deposits were spatially correlated with larger and smaller blood vessels both in three and four cases respectively. Spatial correlations between the Abeta deposits and blood vessels may be more prevalent in cases expressing amyloid precursor protein (APP) than presenilin 1 (PSEN1) mutations. Apolipoprotein E (Apo E) genotype of the patient did not appear to influence the spatial correlation with blood vessel profiles. The data suggest that the larger diameter blood vessels are less important in the pathogenesis of the classic Abeta deposits in early-onset compared with late-onset AD.

Spatial correlations between the vacuolation, prion protein (PrPsc) deposits and the cerebral blood vessels in sporadic Creutzfeldt-Jakob disease

In the variant form of Creutzfeldt-Jakob disease (vCJD), 'florid' deposits of the protease resistant form of prion protein (PrP(sc)) were aggregated around the cerebral blood vessels suggesting the possibility that prions may spread into the brain via the cerebral microcirculation. The objective of the present study was to determine whether the pathology was spatially related to blood vessels in cases of sporadic CJD (sCJD), a disease without an iatrogenic etiology, and therefore, less likely to be caused by hematogenous spread. Hence, the spatial correlations between the vacuolation ('spongiform change'), PrP(sc) deposits, and the blood vessels were studied in immunolabelled sections of the cerebral cortex and cerebellum in eleven cases of the common M/M1 subtype of sCJD. Both the vacuolation and the PrP(sc) deposits were spatially correlated with the blood vessels; the PrP(sc) deposits being more focally distributed around the vessels than the vacuoles. The frequency of positive spatial correlations was similar in the different gyri of the cerebral cortex, in the upper and lower cortical laminae, and in the molecular layer of the cerebellum. It is hypothesized that the spatial correlation is attributable to factors associated with the blood vessels which promote the aggregation of PrP(sc) to form deposits rather than representing the hematogenous spread of the disease. The aggregated form of PrP(sc) then enhances cell death and may encourages the development of vacuolation in the vicinity of the blood vessels.

Spatial correlations between beta-amyloid (A beta) deposits and blood vessels in familial Alzheimer's disease

In sporadic Alzheimer’s disease (SAD), the classic (‘dense-cored’) ß-amyloid (Aß) deposits are aggregated around the larger blood vessels in the upper laminae of the cerebral cortex. To determine whether a similar relationship exists in familial AD (FAD), the spatial correlations between the diffuse, primitive, and classic ß-amyloid (Aß deposits and blood vessels were studied in ten FAD cases including cases linked to amyloid precursor protein (APP) and presenilin (PSEN) gene mutations and expressing apolipoprotein E (apo E) allele E4. Sections of frontal cortex were immunolabelled with antibodies against Aß and with collagen IV to reveal the Aß deposits and blood vessel profiles. In the FAD cases as a whole, Aßdeposits were distributed in clusters. There was a positive spatial correlation between the clusters of the diffuse Aßdeposits and the larger (>10 µm) and smaller diameter (<10 µm) blood vessels in one and three cases respectively. The primitive Aß deposits were spatially correlated with larger and smaller blood vessels each in four cases and the classic deposits in three and four cases respectively. Apo E genotype of the patient did not influence spatial correlation with blood vessels. Hence, spatial correlations between the classic deposits and larger diameter blood vessels were significantly less frequent in FAD compared with SAD. It was concluded that both Aß deposit morphology and AD subtype determine spatial correlations with blood vessels in AD.

Spatial patterns of the pathological changes in the temporal lobe of patients with neuronal intermediate filament inclusion disease

Neuronal intermediate filament inclusion disease (NIFID) is a new neurodegenerative disease characterized histologically by the presence of neuronal cytoplasmic inclusions (NI) immunopositive for intermediate filament proteins, neuronal loss, swollen achromatic neurons (SN), and gliosis. We studied the spatial patterns of these pathological changes parallel to the pia mater in gyri of the temporal lobe in four cases of NIFID. Both the NI and SN occurred in clusters that were regularly distributed parallel to the pia mater, the cluster sizes of the SN being significantly greater than those of the NI. In a significant proportion of areas studied, there was a spatial correlation between the clusters of NI and those of the SN and with the density of the surviving neurons. In addition, the clusters of surviving neurons were negatively correlated (out of phase) with the clusters of glial cell nuclei. The pattern of clustering of these histological features suggests that there is degeneration of the cortico-cortical projections in NIFID leading to the formation of NI and SN within the same vertical columns of cells. The glial cell reaction may be a response to the loss of neurons rather than to the appearance of the NI or SN.

Spatial correlations between the neuronal inclusions, swollen achromatic neurons, and glial cells in neuronal intermediate filament inclusion disease (NIFID)

Neuronal intermediate filament (IF) inclusion disease (NIFID) is characterized by neuronal loss, neuronal cytoplasmic IF-positive inclusions (NI), swollen neurons (SN), and a glial cell reaction. We studied the spatial correlations between the clusters of NI, SN, and glial cells in four gyri of the temporal lobe (superior temporal gyrus, inferior temporal gyrus, lateral occipitotemporal gyrus, and parahippocampal gyrus) in four cases of NIFID. The densities of histological features (per 50x250 μ sample field) were as follows: NI (mean = 0.41, range 0.28-0.68), SN (mean = 1.41, range 0.47-2.65), glial cell nuclei (mean = 5.21, range 3.63-8.17). The NI and the SN were positively correlated in half of the brain regions examined, the correlations being present at the smallest field size (50x250 μm). The NI were also positively or negatively correlated with the glial cell nuclei in different areas, the negative correlations being present at the smallest field size. Glial cell nuclei were positively or negatively correlated with the SN in different brain areas, mainly at the larger field sizes (400x250 and 800x250 μm). The spatial correlation between the clusters of NI and SN in the cortex suggests their development within the same columns of cells. At first, the glial cell reaction is also confined to these columns but later becomes more generally distributed across the cortex. © Springer-Verlag 2004.

The spatial patterns of Pick bodies, Pick cells and Alzheimer’s disease pathology in Pick’s disease

The spatial patterns of Pick bodies (PB), Pick cells (PC), senile plaques (SP) and neurofibrillary tangles (NFT) were studied in the frontal and temporal lobe in nine cases of Pick’s disease (PD). Pick bodies exhibited clustering in 41/44 (93%) of analyses and clusters of PB were regularly distributed parallel to the tissue boundary in 24/41 (58%) of analyses. Pick cells exhibited clustering with regular periodicity of clusters in 14/16 (88%) analyses, SP in three out of four (75%) analyses and NFT in 21/27 (78%) analyses. The largest clusters of PB were observed in the dentate gyrus and PC in the frontal cortex. In 10/17 (59%) brain areas studied, a positive or negative correlation was observed between the densities of PB and PC. The densities of PB and NFT were not significantly correlated in the majority of brain areas but a negative correlation was observed in 7/29 (24%) brain areas. The data suggest that PB and PC in patients with PD exhibit essentially the same spatial patterns as SP and NFT in Alzheimer’s disease (AD) and Lewy bodies (LB) in dementia with Lewy bodies (DLB). In addition, there was a spatial correlation between the clusters of PB and PC, suggesting a pathogenic relationship between the two lesions. However, in the majority of tissues examined there was no spatial correlation between the clusters of PB and NFT, suggesting that the two lesions develop in association with different populations of neurons.

Density and spatial pattern of the neuropathological changes in the cerebellum in the prion disease variant Creutzfeldt-Jakob disease (vCJD)

The objective of this chapter is to quantify the neuropathology of the cerebellar cortex in cases of the prion disease variant Creutzfeldt-Jakob disease (vCJD). Hence, sequential sections of the cerebellum of 15 cases of vCJD were stained with H/E, or immunolabelled with a monoclonal antibody 12F10 against prion protein (PrP) and studied using quantitative techniques and spatial pattern analysis. A significant loss of Purkinje cells was evident in all cases. Densities of the vacuolation and the protease resistant form of prion protein (PrPSc) in the form of diffuse and florid plaques were greater in the granule cell layer (GL) than the molecular layer (ML). In the ML, vacuoles and PrPSc plaques, occurred in clusters which were regularly distributed along the folia, larger clusters of vacuoles and diffuse plaques being present in the GL. There was a negative spatial correlation between the vacuoles and the surviving Purkinje cells in the ML and a positive spatial correlation between the clusters of vacuoles and the diffuse PrPSc plaques in the ML and GL in five and six cases respectively. A canonical variate analysis (CVA) suggested a negative correlation between the densities of the vacuolation in the GL and the diffuse PrPSc plaques in the ML. The data suggest: 1) all laminae of the cerebellar cortex were affected by the pathology of vCJD, the GL more severely than the ML, 2) the pathology was topographically distributed especially in the Purkinje cell layer and GL, 3) pathological spread may occur in relation to a loop of anatomical projections connecting the cerebellum, thalamus, cerebral cortex, and pons, and 4) there are differences in the pathology of the cerebellum in vCJD compared with the M/M1 subtype of sporadic CJD (sCJD).

The carbon dioxide emissions of firms:a spatial analysis

In order to gain a greater understanding of firms' 'environmental behaviour' this paper explores the factors that influence firms' emissions intensities and provides the first analysis of the determinants of firm level carbon dioxide (CO2) emissions. Focussing on Japan, the paper also examines whether firms' CO2 emissions are influenced by the emissions of neighbouring firms and other possible sources of spatial correlation. Results suggest that size, the capital-labour ratio, R&D expenditure, the extent of exports and concern for public profile are the key determinants of CO2 emissions. Local lobbying pressure, as captured by regional community characteristics, does not appear to play a role, however emissions are found to be spatially correlated. This raises implications for the manner in which the environmental performance of firms is modelled in future. © 2013 Elsevier Inc.

Spatial modeling of data with excessive zeros applied to reindeer pellet-group counts

We analyze a real data set pertaining to reindeer fecal pellet-group counts obtained from a survey conducted in a forest area in northern Sweden. In the data set, over 70% of counts are zeros, and there is high spatial correlation. We use conditionally autoregressive random effects for modeling of spatial correlation in a Poisson generalized linear mixed model (GLMM), quasi-Poisson hierarchical generalized linear model (HGLM), zero-inflated Poisson (ZIP), and hurdle models. The quasi-Poisson HGLM allows for both under- and overdispersion with excessive zeros, while the ZIP and hurdle models allow only for overdispersion. In analyzing the real data set, we see that the quasi-Poisson HGLMs can perform better than the other commonly used models, for example, ordinary Poisson HGLMs, spatial ZIP, and spatial hurdle models, and that the underdispersed Poisson HGLMs with spatial correlation fit the reindeer data best. We develop R codes for fitting these models using a unified algorithm for the HGLMs. Spatial count response with an extremely high proportion of zeros, and underdispersion can be successfully modeled using the quasi-Poisson HGLM with spatial random effects.

Spatial Assessment of Indices for Characterizing the Erosive Force of Rainfall in El Salvador Republic

Soil erosion data in El Salvador Republic are scarce and there is no rainfall erosivity map for this region. Considering that rainfall erosivity is an important guide for planning soil erosion control practices, a spatial assessment of indices for characterizing the erosive force of rainfall in El Salvador Republic was carried out. Using pluviometric records from 25 weather stations, we applied two methods: erosivity index equation and the Fournier index. In all study area, the rainiest period is from May to November. Annual values of erosivity index ranged from 7,196 to 17,856 MJ mm ha(-1) h(-1) year(-1) and the Fournier index ranged from 52.9 to 110.0 mm. The erosivity map showed that the study area can be broadly divided into three major erosion risk zones, and the Fournier index map was divided into four zones. Both methods revealed that the erosive force is severe in all study area and presented significant spatial correlation with each other. The erosive force in the country is concentrated mainly from May to November.