17  Tests for patterns

Last updated

October 4, 2026

A range of techniques, most based on statistical tests, can be used to detect deviations from randomness in space. The choice of the methods depends on the question, the nature of the data and scale of observation. Usually, more than one test is applied for the same or different scales of interest depending on how the data are collected.

The several exploratory or inferential methods can be classified based on the spatial scale and type of data (binary, count, etc.) collected, but mainly if the spatial location of the unit is known (mapped) or not known (sampled). Following Madden et al. (2007b), two major groups can be formed. The first group uses intensively mapped data for which the location [x,y] of the sampling unit is known. Examples of data include binary data (plant is infected or not infected) in planting row, point pattern (spatial arrangements of points in a 2-D space) and quadrat data (grids are superimposed on point pattern data). The second group is the sparsely sampled data in the form of count or proportion (incidence) data for which the location is not known or, if known, not taken into account in the analysis.

Figure 17.1: Classification of spatial data and methods used to study spatial patterns of plant disease based on knowledge of the location of the sampling units and nature of the data.

17.1 Intensively mapped

17.1.1 Binary data

In this situation the individual plants are mapped, meaning that their relative positions to one another are known. It is the case when a census is used to map presence/absence data. The status of each unit (usually a plant) is noted as a binary variable. The plant is either diseased (D or 1) or non-diseased or healthy (H or 0). Several statistical tests can be used to detect a deviation from randomness. The most commonly used tests are runs, doublets and join count.

17.1.1.1 Runs test

A run is defined as a succession of one or more diseased (D) or healthy (H) plants, which are followed and preceded by a plant of the other disease status or no plant at all (Madden 1982). There would be few runs if there is an aggregation of diseased or healthy plants and a large number of runs for a random mixing of diseased and healthy plants.

Let’s create a vector of binary (0 = non-diseased; 1 = diseased) data representing a crop row with 32 plants and assign it to y. For plotting purposes, we make a data frame for more complete information.

library(tidyverse) 
y1 <- c(1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,1,1,1,1,0,0,1,1,1,1,0,0,0)
x1 <- c(1:32) # position of each plant
z1 <- 1
row1 <- data.frame(x1, y1, z1) # create a dataframe

We can then visualize the series using ggplot and count the number of runs as 8, aided by the different colors used to identify a run.

library(grid)
runs2 <- row1 |>
  ggplot(aes(x1, z1, label = x1, color = factor(y1))) +
  geom_point(shape = 15, size = 10) +
  scale_x_continuous(breaks = max(z1)) +
  scale_color_manual(values = c("gray70", "darkred")) +
  geom_text(vjust = 0, nudge_y = 0.5) +
  coord_fixed() +
  ylim(0, 2.5) +
  theme_void() +
   theme(legend.position = "top")+
  labs(color = "Status")
ggsave("imgs/runs2.png", runs2, width = 12, height = 2, bg = "white")
Figure 17.2: Sequence of diseased (dark red) or non-diseased (gray) units (plants). The numbers represent the position of the unit

We can obtain the number of runs and related statistics using the oruns_test() function of the {r4pde}.

library(r4pde)
oruns_test(row1$y1)
Ordinary Runs Test of Data Sequence:
 -------------------------------------
 Total Number of Runs (U): 8
 Expected Number of Runs (EU): 16.75
 Standard Deviation of Runs (sU): 2.74
 Z-score: -3.20
 P-value: 0.0007

 Interpretation:
 Based on the Z-score, the sequence exhibits 'aggregation or clustering'.

17.1.1.2 Join count

In join-count statistics, two adjacent plants may be classified by the type of join that links them: D-D, H-H or H-D. The number of joins of the specified type in the orientation(s) of interest is then counted. The question is whether the observed join-count is large (or small) relative to that expected for a random pattern. The join-count statistics provide a basic measure of spatial autocorrelation. The expected number of join counts can be defined under randomness and the corresponding standard errors can be estimated after constants are calculated based on the number of rows and columns of the matrix.

Let’s use the join_count() function of the {r4pde} package to perform a join count test. The formulations in this function apply only for a rectangular array with no missing values, with the “rook” definition of proximity, and they were all presented in the book The Study of Plant Disease Epidemics, page 261 (Madden et al. 2007a).

Let’s create a series of binary data from left to right and top to bottom. The data is displayed in Fig. 9.13 in page 260 of the book (Madden et al. 2007a). In the example, there are 5 rows and 5 columns. This is specified later to compose the matrix which is the data format for analysis.

m1 <- c(1,0,1,1,0,
       1,1,0,0,0,
       1,0,1,0,0,
       1,0,0,1,0,
       0,1,0,1,1)
matrix1 <- matrix(m1, 5, 5, byrow = TRUE)

We can visualize the two-dimensional array by converting to a raster.

# Convert to raster 
mapS2 <- terra::rast(matrix1)
# Convert to data frame
mapS3 <- terra::as.data.frame(mapS2, xy = TRUE)
mapS3 |>
  ggplot(aes(x, y, label = lyr.1, fill = factor(lyr.1))) +
  geom_tile(color = "white", linewidth = 0.5) +
  theme_void() +
  coord_fixed()+
  labs(fill = "Status") +
  scale_fill_manual(values = c("gray70", "darkred"))+
  theme(legend.position = "top")
Figure 17.3: Visualization of a matrix of presence/absence data representing a disease spatial pattern

The function join_count() calculates spatial statistics for a matrix based on the specifications and calculations shown in the book (Madden et al. 2007a). It identifies patterns of aggregation for values in a binary matrix based on join count statistics. The results determine whether the observed spatial arrangement is aggregated or non-aggregated (random) based on a standard normal distribution test statistic Z-score applied separately for HD or DD sequences. For HD, Z-score lower than - 1.64 (more negative) is taken as a basis for rejection of hypothesis of randomness (P = 0.05). For DD sequences, Z-score greater than 1.64 indicates aggregation.

library(r4pde)
join_count(matrix1)
Join Count Analysis of Spatial Patterns of Plant Diseases:
----------------------------------------------------------
1) 'HD' Sequences:
   - Observed Count: 23
   - Expected Count: 19.97
   - Standard Deviation: 3.17
   - Z-score: 0.80
   - Pattern: not aggregated

2) 'DD' Sequences:
   - Observed Count: 7
   - Expected Count: 9.22
   - Standard Deviation: 4.23
   - Z-score: -0.41
   - Pattern: not aggregated
----------------------------------------------------------

In the example above, the observed count of HD sequence was larger than the expected count, so the Z-score suggests non-aggregation, or randomness. The DD sequence observed count was lower than expected count, confirming non-aggregation. Let’s repeat the procedure using the second array of data shown in the book chapter, for which the result is different. In this case, there is evidence of aggregation of diseased plants, because the observed DD is greater than expected and observed HD is lower than expected. The Z-score for HD is less than -1.64 (P< 0.05) and the Z-score for DD is greater than 1.64 (P < 0.05).

