## Package # install.packages("stopp") library(stopp) ## for saving plots: require("grDevices") # dir.create("Figures") #### Codes from the paper # 2. Data types # 2.1. Spatio-temporal point patterns set.seed(2) df <- data.frame(runif(100), runif(100), runif(100)) stp1 <- stp(df) stp1 # pdf("Figures/stp.pdf", width = 12, height = 4) plot(stp1) # dev.off() set.seed(2) df_net <- data.frame(runif(100, 0, 0.85), runif(100, 0, 0.85), runif(100)) stlp1 <- stp(df_net, L = chicagonet) stlp1 # pdf("Figures/stpL.pdf", width = 12, height = 4) plot(stlp1, tcum = FALSE) # dev.off() # 2.2. Marked point processes set.seed(2) dfA <- data.frame(x = runif(100), y = runif(100), t = runif(100), m1 = rnorm(100), m2 = rep(c("C"), times = 100)) dfB <- data.frame(x = runif(50), y = runif(50), t = runif(50), m1 = rnorm(25), m2 = rep(c("D"), times = 50)) stpm2 <- stpm(rbind(dfA, dfB), names = c("continuous", "dichotomous")) # pdf("Figures/mark2.pdf", width = 12, height = 6) plot(stpm2) # dev.off() # 2.3. Spatio-temporal covariates set.seed(2) df <- data.frame(runif(100), runif(100), runif(100), rpois(100, 15)) sim_cov <- stcov(df, interp = FALSE, names = "SimulatedCovariate") interp_cov <- stcov(df, mult = 20, names = "InterpolatedCovariate") # pdf("Figures/simcov.pdf", width = 6, height = 6) plot(sim_cov) # dev.off() # pdf("Figures/cov.pdf", width = 6, height = 6) plot(interp_cov) # dev.off() # 3. Datasets data("greececatalog", package = "stopp") # pdf("Figures/greece.pdf", width = 12, height = 4) plot(greececatalog) # dev.off() data("valenciacrimes", package = "stopp") # pdf("Figures/valencia.pdf", width = 18, height = 12) plot(valenciacrimes) # dev.off() data("chicagonet", package = "stopp") data("valencianet", package = "stopp") # pdf("Figures/nets.pdf", width = 12, height = 6) par(mfrow = c(1, 2)) plot(valencianet, col = "grey"); axis(1); axis(2) plot(chicagonet, col = "grey"); axis(1); axis(2) # dev.off() # 4. Simulations rstpp(500) rstpp(lambda = function(x, y, t, a) exp(a[1] + a[2] * x), par = c(2, 6)) set.seed(95) X <- rETASp(c(0.1293688525, 0.003696, 0.013362, 1.2,0.424466, 1.164793), betacov = 0.5, xmin = 600, xmax = 2200, ymin = 4000, ymax = 5300) # pdf("Figures/etas.pdf", width = 6, height = 6) plot(X) # dev.off() # 5. Exploratory analysis set.seed(2) df_net <- data.frame(runif(25, 0, 0.85), runif(25, 0, 0.85), runif(25)) stlp1 <- stp(df_net, L = chicagonet) lambda <- rep(diff(range(stlp1$df$x)) * diff(range(stlp1$df$y)) * diff(range(stlp1$df$t)) / spatstat.geom::volume(stlp1$L), nrow(stlp1$df)) k <- localSTLKinhom(stlp1, lambda = rep( 0.02708656 , 25), normalize = TRUE) # pdf("Figures/listas.pdf", width = 12, height = 4) plot(k, 1:3) # dev.off() # 5.1 Local test set.seed(2) X <- rstpp(lambda = function(x, y, t, a) {exp(a[1] + a[2]*x)}, par = c(.005, 5)) set.seed(2) Z <- rstpp(lambda = 30) test <- localtest(X, Z, method = "K", k = 3) test # pdf("Figures/localtest.pdf", width = 12, height = 4) plot(test) # dev.off() # 6. Model fitting # 6.1. Homgeneous and inhomogeneous spatio-temporal Poisson point processes set.seed(2) ph <- rstpp(lambda = 200) hom1 <- stppm(ph, formula = ~ 1, seed = 2) hom1 set.seed(2) pin <- rstpp(lambda = function(x, y, t, a) {exp(a[1] + a[2]*x)}, par = c(2, 6)) inh1 <- stppm(pin, formula = ~ x, seed = 2) inh1 inh2 <- stppm(pin, formula = ~ s(x, y, bs = "tp", k = 30)) # pdf("Figures/nonpar.pdf", width = 12, height = 6) plot(inh2) # dev.off() # 6.2. Spatio-temporal Poisson point processes with dependence on external covariates set.seed(2) df1 <- data.frame(runif(100), runif(100), runif(100), rpois(100, 15)) df2 <- data.frame(runif(100), runif(100), runif(100), rpois(100, 15)) obj1 <- stcov(df1, names = "cov1") obj2 <- stcov(df2, names = "cov2") covariates <- list(cov1 = obj1, cov2 = obj2) inh3 <- stppm(pin, formula = ~ x + cov2, covs = covariates, spatial.cov = TRUE, seed = 2) inh3 # 6.3. Multitype spatio-temporal Poisson point processes set.seed(2) dfA <- data.frame(x = runif(100), y = runif(100), t = runif(100), m1 = rep(c("A"), times = 100)) dfB <- data.frame(x = runif(50), y = runif(50), t = runif(50), m1 = rep(c("B"), each = 50)) stpm1 <- stpm(rbind(dfA, dfB)) # pdf("Figures/multi.pdf", width = 6, height = 6) plot(stpm1) # dev.off() inh4 <- stppm(stpm1, formula = ~ x + s(m1, bs = "re"), marked = TRUE, seed = 2) # pdf("Figures/multimodel.pdf", width = 12, height = 6) plot(inh4) #dev.off() # 6.4. Spatio-temporal Poisson point processes with separable intensity crimesub <- stpm(valenciacrimes$df[101:200, ], names = colnames(valenciacrimes$df)[-c(1:3)], L = valencianet) mod1 <- sepstlppm(crimesub, spaceformula = ~x , timeformula = ~ day) # pdf("Figures/sepL.pdf", width = 12, height = 6) plot(mod1) # dev.off() # 6.5 nonsepmod <- stppm(greececatalog, formula = ~ x + y + t + x:y + y:t + I(x^2) + I(y^2) + I(t^2) + I(x^2):I(y^2), seed = 2) summary(nonsepmod) # pdf("Figures/nonsepmod.pdf", width = 12, height = 6) plot(nonsepmod) # dev.off() summary(nonsepmod$mod_global) # 6.6. Log-Gaussian Cox processes catsub <- stp(greececatalog$df[1:200, ]) lgcp1 <- stlgcppm(catsub, seed = 2) lgcp1 # 6.7. Local models # Local spatio-temporal Poisson point processes set.seed(2) inh <- rstpp(lambda = function(x, y, t, a) {exp(a[1] + a[2]*x)}, par = c(0.005, 5)) inh_local <- locstppm(inh, formula = ~ x, seed = 2) summary(inh_local) # pdf("Figures/localPlot.pdf", width = 12, height = 6) localplot(inh_local) # dev.off() # pdf("Figures/localSum.pdf", width = 10, height = 6) localsummary(inh_local) # dev.off() # Local spatio-temporal log-Gaussian Cox processes lgcp2 <- stlgcppm(catsub, second = "local", seed = 2) lgcp2 # pdf("Figures/plotl2.pdf", width = 12, height = 6) plot(lgcp2) # dev.off() # pdf("Figures/localplotstlgcppm.pdf", width = 12, height = 4) localplot(lgcp2) # dev.off() # 7. Diagnostics # 7.1. Global diagnostics set.seed(2) inh <- rstpp(lambda = function(x, y, t, a) {exp(a[1] + a[2]*x)}, par = c(.3, 6)) mod1 <- stppm(inh, formula = ~ 1, seed = 2) mod2 <- stppm(inh, formula = ~ x, seed = 2) (g1 <- globaldiag(inh, mod1$l)) (g2 <- globaldiag(inh, mod2$l)) # pdf("Figures/diag1.pdf", width = 12, height = 4) plot(g1) # dev.off() # pdf("Figures/diag2.pdf", width = 12, height = 4) plot(g2) # dev.off() # 7.2. Local diagnostics res <- localdiag(inh, mod1$l, p = .9) res # pdf("Figures/ldiag.pdf", width = 12, height = 12) plot(res) # dev.off() # pdf("Figures/infl.pdf", width = 12, height = 12) infl(res) # dev.off()()