# Graphical normality test (Myllymaki and Mrkvicka, 2020, Section 3.3.)
#=========================
if(require("fda.usc", quietly=TRUE)) {
data("poblenou")
dat <- poblenou[['nox']][['data']][,'H10']
n <- length(dat)
# The number of simulations
nsim <- nsimsub <- 199
nsim <- nsimsub <- 19
set.seed(200127)
# General setup
#==============
# 1. Fit the model
mu <- mean(dat)
sigma <- sd(dat)
# 2. Simulate a sample from the fitted null model and
# compute the test vectors for data (obs) and each simulation (sim)
r <- seq(min(dat), max(dat), length=100)
obs <- stats::ecdf(dat)(r)
sim <- sapply(1:nsimsub, function(i) {
x <- rnorm(n, mean = mu, sd = sigma)
stats::ecdf(x)(r)
})
cset <- create_curve_set(list(r = r, obs = obs, sim_m = sim))
# 3. Simulate another sample from the fitted null model.
# 4. Fit the null model to each of the patterns,
# simulate a sample from the null model,
# and compute the test vectors for all.
cset.ls <- vector("list", nsim)
for(rep in 1:nsim) {
x <- rnorm(n, mean = mu, sd = sigma)
mu2 <- mean(x)
sigma2 <- sd(x)
obs2 <- stats::ecdf(x)(r)
sim2 <- sapply(1:nsimsub, function(i) {
x2 <- rnorm(n, mean = mu2, sd = sigma2)
stats::ecdf(x2)(r)
})
cset.ls[[rep]] <- create_curve_set(list(r = r, obs = obs2, sim_m = sim2))
}
# Perform the adjusted test
res <- GET.composite(X = cset, X.ls = cset.ls, type = 'erl')
plot(res) + ggplot2::labs(x = "NOx", y = "Ecdf")
}
# Example of a point pattern data
#================================
# Test the fit of a Matern cluster process.
# \donttest{
if(require("spatstat.model", quietly=TRUE)) {
data("saplings")
saplings <- as.ppp(saplings, W = square(75))
# First choose the r-distances
rmin <- 0.3; rmax <- 10; rstep <- (rmax-rmin)/500
r <- seq(0, rmax, by = rstep)
# Number of simulations
nsim <- 19 # Increase nsim for serious analysis!
# Option 1: Give the fitted model object to GET.composite
#--------------------------------------------------------
# This can be done and is preferable when the model is
# a point process model of spatstat.
# 1. Fit the Matern cluster process to the pattern
# (using minimum contrast estimation with the K-function)
M1 <- kppm(saplings~1, clusters = "MatClust", statistic = "K")
summary(M1)
# 2. Make the adjusted global area rank envelope test using the L(r)-r function
adjenvL <- GET.composite(X = M1, nsim = nsim,
testfuns = list(L =list(fun="Lest", correction="translate",
transform=expression(.-r), r=r)), # passed to envelope
type = "area", r_min = rmin, r_max = rmax)
# Plot the test result
plot(adjenvL)
# Option 2: Generate the simulations "by yourself"
#-------------------------------------------------
# and provide them as curve_set or envelope objects
# Preferable when you want to have a control
# on the simulations yourself.
# 1. Fit the model
M1 <- kppm(saplings~1, clusters = "MatClust", statistic = "K")
# 2. Generate nsim simulations by the given function using the fitted model
X <- envelope(M1, nsim = nsim, savefuns = TRUE,
fun = "Lest", correction = "translate",
transform = expression(.-r), r = r)
plot(X)
# 3. Create another set of simulations to be used to estimate
# the second-state p-value (as proposed by Baddeley et al., 2017).
simpatterns2 <- simulate(M1, nsim = nsim)
# 4. Calculate the functions for each pattern
simf <- function(rep) {
# Fit the model to the simulated pattern Xsims[[rep]]
sim_fit <- kppm(simpatterns2[[rep]], clusters = "MatClust", statistic = "K")
# Generate nsim simulations from the fitted model
envelope(sim_fit, nsim = nsim, savefuns = TRUE,
fun = "Lest", correction = "translate",
transform = expression(.-r), r = r)
}
X.ls <- parallel::mclapply(X = 1:nsim, FUN = simf, mc.cores = 1) # list of envelope objects
# 5. Perform the adjusted test
res <- GET.composite(X = X, X.ls = X.ls, type = "area",
r_min = rmin, r_max = rmax)
plot(res)
}# }
Run the code above in your browser using DataLab