Modulo 2 - Unidad 2
David Arango-Londoño
11/04/2023
Patrones Puntuales con Datos de Sismos para Colombia
A continuación se realiza una aplicación de patrones puntuales con base en los datos de sismos registrados en la pagina del USGS de los estados unidos: https://earthquake.usgs.gov/earthquakes/map/
#cargando librerias
require(spatstat)
require(rgdal)
require(maptools)
require(sp)
require(raster)
require(stpp)
require(plot3D)
require(KernSmooth)
library(sf) # Development version in the ppp branch
library(tmap)
library(here)
#leyendo shape, estableciendo ventana de trabajo
global <- readShapePoly("~/Documents//Univalle 2019-1/Tratamientos Datos Espaciales/shape_global/g2008_0.shp")
Colombia=global[global$ADM0_NAME=="Colombia",1]
global <- st_read(dsn = here("~/Documents/Univalle 2019-1/Tratamientos Datos Espaciales/shape_global/g2008_0.shp"))## Reading layer `g2008_0' from data source `/Users/davidarango/Documents/Univalle 2019-1/Tratamientos Datos Espaciales/shape_global/g2008_0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 534 features and 11 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -180 ymin: -55.72333 xmax: 180 ymax: 83.62742
## geographic CRS: WGS 84Colombia=global[global$ADM0_NAME=="Colombia",1]
Colombia = st_transform(Colombia, crs = 6345)
colombia_win=as.owin(Colombia)
colombia_win=simplify.owin(colombia_win,0.5)#analisis espacial y temporal de las propiedades de primer orden
#cargando csv - estableciendo columnas a trabajar x,y,t
sismos=st_read("~/Documents/Univalle 2018-2/Patrones Puntuales/avance/sismos.shp")## Reading layer `sismos' from data source `/Users/davidarango/Documents/Univalle 2018-2/Patrones Puntuales/avance/sismos.shp' using driver `ESRI Shapefile'
## Simple feature collection with 232 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -78.388 ymin: 1.135 xmax: -71.21 ymax: 11.6675
## geographic CRS: GCS_unknown#sismos = st_transform(sismos, crs = 6345)
terremotos <- as.3dpoints(cbind(st_coordinates(sismos),sismos$year))
#Graficos espacial y temporal
par(mfrow=c(1,1))
plot(terremotos,window=Colombia)
Como se puede observar en la grafica los datos de los sismos se agregan principalmente en el centro del pais, sin embargo debemos tener en cuenta que esta información es espacio temporal ya que tenemos las localizaciones en varios momentos del tiempo (1990 a 2017).
#Grafico 3D espacio-temporal (??A que hacen referencia los parametros theta,phi,etc?)
par(mfrow=c(1,1))
cubo <- scatter3D(terremotos[,1],terremotos[,2],terremotos[,3],theta=-45,phi=30,xlab="\n x",ylab="\n y",zlab="\n t",ticktype="detailed",col="black",pch=".",cex=5,nticks=5,cex.axis=0.6,cex.lab=1)
A continuación se presenta la grafica del patron espacio temporal.
#Datos a patron puntual
require(leaflet)
leaflet() %>% addTiles() %>% addCircleMarkers(lng = terremotos[,1],lat =terremotos[,2] )summary(terremotos[,1:2])## x y
## Min. :-78.39 Min. : 1.135
## 1st Qu.:-76.50 1st Qu.: 5.041
## Median :-74.49 Median : 6.748
## Mean :-74.67 Mean : 6.110
## 3rd Qu.:-72.98 3rd Qu.: 6.884
## Max. :-71.21 Max. :11.668Colombia## Simple feature collection with 1 feature and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 1387465 ymin: -489162.4 xmax: 2787315 ymax: 1427343
## projected CRS: NAD83(2011) / UTM zone 16N
## AREA geometry
## 370 92.77901 POLYGON ((2227113 1359153, ...colombia_win=owin(xrange = c(-79.04723,-66.86983),yrange = c(-4.230484,12.46331))
terremotospp <- ppp(x=terremotos[,1],y=terremotos[,2],window=colombia_win)
#eliminando datos duplicados
terremotospp <- unique(terremotospp)
plot(terremotospp)
#
#
terremotospp2=ripras(terremotospp)
#
terremotospp <- ppp(x=terremotos[,1],y=terremotos[,2],window=terremotospp2)
#
plot(terremotospp)
Prueba Chi2
Se aplica la prueba Chi-cuadrado para evalaur si la intensidad tiene un comportamiento homogeneo o no, bajo la hipótesis nula de aleatoriedad. Esta prueba se calcula bajo el siguiente estadístico:
\(\chi^{2}=\sum\frac{(observado-Esperado)^{2}}{Esperado}\)
En este caso se puede observar que de acuerdo con el valor p se rechaza la hipotesis nula de de aleatoriedad. Sin embargo esta prueba no logra definir si el patron es agregado o regular es por eso que se proponen otras pruebas como la K.
#prueba chi-cuadrado
qc.loc <- quadratcount(terremotospp, nx=5, ny=5)
plot(terremotospp, pch=3, cex=0.6)
plot(qc.loc, add=T, textargs = list(col='red'))
quadrat.test(qc.loc)## Warning: Some expected counts are small; chi^2 approximation may be inaccurate##
## Chi-squared test of CSR using quadrat counts
##
## data:
## X2 = 452.98, df = 21, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 22 tiles (irregular windows)plot(Kest(terremotospp))
Para la estimación de la intensidad se utilizó el método no parámetrico kernel, para los diferentes modelos: “gaussiano” “epanechnikov” y “quartic”. En todos se incluyó la corrección de borde de Diggle. También se variaron los anchos de banda para el modelo “gaussiano” teniendo en cuenta la opción por defecto que usa el software R a través de la función density.ppp, el método diggle, CvL, Scott y ppl.
El kernel corregido de Diggle se estima a través de la siguiente fórmula (4):
\(\widetilde{\lambda}^{(D)}(u)=\sum_i^n\frac{1}{e(x_{i})}\kappa(u-x_{i})\)
Se debe tener en cuenta que para el primer caso la determinación del ancho de banda se realiza de la siguiente forma (5):
Métodos por máxima verosimilitud para modelar la Intensidad
Para un proceso poisson homogéneo con intensidad \(\lambda\) el logaritmo de la verosimilitud es:
\(log L(\lambda;x)=n(x)log(\lambda)-\lambda area(W)\)
Para un proceso poisson inhomogéneo con función de intensidad \(\lambda_\theta(x_i)x\) dependiendo del parámetro \(\theta\). El logaritmo de la verosimilitud para \(\theta\) es:
\(log L(\lambda;x)=\sum_{i=1}^nlog\lambda_{\theta}(x_i)-\int_{W}\lambda_{\theta}(\mu)d\mu\)
donde \(\int_{W}\lambda_{\theta}(\mu)d\mu\) es el número esperado de casos del proceso poisson inhomogénero con intensidad \(\lambda_{\theta}(\mu)\) en la región W.
Kernel Adaptativo
El kernel adaptativo es una prpuesta de estimación de la intensidad tradicional la cual posee un ancho de banda “bandwidth” fijo y el cual para el caso del kernel adaptativo se considera como variable y se estima con la siguiente \[ \widehat{\lambda}_{h0} (x)=\sum_{y\in{Y}} h(y)^{-2} K(h(y)^{-1} (x-y)) =\sum_{y\in{Y}} K_{h(y)}(x-y) \]
en donde el ancho de banda variable se muestra en la siguiente expresión:
\[h(y)=h_0 n^{1/2}\lambda(y)^{-1/2}\]
##intensidad proceso aleatorio##
#intensidad espacial
intensidad_terremotos = density.ppp(terremotospp,diggle=TRUE,sigma=bw.diggle(terremotospp))
intensidad_terremotos1 = density.ppp(terremotospp,diggle=TRUE,sigma=bw.ppl(terremotospp))
intensidad_terremotos2 = density.ppp(terremotospp,diggle=TRUE,sigma=bw.ppl(terremotospp), kernel="epanechnikov")
intensidad_terremotos3 = density.ppp(terremotospp,diggle=TRUE,sigma=bw.ppl(terremotospp), kernel="quartic")
intensidad_terremotos4 = density.ppp(terremotospp,diggle=TRUE,sigma=bw.scott(terremotospp))
intensidad_terremotos5=densityAdaptiveKernel.ppp(terremotospp)En la siguiente figura se muestran las 6 diferentes formas que se presentaron para estimar la función de intensidad en donde podemos observar que algunas son muy similares y reflejan mejor los puntos observados de terremostos, dentro de las cuales se destaca la intensidad estimada por el maxima verosimulitud y scoot que reflejan bien dos nichos de sismos la meseta de todos los santos en el centro del pais y otra en el pacifico colombiano que son tradicionales de alta sismisidad. Los metodos de diggle inicial y el kernel adaptativo presentan una intensidad mas amplia y muestra mayores zonas de riesgo que en este caso no se observan muy marcados en los datos obervados.
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda (método Diggle)
par(mfrow=c(1,1))
plot(intensidad_terremotos,xlab="x",ylab="y",main=paste("Modelo Gaussiano sigma=",round(attr(intensidad_terremotos,"sigma"),5)))
contour(density(terremotospp),add=TRUE)
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda (método verosimilitud)
plot(intensidad_terremotos1,xlab="x",ylab="y",main=paste("Modelo Gaussiano sigma=",round(attr(intensidad_terremotos1,"sigma"),5)))
contour(density(terremotospp),add=TRUE)
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda y método epanechnikov
plot(intensidad_terremotos2,xlab="x",ylab="y",main=paste("Modelo epanechnikov sigma=",round(attr(intensidad_terremotos2,"sigma"),5)))
contour(density(terremotospp),add=TRUE)
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda y método quartic
plot(intensidad_terremotos3,xlab="x",ylab="y",main=paste("Modelo método quartic sigma=",round(attr(intensidad_terremotos3,"sigma")),5))
contour(density(terremotospp),add=TRUE)
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda y método scott
plot(intensidad_terremotos4,xlab="x",ylab="y",main=paste("Modelo método scott sigmas=",round(bw.scott(terremotospp),5)))
contour(density(terremotospp),add=TRUE)
#Gráfico espacial de la función de intensidad usando lkernel adaptativo
plot(intensidad_terremotos5,xlab="x",ylab="y",main="Kernel Adaptativo")
contour(density(terremotospp),add=TRUE)
Sin embargo es el metodo se scoot el que logra con base en los puntos observados reflejar mejor el patron.
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda y método scott
plot(intensidad_terremotos4,xlab="x",ylab="y",main=paste("Modelo método scott sigmas=",round(bw.scott(terremotospp),5)))
contour(density(terremotospp),add=TRUE)
points(terremotospp,col="red",pch=16)
Se observa en la intensidad temporal dos nichos muy marcados uno con en el 2000 el mas alto y otro en el 2010 periodos con gran cantidad de actividad sismica.
# Intensidad temporal (??A que hace referencia n=400,"epanechnikov" y los parametros de la segunda linea?
intensidad_temporal <- density(terremotos[,3], kernel="epanechnikov", n=400)
intensidad_temporal1 <- intensidad_temporal$y[findInterval(terremotos[,3],intensidad_temporal$x)]*dim(terremotos)[1]
#Graficos espacial y temporal de la intensidad
par(mfrow=c(1,2))
plot(intensidad_temporal,type="l",xlab="\n t = time",ylab="",main="Intensidad temporal")
plot(intensidad_terremotos4,xlab="x",ylab="y",main="Intensidades espacial")
Ahora realizamos una estimación de las intensidades para cada año y se muestra en la figura de los mapas la evolución de los ultomos 6 años en donde logramos ver como ciertos años tienen mayor intensidad de sismos.
##replicar el grafico con plotly (dinamico)
intensidad_terremotos4## real-valued pixel image
## 128 x 128 pixel array (ny, nx)
## enclosing rectangle: [-78.466, -71.128] x [1.0231, 11.79] unitsintensidad_raster=rasterFromXYZ(data.frame(intensidad_terremotos4)[,1:3])
plot(intensidad_raster)
intensidades_raster=list()
for(i in 2010:2016){
terremotos2=terremotos[terremotos[,3]==i,]
terremotospp_sub <- ppp(x=terremotos2[,1],y=terremotos2[,2],window=colombia_win)
intensidad_terremotos<-density.ppp(terremotospp_sub,diggle=TRUE,sigma=bw.scott(terremotospp_sub))
intensidad_raster=rasterFromXYZ(data.frame(intensidad_terremotos)[,1:3])
intensidades_raster[[i-2009]]=intensidad_raster
print(i-2009)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7par(mfrow=c(1,1))
##Adicionar el shape
require(rasterVis)
intensidades_raster=stack(intensidades_raster)
names(intensidades_raster)=2010:2016
levelplot(intensidades_raster[[1:6]])
A continuación realizamos la estimación de la intensidad y comparamos el kernel tradicional con el adaptativo, en donde podemos observar una diferencia muy fuerte entre ambas estimaciones ya que el adaptativo se enfoca en un nicho muy particular conocido como la mesa de todos los santos que tiene una muy alta actividad sismica en el pais.
Se ajustaron dos modelos paramétricos para estimar la intensidad, un modelo log-lineal (1) y un log-cuadrático (2). Denotados por las siguientes expresiones:
(1)\(\lambda_{\theta}(x,y)=exp(\theta_{0}+\theta_{1}x+\theta_{2}y)\)
(2)\(\lambda_{\theta}(x,y)=exp(\theta_{0}+\theta_{1}x+\theta_{2}y+\theta_{3}x^{2}+\theta_{4}y^{2})\)
#modelo log-lineal
loglin.loc <- ppm(terremotospp, ~x+y)
#modelo log-cuadrático
logquad.loc <- ppm(terremotospp, ~polynom(x, y, 2))
#predicción del modelo
par(mfrow=c(1,2),mar=c(0,0,1,2))
plot(predict(loglin.loc), main='log-linear')
points(terremotospp, pch=16, cex=0.2)
plot(predict(logquad.loc), main='log-quadratic')
points(terremotospp, pch=16, cex=0.2)
Posterior a la estimación de los modelos se revisó la bondad de ajuste de los mismos a través del análisis de residuales suavizados usando un método kernel.
En conclusión respecto a los modelos parametricos ambos presentan solo un nicho de sismos lo cual no se refleja muy bien en los datos ya que en los datos se observan al menos dos nichos. Esta puede ser una limitación de los metodos parametricos y una ventaja del kernel. Sin embargo los metodos parametricos permiten la interpretación e inferencia de sus paramtros que en una aplicación pueden llegar a ser interesantes.
par(mfrow=c(1,2),mar=c(1,1,2,2))
diagnose.ppm(loglin.loc, which = "smooth")## Model diagnostics (raw residuals)
## Diagnostics available:
## smoothed residual field
## range of smoothed field = [-4.34, 6.777]points(terremotospp, pch=16, cex=0.2)
diagnose.ppm(logquad.loc, which = "smooth")## Model diagnostics (raw residuals)
## Diagnostics available:
## smoothed residual field
## range of smoothed field = [-2.288, 4.837]points(terremotospp, pch=16, cex=0.2)
Función K
Se estimó la función K, la cual se calcula a través de la siguiente formula en el caso homogéneo (1):
\(\widehat{K}(r)=\frac{\mid W\mid}{n(n-1)}\sum_i^n\sum_j^n1 \left\{d_{ij}\leq r\right\}e_{ij}(r)\)
Para el caso inhomogéneo su cálculo se realiza a través de la siguiente función (2):
\(\widehat{K}_{inhom}(t)=\frac{1}{\mid W\mid}\sum_i^n\sum_j^n\frac{w_{y_{i}y_{j}}1 \left\{\parallel y_{i}-y_{j}\parallel\leq t\right\}}{\lambda(y_{i})\lambda(y_{j})}\)
Donde \(w_{y_{i}y_{j}}\) es el factor de corrección de borde de Ripley´s.
En las siguientes graficas se presenta la función K homogenea e Inhomogenea en donde se puede observar que existe un patron agregado ya que la linea se encuentra por encima del poisson, esto en ambos casos. Sin embargo vale la pena destacar que en el caso de patron homegeneo que es el mas apropiada la linea corta la poisson lo cual nos da información adicional que indica el tamaño promedio de los radios de los clusters. Las lineas de correción de borde presentan un comportamiento inestable sin embargo al inicio tambien muestran agregación.
par(mfrow=c(1,3))
plot(Kest(terremotospp,correction = "Ripley"), main = "K-estimada (Homogenea)",)
plot(Kinhom(terremotospp,correction = "Ripley"), main = "K-estimada (Inhomogenea)")
intensidad=density(terremotospp)
plot(Kinhom(terremotospp,lambda = intensidad_terremotos4), main = "K-estimada (Inhomogenea) - Basada en Intensidad")
Algo similar ocurre con la función de correlación por pares, la cual se calcula a través de la siguiente fórmula (1):
\(g(\xi,\eta)=\frac{\rho^{(2)}(\xi,\eta)}{\rho(\xi) \rho(\eta)}\)
Donde \(g(\xi,\eta)>1\) indica agregación o aglomeración a distancias r, si \(g(\xi,\eta)<1\) indica regularidad (1).
La función pair-correlation g tambien muestra un comportamiento agregado en ambos casos: Homogeneo e Inhomogeneo ya que la linea esta por encima de cero.
par(mfrow=c(1,3))
plot(pcf(terremotospp), main = "pair-correlation estimada (Homogenea)")
plot(pcfinhom(terremotospp), main = "pair-correlation estimada (Inhomogenea)")
plot(pcfinhom(terremotospp,lambda = intensidad_terremotos4), main = "pair-correlation estimada (Inhomogenea)- basada en intensidad")
A continuación se presentan los envelops para ambas funciones K y g con un total de 10 simulaciones y em terminos generales los resultados son iguales es decir se valida que el proceso es agregado.
# Pointwise envelopes
#bandas de bondad
opa <- par(mfrow=c(1,2)) #Graficos
ksismos <- envelope(terremotospp, Kinhom, nsim = 99, fix.n = TRUE) #99 simulaciones## Generating 99 simulations of CSR with fixed number of points ...
## 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.gsismos <- envelope(terremotospp, pcfinhom, nsim = 99, fix.n = TRUE) #99 simulaciones## Generating 99 simulations of CSR with fixed number of points ...
## 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.plot(ksismos, main = "Envelopes sismos usando K") #agregado
plot(gsismos, main = "Envelopes sismos usando g")
También se estimó la función J, esta función en el caso homogéneo se calcula de la siguiente forma (1):
\(J(r)=(1-G(r))/(1-F(r))\)
Donde \(G(r)\) es la función de la distancia del vecino más cercano y \(F(r)\) es la función de espacio vacío. Para el caso inhomogéneo, su cálculo se realiza de la siguiente forma:
\(J(r)_{inhom}=(1-G(r)_{inhom})/(1-F(r)_{inhom})\)
Donde el numerador y denominador de acuerdo con Van Lieshout corresponden con (3):
Para la función J, se observa de igual forma que la función homogénea e inhomogénea están por debajo de 1, indicando agregación.
par(mfrow=c(1,2))
plot(Jest(terremotospp), main = "J-estimada (Homogénea)")
plot(Jinhom(terremotospp), main = "J-estimada (Inhomogénea)")
Ajuste de Modelos
##Modelo Thomas
mod_thomas=kppm(terremotospp~1,"Thomas",method ="clik2")
par(mfrow=c(1,2))
plot(mod_thomas,pause=F)
summary(mod_thomas)## Stationary cluster point process model
## Fitted to point pattern dataset 'terremotospp'
## Fitted by maximum second order composite likelihood
## rmax = 1.83447149519056
## weight function: Indicator(distance <= 0.917235747595278)
## Converged successfully after 103 function evaluations
##
## ----------- TREND MODEL -----## Warning in print.summary.ppm(x$trend, ...): Model was fitted by an older version
## of spatstat## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## nd = nd, eps = eps)
## ** Executed by old spatstat version 1.64.1 **
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 232 points
## Average intensity 4.97 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 0.1146545 x 0.1682360 units
##
## Original dummy parameters: =
## Planar point pattern: 2482 points
## Average intensity 53.2 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [0.000689, 0.0193] total: 46.6
## Weights on data points:
## range: [0.000689, 0.00964] total: 1.38
## Weights on dummy points:
## range: [0.000689, 0.0193] total: 45.2
## --------------------------------------------------------------------------------
## FITTED MODEL:
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 4.978604
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.06565322 1.476472 1.733827 *** 24.44891
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept)
## 1.605149
##
## Fitted exp(theta):
## (Intercept)
## 4.978604
##
## ----------- CLUSTER MODEL -----------
## Model: Thomas process
##
## Fitted cluster parameters:
## kappa scale
## 0.25409445 0.06351616
## Mean cluster size: 19.56554 points
##
## Final standard error and CI
## (allowing for correlation of cluster process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.2965409 1.02394 2.186359 *** 5.412912
##
## ----------- cluster strength indices ----------
## Sibling probability 0.9872821
## Count overdispersion index (on original window): 20.82303
##Modelo Mattern Cluster
mod_matclus=kppm(terremotospp~1,"MatClust",method ="clik2")
par(mfrow=c(1,2))
plot(mod_matclus,pause=F)
summary(mod_matclus)## Warning in info$range(..., thresh = thresh): Argument 'thresh' is ignored for
## Matern Cluster model## Stationary cluster point process model
## Fitted to point pattern dataset 'terremotospp'
## Fitted by maximum second order composite likelihood
## rmax = 1.83447149519056
## weight function: Indicator(distance <= 0.917235747595278)
## Converged successfully after 99 function evaluations
##
## ----------- TREND MODEL -----## Warning in print.summary.ppm(x$trend, ...): Model was fitted by an older version
## of spatstat## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## nd = nd, eps = eps)
## ** Executed by old spatstat version 1.64.1 **
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 232 points
## Average intensity 4.97 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 0.1146545 x 0.1682360 units
##
## Original dummy parameters: =
## Planar point pattern: 2482 points
## Average intensity 53.2 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [0.000689, 0.0193] total: 46.6
## Weights on data points:
## range: [0.000689, 0.00964] total: 1.38
## Weights on dummy points:
## range: [0.000689, 0.0193] total: 45.2
## --------------------------------------------------------------------------------
## FITTED MODEL:
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 4.978604
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.06565322 1.476472 1.733827 *** 24.44891
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept)
## 1.605149
##
## Fitted exp(theta):
## (Intercept)
## 4.978604
##
## ----------- CLUSTER MODEL -----------
## Model: Matern cluster process
##
## Fitted cluster parameters:
## kappa scale
## 0.2734631 0.1284259
## Mean cluster size: 18.17977 points
##
## Final standard error and CI
## (allowing for correlation of cluster process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.2946705 1.027606 2.182693 *** 5.447268
##
## ----------- cluster strength indices ----------
## Sibling probability 0.9860285
## Count overdispersion index (on original window): 19.17551
##Modelo Log Gaussian Cox
mod_logauscox=kppm(terremotospp~1,clusters ="LGCP",method ="clik2")
summary(mod_matclus)## Warning in info$range(..., thresh = thresh): Argument 'thresh' is ignored for
## Matern Cluster model## Stationary cluster point process model
## Fitted to point pattern dataset 'terremotospp'
## Fitted by maximum second order composite likelihood
## rmax = 1.83447149519056
## weight function: Indicator(distance <= 0.917235747595278)
## Converged successfully after 99 function evaluations
##
## ----------- TREND MODEL -----## Warning in print.summary.ppm(x$trend, ...): Model was fitted by an older version
## of spatstat## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## nd = nd, eps = eps)
## ** Executed by old spatstat version 1.64.1 **
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 232 points
## Average intensity 4.97 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 0.1146545 x 0.1682360 units
##
## Original dummy parameters: =
## Planar point pattern: 2482 points
## Average intensity 53.2 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [-78.46637, -71.12849] x [1.023142, 11.790247] units
## (7.338 x 10.77 units)
## Window area = 46.666 square units
## Fraction of frame area: 0.591
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [0.000689, 0.0193] total: 46.6
## Weights on data points:
## range: [0.000689, 0.00964] total: 1.38
## Weights on dummy points:
## range: [0.000689, 0.0193] total: 45.2
## --------------------------------------------------------------------------------
## FITTED MODEL:
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 4.978604
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.06565322 1.476472 1.733827 *** 24.44891
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept)
## 1.605149
##
## Fitted exp(theta):
## (Intercept)
## 4.978604
##
## ----------- CLUSTER MODEL -----------
## Model: Matern cluster process
##
## Fitted cluster parameters:
## kappa scale
## 0.2734631 0.1284259
## Mean cluster size: 18.17977 points
##
## Final standard error and CI
## (allowing for correlation of cluster process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 1.605149 0.2946705 1.027606 2.182693 *** 5.447268
##
## ----------- cluster strength indices ----------
## Sibling probability 0.9860285
## Count overdispersion index (on original window): 19.17551- Moller, J., & Waagepetersen, R. P. (2003). Statistical inference and simulation for spatial point processes. CRC Press.
- Baddeley, A. J., Moller, J., & Waagepetersen, R. (2000). Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54(3), 329–350. doi:10.1111/1467-9574.00144.
- Van Lieshout, M. N. M. (2011). A J-function for inhomogeneous point processes. Statistica Neerlandica, 65(2), 183–201. doi:10.1111/j.1467-9574.2011.00482.
- Baddeley, A., Rubak, E., & Turner, R. (2015). Spatial point patterns: methodology and applications with R. CRC press. 5.Cronie, O., & Van Lieshout, M. N. M. (2018). A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions. Biometrika, 105(2), 455-462.
- Davies, T. M., & Baddeley, A. (2018). Fast computation of spatially adaptive kernel estimates. Statistics and Computing, 28(4), 937-956.