m2 <- c(1,1,1,0,0,
       1,1,1,0,0,
       1,1,1,0,0,
       1,1,1,0,0,
       0,0,0,0,0)
matrix2 <- matrix(m2, 5, 5, byrow = TRUE)

join_count(matrix2)
Join Count Analysis of Spatial Patterns of Plant Diseases:
----------------------------------------------------------
1) 'HD' Sequences:
   - Observed Count: 7
   - Expected Count: 19.97
   - Standard Deviation: 3.17
   - Z-score: -4.24
   - Pattern: aggregated

2) 'DD' Sequences:
   - Observed Count: 17
   - Expected Count: 9.22
   - Standard Deviation: 4.23
   - Z-score: 1.96
   - Pattern: aggregated
----------------------------------------------------------

We can apply these tests for a real example epidemic data provided by the {r4pde} R package. Let’s work with a portion of the intensively mapped data on the occurrence of gummy stem blight (GSB) of watermelon (Café-Filho et al. 2010).

library(r4pde)
gsb <- DidymellaWatermelon
gsb 
# A tibble: 1,344 × 4
     dap NS_col EW_row severity
   <int>  <int>  <int>    <int>
 1    50      1      1        0
 2    50      1      2        0
 3    50      1      3        0
 4    50      1      4        0
 5    50      1      5        0
 6    50      1      6        0
 7    50      1      7        0
 8    50      1      8        0
 9    50      1      9        0
10    50      1     10        0
# ℹ 1,334 more rows

The inspection of the data frame shows four variables where dap is the number of days after planting, NS_col is the north-south direction, EW_row is the east-west direction and severity is the severity in ordinal score (0-4). Let’s produce a map for the 65 and 74 dap, but first we need to create the incidence variable based on severity.

gsb2 <- gsb |>
  filter(dap %in% c(65, 74)) |> 
  mutate(incidence = case_when(severity > 0 ~ 1,
                               TRUE ~ 0))

gsb2 |> 
  ggplot(aes(NS_col, EW_row, fill = factor(incidence)))+
  geom_tile(color = "white") +
  coord_fixed() +
  scale_fill_manual(values = c("gray70", "darkred")) +
  facet_wrap( ~ dap) +
  labs(fill = "Status")+
  theme_void()+
  theme(legend.position = "top")
Figure 17.4: Incidence maps for gummy stem blight of watermelon

Now we should check the number of rows (y) and columns (x) for further preparing the matrix for the join count statistics. The m65 object should be in the matrix format before running the join count test.

max(gsb2$NS_col)
[1] 12
max(gsb2$EW_row)
[1] 16
m65 <- gsb2 |> 
  filter(dap == 65) |> 
  pull(incidence) |> 
  matrix(12, 16, byrow = TRUE)

join_count(m65)
Join Count Analysis of Spatial Patterns of Plant Diseases:
----------------------------------------------------------
1) 'HD' Sequences:
   - Observed Count: 104
   - Expected Count: 108.47
   - Standard Deviation: 13.89
   - Z-score: -0.36
   - Pattern: not aggregated

2) 'DD' Sequences:
   - Observed Count: 15
   - Expected Count: 12.52
   - Standard Deviation: 4.76
   - Z-score: 0.63
   - Pattern: not aggregated
----------------------------------------------------------

We can apply the join count test for 74 dap. The result shows that the pattern is now aggregated, differing from 65 dap.

# Pull the binary sequence of time 2
m74 <- gsb2 |> 
  filter(dap == 74) |> 
  pull(incidence) |> 
  matrix(12, 16, byrow = TRUE)

join_count(m74)
Join Count Analysis of Spatial Patterns of Plant Diseases:
----------------------------------------------------------
1) 'HD' Sequences:
   - Observed Count: 110
   - Expected Count: 155.67
   - Standard Deviation: 11.91
   - Z-score: -3.88
   - Pattern: aggregated

2) 'DD' Sequences:
   - Observed Count: 58
   - Expected Count: 37.12
   - Standard Deviation: 8.85
   - Z-score: 2.42
   - Pattern: aggregated
----------------------------------------------------------

17.1.1.3 Foci analysis

The Analysis of Foci Structure and Dynamics (AFSD), introduced by Nelson (1996) and further expanded by Laranjeira et al. (1998), was used in several studies on citrus diseases in Brazil. In this analysis, the data come from incidence maps where both the diseased and non-diseased trees are mapped in the 2D plane (Jesus Junior and Bassanezi 2004; Laranjeira et al. 2004).

Here is an example of an incidence map with four foci (adapted from Laranjeira et al. (1998)). The data is organized in the wide format where the first column x is the index for the row and each column is the position of the plant within the row. The 0 and 1 represent the non-diseased and diseased plant, respectively.

foci <- tibble::tribble(
           ~x, ~`1`, ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, ~`7`, ~`8`, ~`9`,
           1,   0,   0,   0,   0,   0,   0,   0,   0,   0,
           2,   1,   1,   1,   0,   0,   0,   0,   1,   0,
           3,   1,   1,   1,   0,   0,   0,   1,   1,   1,
           4,   0,   1,   1,   0,   0,   0,   0,   1,   0,
           5,   0,   1,   1,   0,   0,   0,   0,   0,   0,
           6,   0,   0,   0,   1,   0,   0,   0,   0,   0,
           7,   0,   0,   0,   0,   0,   0,   0,   0,   0,
           8,   0,   0,   0,   0,   0,   0,   0,   0,   0,
           9,   0,   0,   0,   0,   0,   1,   0,   1,   0,
          10,   0,   0,   0,   0,   0,   0,   1,   0,   0,
          11,   0,   1,   0,   0,   0,   1,   0,   1,   0,
          12,   0,   0,   0,   0,   0,   0,   0,   0,   0
          )

Since the data frame is in the wide format, we need to reshape it to the long format using pivot_longer function of the tidyr package before plotting using ggplot2 package.

library(tidyr)

foci2 <- foci |> 
  pivot_longer(2:10, names_to = "y", values_to = "i")
foci2
# A tibble: 108 × 3
       x y         i
   <dbl> <chr> <dbl>
 1     1 1         0
 2     1 2         0
 3     1 3         0
 4     1 4         0
 5     1 5         0
 6     1 6         0
 7     1 7         0
 8     1 8         0
 9     1 9         0
10     2 1         1
# ℹ 98 more rows

Now we can make the plot.

library(ggplot2)
foci2 |> 
  ggplot(aes(x, y, fill = factor(i)))+
  geom_tile(color = "black")+
  scale_fill_manual(values = c("grey70", "darkred"))+
  theme_void()+
  coord_fixed()+
  theme(legend.position = "none")
Figure 17.5: Examples of foci of plant diseases - see text for description

