Análisis exploratorio de datos - PPP

library(spatstat)

Intensidad

Intensidad homogénea

Cuando en cualquier sub-región B bidimensional el número esperado de puntos de X que caen en B es proporcional al área de B.

El estimador de la intensidad puede ser computado usando las funciones “npoints” y “area.owin” o usando la función “intensity.ppp”:

Ejemplo 1:

La intensidad estimada es 0.0074 árboles por decímetro cuadrado.

data(swedishpines)
intensity(swedishpines)

[1] 0.007395833

unitname(swedishpines)

0.1 metres 

intensity(rescale(swedishpines))

[1] 0.7395833

plot(swedishpines)

Ejemplo 2

Si el patrón de puntos es múltiple entonces la función “intensity.ppp” retorna un vector de intensidades, una para cada tipo.

data(amacrine)
intensity(amacrine)

     off       on 
88.68302 94.92830 

unitname(amacrine)

662 microns 

summary(amacrine)

Marked planar point pattern:  294 points
Average intensity 183.6113 points per square unit (one unit = 
662 microns)

Coordinates are given to 4 decimal places

Multitype:
    frequency proportion intensity
off       142  0.4829932  88.68302
on        152  0.5170068  94.92830

Window: rectangle = [0, 1.6012085] x [0, 1] units
Window area = 1.60121 square units
Unit of length: 662 microns

sum(intensity(amacrine))

[1] 183.6113

plot(amacrine)

Conteo de cuadrantes

Una forma sencilla de verificar la falta de homogeneidad es probar si las regiones de igual área contienen aproximadamente el mismo número de puntos, como debe esperarse en un proceso puntual homogéneo.

swp = rescale(swedishpines)
Q3 = quadratcount(swp, nx=3, ny=4)
Q3

          x
y          [0,3.2) [3.2,6.4) [6.4,9.6]
  [7.5,10]       5         4         6
  [5,7.5)        8         8         6
  [2.5,5)        5         6         6
  [0,2.5)        3         5         9

plot(intensity (Q3, image=T))
plot(Q3, add=T)

plot(swp, add=T)

Test de homogeneidad

# Prueba Chi-cuadrada
#**********************************
# H0: La intensidad es homogenea
# Ha: La intensidad es no homogenea
# Si p<a entonces rechazamos H0.
#**********************************
tS = quadrat.test(swp, 3, 3)
tS

    Chi-squared test of CSR using quadrat counts

data:  swp
X2 = 4.6761, df = 8, p-value = 0.4169
alternative hypothesis: two.sided

Quadrats: 3 by 3 grid of tiles

plot(swp)
plot(tS, add=T)

• Dado que 0.41>0.05 No rechazamos la hipótesis nula, y la intensidad es homogénea.

Estimación usando kernel

La función intensidad puede ser estimada usando una función kernel. Con spatstat utilzamos la función density, por defecto, utilizamos un kernel gaussiano, con el ancho de banda de suavisado determinado por el argumento sigma. Por defecto se utiliza la corrección uniforme, y la salida es una imagen im.

D = density(swedishpines)
plot(D)
plot(swedishpines, add=T)

contour(D, add=T)

persp(D, theta= 30, phi = 30, col = "lightblue")

Sigma controla el grado del suavisado:

D05 = density(swedishpines, sigma=0.5)
D1 = density(swedishpines, sigma=1)
D15 = density(swedishpines, sigma=1.5)

Seleccionando ancho de banda bw

# Seleccionando ancho de banda bw
bw = bw.ppl (rescale(swedishpines))

Likelihood Cross-Validation criterion was minimised at right-hand end of interval [0.224, 6.93]; use argument 㤼㸱srange㤼㸲 to specify a wider interval for bandwidth 㤼㸱sigma㤼㸲

bw

   sigma 
6.931089 

plot(bw)

Cuadrantes escogidos por una covariable

Nuevos datos: bei
Un patrón de puntos que muestra la ubicación de 3605 árboles en una selva tropical. Acompañado de datos de covariables que dan la elevación (altitud) y pendiente de elevación en la región de estudio

data(bei)
class(bei)
bei

Intensidad usando kernel

data(bei)
class(bei)

[1] "ppp"

bei

Planar point pattern: 3604 points
window: rectangle = [0, 1000] x [0, 500] metres

# Intensidad usando kernel
intensity(bei)

[1] 0.007208

D = density(bei)
plot(D)
plot(bei, add=T)

contour(D, add=T, col="white")

Covariable elevación

# Covariable
elev = bei.extra$elev
class(elev)

[1] "im"

plot(elev)
plot(bei, add=T)

Dividiendo la covariable en 4 cuantiles iguales

# Dividiendo la covariable en 4 cuantiles iguales
b = quantile(elev, probs = (0:4)/4, type = 2)
Zcut = cut(elev, breaks = b, labels = 1:4)
V = tess(image=Zcut)
textureplot (V)

plot(V)

Conteo por cuadrante

# Conteo por cuadrante

There were 12 warnings (use warnings() to see them)

qb = quadratcount(bei, tess = V)
qb

tile
   1    2    3    4 
 714  883 1344  663 

Ejercicio

Realizar lo anterior con los datos de Medellín.

Planar point pattern: 1702 points
window: polygonal boundary
enclosing rectangle: [827.7456, 840.7194] x [1174.8959, 
1190.0705] units
*** 4 illegal points stored in attr(,“rejects”) ***

1. Conteo de cuadrantes

# Conteo de cuadrantes:
Q3S = quadratcount(S, nx=3, ny=4)
Q3S

             x
y             [828,832) [832,836) [836,841]
  [1186,1190]        31       344        99
  [1182,1186)       308       205       189
  [1179,1182)        98       191       121
  [1175,1179)         3       109         4

4 illegal points also plotted

# Test de homogeneidad
tS1 = quadrat.test(S, 3, 4)
tS1 

    Chi-squared test of CSR using quadrat counts

data:  S
X2 = 271.37, df = 11, p-value < 2.2e-16
alternative hypothesis: two.sided

Quadrats: 12 tiles (irregular windows)

• Dado que p es menor que 0.05, rechazamos la hipótesis nula, y por ende, la intensidad es no homogénea.

2. Estimación usando kernel

D1 = density(S)
plot(D1)
plot(S, add=T)

4 illegal points also plotted

contour(D1, add=T, col="white")