R语言spatstat包中mppm对象包络线绘制时nrank参数的调整方法及无效问题解决

Great question—this is a common gotcha when working with mppm objects and pooled envelopes in spatstat! Let’s unpack what’s happening here and how to choose the right nrank value.

Why nrank Seems to Have No Effect #

First, let’s clarify why your nrank parameter isn’t showing obvious changes:

• When you use pre-generated simulations via simulate=Sims, nrank defines how many extreme simulation results to exclude from each tail of the distribution to form the envelope. If your nsims value is small (like ~58 in your code), shifting nrank by 1 only changes the envelope from absolute min/max to the second min/max—this difference can be subtle, especially with noisy spatial data.
• The VARIANCE=TRUE flag switches to a parametric envelope (mean ± standard deviation of simulations) instead of the nonparametric rank-based method. This is an entirely different calculation, so it’s expected to produce distinct results.

How to Choose nrank for mppm Objects #

Since you’re working with multiple independent datasets (the foundation of your mppm model), you need to account for multiple testing to maintain your desired overall confidence level. Here’s a step-by-step approach:

1. Define Your Confidence Target #

Start with your overall desired confidence level (e.g., 95%). For k independent datasets, use the Bonferroni correction to adjust the significance level for each individual dataset:

• Overall confidence: alpha_total = 0.95
• Per-dataset confidence: alpha_single = alpha_total^(1/k) (this ensures the combined confidence stays near your target)

2. Calculate nsims and nrank #

Your initial gamma calculation already uses this correction for k=3 datasets. To map this to nrank:

• Compute nsims as you did: nsims = round(1/gamma - 1)
• Calculate the per-dataset significance threshold: alpha_per_test = 1 - alpha_single
• Use this formula to get the correct nrank (the number of extreme simulations to exclude from each tail):

  nrank_val <- floor( (nsims + 1) * alpha_per_test / 2 )
  

  This ensures the proportion of simulations outside the envelope matches your corrected significance level.

3. Refined Code Example #

Here’s a cleaned-up version of your code that ties this all together explicitly:

# Get number of independent datasets
k <- nrow(data)
# Overall desired confidence level
alpha_total <- 0.95
# Bonferroni-corrected per-dataset confidence
alpha_single <- alpha_total^(1/k)
gamma <- 1 - alpha_single
nsims <- round(1/gamma - 1)

# Generate simulations for each dataset
sims <- simulate(model, nsim = nsims)
SIMS <- lapply(1:nrow(sims), function(i) as.solist(sims[i,, drop=TRUE]))
Hplus <- cbind(data, hyperframe(Sims = SIMS))

# Calculate appropriate nrank
nrank_val <- floor( (nsims + 1) * (1 - alpha_single) / 2 )

# Compute envelopes with correct nrank
EE <- with(Hplus, envelope(Points, Kcross.inhom,
                           funargs = list("A", "A"),
                           nsim = nsims,
                           simulate = Sims,
                           savefuns = TRUE,
                           nrank = nrank_val))

# Plot pooled envelope
plot(pool(EE), main = "A to A interactions")

Key Tips #

• Make nrank changes visible: If you want to see clear differences from adjusting nrank, increase nsims (e.g., to 199). Larger simulation sets make the shift between rank thresholds more noticeable.
• Validate simulation structure: Double-check that SIMS is a list where each element is a solist containing exactly nsims simulated point patterns. Use str(SIMS[[1]]) to confirm.
• Avoid over-reliance on VARIANCE=TRUE: This parametric approach assumes simulation results are normally distributed, which isn’t always true for spatial summary functions. Stick to rank-based envelopes unless you have a strong reason to use parametric ones.

内容的提问来源于stack exchange，提问作者Camille Gontier