In the above plot, the upper left focus is composed of four diseased plants with a pattern of vertical and horizontal proximity to the central unit (or the Rook’s case). The upper right focus, also with four diseased plants denotes a pattern of diagonal proximity to the central unit (or the Bishop’s case). The lower left focus is composed of 11 diseased plants with 5 rows and 4 columns occupied by the focus; the shape index of the focus (SIF) is 1.25 and the compactness index of the focus (CIF) is 0.55. The lower right is a single-unit focus.

In this analysis, several statistics can be summarized, both at the single focus and averaging across all foci in the area, including:

  • Number of foci (NF) and number of single-plant foci (NSF)

  • To compare maps with different number of plants, NF and NSF can be normalized to 1000 plants as NF1000 and NSF1000

  • Number of plants in each focus i (NPFi)

  • Maximum number of rows of the focus i (rfi) and maximum number of columns of the focus i (cfi)

  • Mean shape index of foci (meanSIF = [∑(rfi / cfi)]/NF), where SIF values equal to 1.0 indicate isodiametrical foci; values greater than 1.0 indicate foci with greater length in the direction between the planting rows and values less than 1 indicate foci with greater length in the direction of the planting row.

  • Mean compactness index of foci (meanCIF = [∑(NPFi/rfi*cfi)]/NF), where CIF values close to 1.0 indicate more compact foci, that is, greater aggregation and proximity among all the plants belonging to the focus

We can obtain the above-mentioned foci statistics using the AFSD function of the r4pde package. Let’s calculate for the foci2 dataset already loaded, but first we need to check whether all variables are numeric or integer.

str(foci2) # y was not numeric
tibble [108 × 3] (S3: tbl_df/tbl/data.frame)
 $ x: num [1:108] 1 1 1 1 1 1 1 1 1 2 ...
 $ y: chr [1:108] "1" "2" "3" "4" ...
 $ i: num [1:108] 0 0 0 0 0 0 0 0 0 1 ...
foci2$y <- as.integer(foci2$y) # transform to numeric

library(r4pde)
result_foci <- AFSD(foci2)

The AFSD function returns a list of three data frames. The first is summary statistics of this analysis, together with the disease incidence (DIS_INC), for the data frame in analysis.

knitr::kable(result_foci[[1]])
stats value
NF 4.0000000
NF1000 37.0370370
NSF 1.0000000
NSF1000 9.2592593
DIS_INC 0.2037037
mean_SIF 1.0625000
mean_CIF 0.6652778

The second object in the list is a data frame with statistics at the focus level, including the number of rows and columns occupied by each focus as well as the two indices for each focus: shape and compactness.

knitr::kable(result_foci[[2]])
focus_id size rows cols SIF CIF
1 11 4 5 0.8 0.5500000
2 5 3 3 1.0 0.5555556
3 5 3 3 1.0 0.5555556
4 1 1 1 1.0 1.0000000

The third object is the original data frame amended with the id for each focus which can be plotted and labelled (the focus ID) using the plot_AFSD() function.

foci_data <- result_foci[[3]]
foci_data
# A tibble: 22 × 4
       x     y     i focus_id
   <dbl> <int> <dbl>    <dbl>
 1     2     1     1        1
 2     2     2     1        1
 3     2     3     1        1
 4     2     8     1        2
 5     3     1     1        1
 6     3     2     1        1
 7     3     3     1        1
 8     3     7     1        2
 9     3     8     1        2
10     3     9     1        2
# ℹ 12 more rows

The plot shows the ID for each focus.

plot_AFSD(foci_data)+
  theme_void()

We will now analyze the gummy stem blight dataset again focusing on 65 and 74 dap. We first need to prepare the data for the analysis by getting the x, y and i vectors in the dataframe.

library(r4pde)
gsb65 <- DidymellaWatermelon |> 
  filter(dap %in% c(65)) |> 
  mutate(incidence = case_when(severity > 0 ~ 1,
                               TRUE ~ 0)) |> 
  select(EW_row, NS_col, incidence) |> 
  rename(x = EW_row, y = NS_col, i = incidence)

Now we can run the AFSD function and obtain the statistics.

result_df1 <- AFSD(gsb65)

knitr::kable(result_df1[[1]])
stats value
NF 12.0000000
NF1000 62.5000000
NSF 2.0000000
NSF1000 10.4166667
DIS_INC 0.1875000
mean_SIF 1.0722222
mean_CIF 0.7162037

This analysis is usually applied to multiple maps and the statistics are visually related to the incidence in the area in a scatter plot. Let’s calculate the statistics for the two dap. We can split it by dap before applying the function. We can do it using the map function of the purrr package.

library(purrr)
gsb65_74 <- DidymellaWatermelon |> 
  filter(dap %in% c(65, 74)) |> 
  mutate(incidence = case_when(severity > 0 ~ 1,
                               TRUE ~ 0)) |> 
  select(dap, EW_row, NS_col, incidence) |> 
  rename(dap = dap, x = EW_row, y = NS_col, i = incidence)

# Split the dataframe by 'time'
df_split <- split(gsb65_74, gsb65_74$dap)

# Apply the AFSD function to each split dataframe
results <- map(df_split, AFSD)

We can check the summary results for time 2 and time 3.

time65 <- data.frame(results[[1]][1])
knitr::kable(time65)
stats value
NF 12.0000000
NF1000 62.5000000
NSF 2.0000000
NSF1000 10.4166667
DIS_INC 0.1875000
mean_SIF 1.0722222
mean_CIF 0.7162037 
time74 <- data.frame(results[[2]][1])
knitr::kable(time74)
stats value
NF 10.0000000
NF1000 52.0833333
NSF 2.0000000
NSF1000 10.4166667
DIS_INC 0.3229167
mean_SIF 1.1338095
mean_CIF 0.7525022 
# Plot the results to see the two foci in time 1
plot_AFSD(results[[1]][[3]])+
  theme_void()+
  coord_fixed()

# Plot time 2
plot_AFSD(results[[2]][[3]])+
  theme_void()+
  coord_fixed()

17.1.2 Point pattern analysis

Point pattern analysis studies the arrangement of points in a two-dimensional space. In the simplest case each point is the location of a diseased plant, and the pattern is a map of these locations. A map (a scatter plot with equally scaled axes) gives a first impression, but we need quantitative methods to decide whether the points are random, clustered (closer together than expected by chance) or regular (more evenly spaced than expected by chance). Each of these patterns points to different processes behind the epidemic.

Let’s work with two simulated datasets that were originally generated to produce a random or an aggregated (clustered) pattern.

library(r4pde)
rand <- SpatialRandom
aggr <- SpatialAggregated

Let’s produce a 2-D map for each data frame.

