4  Simulation

The examples in Chapter 4 require that the search path contains the following namespaces,

library(groupedHyperframe.random)
# Loading required package: groupedHyperframe
search path on author’s computer running RStudio (Posit Team 2025)
search()
#  [1] ".GlobalEnv"                       "package:groupedHyperframe.random" "package:groupedHyperframe"       
#  [4] "package:stats"                    "package:graphics"                 "package:grDevices"               
#  [7] "package:utils"                    "package:datasets"                 "package:methods"                 
# [10] "Autoloads"                        "package:base"

4.1 Simulated Point-Pattern

Function .rppp() simulates superimposed point-pattern objects (ppp.object, Chapter 23) via vectorized parameterization (Section 33.7), which is supported in the parameter-specification of both the random point-pattern and the random marks.

Listing 4.1 simulates an \(x\)- and \(y\)-coordinates only, two superimposed Matérn (1986)’s cluster processes \((\kappa, \mu, s)\) with parameters \((\kappa_1, \mu_1, s_1) = (10,8,.15)\) and \((\kappa_2, \mu_2, s_2) = (5,4,.06)\) using a vectorized parameterization of kappa \(\kappa= (\kappa_1, \kappa_2)\), mu \(\mu= (\mu_1, \mu_2)\) and scale \(s= (s_1, s_2)\).

Listing 4.1: Superimposed Matérn (1986)’s cluster processes
set.seed(12); r = .rppp(rMatClust(kappa = c(10, 5), mu = c(8, 4), scale = c(.15, .06)))
# Point-pattern simulated by `spatstat.random::rMatClust()`
# 

Listing 4.2 simulates two superimposed marked point-pattern objects,

  1. one marked point-pattern object with
    • \(x\)- and \(y\)-coordinates by Matérn (1986)’s cluster process with parameters \((\kappa_1,\mu_1,s_1) = (10,8,.15)\);
    • log-normal \((\log\mu, \log\sigma)\) marks with parameters \((\log\mu_1,\log\sigma_1)=(3,.4)\);
    • negative-binomial \((r, p)\) marks with parameters \((r,p)=(4,.3)\).
  2. another marked point-pattern object with
    • \(x\)- and \(y\)-coordinates by Matérn (1986)’s cluster process with parameters \((\kappa_2,\mu_2,s_2) = (5,4,.06)\);
    • log-normal marks with parameters \((\log\mu_2,\log\sigma_2)=(5,.2)\)
    • negative-binomial marks with parameters \((r,p)=(4,.3)\).

using a vectorized parameterization, in addition to that of Listing 4.1, of meanlog \(\log\mu= (\log\mu_1, \log\mu_2)\), sdlog \(\log\sigma= (\log\sigma_1, \log\sigma_2)\), size \(r\) and prob \(p\).

Listing 4.2: Superimposed Matérn (1986)’s cluster processes with log-normal and negative-binomial marks
set.seed(25); r1 = .rppp(
  rMatClust(kappa = c(10, 5), mu = c(8, 4), scale = c(.15, .06)), 
  rlnorm(meanlog = c(3, 5), sdlog = c(.4, .2)),
  rnbinom(size = 4, prob = .3) # shorter arguments recycled
)
# Point-pattern simulated by `spatstat.random::rMatClust()`
# Mark simulated by `stats::rlnorm()`
# Mark simulated by `stats::rnbinom()`

Listing 4.3 simulates two superimposed point-pattern objects,

  1. one marked point-pattern object with
    • \(x\)- and \(y\)-coordinates by Poisson point-pattern \(\lambda\) with parameter \(\lambda_1=3\),
    • log-normal marks with parameters \((\log\mu_1,\log\sigma_1)=(3,.4)\);
    • negative-binomial marks with parameters \((r_1,p_1)=(4,.3)\).
  2. another marked point-pattern object with
    • \(x\)- and \(y\)-coordinates by Poisson point-pattern with parameters \(\lambda_2=6\)
    • log-normal marks with parameters \((\log\mu_2,\log\sigma_2)=(5,.2)\)
    • negative-binomial marks with parameters \((r_2,p_2)=(6,.1)\).

using a vectorized parameterization, in addition to that of Listing 4.2, of lambda \(\lambda= (\lambda_1, \lambda_2)\).

Listing 4.3: Superimposed Poisson point-patterns with log-normal and negative-binomial marks
set.seed(62); r2 = .rppp(
  rpoispp(lambda = c(3, 6)),
  rlnorm(meanlog = c(3, 5), sdlog = c(.4, .2)),
  rnbinom(size = c(4, 6), prob = c(.3, .1))
)
# Point-pattern simulated by `spatstat.random::rpoispp()`
# Mark simulated by `stats::rlnorm()`
# Mark simulated by `stats::rnbinom()`

Listing 4.4 superimposes multiple simulated (marked) point-pattern objects.

Listing 4.4: Review: Superimpose simulated point-patterns
spatstat.geom::superimpose(r1, r2)
# Marked planar point pattern: 140 points
# Mark variables: m1, m2 
# window: rectangle = [0, 1] x [0, 1] units