prand <- rand |> 
  ggplot(aes(x, y))+
  geom_point()+
  scale_color_manual(values = c("black", NA))+
  coord_fixed()+
  coord_flip()+
   theme_r4pde(font_size = 12)+
  theme(legend.position = "none")+
  labs (title = "Random", 
        x = "Latitude",
        y = "Longitude",
        caption = "Source: r4pde R package")

paggr <- aggr |> 
  ggplot(aes(x, y))+
  geom_point()+
  scale_color_manual(values = c("black", NA))+
  coord_fixed()+
  coord_flip()+
  theme_r4pde(font_size = 12)+
  theme(legend.position = "none")+
  labs (title = "Aggregated", 
        x = "Latitude",
        y = "Longitude",
        caption = "Source: r4pde R package")

library(patchwork)
prand | paggr

17.1.2.1 Quadrat count

In the Quadrat Count method, the study region is divided into a regular grid of smaller, equally-sized rectangular or square subregions known as quadrats. For each quadrat, the number of points falling inside it is counted. If the points are uniformly and independently distributed across the region (i.e., random), then the number of points in each quadrat should follow a Poisson distribution. If the variance of the counts is roughly equal to the mean of the counts, then the pattern is considered to be random. If the variance is greater than the mean, it suggests that the pattern is aggregated or clumped.

Using the {spatstat} package, we first need to create a ppp object which represents a point pattern. This is the primary object type in spatstat for point patterns.

library(spatstat)

### Create rectangular window around the points
window_rand <- bounding.box.xy(rand$x, rand$y)

# create the point pattern object
ppp_rand <- ppp(rand$x, rand$y, window_rand)
plot(ppp_rand)

Using the quadratcount function, we can divide the study region into a grid and count the number of points in each cell:

## Quadrat count 8 x 8
qq <- quadratcount(ppp_rand,8,8, keepempty = TRUE) 

# plot the quadrat count
plot(qq)

To determine whether the observed distribution of points is consistent with a random Poisson process, we can use the quadrat.test function:

# Quadrat test
qt <- quadrat.test(qq, alternative="clustered", method="M")
qt

    Conditional Monte Carlo test of CSR using quadrat counts
    Test statistic: Pearson X2 statistic

data:  
X2 = 47.04, p-value = 0.9375
alternative hypothesis: clustered

Quadrats: 8 by 8 grid of tiles

The result will give a Chi-squared statistic and a p-value. If the p-value is very low, then the pattern is likely not random. Keep in mind that the choice of the number and size of quadrats can affect the results. It’s often helpful to try a few different configurations to ensure robust conclusions.

Let’s repeat this procedure for the situation of aggregated data.

### Create a rectangular window around the points 
window_aggr <- bounding.box.xy(aggr$x, aggr$y)

# create the point pattern object
ppp_aggr <- ppp(aggr$x, aggr$y, window_aggr)
plot(ppp_aggr)

## Quadrat count 8 x 8
qq_aggr <- quadratcount(ppp_aggr,8,8, keepempty=TRUE) 

# plot the quadrat count
plot(qq_aggr)

# Quadrat test
qt_aggr <- quadrat.test(qq_aggr, alternative="clustered", method="M")
qt_aggr

    Conditional Monte Carlo test of CSR using quadrat counts
    Test statistic: Pearson X2 statistic

data:  
X2 = 244.92, p-value = 5e-04
alternative hypothesis: clustered

Quadrats: 8 by 8 grid of tiles

17.1.2.2 Spatial KS test

The spatial Kolmogorov-Smirnov (KS) test compares the distribution of a spatial covariate at the observed points with the distribution expected if the points were completely spatially random (CSR) (Baddeley et al. 2005). If the points tend to lie where the covariate takes particular values (for example, low elevation or close to a field edge), the test rejects CSR.

The covariate can be any spatially varying quantity that might influence the process, such as elevation, soil quality or distance to a feature. Here we use the coordinates themselves as covariates, so the test asks whether the distribution of the x-coordinates, the y-coordinates, or a combination of both, differs from what is expected under CSR.

Let’s test for the aggregated data.

# y as covariate
ks_y <- cdf.test(ppp_aggr, test="ks", "y", jitter=FALSE)
ks_y

    Spatial Kolmogorov-Smirnov test of CSR in two dimensions

data:  covariate 'y' evaluated at points of 'ppp_aggr' 
     and transformed to uniform distribution under CSR
D = 0.070995, p-value = 0.2453
alternative hypothesis: two-sided
plot(ks_y)

# x as covariate
ks_x <- cdf.test(ppp_aggr, test="ks", "x", jitter=FALSE)
ks_x

    Spatial Kolmogorov-Smirnov test of CSR in two dimensions

data:  covariate 'x' evaluated at points of 'ppp_aggr' 
     and transformed to uniform distribution under CSR
D = 0.093863, p-value = 0.0512
alternative hypothesis: two-sided
plot(ks_x)

# x and y as covariates
fun <- function(x,y){2* x + y}
ks_xy <- cdf.test(ppp_aggr, test="ks", fun, jitter=FALSE)
ks_xy

    Spatial Kolmogorov-Smirnov test of CSR in two dimensions

data:  covariate 'fun' evaluated at points of 'ppp_aggr' 
     and transformed to uniform distribution under CSR
D = 0.12687, p-value = 0.002471
alternative hypothesis: two-sided
plot(ks_xy)

As shown above, we have sufficient evidence to reject the null hypothesis of complete spatial randomness.

17.1.2.3 Distance based

17.1.2.3.1 Ripley’s K

A spatial point process is a set of irregularly distributed locations in a region, assumed to be generated by a stochastic mechanism. Ripley’s K-function, \(K(r)\), describes the process at all distances at once. Assuming a constant intensity of points, \(K(r)\) is the expected number of further points within distance \(r\) of an arbitrary point, divided by the intensity. For diseased trees in a plantation, it answers the question: how many other diseased trees do we expect within distance \(r\) of a diseased tree chosen at random?

To interpret the estimated function:

  1. Random pattern: \(K(r)\) follows the theoretical curve for complete spatial randomness, \(K(r) = \pi r^2\).

  2. Clustered pattern: \(K(r)\) lies above the theoretical curve, because points have more neighbors than expected by chance.

  3. Regular pattern: \(K(r)\) lies below the theoretical curve, because points have fewer neighbors than expected by chance.

Let’s use the Kest function of the spatstat package to obtain K(r).

k_rand <- Kest(ppp_rand)
plot(k_rand)

k_aggr <- Kest(ppp_aggr)
plot(k_aggr)

The envelope function performs simulations and computes envelopes of a summary statistic based on the simulations. The envelope can be used to assess the goodness-of-fit of a point process model to point pattern data (Baddeley et al. 2014). Let’s simulate the envelope and plot the values using ggplot. Because observed K(r) (solid line) lies outside the simulation envelope, aggregation was detected.

ke <- envelope(ppp_aggr, fun = Kest)
Generating 99 simulated realisations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
data.frame(ke) |> 
  ggplot(aes(r, theo))+
  geom_line(linetype =2)+
  geom_line(aes(r, obs))+
  geom_ribbon(aes(ymin = lo, ymax = hi),
              fill = "steelblue", alpha = 0.5)+
  labs(y = "K(r)", x = "r")+
  theme_r4pde(font_size = 16)
Figure 17.6: Estimated K(r) of the aggregated pattern (solid line) and the simulation envelope under complete spatial randomness (shaded area); the dashed line is the theoretical K(r)

mad.test performs the ‘global’ or ‘Maximum Absolute Deviation’ test described by Ripley (1977, 1981). See (Baddeley et al. 2014). This performs hypothesis tests for goodness-of-fit of a point pattern data set to a point process model, based on Monte Carlo simulation from the model.

# Maximum absolute deviation test
mad.test(ppp_aggr, Kest)
Generating 99 simulated realisations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

    Maximum absolute deviation test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: K(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 49.477517555]
    Test statistic: Maximum absolute deviation
    Deviation = observed minus theoretical

data:  ppp_aggr
mad = 2234.7, rank = 1, p-value = 0.01
mad.test(ppp_rand, Kest)
Generating 99 simulated realisations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

    Maximum absolute deviation test of CSR
    Monte Carlo test based on 99 simulations
    Summary function: K(r)
    Reference function: theoretical
    Alternative: two.sided
    Interval of distance values: [0, 44.352774865]
    Test statistic: Maximum absolute deviation
    Deviation = observed minus theoretical

data:  ppp_rand
mad = 106.7, rank = 87, p-value = 0.87
17.1.2.3.2 O-ring statistics

The O-ring statistic estimates the average density of points in a ring of radius \(r\) around the points of a pattern, so it shows at which distances points are more or less frequent than expected (Wiegand and A. Moloney 2004). In the univariate form used here the rings are placed around every point of a single pattern. In the bivariate form the rings are placed around points of one type and the points of a second type are counted.

The statistic is interpreted relative to the intensity of the pattern, which is the value expected under CSR:

  1. Random pattern: \(O(r)\) is flat and close to the intensity at all distances.

  2. Aggregation: \(O(r)\) is above the intensity at short distances and decreases toward it as \(r\) increases, because points are surrounded by more neighbors than expected.

  3. Regularity: \(O(r)\) is below the intensity at short distances and increases toward it as \(r\) increases.

As with \(K(r)\), a simulation envelope under CSR is used for inference: if the observed statistic lies inside the envelope, the pattern is consistent with randomness, and if it lies outside, it is not.

In R, the O-ring statistic is \(O(r) = \lambda \, g(r)\), the intensity \(\lambda\) of the pattern multiplied by the pair correlation function \(g(r)\). The onpoint package provides a fast estimate of \(g(r)\) with pcf_fast() (older versions of the package had a function called estimate_o_ring(), which we define here as a small wrapper). We will use the point pattern objects ppp_rand and ppp_aggr created in the previous examples

library(onpoint)

# O-ring statistic: intensity x pair correlation function
estimate_o_ring <- function(X, ...) {
  intensity(X) * pcf_fast(X, ...)
}

plot(estimate_o_ring(ppp_rand))

plot(estimate_o_ring(ppp_aggr))

The function can be used in combination with spatstat’s envelope() function.

oring_envelope <- envelope(ppp_aggr, fun = estimate_o_ring, nsim = 199, verbose = FALSE)
plot(oring_envelope)

To plot simulation envelopes using quantum plots (Esser et al. 2014), just pass an envelope object as input to plot_quantums().

plot_quantums(oring_envelope, ylab = "O-ring")+
  theme_r4pde()

17.1.3 Grouped data

If the data are intensively mapped, meaning that the spatial locations of the sampling units are known, we are not limited to analyze presence/absence (incidence) only data at the unit level. The sampling units may be quadrats where the total number of plants and the number of disease plants (or number of pathogen propagules) are known. Alternatively, it could be a continuous measure of severity. The question here, similar to the previous section, is whether a plant being diseased makes it more (or less) likely that neighboring plants will be diseased. If that is the case, diseased plants are exhibiting spatial autocorrelation. The most common methods are:

  • Autocorrelation (known as Moran’s I)

  • Semivariance

  • SADIE (an alternative approach to autocorrelation.)

17.1.3.1 Autocorrelation

Spatial autocorrelation analysis provides a quantitative assessment of whether a large value of disease intensity in a sampling unit makes it more (positive autocorrelation) or less (negative autocorrelation) likely that neighboring sampling units tend to have a large value of disease intensity (Madden et al. 2007a).

We will illustrate the basic concept by reproducing the example provided in page 264 of the chapter on spatial analysis (Madden et al. 2007a), which was extracted from table 11.3 of Campbell and Madden. L. (1990). The data represent a single transect with the number of Macrophomina phaseolina propagules per 10 g air-dry soil recorded in 16 contiguous quadrats across a field.

mp <- data.frame(
  i = c(1:16),
  y = c(41, 60, 81, 22, 8, 20, 28, 2, 0, 2, 2, 8, 0, 43, 61, 50)
)
mp
    i  y
1   1 41
2   2 60
3   3 81
4   4 22
5   5  8
6   6 20
7   7 28
8   8  2
9   9  0
10 10  2
11 11  2
12 12  8
13 13  0
14 14 43
15 15 61
16 16 50

We can produce a plot to visualize the number of propagules across the transect.

mp |>
  ggplot(aes(i, y)) +
  geom_col(fill = "darkred") +
  theme_r4pde()+
  labs(
    x = "Relative position within a transect",
    y = "Number of propagules",
    caption = "Source: Campbell and Madden (1990)"
  )
Figure 17.7: Number of propagules of Macrophomina phaseolina in the soil at various positions within a transect

To calculate the autocorrelation coefficient in R, we can use the acf() function of the stats package (part of base R).

ac_mp <- acf(mp$y, lag = 5, pl = FALSE)
ac_mp

Autocorrelations of series 'mp$y', by lag

     0      1      2      3      4      5 
 1.000  0.586  0.126 -0.033 -0.017 -0.181

Let’s store the results in a data frame to facilitate visualization.

ac_mp_dat <- data.frame(index = ac_mp$lag, ac_mp$acf)
ac_mp_dat
  index   ac_mp.acf
1     0  1.00000000
2     1  0.58579374
3     2  0.12636306
4     3 -0.03307249
5     4 -0.01701392
6     5 -0.18092810

And now the plot known as autocorrelogram.

ac_mp_dat |>
  ggplot(aes(index, ac_mp.acf, label = round(ac_mp.acf, 3))) +
  geom_col(fill = "darkred") +
  theme_r4pde()+
  geom_text(vjust = 0, nudge_y = 0.05) +
  scale_x_continuous(n.breaks = 6) +
  geom_hline(yintercept = 0) +
  labs(x = "Distance lag", y = "Autocorrelation coefficient")