4.2 Simulated Grouped Hyper Data Frame

Function grouped_rppp() simulates a grouped hyper data frame containing one-and-only-one point-pattern (ppp) hypercolumn (Chapter 3) via matrix parameterization (Section 33.3), which is supported in the subject-specific parameter-specification of both the random point-pattern and the random marks. The full process consists of three steps:

  1. Simulate the subject-specific parameters from a multivariate normal distribution using the function mvrnorm2(), e.g., Section 4.2.1.1.
  2. Simulate the number of (marked) point-patterns per subject using the function base::sample(), e.g., Section 4.2.1.2;
  3. Simulate the grouped hyper data frame using the function grouped_rppp(), e.g., Section 4.2.1.3.

4.2.1 with Superimposed Matérn (1986)’s Cluster Processes

4.2.1.1 Simulate Subject-Specific Parameters

Consider two superimposed Matérn (1986)’s cluster processes, each with one log-normal mark. The population parameters are

  • \(x\)- and \(y\)-coordinates by Matérn (1986)’s cluster process \((\kappa_1,\mu_1,s_1) = (3,10,.4)\), with a log-normal mark \((\log\mu_1,\log\sigma_1)=(3,.4)\)
  • \(x\)- and \(y\)-coordinates by Matérn (1986)’s cluster process \((\kappa_2,\mu_2,s_2) = (2,5,.2)\), with a log-normal mark \((\log\mu_2,\log\sigma_2)=(5,.2)\)

Listing 4.5 simulates for three (3L) subjects (e.g., patients). The subject-specific parameters deviate from the population parameters under a multivariate normal distribution specified by the standard deviation(s) \(\sigma\) of the subject-specific random effects, i.e., \(\sigma_\kappa=.2\), \(\sigma_\mu=.5\), \(\sigma_s=.05\), \(\sigma_{\log\mu}=.1\) and \(\sigma_{\log\sigma}=.01\) (Section 33.4). Listing 4.5 makes sure that all subject-specific parameters satisfy that kappa \(\kappa>1\), mu \(\mu>1\), scale \(s>0\) for Matérn (1986)’s cluster processes, and that sdlog \(\log\sigma>0\) for log-normal distribution.

Listing 4.5: Example: Create subject-specific parameters p_Matern
set.seed(39); p_Matern = mapply(
  FUN = mvrnorm2, 
  mu = list(kappa = c(3,2), mu = c(10,5), scale = c(.4,.2), meanlog = c(3,5), sdlog = c(.4,.2)), 
  sd = list(kappa = .2, mu = .5, scale = .05, meanlog = .1, sdlog = .01), 
  MoreArgs = list(n = 3L), 
  SIMPLIFY = FALSE
) |>
  within.list(expr = {
    kappa = pmax(kappa, 1 + .Machine$double.eps)
    mu = pmax(mu, 1 + .Machine$double.eps)
    scale = pmax(scale, .Machine$double.eps)
    sdlog = pmax(sdlog, .Machine$double.eps)
  })

Listing 4.6 shows the subject-specific parameters p_Matern as a list of matrices, where each matrix represents one parameter used in generating the random marked point-patterns. In these matrices,

  • Each row corresponds to a subject; i.e., each matrix in p_Matern contains three (3L) rows, one for each of the three (3L) subjects;
  • Each column corresponds to a random generation parameter for each superimposed marked point-pattern; i.e., each matrix in p_Matern contains two (2L) columns, one for each of the two (2L) Matérn (1986)’s cluster processes to be superimposed.
Listing 4.6: Example: Subject-specific parameters p_Matern
Code
p_Matern
# $kappa
#          [,1]     [,2]
# [1,] 3.119196 1.962887
# [2,] 2.906535 1.754151
# [3,] 2.915672 1.914559
# 
# $mu
#           [,1]     [,2]
# [1,]  9.519187 4.489174
# [2,]  9.947885 4.688761
# [3,] 10.029897 5.418501
# 
# $scale
#           [,1]      [,2]
# [1,] 0.4227306 0.2351866
# [2,] 0.4058303 0.1622515
# [3,] 0.3726093 0.1727232
# 
# $meanlog
#          [,1]     [,2]
# [1,] 2.982251 5.007964
# [2,] 2.878996 4.976754
# [3,] 3.034019 4.852771
# 
# $sdlog
#           [,1]      [,2]
# [1,] 0.4055664 0.1975062
# [2,] 0.4172841 0.2115214
# [3,] 0.4076479 0.1970090

4.2.1.2 Simulate Number of Point-Patterns per Subject

Listing 4.7 simulates one (1L) to four (4L) point-patterns (e.g., medical images) for each of the three (3L) subjects (e.g., patients).

Listing 4.7: Example: Simulate number of point-patterns per subject
set.seed(37); (n = sample(x = 1:4, size = 3L, replace = TRUE)) 
# [1] 2 3 4

4.2.1.3 Simulate Grouped Hyper Data Frame