Figure 17.8: Autocorrelogram for the spatial distribution of Macrophomina phaseolina in soil

The values we obtained here are not the same but quite close to the values reported in Madden et al. (2007b). For the transect data, the calculated coefficients in the book example for lags 1, 2 and 3 are 0.625, 0.144, and - 0.041. The conclusion is the same, the smaller the distance between sampling units, the stronger is the correlation between the count values.

17.1.3.1.1 Moran’s I

The autocorrelation coefficient above is closely related to Moran’s I (Moran 1950), a widely used statistic of spatial autocorrelation. Positive autocorrelation means that nearby locations have more similar values than expected under a random arrangement, and negative autocorrelation means that they have more dissimilar values.

Let’s use another example dataset from the book to calculate the Moran’s I in R. The data is shown in page 269 of the book. The data represent the number of diseased plants per quadrat (out of a total of 100 plants in each) in 144 quadrats. It was based on an epidemic generated using the stochastic simulator of Xu and Madden (2004). The data is stored in a CSV file.

epi <- read_csv("https://raw.githubusercontent.com/emdelponte/epidemiology-R/main/data/xu-madden-simulated.csv")
epi1 <- epi |>
  pivot_longer(2:13,
               names_to = "y",
               values_to = "n") |>
  pull(n)

The {spdep} package provides the function moran() to calculate Moran’s I. It needs two inputs: a numeric vector with the values (here, the number of diseased plants in each quadrat) and a spatial weights list that defines which quadrats are neighbors. Let’s load the library and create the weights.

set.seed(100)
library(spdep)

The cell2nb() function creates the neighbor list with 12 rows and 12 columns, which is how the 144 quadrats are arranged.

nb <- cell2nb(12, 12, type = "queen", torus = FALSE)

The nb2listw() function supplements a neighbors list with spatial weights for the chosen coding scheme. We use the default W, which is the row standardized (sums over all links to n). We then create the col.W neighbor list.

col.W <- nb2listw(nb, style = "W")

The Moran’s I statistic is given by the moran() function

moran(x = epi1, # numeric vector
      listw = col.W, # the nb list
      n = length(epi1), # number of zones
      S0 = Szero(col.W)) # global sum of weights
$I
[1] 0.6982314

$K
[1] 2.878088

$I is Moran’s I and $K is the sample kurtosis of the values.

A positive Moran’s I indicates positive spatial autocorrelation (nearby quadrats have similar values), a negative value indicates negative autocorrelation (neighboring quadrats tend to be dissimilar), and a value close to zero indicates no spatial autocorrelation.

17.1.3.1.2 Moran’s test

To test whether the autocorrelation is statistically significant, we use moran.test() of the {spdep} package. Its two main inputs are x, the numeric vector of values, and listw, the spatial weights list created above.

moran.test(x = epi1, 
           listw = col.W)

    Moran I test under randomisation

data:  epi1  
weights: col.W    

Moran I statistic standard deviate = 15.919, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.698231416      -0.006993007       0.001962596

The function returns Moran’s I, its expected value under the null hypothesis of no spatial autocorrelation, its variance and a P-value. A significant positive value indicates that similar values cluster together, a significant negative value indicates that dissimilar values are neighbors, and a non-significant value gives no evidence of spatial autocorrelation.

As before, we can construct a correlogram using the output of the sp.correlogram() function. Note that the figure below is very similar to the one shown in Figure 91.5 in page 269 of the book chapter (Madden et al. 2007a). Let’s store the results in a data frame.

correl_I <- sp.correlogram(nb, epi1, 
                           order = 10,
                           method = "I",  
                           zero.policy = TRUE)
Warning in nblag(neighbours, maxlag = order): lag 7 neighbour object has 5
sub-graphs
Warning in nblag(neighbours, maxlag = order): lag 8 neighbour object has 17
sub-graphs
Warning in nblag(neighbours, maxlag = order): lag 9 neighbour object has 37
sub-graphs
Warning in nblag(neighbours, maxlag = order): lag 10 neighbour object has 65
sub-graphs
df_correl <- data.frame(correl_I$res) |> 
  mutate(lag = c(1:10))
# Show the spatial autocorrelation for 10 distance lags
round(df_correl$X1,3)
 [1]  0.698  0.340  0.086 -0.002 -0.009 -0.024 -0.090 -0.180 -0.217 -0.124

Then, we can generate the plot using ggplot.

df_correl |>
  ggplot(aes(lag, X1)) +
  geom_col(fill = "darkred") +
  theme_r4pde()+
  scale_x_continuous(n.breaks = 10) +
  labs(x = "Distance lag", y = "Spatial autocorrelation")
Figure 17.9: Autocorrelogram for the spatial distribution of simulated epidemics

17.1.3.2 Semivariance

Semi-variance is a key quantity in geostatistics. This differs from spatial autocorrelation because distances are usually measured in discrete spatial lags. The semi-variance can be defined as half the variance of the differences between all possible points spaced a constant distance apart.

The semi-variance at a distance d = 0 will be zero, because there are no differences between points that are compared to themselves. However, as points are compared to increasingly distant points, the semi-variance increases. At some distance, called the Range, the semi-variance will become approximately equal to the variance of the whole surface itself. This is the greatest distance over which the value at a point on the surface is related to the value at another point. In fact, when the distance between two sampling units is small, the sampling units are close together and, usually, variability is low. As the distance increases, so (usually) does the variability.

Results of semi-variance analysis are normally presented as a graphical plot of semi-variance against distance, which is referred to as a semi-variogram. The main characteristics of the semi-variogram of interest are the nugget, the range and the sill, and their estimations are usually based on an appropriate (non-linear) model fitted to the data points representing the semi-variogram.

For the semi-variance, we will use the variog() function of the geoR package. We need the data in the long format (x, y and z). Let’s reshape the data to the long format and store it in epi2 dataframe.

epi2 <- epi |>
  pivot_longer(2:13,
               names_to = "y",
               values_to = "n") |>
  mutate(y = as.numeric(y))

head(epi2)
# A tibble: 6 × 3
      x     y     n
  <dbl> <dbl> <dbl>
1     1     1     2
2     1     2     2
3     1     3     3
4     1     4    33
5     1     5     4
6     1     6     0
library(geoR)
# the coordinates are x and y and the data is the n
v1 <- variog(coords = epi2[,1:2], data = epi2[,3])
variog: computing omnidirectional variogram
v2 <- variofit(v1, ini.cov.pars = c(1200, 12), 
               cov.model = "exponential", 
               fix.nugget = F)
variofit: covariance model used is exponential 
variofit: weights used: npairs 
variofit: minimisation function used: optim 
# Plotting 
plot(v1, xlim = c(0,15))
lines(v2, lty = 1, lwd = 2)
Figure 17.10: Semivariance plot for the spatial distribution simulated epidemic

17.1.3.3 SADIE

SADIE (spatial analysis by distance indices) is an alternative to autocorrelation and semi-variance methods described previously, which has found use in plant pathology (Madden et al. 2007a; Xu and Madden 2004; Li et al. 2011). Similar to those methods, the spatial coordinates for the disease intensity (count of diseased individuals) or pathogen propagules values should be provided.

SADIE quantifies spatial pattern by calculating the minimum total distance to regularity. That is, the distance that individuals must be moved from the starting point defined by the observed counts to the end point at which there is the same number of individuals in each sampling unit. Therefore, if the data are highly aggregated, the distance to regularity will be large, but if the data are close to regular to start with, the distance to regularity will be smaller.

The null hypothesis to test is that the observed pattern is random. SADIE calculates an index of aggregation (Ia). When this is equal to 1, the pattern is random. If this is greater than 1, the pattern is aggregated. Hypothesis testing is based on the randomization procedure. The null hypothesis of randomness, with an alternative hypothesis of aggregation.

An extension was made to quantify the contribution of each sampling unit count to the observed pattern. Regions with large counts are defined as patches and regions with small counts are defined as gaps. For each sampling unit, a clustering index is calculated and can be mapped.

In R, we can use the sadie() function of the epiphy package (Gigot 2018). The function computes the different indices and probabilities based on the distance to regularity for the observed spatial pattern and a specified number of random permutations of this pattern. To run the analysis, the dataframe should have only three columns: the first two must be the x and y coordinates and the third one the observations. Let’s continue working with the simulated epidemic dataset named epi2. We can map the original data as follows:

epi2 |>
  ggplot(aes(x, y, label = n, fill = n)) +
  geom_tile() +
  geom_text(size = 5, color = "white") +
  theme_void() +
  coord_fixed() +
  scale_fill_gradient(low = "gray70", high = "darkred")
Figure 17.11: Spatial map for the number of diseased plants per quadrat (n = 144) in simulated epidemic 
library(epiphy)
sadie_epi2 <- sadie(epi2)
Computation of Perry's indices:
sadie_epi2
Spatial Analysis by Distance IndicEs (sadie)

Call:
sadie.data.frame(data = epi2)

Ia: 2.4622 (Pa = < 2.22e-16)

The simple output shows the Ia value and associated P-value. As suggested by the low value of the P-value, the pattern is highly aggregated. The summary() function provides more complete information such as the overall inflow and outflow measures. A data frame with the clustering index for each sampling unit is also provided using the summary() function.

summary(sadie_epi2)

Call:
sadie.data.frame(data = epi2)

First 6 rows of clustering indices:
  x y  i cost_flows      idx_P idx_LMX prob
1 1 1  2 -11.382725 -7.2242617      NA   NA
2 1 2  2  -9.461212 -6.2258877      NA   NA
3 1 3  3  -7.299482 -5.3390880      NA   NA
4 1 4 33   1.000000  0.8708407      NA   NA
5 1 5  4  -5.830952 -3.6534511      NA   NA
6 1 6  0  -5.301329 -2.9627172      NA   NA

Summary indices:
                      overall    inflow  outflow
Perry's index        2.495346 -2.811023 2.393399
Li-Madden-Xu's index       NA        NA       NA

Main outputs:
Ia: 2.4622 (Pa = < 2.22e-16)

'Total cost': 201.6062
Number of permutations: 100

The plot() function allows us to map the clustering indices and so to identify regions of patches (red, outflow) and gaps (blue, inflow).

plot(sadie_epi2) 
Warning: `aes_()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`
ℹ The deprecated feature was likely used in the epiphy package.
  Please report the issue at <https://github.com/chgigot/epiphy/issues>.
Figure 17.12: Map of the SADIE clustering indices where red identify patches (outflow) and blue identify gaps (inflow)

An isocline plot can be obtained by setting the isoclines argument as TRUE.

plot(sadie_epi2, isoclines = TRUE)
Figure 17.13: Map of the SADIE clustering indices

17.2 Sparsely sampled data

In sparsely sampled data, the locations of the sampling units are unknown or are ignored in the analysis (Madden et al. 2007a). The aim is to describe the variability of disease intensity among sampling units, and there are two approaches: (1) testing the goodness of fit of the data to a statistical distribution, and (2) calculating indices of aggregation. The distributions and indices depend on whether the data are counts or incidences (proportions), so we treat the two cases separately.

17.2.1 Count data

17.2.1.1 Fit to distributions

Two statistical distributions can be adopted as reference for the description of random or aggregated patterns of disease data in the form of counts of infection within sampling units. Take the count of lesions on a leaf, or the count of diseased plants on a quadrat, as an example. If the presence of a lesion/diseased plant does not increase or decrease the chance that other lesions/diseased plants will occur, the Poisson distribution describes the distribution of lesions on the leaf. Otherwise, the negative binomial provides a better description.

Let’s work with the previous simulation data of 144 quadrats with a variable count of diseased plants per quadrat (in a maximum of 100). Notice that we won’t consider the location of each quadrat as in the previous analyses of intensively mapped data. We only need the vector with the number of infected units per sampling unit.

The epiphy package provides a function called fit_two_distr(), which allows fitting these two distributions for count data. In this case, either randomness assumption (Poisson distributions) or aggregation assumption (negative binomial) are made, and then, a goodness-of-fit comparison of both distributions is performed using a log-likelihood ratio test. The function requires a data frame created using the count() function where the number of infection units is designated as i. It won’t work with a single vector of numbers. We create the data frame using:

data_count <- epi2 |> 
  mutate(i = n) |>  # create i vector
  epiphy::count()   # create the map object of count class

We can now run the function that will look for the vector i. The function returns a list of four components including the outputs of the fitting process for both distribution and the result of the log-likelihood ratio test, the llr.

fit_data_count <- fit_two_distr(data_count)
summary(fit_data_count)
Fitting of two distributions by maximum likelihood
for 'count' data.
Parameter estimates:

(1) Poisson (random):
       Estimate  Std.Err Z value    Pr(>z)    
lambda 27.85417  0.43981  63.333 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(2) Negative binomial (aggregated):
       Estimate    Std.Err Z value    Pr(>z)    
k     0.6327452  0.0707846  8.9390 < 2.2e-16 ***
mu   27.8541667  2.9510195  9.4388 < 2.2e-16 ***
prob  0.0222118  0.0033463  6.6378 3.184e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit_data_count$llr
Likelihood ratio test

               LogLik Df  Chisq Pr(>Chisq)    