Function grouped_rppp() takes (each matrix in) the matrix parameterization p_Matern (Listing 4.6) and simulates a grouped hyper data frame with (Listing 4.8)

  • one point-pattern (ppp) hypercolumn, and
  • one-or-more regular columns representing the grouping structure.
Listing 4.8: Example: Simulate grouped hyper data frame r_Matern
set.seed(76); r_Matern = p_Matern |> 
  with.default(expr = {
    grouped_rppp(
      rMatClust(kappa = kappa, scale = scale, mu = mu), 
      rlnorm(meanlog = meanlog, sdlog = sdlog),
      n = n
    )
  })
Simulated grouped hyper data frame r_Matern
r_Matern
# Grouped Hyperframe: ~g1/g2
# 
# 9 g2 nested in
# 3 g1
# 
# Preview of first 10 (or less) rows:
# 
#     ppp g1 g2
# 1 (ppp)  1  1
# 2 (ppp)  1  2
# 3 (ppp)  2  1
# 4 (ppp)  2  2
# 5 (ppp)  2  3
# 6 (ppp)  3  1
# 7 (ppp)  3  2
# 8 (ppp)  3  3
# 9 (ppp)  3  4

4.2.2 with Superimposed Poisson Point-Pattern

Consider two superimposed Poisson point-pattern, each with one negative-binomial mark. The population parameters are

  • \(x\)- and \(y\)-coords by Poisson point-pattern \(\lambda_1=3\), with a negative-binomial mark \((r_1,p_1)=(4,.3)\).
  • \(x\)- and \(y\)-coords by Poisson point-pattern \(\lambda_2=6\), with a negative-binomial mark \((r_2,p_2)=(6,.1)\).

Listing 4.9 also simulate for three (3L) subjects (e.g., patients). The subject-specific parameters deviate from the population parameters under a multivariate normal distribution specified by the standard deviation(s) \(\sigma\) of the subject-specific random effects, i.e., \(\sigma_\lambda=.1\), \(\sigma_r=.2\) and \(\sigma_p=.05\) (Section 33.4). Listing 4.9 makes sure that all subject-specific parameters satisfy that lambda \(\lambda>0\) for Poisson point-pattern, and that size \(r>0\), prob \(0<p<1\) for negative-binomial distribution.

Listing 4.9: Example: Create subject-specific parameters p_Poisson
set.seed(39); p_Poisson = mapply(
  FUN = mvrnorm2, 
  mu = list(lambda = c(3, 6), size = c(4, 6), prob = c(.3, .1)), 
  sd = list(lambda = .1, size = .2, prob = .05), 
  MoreArgs = list(n = 3L), 
  SIMPLIFY = FALSE
) |>
  within.list(expr = {
    lambda = pmax(lambda, .Machine$double.eps)
    size = pmax(size, .Machine$double.eps)
    prob = pmin(pmax(prob, .Machine$double.eps), 1 - .Machine$double.eps)
  })
Subject-specific parameters p_Poisson
p_Poisson
# $lambda
#          [,1]     [,2]
# [1,] 3.059598 5.981443
# [2,] 2.953268 5.877076
# [3,] 2.957836 5.957280
# 
# $size
#          [,1]     [,2]
# [1,] 3.807675 5.795670
# [2,] 3.979154 5.875505
# [3,] 4.011959 6.167400
# 
# $prob
#           [,1]       [,2]
# [1,] 0.3227306 0.13518657
# [2,] 0.3058303 0.06225150
# [3,] 0.2726093 0.07272321

Listing 4.10 simulates a grouped hyper data frame using the same grouping size as in Listing 4.7.

Listing 4.10: Example: Simulate grouped hyper data frame r_Poisson
set.seed(76); r_Poisson = p_Poisson |> 
  with.default(expr = {
    grouped_rppp(
      rpoispp(lambda = lambda),
      rnbinom(size = size, prob = prob),
      n = n
    )
  })
Simulated grouped hyper data frame r_Poisson
r_Poisson
# Grouped Hyperframe: ~g1/g2
# 
# 9 g2 nested in
# 3 g1
# 
# Preview of first 10 (or less) rows:
# 
#     ppp g1 g2
# 1 (ppp)  1  1
# 2 (ppp)  1  2
# 3 (ppp)  2  1
# 4 (ppp)  2  2
# 5 (ppp)  2  3
# 6 (ppp)  3  1
# 7 (ppp)  3  2
# 8 (ppp)  3  3
# 9 (ppp)  3  4

4.2.3 Superimpose Grouped Hyper Data Frames

Listing 4.11 (Table 16.2, Section 17.8) superimposes multiple simulated grouped hyper data frames.

Listing 4.11: Example: Superimpose simulated grouped hyper data frame
r = superimpose.groupedHyperframe(r_Matern, r_Poisson)

Listing 4.12 checks the return of Listing 4.11 (note that the marks are plotted on different scales).

Listing 4.12: Example: a quick check
Code
spatstat.geom::solist(
  Mattern = r_Matern$ppp[[1L]],
  Poisson = r_Poisson$ppp[[1L]],
  Superimposed = r$ppp[[1L]]
) |>
  spatstat.geom::plot.solist(main = '')