random :     -2654.71                         
aggregated :  -616.51  1 4076.4  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The very low value of the P-value of the LLR test suggests that the negative binomial provides a better fit to the data. The plot() function allows for visualizing the expected random and aggregated frequencies together with the observed frequencies. The number of breaks can be adjusted as indicated.

plot(fit_data_count, breaks = 5)
Figure 17.14: Frequencies of the observed and expected aggregated and random distributions

See below another way to plot by extracting the frequency data (and pivoting from wide to long format) from the generated list and using ggplot. Clearly, the negative binomial is a better description for the observed count data.

df <- fit_data_count$freq |>
  pivot_longer(2:4, names_to = "pattern", values_to = "value")

df |>
  ggplot(aes(category, value, fill = pattern)) +
  geom_col(position = "dodge", width = 2) +
  scale_fill_manual(values = c("gray70", "darkred", "steelblue")) +
  theme_r4pde()+
  theme(legend.position = "top")
Figure 17.15: Frequencies of the observed and expected aggregated and random distributions

17.2.1.2 Aggregation indices

Several indices summarize the degree of aggregation of counts (Fisher, Lloyd and Morisita).

idx <- agg_index(data_count, method = "fisher")
idx
Fisher's index of dispersion:
(Version for count data)
34.25
chisq.test(idx)

    Chi-squared test for (N - 1)*index following a chi-squared distribution
    (df = N - 1)

data:  idx
X-squared = 4897.2, df = 143, p-value < 2.2e-16
z.test(idx)

    One-sample z-test

data:  idx
z = 82.085, p-value < 2.2e-16
alternative hypothesis: two.sided
# Lloyd index

idx_lloyd <- agg_index(data_count, method = "lloyd")
idx_lloyd
Lloyd's index of patchiness:
2.194
idx_mori <- agg_index(data_count, method = "morisita")
idx_mori
Morisita's coefficient of dispersion:
(Version for count data)
2.186
# Using the vegan package
library(vegan)
z <- data_count$data$i
mor <- dispindmorisita(z)
mor
      imor     mclu      muni      imst pchisq
1 2.185591 1.008728 0.9922162 0.5041152      0

17.2.2 Incidence data

17.2.2.1 Fit to distributions

tas <-
  read.csv(
    "data/tasmania_test_1.txt",
    sep = ""
  )
head(tas,10)
   quad group_size count
1     1          6     4
2     2          6     6
3     3          6     6
4     4          6     6
5     5          6     6
6     6          6     6
7     7          6     6
8     8          6     6
9     9          6     4
10   10          6     6
# Create incidence object for epiphy
dat_tas <- tas |>
  mutate(n = group_size, i = count) |>
  epiphy::incidence()

## Fit to two distributions
fit_tas <- fit_two_distr(dat_tas)
summary(fit_tas)
Fitting of two distributions by maximum likelihood
for 'incidence' data.
Parameter estimates:

(1) Binomial (random):
     Estimate Std.Err Z value    Pr(>z)    
prob  0.90860 0.01494  60.819 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(2) Beta-binomial (aggregated):
      Estimate  Std.Err Z value    Pr(>z)    
alpha 1.923479 0.869621  2.2119  0.026976 *  
beta  0.181337 0.075641  2.3973  0.016514 *  
prob  0.913847 0.023139 39.4943 < 2.2e-16 ***
rho   0.322080 0.096414  3.3406  0.000836 ***
theta 0.475101 0.209789  2.2647  0.023534 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit_tas$llr
Likelihood ratio test

              LogLik Df  Chisq Pr(>Chisq)    
random :     -75.061                         
aggregated : -57.430  1 35.263   2.88e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fit_tas)

17.2.2.2 Aggregation indices

Overdispersion can first be checked by fitting an intercept-only binomial GLM and testing for overdispersion.

binom.tas = glm(cbind(count, group_size - count) ~ 1,
                family = binomial,
                data = tas)
summary(binom.tas)

Call:
glm(formula = cbind(count, group_size - count) ~ 1, family = binomial, 
    data = tas)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.2967     0.1799   12.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 117.76  on 61  degrees of freedom
Residual deviance: 117.76  on 61  degrees of freedom
AIC: 152.12

Number of Fisher Scoring iterations: 5
library(performance)
check_overdispersion(binom.tas)
# Overdispersion test

 dispersion ratio =   2.412
          p-value = < 0.001

The C(alpha) test from the epiphy package tests the binomial null hypothesis directly.

library(epiphy)
tas2 <- tas |>
  mutate(i = count,
         n = group_size) |>  # create i vector
  epiphy::incidence()

t <- agg_index(tas2, flavor = "incidence")
t
Fisher's index of dispersion:
(Version for incidence data)
2.348
calpha.test(t)

    C(alpha) test

data:  t
z = 7.9886, p-value = 1.365e-15

17.2.2.3 Binary power law

The binary power law (BPL), introduced by Hughes and Madden (1992), relates the observed variance of disease incidence among sampling units to the variance expected for a binomial distribution (a random pattern). It is widely used to describe spatial heterogeneity in incidence data (Madden et al. 2007a; Madden et al. 2018). In its general form, \(v = A_p \, v_{bin}^{\beta}\), where \(A_p\) and \(\beta\) are parameters. For a random (binomial) pattern, \(A_p = 1\) and \(\beta = 1\).

The BPL relates two variances:

  • Observed Variance (V): \(V = \text{var}[X] = E[(X - np)^2]\)

  • Expected Binomial Variance (Vbin): \(V_{bin} = np(1 - p)\)

For proportions:

  • Observed Variance (v): \(v = \text{var}[x] = E[(x - p)^2]\)

  • Expected Binomial Variance (v_bin): \(v_{bin} = \frac{p(1 - p)}{n}\)

The binomial distribution assumes independent observations of \(X\), often interpreted as a random spatial pattern. If \(X\) follows a binomial distribution, then \(\beta = 1\) in the BPL equations. Deviations from this indicate spatial heterogeneity or overdispersion. The log-log transformation is often used in statistical analyses due to the linear relationship between the log-transformed variables. Interpretation is as follows: when \(\beta = 1\): Spatial heterogeneity is constant, indicating fixed overdispersion; when \(\beta > 1\): Overdispersion varies systematically with incidence, showing a change in spatial heterogeneity with disease incidence levels. In practice, researchers fit the log-log model to empirical data to estimate \(\beta\) and \(A_p\). This involves ordinary least squares (OLS) regression:

\[ \ln(v_j) = \ln(A_p) + \beta \ln\left[\frac{\hat{p}_j(1 - \hat{p}_j)}{n_j}\right] + \epsilon_j \]

where: - \(j\) refers to specific observations, - \(\hat{p}\) is the estimated proportion, - \(\epsilon_j\) is the residual error.

The intercept of the fitted model (\(\ln(A_p)\)) provides the scale parameter, and the slope (\(\beta\)) gives insight into the nature of spatial heterogeneity and overdispersion in the data.