Revisión del material "Analysing spatial point patterns in R" de Adrian Baddleley, 2008.
#######################################################
### Baddley 2008. R-MPP. Procesos puntuales ###########
#######################################################
######## parte I. Perspectiva general #######
### 4. Introducción al paquete spatstat de R
# 4.4 Inspección de los datos
data(swedishpines)
X <- swedishpines
plot(X) #patrón de puntos
X #información básica
summary(X) #estudia la intensidad (densidad de puntos) en el patrón.
plot(density(X, 10)) #nos permite ver la variación local de la intensidad, graficando la estimación kernel de la intensidad. El 10 es el valor elegido para la desviación estándar del "Gaussian smoothing kernel"
contour(density(X, 10), axes = FALSE) #otra opción para el mismo gráfico, que obtiene líneas de contorno nombradas en la unidad de densidad dada
# 4.5 Análisis de datos exploratorio (EDA)
#conteo por cuadrante "quadrat counting"
Q<-quadratcount(X, nx = 4, ny = 3) #la región se divide en rectáncilos (cuadrantes) del mismo tamaño, y se cuenta el número de puntos en cada rectángulo.
plot(X); plot(Q, add = TRUE, cex = 2)
#función K de Ripley
K <- Kest(X) #
plot(K)
# 4.6 Patrón de puntos múltiple (o marcado, donde las marcas son variables categóricas)
data(lansing)
lansing
summary(lansing)
plot(lansing) #devuelve una linea donde explica el código para cada marca
plot(split(lansing)) # una forma alternativa de graficar los datos, dividiéndolo en tantos patrones puntuales como tipos de marcas
hick <- split(lansing)$hickory; plot(hick) #para extraer uno de los patrones
# 4.7 Bases de datos
# 4.8 Interfaz gráfica para dibujar un patron de puntos en el screen
X <- clickppp(10)
plot(X)
######## parte II. Tipos de datos ######
### 6 Patrones puntuales en spatstat
# 6.1 para crear un objeto de clase "ppp":
#1. guardar las coordenadas x e y de los puntos en dos vectores x e y
#2. si existen marcas asociadas a los puntos, guardarlas en un vector m.
#3. crear un objeto de patrón de puntos mediante:
# si la ventana es un rectángulo, incluir las dimensiones de la ventana rectangular (el rango de las coordenadas xrange e yrange)
ppp(x, y, xrange, yrange) #sin marcas y
ppp(x, y, xrange, yrange, marks = m) #con marcas (si las marcas son categóricas utilizar "factor(m)")
#si la ventana no es rectangular, indicamos por W el objeto ventana
ppp(x, y, window = W)
# 6.2 chequear los datos
print(X) #errores en las entradas y coordenadas
Y <- unique(X) #remover puntos duplicados
plot(X) #observar patrones no esperados y puntos fuera de la ventana
identify(X) #identificar un punto individual
# 6.3 unidades
unitname(P) <- c("metre", "metres") #si las coordenadas x e y están medidas en metros
# 6.4 otras formas de obtener procesos puntuales
#crear: "ppp"; convertir: "as.ppp"; seleccionar: "clickppp"; leer: "scanpp"; transformar:; generar por rutinas de simulación o utilizar las bases de datos.
### 7 Ventanas en spatstat
# 7.1 construir ventanas
W <- owin(xrange, yrange) #ventana rectangular
W <- disc(radius = 3, centre = c(0, 0)) #ventana circular
W <- owin(poly = p); Z <- owin(poly = list(x = c(0, 1, 0), y = c(0, 0, 1))) #ventana poligonal
Z <- owin(poly = list(list(x = c(0, 8, 0), y = c(0, 0, 8)), list(x = c(2, 2, 3, 3), y = c(2, 3, 3, 2)))) #ventana de varios polígonos
owin(mask=m, xrange, yrange) #ventana definida por una aproximación de pixel discreta
Z <- owin(poly = list(x = c(0, 1, 0), y = c(0, 0, 1))); W <- as.mask(Z) #para convertir una ventana un una máscara binaria
# 7.2 Conversión desde formato GIS
# leer los datos utilizando un paquete diseñado específicamente para el tipo de formato deseado (e.g. shapefiles para ESRI shapefiles)
# convertir datos en un formato genérico usado por "sp"
# convertir el formato genérico "sp" al formato requerido por "spatstat" usando "sp"
library(maptools); library(sp)
S <- readShapePoly("myfile.shp") #de formato "shp" a "SpatialPolygonsDataFrame"
SP <- as(S, "SpatialPolygons") #de formato "SpatialPolygonsDataFrame" a "SpatialPolygons"
W <- as(SP, "owin") #de "SpatialPolygons" a un objeto ventana.
S <- readShapePoints("myfile.shp")
SP <- as(S, "SpatialPoints")
P <- as(SP, "ppp")
# 7.3 y 7.4 ventanas
# 7.5 crear un patrón de puntos en una ventana
ppp(x, y, xrange, yrange) #point pattern in rectangle
ppp(x, y, poly=p) #point pattern in polygonal window
ppp(x, y, poly=p, xrange, yrange) #point pattern in polygonal window
ppp(x, y, mask=m, xrange, yrange) #point pattern in binary mask window
ppp(x, y, window = W) #if W is a window object (class "owin")
### 8 Manipulación de patrones puntuales
#el patrón puntual P tiene los siguientes componentes: P$n (número de puntos), P$x o P$y (vector con coordenadas x o y de los puntos), P$n largo del vector de coordenadas, P$marks (marcas), P$window (ventana de observacion)
#split(X) #divide el patrón en sub-patrones; by(X) #para hacer una operación en cada sub-patrón
#si X es un Patrón de Puntos Marcado con factores como marcas, utilizamos split(X); si Z es una imágen pixelada con factores como valores, o si Z es una teselado (mosaicos), usamos split(X,Z)
#ejemplo: para aplicar una estimación adaptativa de la intensidad en cada especie de árbol de los datos Lansing Woods
data(lansing)
V <- split(lansing); A <- lapply(V, adaptive.density)
plot(as.listof(A))
#otra forma:
data(lansing)
A <- by(lansing, FUN = adaptive.density)
plot(A)
#marcas:
m<-marks(X) #extraer las marcas de un patrón.And unmark to remove marks and the binary operator %mark% to add marks.
Y <- cut(X, breaks = 3) #divide the range of mark values into several discrete bands, yielding a point pattern with categorical marks:
#escala:
data(swedishpines); #problemas de escala
X <- affine(swedishpines, mat = diag(c(1/10, 1/10))); unitname(X) <- c("metre", "metres"); X #aumento de la escala
X <- rescale(swedishpines, 10); X #aumento de la escala
#transformaciones geométricas: rotate() #rotación 2D; shidt() #cambio a vector;
#affine() #transformaciones geométricas; rshift() #el mismo cambio aleatorio a cada punto; rjitter() #diferente cambio aleatorio a cada punto; quadratresample() #muestreo por bloque
#ejemplo:
data(nztrees)
nztrees
plot(nztrees); contour(density(nztrees, 10), axes = FALSE)
hist(nztrees$x, nclass = 25) #observamos un límite al lado derecho y un cluster encima a la derecha de la región
chopped <- owin(c(0, 148), c(0, 95)) #decidimos cortar el margen derecho a 5-pies y dedicarnos al estudio del resto de la zona
win <- nztrees$window; chopped <- trim.rectangle(win, xmargin = c(0, 5), ymargin = 0); chopped #otra forma de crear una nueva región de estudio
nzchop <- nztrees[chopped] #extraemos el patrón de puntos original que cae en nuestra nueva ventana
summary(nzchop) #estudiamos el patrón "chopped"
plot(density(nzchop, 10)); plot(nzchop, add = TRUE) #observamos un patrón más uniforme luego de excluir el margen derecho
#combinar patrones de puntos:
X <- runifpoint(20); Y <- runifpoint(10)
superimpose(X, Y)
superimpose(Hooray = X, Boo = Y) #para designar una marca diferente a cada componentes del patrón
superimpose(X, Y, W = square(2)) #para especificar una vetana para el patrón combinado
### 9. Imágenes de pixel en spatstat. clase "im"
#creación
im(mat, xcol, yrow) #a partir de datos fila, donde mat es una matriz con los valores de los pixeles, xcol e ycol las coordenadas.
#cuando la matriz mat es una función f(x,y) debemos convertirla a imagen pixelada.
f <- function(x, y) { x^2 + y^2 }
w <- owin(c(-1, 1), c(-1, 1)); Z <- as.im(f, w)
#inspección
opa <- par(mfrow = c(1, 3))
data(redwood); D <- density(redwood)
plot(D); persp(D); contour(D); par(opa)
#análisis exploratorio
#as.matrix extract matrix of pixel values from image
#cut.im convert numeric image to factor image
#hist.im histogram of pixel values
#manipulación
# X[S, drop=TRUE] The subset to be extracted is determined by the index argument S. The logical argument drop determines whether pixel values that are undefined are omitted (drop = TRUE) or returned as the value NA (drop=FALSE).
#eval.im allows us to perform pixel-by-pixel calculations on an image or on several compatible images.
#shift.im vector shift of an image
#cut.im convert numeric image to factor image
#split.im divide pixel image into sub-images
#by.im apply function to subsets of pixel image
#interp.im spatially interpolate an image
#levelset threshold an image (produces a window)
#solutionset find the region where a statement is true (produces a window)
### 10. Teselas (mosaicos). Clase "tess": división del espacio entre regiones ("tiles") que no se solapan
#tipos y creación:
tess(xgrid = xg, ygrid = yg); quadrats(W, nx, ny) #rectangular tessellations in which the tiles are rectangles with sides parallel to the coordinate axes;
tess(tiles = z) #tile lists, tessellations consisting of a list of windows, usually polygonal windows;
tess(image = Z) #pixellated tessellations, in which space is divided into pixels and each tile occupies a subset of the pixel grid.
as.tess() #para convertir desde otros tipos de datos a una tesela
#convertir desde un patrón de puntos a una tesela
X <- runifpoint(42)
plot(dirichlet(X)) #calcula la teselas de Dirichlet o la teselas de Voroni del patrón de puntos X. La tesela asociada a un punto dado del patrón X es la región del espacio que está más cercana a ese punto, respecto a otros puntos de X.
plot(delaunay(X)) #calcula la triangulación de Delaunay del patrón de puntos X.Esto realmente sería una gráfico de redes. Los puntos de X se unen si la tesela de Dirichlet tiene un borde común. Se contruye así un conjunto de triángulos no solapados que cubren todo X más que toda la ventana de X (caso anterior).
#operaciones:
X <- runifpoint(10)
V <- dirichlet(X)
U <- tiles(V) #extrae una lista de teselas
unlist(lapply(U, area.owin)) #construye una lista de ventanas
#manipulación
# cut(X,V) attaches marks to the points of X identifying which tile of V each point falls into;
# split(X,V) divides the point pattern into sub-patterns according to the tiles of V, and returns a list of the sub-patterns;
# by(X,V,FUN) divides the point pattern into sub-patterns according to the tiles of V, applies the function FUN to each sub-pattern, and returns the results as a list.
par(mfrow = c(1, 3))
X <- runifpoint(100)
plot(X)
Z <- dirichlet(runifpoint(16))
plot(Z); plot(cut(X, Z)); plot(split(X, Z))
par(mfrow = c(1, 1))
#la "intersección" entre dos teselados X e Y es el teselado cuyas teselas son las intersecciones entre las teselas de X e Y.
opa <- par(mfrow = c(1, 3))
plot(X); plot(Y); plot(intersect.tess(X, Y))
par(opa)
######## parte III: intensidad y aleatoriedad #######
### 11. Método 1: investigar la intensidad
#la intensidad es la densidad media de puntos (número esperado de puntos por unidad de área). Puede ser constante (uniforme u homogénea) o puede variar de localización en localización (no-homogénea).
#1. intensidad uniforme. aquí la densidad empírica de puntos es lambda=n(x)/area(W)
data(bei); summary(bei) #calcula la intensidad
lamb <- summary(bei)$intensity #extrae la intensidad
#2. intensidad no-homogénea. aquí estimamos la "función intensidad" lambda(u) o la "medida de densidad" delta(B) (cuando la función intensidad puede que no exista por concentraciones singulares).
#podemos estimar ambas con métodos paramétricos (apartado 13) o no-paramétricos (como el conteo por cuadrantes o el suavizado con una función núcleo -kernel smoothing-)
Q <- quadratcount(bei, nx = 6, ny = 3) #método de conteo por cuadrante
plot(bei, cex = 0.5, pch = "+"); plot(Q, add = TRUE, cex = 2)
den <- density(bei, sigma = 70) #método de estimación de la densidad de núcleo (o intensidad) mediante un núcleo gaussiano isotrópico
plot(den); plot(bei, add = TRUE, cex = 0.5)
persp(den); contour(den)
aden <- adaptive.density(bei, f = 0.01, nrep = 10) #estimador adaptativo de la intensidad, que usa la fracción f de los datos para construir un teselado de Dirichlet, y forma una estimación de intensidad que es constante en cada tesela del teselado.
plot(aden, main = "Adaptive intensity"); plot(bei, add = TRUE, cex = 0.5)
#3. cuadrantes determinados por una covariable
#el conteo por cuadrante es más útil si elegimos los cuadrantes de manera apropiada, por ejemplo, usando la información de covariables.
data(bei); Z <- bei.extra$grad #ej: partimos la región de estudio en distintas subregiones según la pendiente del terreno
b <- quantile(Z, probs = (0:4)/4)#nos da los cuantiles de los valores de pendiente.
Zcut <- cut(Z, breaks = b, labels = 1:4) #dividimos la región en 4 zonas de igual área según la pendiente del terreno.
V <- tess(image = Zcut)
plot(V); plot(bei, add = TRUE, pch = "+")
qb <- quadratcount(bei, tess = V); qb; plot(qb) #usamos el teselado para estudiar el patrón de puntos bei. Calculamos los puntos según los cuadrantes tesela.
### 12. Método 2: Probar la aleatoriedad espacial completa (CSR)
#el modelo básico de "referencia" (o modelo nulo) es el proceso puntual de Poisson uniforme (o CSR) en el plano, con intensidad lambda, que tiene las siguientes propiedades:
#el número de puntos en la región A tiene distribución Poisson con media lambda*area(A)
#dado que existen n puntos en la región A , las localizaciones de los puntos son iid y uniformemente distribuidas dentro de A
#los contenidos de dos regiones disjuntas A y B son independientes.
plot(rpoispp(100)) #simula un CSR. También se puede establecer una ventana: plot(rpoispp(100,win=W)); o simularlo condicional a un número fijo de puntos: runifpoint(100).
#test de chi-cuadrado por conteo en cuadrantes para evaluar CSR: bajo la hipótesisi nula de CSR, los nj son iid de variables aleatorias de poisson con el mismo valor esperado.
M<-quadrat.test(nzchop, nx = 3, ny = 2); M #test de Pearson chi-cuadrado. En este caso el test no refuta la hipótesis nula para los datos de árboles de Nueva Zelanda (chopped)
plot(nzchop); plot(M, add = TRUE, cex = 2) #observamos para cada cuadrante el conteo observado (arriba izquierda) y esperado (arriba derecha), así como los residuos de Pearson. (abajo)
#crítica: falta de información (existen muchos tipos de desviaciones del CSR).
#test Kolmogorov-Smirnov para evaluar CSR: comparamos las distribuciones observada y esperada de los valores de alguna función T(x,y).
KS<-kstest(nzchop, fun=function(x, y) { x }); KS #donde fun es una función(x,y). Por ej. para los datos nzchop, elegimos la función T como la coordenada x (T(X,y)=x)
plot(KS)
pval <- KS$p.value; pval
#usando datos covariables: cuando queremos saber si la intensidad de un proceso puntial depende de la covariable.
data(bei); Z <- bei.extra$grad #queremos probar formalmente si el bosque lluvioso depende de la pendiente del terreno.
b <- quantile(Z, probs = (0:4)/4) #dividimos la región en cuadrantes regulares según la pendiente.
Zcut <- cut(Z, breaks = b, labels = 1:4)
V <- tess(image = Zcut)
quadrat.test(bei, tess = V) #aplicamos el test chi-cuadrado: dado el gran número en estas regiones, podemos ignorar el tema de la independencia y concluir que los árboles no son unifomes en su intensidad.
KS <- kstest(bei, Z); KS #aplicamos el test de KS: es efectivo para covariables continuas (no para covariables factores o discretas)
plot(KS)
### Método 3: máxima probabilidad (ML) para procesos de Poisson
#si asumimos (tentativamente) que los puntos son independientes, podemos aplicar algunos métodos estadísticos para investigar la intensidad
#IPP con función de intensidad lambda(u)
lambda <- function(x, y) { 100 * (x + y) } #simulación de un IPP
plot(rpoispp(lambda))
#métodos de probabilidad:
#el log-likelihood para un HPP con intensidad lambda es logL(lambda,x)=n(x)*log(lambda)-lambda*area(W), donde n(x) es el número de puntos en el conjunto de datos x. Entonces el estimados ML de lambda es: lambda_hat=n(x)/area(W) (estimador insesgado).
#el log-likelihood de un IPP con función de intensidad lambda_theta(u) (depende del parámetro theta) es logL(theta,x)=sum(log(lambda_theta(xi)))-integral(lambda_theta(u)*du)
#si queremos probar CSR el estadístico de proporción de probabilidad (likelihood ratio test) es: R=q*log(L(theta)/L(lambda)) que es asintóticamente chi-cuadrado bajo CSR.
#Ajustar un proceso Poisson con spatstat:
#el algoritmo de Berman-Turner obtiene el MLE de theta explotando la similaridad entre el log-likelihood de Poisson y aquel para la regresión de Poisson no-lineal
#la función intensidad lambda_theta(u) debe ser log-lineal en el parámetro theta: log(lambda_theta(u))=theta*S(u) donde S(u) es uns función de valor real o vector de la localizaciónd e u.
ppm(X, ~trend) #la función de ajuste se llama "ppm" (point ptocess model) y es análoga a lm o glm. El estad´sitico S(u) se especifica por una fórmula, análoga a los modelos lm o glm.
#este corresponde a un modelo con log link, donde la fórmula ~trend especifica la forma del logaritmo de la funciónd e intensidad
ppm(bei, ~1) #HPP con intensidad uniforme
ppm(bei, ~x + y) #IPP con intensidad que es log-lineal en las coordenadas cartesianas: lambda_theta((x,y))=exp(theta0+theta1*x+theta2*y)
#las salidas contienen los coeficientes ajustados que son estimaciones por ML de los theta, los coeficientes del predictor lineal.
ppm(bei, ~polynom(x, y, 2)) #IPP con intensidad que es log-cuadrática en las coordenadas cartesianas: log(lambda_theta((x,y)))
side <- function(z) factor(ifelse(z < 500, "left", "right")); ppm(bei, ~side(x)) #para ajustar un modelo con constante pero intensidad desigual en cada sitio de una línea vertical x=500. S(u) es un factor con dos niveles (left y rigth).
#cuando se incluyen factores, la interpretación de los coeficientes depende en el "contraste" (por default se asume el "treatment contrasts", es decir, el efecto del tratamiento es cero para el primer nivel del factor, y los efectos de tratamiento estimados para otros niveles son estimaciones de las diferencias con el primer nivel. en este caso, "left" está alfabéticamente después de "rigth" por lo que el primer nivel es "left" por default.
#el modelo ajustado es: lambda_theta((x,y))= 1) exp(intercept) si x<500 o 2) exp(intercept+side(x)rigth) si x>=500.
#modelo con covariables espaciales:
#ajustar un IPP con una función intensidad que dependa de una covariable observada Z(u). Entonces Z(u) o cualquier transformación de él, puede servir como estadístico S(u) en la forma paramétrica de la función de intensidad.
data(bei); grad <- bei.extra$grad; plot(grad)
ppm(bei, ~slope, covariates = list(slope = grad)) #ajustamos un IPP con intensidad que es una función log-lineal de la pendiente: lambda(u)=exp(beta0+beta1*Z(u)), donde beta0 (intercept) y beta1 (slope) son los parámetros y Z(u) la pendiente en la localización u.
ppm(bei, ~offset(log(slope)), covariates = list(slope = grad)) #ajjustamos un IPP con intensidad que es proporcional a la pendiente: lambda(u)=beta*Z(u), con Z(u) la pendiente en u. O equivalentemente: log(lambda(u))=log(beta)+log(Z(u)). Cuando no hay coeficiente delante de log(Z(u)) el término es "offset". Entonces lambda(u)=exp(intercept)*Z(u).
#modelos ajustados
#1) operaciones estándar
#print print basic information
#summary print detailed summary information
#plot plot the fitted intensity
#predict compute the fitted intensity
#fitted compute the fitted intensity at data points
#update re-fit the model
#coef extract the fitted coefficient vector b
#vcov variance-covariance matrix of b
#anova analysis of deviance
#logLik log-likelihood value
#formula extract the model formula
#terms extract the terms in the model
#model.matrix compute the design matrix
#step stepwise model selection
#drop1 one step model deletion
#AIC Akaike Information Criterion
#ejemplo:
fit <- ppm(bei, ~x + y); fit #el modelo ajustado tiene intensidad: lambda_theta((x,y))=exp(theta0+theta1*x+theta2*y)
plot(fit, how = "image", se = FALSE)
predict(fit, type = "trend")
predict(fit, type = "cif", ngrid = 256)
coef(fit) #thetai
vcov(fit) #la diagonal nos da la var(thetai)
sqrt(diag(vcov(fit))) #desviaciones estándar
round(vcov(fit, what = "corr"), 2) #correlaciones
SE <- predict(fit, type = "se"); plot(SE, main = "standard error of fitted intensity") #error estándar de la intensidad ajustada en cada localización u. Use predict(fit, type="se") or plot(fit, se=TRUE).
fit <- ppm(bei, ~sqrt(slope) + x, covariates = list(slope = grad)) #si el modelo incluye transformaciones de las covariables originales
mo <- model.images(fit); mo; plot(mo[[2]])
fit <- ppm(bei, ~slope, covariates = list(slope = grad)) #si queremos graficar el efecto de una covariable en el modelo.
plot(effectfun(fit, "slope"))
#selección del modelo: análisis de devianza para el PP Poisson anidado
fit <- ppm(bei, ~slope, covariates = list(slope = grad)) #el primer modelo debe ser un sub-modelo del segundo
fitnull <- update(fit, ~1)
anova(fitnull, fit, test = "Chi") #realiza un likelihood ratio test con la hipótesisi nula de CSR vs. la hipótesis alternativa de IPP cn intensidad que es función log-lineal de la covariable pendiente.
#aquí el p-valor es muy pequeño, por lo que se refuta la CSR
#Como el análisis de devianza estánder no permite requiere que la hipótesis nila sea un submodelo de la alternatica, el modelo lamda(u) = beta*Z(u) (intensidad proporcional a la pendiente) no incluye un HPP como caso especial, y no podremos utilizar este análisis para testear la hipótesis nula de HPP vs. IPP.
#una solución es usar el modo de selección AIC
fitprop <- ppm(bei, ~offset(log(slope)), covariates = list(slope = grad))
fitnull <- ppm(bei, ~1)
AIC(fitprop) #un pequeño valor de AIC favorece el modelo lambda(u)=beta*Z(u)
#selección del modelo automática
X <- rpoispp(100)
fit <- ppm(X, ~x + y)
step(fit)
#simulando el modelo ajustado;
X <- rmh(fitprop)
plot(X, main = "realisation of fitted model")
### método 4: evaluar el modelo de Poisson ajustado (goodness-of-fit y validation)
#técnicas formales: basadas en supuestos probabilisticos detallados sobre los datos, permite hacer afirmaciones probabilisticas sobre los resultados: test chi-cuadrado, test goodness-of-fit, test Monte Carlo, y selección bayesiana del modelo
#técnicas informales: no imponen supuestos sobre los datos y sus interpretaciones dependen del juicio personal: análisis resudial.
#goodness-of-fit: hipótesisi nula de que el modelo es verdadero vs. hipótesisi alternativa de que el modelo es incorrecto.
#1)aplicación del test chi-cuadrado de G-O-F basada en conteo por cuadrantes al modelo ajustado de Poisson (HPP o IPP). H0: conteo por cuadrante son variabels independientes de Poisson con diferentes medias (estimadas por el modelo de ajuste)
data(bei)
fit <- ppm(bei, ~x)
M <- quadrat.test(fit, nx = 4, ny = 2)
M #aquí el test refuta el modelo ajustado
plot(bei, pch = "."); plot(M, add = TRUE, cex = 1.5, col = "red") #obtenemos una impresión informal del tipo de desvío del modelo propuesto para formular un mejor modelo.
#si el modelo huboera sido de Poisson la transformación de los residuos de Pearson (abajo) estandariza los residuos por lo cual la media=0 y varianza=1. Observamos residuos de Pearson de -14, un gran desvío del modelo ajustado.
#2)el test de KS
kstest(fit, function(x, y) { y })
#validación con residuales: si el modelo es correcto, los residuales negativos/positivos se deben cancelar
data(bei); fit <- ppm(bei, ~x + y)
plot(predict(fit)); plot(bei, add = TRUE, pch = "+") #para observar los residuales (la imagen a color indica la densidad de carga negativa)
#una forma más útil de visualizar los residuales es suavizarlos:
fitx <- ppm(bei, ~x)
diagnose.ppm(fitx, which = "smooth") #observamos una tendencia visible que sugiere que el modelo es inapropiado
#si existe una covariable espacial Z(u) podemos graficar la variable lurking de los residuales vs. Z: C(z)=R(B(z))
grad <- bei.extra$grad
lurking(fitx, grad, type = "raw") #observamos que la variable lurking comienza y termina en el eje x, mientras que para cualquier modelo con intercepto, los residuos totales para la ventana W deben ser 0. Aquí vemos que el modelo es inadecuado (la función de intensidad no captura la dependencia de la intensidad en la pendiene)
lurking(fitx, grad, type = "raw", cumulative = FALSE) #calcula la derivada de C(z) que indica qué valores de z están asociados a la falta de ajuste
#gráfico de 4 paneles: si no hay covariables espaciales, usamos el comando diagnose.ppm para graficar los residuales
data(japanesepines)
fit <- ppm(japanesepines, ~x + y)
diagnose.ppm(fit) #observamos una falta de ajuste en y=.15 y un exceso de residuales positivos cerca de x=.8 e y=.15, indicando que el modelo subestima la verdadera intensidad de los puntos en su vecindad.
####### parte IV. Interacción: how to investigate dependence between the points in a point pattern. ###########
### 15. Modelos simples de patrones que no son de Poisson: es decir, que exhiben interacción o dependencia entre puntos.
#15.1 Proceso de Poisson agregado (cluster): Y es un proceso Poisson de puntos "parentañes", donde cada uno yi causa un conjunto finito Zi de puntos progenie según cierto mecanismo estocástico. El conjunto de puntos progenie forma un proceso puntial X.
#un ejemplo de proceso cluster Matern en el cual los puntos parentales provienen de un HPP con intensidad k, y cada parental tiene un número de hijos Poison (mu), independiente y uniformemente distribuidos en un disco de radio r centrado en el parental.
plot(rMatClust(kappa = 10, r = 0.1, mu = 5))
#otros procesos Poisson cluster:
#rThomas: the Thomas process, in which each cluster consists of a Poisson(ì) number of random points, each having an isotropic Gaussian N(0, 2I) displacement from its parent.
#rGaussPoisson: the Gauss-Poisson process in which each cluster is either a single point or a pair of points.
#rNeymanScott: the general Neyman-Scott cluster process in which the cluster mechanism is arbitrary.
#15.2 Proceso Cox: es un proceso Poison con función de intensidad aleatoria delta(u).
#A trivial example is the “mixed Poisson” process in which we generate a random variable delta and, conditional on delta, generate a uniform Poisson process with intensity delta. Following are three different realisations of this process:
par(mfrow = c(1, 3))
for (i in 1:3) { lambda <- rexp(1, 1/100); X <- rpoispp(lambda); plot(X) }
par(mfrow = c(1, 1))
#el modelo Cox es el análogo al modelo de efectos aleatorios.
#A particularly interesting and useful class is that of log-Gaussian Cox processes (LGCP) in which log delta(u) is a Gaussian random function
#The Mat´ern Cluster process and the Thomas process are both Cox processes.
#Currently there are no functions in spatstat for generating the general Cox process, but if you have a way of generating realisations of a random function of interest, then you can use rpoispp to generate the Cox process. The intensity argument lambda of rpoispp can be a function(x,y) or a pixel image.
#15.3 Procesos Thinned: "thinning" significa borrar algunos puntos del patrón. Bajo "independent thinning" el destino de cada punto es independiente de otros puntos, y el proceso resultante es de Poison. Pero si se utiliza un mecanismo "dependen thinning" el proceso obtenido es no-Poison.
#In Mat´ern’s Model I, a homogeneous Poisson process Y is first generated. Any point in Y that lies closer than a distance r from the nearest other point of Y, is deleted. Thus, pairs of close neighbours annihilate each other.
plot(rMaternI(70, 0.05))
#In Mat´ern’s Model II, the points of the homogeneous Poisson process Y are marked by ‘arrival times’ ti which are independent and uniformly distributed in [0, 1]. Any point in Y that lies closer than distance r from another point that has an earlier arrival time, is deleted.
plot(rMaternII(70, 0.05))
#15.4 Modelos secuenciales: comenzamos con una ventana vacía y lso puntos son localizadas dentro de la ventana uno a la vez según algún criterio.
#In Simple Sequential Inhibition, each new point is generated uniformly in the window and independently of preceding points. If the new point lies closer than r units from an existing point, then it is rejected and another random point is generated. The process terminates when no further points can be added.
plot(rSSI(0.05, 200))
#Sequential point processes are the easiest way to generate highly ordered patterns with high intensity.
### 16. Métodos 5: métodos de distancia para patrones puntuales
#si un patrón puntual parede tener intensidad constante y queremos evaluar si es un patrón de Poisson (puntos independientes) vs. puntos dependientes (con interacción).
#def: ‘independence’ (the Poisson process), ‘regularity’ (where points tend to avoid each other), and ‘clustering’ (where points tend to be close together).
#16.1 distancias: método clásico para investigar interacción entre puntos, basados en las distancia entre ellos.
# pairdist(X) pairwise distances sij = ||xi - xj || between all distinct pairs of points xi and xj (i 6= j) in the pattern;
# nndist(X) nearest neighbour distances ti = minj6=i sij , the distance from each point xi to its nearest neighbour;
# distmap(X) empty space distances d(u) = mini ||u-xi||, the distance from a fixed reference location u in the window to the nearest data point.
data(cells)
emp <- distmap(cells)
plot(emp, main = "Empty space distances")
plot(cells, add = TRUE)
#Quite a useful exploratory tool is the Stienen diagram obtained by drawing a circle around each data point of diameter equal to its nearest neighbour distance:
plot(X %mark% (nndist(X)/2), markscale = 1, main = "Stienen diagram")
#16.2 distancias de espacio vacío: efectos de borde y funciónd e espacio vacío F
#la distribución empírica de las distancias del espacio vacío depende de la geometría de la ventana W y de las características del proceso puntual X. =problemas de borde
#ingnorando el problema de borde, sea X estacionaria, definimos la función de distribución acumulativa de la distancia de espacio vacío (ecdf): F(r)=P(d(u,X)<=r).
#utilizando correcciones del efecto de borde, obtenemos una versión ponderada del ecdf: F(r)=sum(e(uj,r)*1(d(uj,x)z=r)) donde e(uj,r) es el peso de la corrección del borde. Para un proceso Poisson F(r)=1.exp(-lambda*phi*r^2)
data(cells)
plot(cells)
Fc <- Fest(cells)
Fc
par(pty = "s"); plot(Fest(cells));
plot(Fest(cells), hazard ~ r, main = "Hazard rate of F") #grafica la tasa hazard vs. r
plot(Fest(cells), cbind(km, rs, raw, theo) ~ r) #grafica la estimación de F(r) incluyendo la estimación no corregida
plot(Fest(cells), cbind(km, rs, theo) ~ theo) #grafica la estimación de F(r) vs. valor de Poisson
plot(Fest(cells), . - theo ~ r) #restando el valor teórico de Poisson a las estimaciones
plot(Fest(cells), asin(sqrt(.)) ~ r) #aplicando la varianza de Fisher estableciendo la transformación phi(F(t))=sin^-1(sqrt(F(t)))
#16.3 distancias al vecino más cercano: función G
#asumiendo que X es estacionario podemos definir la cdf de la distancia al vecino más cercano para un punto: G(r)=P(d(u,X/{u})<=r). La df empírica es G*(r)=1/n(x)*sum(1[ti<=r]).
#para correcciones de borde, tenemos la versión sopesada de ecdf: G(r)=sum(e(xi,r)*1{ti<=r}), donde para HPP de intensidad lambda, G(r)=1-exp(-lambda*phi*r^2)
#interpretación (G(r) se interpreta al revés que F(r)): G(r)>G(r) Poisson sugiere que las distancias al vecino más cercano en el punto son menores que las del PP, por lo cual se sugiere un patrón agrupado, mientras que G(r)<G(r) Poisson sugiere un patrón regular (inhibición).
Gc <- Gest(cells); Gc
par(pty = "s"); plot(Gest(cells)) #aquí se sugiere un patrón regular y que no existen distancias al vecino más cercano para distancias menores de 0.07.
#otros:
#G(r) and Gpois(r) plotted against r plot(Gest(X))
#bG(r) - Gpois(r) plotted against r plot(Gest(X), . - theo ~ r)
#bG(r) plotted against Gpois(r) in P–P style plot(Gest(X), . ~ theo)
fisher <- function(x) { asin(sqrt(x)) } #transformación con la varianza de Fisher
plot(Gest(cells), fisher(.) ~ fisher(theo))
#16.4 distancias pareadas y la función K
#Ripley definió la función K para el PP estacionario tal que lambda*K(r) es el número esperado de otros puntos del proceso dentro de la r a un punto. Para un HPP K(r)=phi*r^2. Hemos de comparar la estimación de K(r) con la función K de Poisson. Valores K(r)>phi*r^2 sugieren agrupamiento mientras que K(r)<phi*r^2 sigieren àtrón regular.
Gc <- Kest(cells); Gc
par(pty = "s"); plot(Kest(cells))
#para nuestro ejemplo, la interpretación de los 3 estadísticos de resumen F, G y K es la misma: evidencia de patrón regular.
#una transformación de K que es muy utilizada es la función L: L(r)=qsqrt(K(r)/phi), que transforma la función K de Poisson a una línea recta L(r)=r.
L <- Lest(cells)
plot(L, main = "L function")
#otra función de resumen es la funciónd e correlación pareada: g(r)=K`(r)/2*phi*r
plot(pcf(cells))
#16.5 función J: una combinación útil de F y G: J(r)=(1-G(r))/(1-F(r)). Para HPP FPois=Gpois, por lo que JPois(r)=1, valores J(r)>1 sugieren regularidad y <1 agrupamiento.
plot(allstats(cells)) #calcula F,G, J y K
#16.6 manipulación y gráfico de funciones de resumen: with, eval.fv
data(redwood)
K <- Kest(redwood)
y <- with(K, iso - theo) #contiene la diferencia entre las columnas iso (estimación de corrección isotrópica) y theo (valor teórico para CSR) de la estimación función K
x <- with(K, r) #contiene los valores del argumento distancia r de la estimación función K
K <- Kest(redwood)
L <- eval.fv(sqrt(K/pi))
#16.7 advertencia (Caveats)
#1. las funciones F, G y K se definen y estiman bajo el supuesto de que el proceso puntual es estacionario (homogéneo)
#2. estas funciones de resumen no caracterizan totalmente el proceso
#3. si el proceso no es estacionario, las desviaciones entre las funciones empírica y teórica (e.g. K y Kpois) no son necesariamente una evidencia de la interacción entre puntos, sino que pueden ser atribuibles a las variaciones en la intensidad.
### 17. métodos 6: simulación envolvente (simulation envelopes) y test de bondad de ajuste (goodness-of-fit)
#también es posible usar estadísticos de resumen como base para la inferencia estadística.
#17.1 test envolvente y test Monte Carlo
#debido a variabilidad aleatoria nunca obtendremos una concordancia perfecta entre K y Kpois.
#para concluir formalmente que existe una diferencia significativa entre K y Kpois usamos el test de hipótesis: H0: los datos son realización de CSR, H1: los datos son una realización de otro proceso puntoal no-especificado.
#el test Monte Carlo es un test basado en las simulaciones de la Ho.
#supongamos la referencia de una función K teórica para el CSR, generamos M simulaciones independientes del CSR dentro del a región de estudio W. Calculamos las funciónes K estimadas para cada realización y obtenemos: L(r)=min(Kj(r)), U(r)=max(Kj(r)); j=1,...,M.
#el test que refuta la Ho de un PP uniforme cuando K(r)cae fuera de [L(r),U(r)], tiene nivel de significación alfa=2/(M+1).
#In spatstat the function envelope computes the pointwise envelopes.
data(cells)
E <- envelope(cells, Kest, nsim = 39, rank = 1); E
plot(E, main = "pointwise envelopes") #por ej. si R=.10 fijo refutamos la Ho de CSR al nivel 5%. El valor M=39 es el menor para producir un test two-sided con nivel significativo al 5%.
#we can construct simultaneous critical bands which have the property that, under H0, the probability that bK ever wanders outside the critical bands is exactly 5%.
#compute, for each estimate bK(r), its maximum deviation from the Poisson K function, D = maxr | bK(r) - Kpois(r)|: L(r)=phi*r^2-Dmax; U(u)=phi*r^2+Dmax.
#The estimated K function for the data transgresses these limits if and only if the D-value for the data exceeds Dmax. Under H0 this occurs with probability 1/(M + 1). Thus, a test of size 5% is obtained by taking M = 19.
E <- envelope(cells, Kest, nsim = 19, rank = 1, global = TRUE)
plot(E, main = "global envelopes")
#A more powerful test is obtained if we (approximately) stabilise the variance, by using the L function in place of K.
E <- envelope(cells, Lest, nsim = 19, rank = 1, global = TRUE)
plot(E, main = "global envelopes of L(r)")
#Monte Carlo testing rationale can be applied to any point process model serving as a null hypothesis. We simply have to generate simulated realisations from the null hypothesis, and compute the summary function for each simulated realisation.
data(bei)
fit <- ppm(bei, ~elev + grad, covariates = bei.extra)
E <- envelope(fit, Lest, nsim = 19, global = TRUE, correction = "border")
plot(E, main = "envelope for inhomogeneous Poisson")
#Envelopes can also be computed using any user-specified procedure to generate the simulated realisations. This allows us to perform randomisation tests, for example.
sim <- expression(rpoispp(100))
#This expression should be passed to the envelope function as the argument simulate. The following code generates simulation envelopes for the L function based on simulations of CSR which have the same number of points as the data pattern.
data(cells)
e <- expression(runifpoint(cells$n, cells$window))
E <- envelope(cells, Lest, nsim = 19, global = TRUE, simulate = e)
plot(E, main = "envelope with fixed n")
#Envelopes can also be computed from a user-supplied list of point patterns, instead of the simulated point patterns generated by a chosen simulation procedure. This improves efficiency and consistency if, for example, we are going to calculate the envelopes of several different summary statistics.
data(cells)
SimPatList <- list()
for (i in 1:1000) SimPatList[[i]] <- runifpoint(cells$n)
EK <- envelope(cells, Kest, simulate = SimPatList, nsim = 1000)
Ep <- envelope(cells, pcf, simulate = SimPatList, nsim = 1000)
### 18. Método 7: Ajuste de modelo utilizando resúmenes estadísticos
#los estadísticos de resumen pueden ser utilizados para ajustar los modelos a los datos.
#ajustamos asemejando el estadístico de resumen (eg la función K) a su valor teórico bajo el modelo.
#cluster processes: ajusta un PP cluster (modelo Thomas) por el método de mínimo contraste.
data(redwood)
fit <- kppm(redwood, ~1, "Thomas"); fit
plot(fit) #The plot shows the theoretical K function of the fitted Thomas process (fit), three nonparametric estimates of the K function (iso, trans, border) and the Poisson K function (theo).
fitM <- kppm(redwood, ~1, "MatClust") #fit the Mat´ern cluster process to the redwood data,
plot(simulate(fit, nsim = 4))
#otros modelos con función K conocida: también disponen de métodos de contrastes mínimos.
fit <- lgcp.estK(redwood, c(sigma2 = 0.1, alpha = 1)); fit #ajusta un poroceso Cox log-Gaussiano con función de covarianza exponencial
#algoritmos genéricos para mínimos contrastes
mincontrast(observed, theoretical, starpar)
### 19. Métodos 8: ajustando para la no-homogeneidad
#si se conoce o sospecha que un patrón puntual sea espacialmente no-homogéneo, entonces nuestro análisis estadístico del patrón debe tener en cuenta esta no-homogeneidad.
#19.1 función K no-homogénea: para un IPP con función de intensidad lambda(u), la función K no-homogénea es Kinhom,pois(r)=phi*r^2.
data(bei)
fit <- ppm(bei, ~elev + grad, covariates = bei.extra)
lam <- predict(fit, locations = bei)
Ki <- Kinhom(bei, lam) #obtenemos la estimación de intensidad lambda(u) ajustando un modelo paramétrico para eliminar el sobreajuste.
plot(Ki, main = "Inhomogeneous K function") #luego de dat cuenta de la dependencia de la altitud y pendiente, los árboles aún parecen estar agrupados (clustered)
Ki2 <- Kinhom(bei) #la función de intensidad lambda(u) también puede ser estimada por kernel smoothing del patrón de puntos con el punto xi eliminado.
plot(Ki2, main = "Kinhom using leave-one-out")
g <- pcf(Kinhom(bei)) # calculamos la función de correlación pareada no-homogénea
plot(g)
#19.2 modelos cluster no-homogéneos: podemos incluir la no-homogeneidad en el proceso parental o el proceso progenie.
Z <- as.im(function(x, y) { 6 * exp(2 * x - 1) }, owin())
plot(rMatClust(10, 0.05, Z))
#19.3 ajuste de modelo no-homogéneo por cntraste mínimo
data(bei)
fit <- kppm(bei, ~elev + grad, "Thomas", covariates = bei.extra); fit
### 20. Modelos de Gibbs: One way to construct a statistical model (in any field of statistics) is to write down its probability density.
#p132/199
### 21. Métodos 9: ajuste de modelos Gibbs
### 22. Métodos 10. validazión de los modelos Gibbs
######## parte V. Patrones puntuados marcados ##################
### 23. Patrones puntuales marcados: las marcas pueden ser variables continuas o categóricas.
#los datos deben ser tratados como MPP? en un MPP los puntos son aleatorios. Tratar los datos como PP es inapropiado si las localizaciones son fijas, o si las localizaciones no son parte de la respuesta.
#cuando existen marcas, hay más elecciones para el análisis. un modelo estad´sitico para el MPP puede ser formulado en distintas formas: condicional en localizaciones, condicional en marcas, conjunto (joint).
#estas aproximaciones conducen a diferentes modelos estadísticos y tienen interpretaciones inferenciales diferentes; y existen diferentes hipótesisi nulas que pueden ser testadas:
# random labelling: given the locations X, the marks are conditionally independent and identically distributed;
# independence of components: the sub-processes Xm of points of each mark m, are independent point processes;
# complete spatial randomness and independence (CSRI): the locations X are a uniform Poisson point process, and the marks are independent and identically distributed. (This implies both random labelling and independence of components).
### 24. Manejo de datos MPP
#creación:
#A marked point pattern dataset can be created using any of the following tools:
#ppp create point pattern dataset
#as.ppp convert other data to point pattern
#superimpose combine several point patterns into a marked point pattern
#marks extract marks from a point pattern
#marks<- attach marks to a point pattern
#%mark% attach marks to a point pattern
#unmark delete marks from a point pattern
#scanpp read point pattern data from text file
#clickppp create a pattern using point-and-click on the screen
#inspección:
data(amacrine) ; amacrine
summary(amacrine)
plot(amacrine)
as.data.frame(amacrine) #convertir el MPP a un data frame
data(longleaf)
m <- marks(longleaf) #extracción de marcas
#manipulación:
data(lansing); d <- marks(lansing) #convierte marcas desde diámetros a áreas circulares
a <- (pi/4) * d^2
marks(lansing) <- a
data(amacrine); Y <- split(amacrine) #separar puntos de diferentes tipos
plot(split(amacrine))
#cortando la escala numérica en bandas
data(longleaf); longleaf #convierte marcas numéricas a factores
X <- cut(longleaf, breaks = c(0, 30, 80), labels = c("juvenile", "adult")); X
par(mfrow = c(1, 2)); plot(longleaf); plot(X, main = "cut(longleaf)"); par(mfrow = c(1, 1))
### 25. Métodos 11: herramientas exploratorias para MPP (para PP múltiples, i.e. donde las marcas son categóricas)
#25.1 Intensidad
data(lansing); summary(lansing)
plot(split(lansing)) #examina sub-patrones de diferentes tipos
plot(density(split(lansing)), ribbon = FALSE) #grafica estimaciones de la intensidad para cada tipo de manera separada.
Y <- density(split(lansing))
attach(Y)
pBlackoak <- eval.im(blackoak/(blackoak + hickory + maple + misc + redoak + whiteoak)) #proporciones relativas de la intensidad
plot(pBlackoak)
detach(Y)
#25.2 Marcas numéricas: distribución y tendencia.
data(longleaf)
hist(marks(longleaf))
plot(smooth.ppp(longleaf)) #una forma de evaluar la tendencia espacial en las marcas es formar el kernel regression smoother: m(u)=sum(mi*k(u-xi))/sum(k(u-xi)), donde k es el smoothing kernel y mi es la marca en el punto i.
data(spruces)
plot(split(cut(spruces, breaks = 3))) #miramos la no-homogeneidad espacial para las marcas
#25.3 Resumenes simples de marcas en el vecindario.
data(amacrine)
M <- marktable(amacrine, R = 0.1); M[1:10, ] #compila la tabla de contingencia de las marcas para todos los ptos dentro del radio dado a cierto punto.
md <- markstat(longleaf, mean, N = 5); md[1:10] #calcula el diámetro promedio de los 5 vecinso más cercanos a cada árbol.
#25.4 Funciones de resumen: F, G, J y K (y L)
#Assume the multitype point process X is stationary. Let Xj denote the sub-pattern of points of type j, with intensity j . Then for any pair of types i and j,
# Fj(r) is the empty space function for Xj .
# Gij(r) is the distribution function of the distance from a point of type i to the nearest point of type j
# Kij(r) is 1/ j times the expected number of points of type j within a distance r of a typical point of type i.
# Lij(r) is the corresponding L-function
# gij(r) is the corresponding analogue of the pair correlation function
# Jij(r)
#The functions Gij ,Kij ,Lij , gij , Jij are called “cross-type” or “i-to-j” summary functions.
data(amacrine); amacrine
plot(Gcross(amacrine, "on", "off")) #para las demás usar Gcross, Kcross, Lcross, pcfcross y Jcros
#The interpretation of the cross-type summary functions is similar, but not identical, to that of the original functions F, G, K etc:
#• if Xj is a uniform Poisson process (CSR), then Fj(r) = 1 - exp(- j r2).
#• if Xj is a uniform Poisson process (CSR) and is independent of Xi, then Gij(r) = 1 -exp(- j r2).
#• if Xi and Xj are independent, then Kij(r) = r2 and so Lij(r) = r and gij(r) = 1.
#• if Xi and Xj are independent, then Jij(r) = 1.
#Here ‘independent’ means that the two point processes are probabilistically independent.
#The command alltypes enables the user to compute the cross-type summary functions between all pairs of types simultaneously.
plot(alltypes(amacrine, "G"))
data(lansing)
a <- alltypes(lansing, "G")
plot(a[2:3, 2:3])
#Also defined are the “i-to-any” summaries
# Gi•(r), the distribution function of the distance from a point of type i to the nearest other point of any type;
# Ki•(r) is 1/ times the expected number of points of any type within a distance r of a typical point of type i. Here = Pj j is the intensity of the entire process X.
# Li•(r) is the corresponding L-function
# Ji•
plot(Gdot(amacrine, "on")) #o Gdot, Kdot, Ldot y Jdot respectivamente,
plot(alltypes(amacrine, "Gdot"))
#gráfico
aG <- alltypes(amacrine, "G")
fisher <- function(x) asin(sqrt(x))
plot(aG, fisher(.) ~ fisher(theo))
### 25.5 Función de correlación de marcas: rho_j(r). Para un MPP estacionario Y, rho(r) es una medida de la dependencia entre las marcas de dos puntos del proceso que están a distancia r: rho_j(r) =E[f(M1,M2)]/E[f(M,M´)], donde f es una función que se elige según:
#for continuous real-valued marks, f(m1,m2) = m1m2;
# for categorical marks (multitype point patterns), f(m1,m2) = 1 {m1 = m2};
# for marks taking values in [0, 2 ], f(m1,m2) = sin(m1 - m2).
#The value 1 suggests glack of correlationh: under random labelling, f (r) ß 1. The interpretation of values larger or smaller than 1 depends on the choice of function f.
data(amacrine)
eqfun <- function(m1, m2) { m1 == m2 }
M <- markcorr(amacrine, eqfun, correction = "translate", method = "density", kernel = "epanechnikov")
plot(M)
### 25.6 test de aleatorización
#Ho Poisson: la Ho de HMPP puede ser probada por simulación directa.
data(amacrine)
E <- envelope(amacrine, Kcross, nsim = 39, i = "on", j = "off") #utiliza la función K de tipo cruzado como el test estadístico
plot(E, main = "test of marked Poisson model")
#independencia de componentes: otras Ho pueden ponerse a prueba por test de aleatorización.
#random labelling: given the locations X, the marks are conditionally independent and identically distributed;
#independence of components: the sub-processes Xm of points of each mark m, are independent point processes.
#In a randomisation test of the independence-of-components hypothesis, the simulated patterns X are generated from the dataset by splitting the data into sub-patterns of points of one type, and randomly shifting these sub-patterns, independently of each other.
E <- envelope(amacrine, Kcross, nsim = 39, i = "on", j = "off", simulate = expression(rshift(amacrine, radius = 0.25)))
plot(E, main = "test of independent components")
#Under the independence hypothesis: Kij(r) = phi*r^2; Gij(r) = Fj(r); Jij(r)=1.
#"random labelling":
#the simulated patterns X are generated from the dataset by holding the point locations fixed, and randomly resampling the marks, either with replacement (independent random sampling) or without replacement (randomly permuting the marks).
#Under random labelling: Ji•(r) = J(r);Ki•(r) = K(r);Gi•(r) = G(r). we would normally use something like Ki•(r) - K(r) to construct a test statistic for random labelling.
Jdif <- function(X, ..., i) {
Jidot <- Jdot(X, ..., i = i)
J <- Jest(X, ...)
dif <- eval.fv(Jidot - J)
return(dif)
}
E <- envelope(amacrine, Jdif, nsim = 39, i = "on", simulate = expression(rlabel(amacrine)))
plot(E, main = "test of random labelling") #parece aceptarse la H de random labelling.
#"arrays of evelopes": To compute a simulation envelope for the function Kij for each pair of types i and j
aE <- alltypes(amacrine, Kcross, nsim = 39, envelope = TRUE)
plot(aE, sqrt(./pi) - r ~ r, ylab = "L(r)-r")
### 26. Métodos 16: modelos múltiples de Poisson
### 26.1 Teoría:
#HPP marcado (o CSRI): randomly marked Poisson process (Poisson [X], iid [M|X]) o superposition of independent Poisson processes (iid [M], Poisson [X|M]) o Poisson marked point process (jointly Poisson [X,M]).
#IPP marcado Y con intensidad conjunta lambda(u,m) para localización u y marca m, es un IPP en R^2xM con itensidad lambda(u,m). Cuando las marcas son categóricas, el proceso tiene las siguientes propiedades:
#localizaciones X (sin marcas) es IPP con intensidad beta(u)
#condicional a las localizaciones X, la smarcas son independientes.
#el sub-proceso Xm de puntos con marca m es IPP con intensidad betam(u)=lambda(u,m).
#el su-proceso Xm de puntosa con diferentes marcas m son procesos independientes.
### 26.2 Simulación:
par(mfrow = c(1, 2))
Xunif <- rmpoispp(100, types = c("A", "B"), win = square(1))
plot(Xunif, main = "CSRI, intensity A=100, B=100")
Xunif <- rmpoispp(c(100, 20), types = c("A", "B"), win = square(1))
plot(Xunif, main = "CSRI, intensity A=100, B=20")
par(mfrow = c(1, 1))
X1 <- rmpoispp(function(x, y, m) { 300 * exp(-3 * x) }, types = c("A", "B"))
lamb <- function(x, y, m) { ifelse(m == "A", 300 * exp(-4 * x), 300 * exp(-4 * (1 - x))) }
X2 <- rmpoispp(lamb, types = c("A", "B"))
par(mfrow = c(1, 2))
plot(X1, main = "")
plot(X2, main = "")
par(mfrow = c(1, 1))
### 26.3 Ajuste de modelos de Poisson
#la densidad de probabilidad de un MPP es una función f(y) definida para todo MPP "y" incluido el patrón evacío.
#si f(y)=1 es un UMPP con intensidad 1 para cada marca.
#si f(y)=exp(sum(1-betam)|W|)prod(betami) es un UMPP con intensidad lambda(u,m)=betam
#si f(y)=exp(sum(integral(1-lambda(u,m)du)))prod(lambda(xi,mi)) es un IMPP con función de intensidad lambda(u,m)
#"Maximum likelihood": logL=sum(log(lambda(xi,mi)))-sum(integral(lambda(u,m)du)). Es equivalente al logL de una regersión Poisson log-lineal, por lo que el algoritmo de Berman-Turner puede ser utilizado para maximizar el logL.
#ajuste:
ppm(X, ~marks) #HMPP (CSRI). La tendencia depende solo de las marcas, no de la localización espacial. Como las marcas es un factor, la tendencia tiene valores constantes separados para cada nivel de la marca.
ppm(X, ~1) #CSRI donde la intensidad betam=alfa, es decir, todas las marcas tienen igual probabilidad.
ppm(lansing, ~marks) #HMPP con marcas como factores, con intensidad betam para puntos de marca m.
ppm(lansing, ~marks + x) #MPP con función intensidad que satisface: log(lambda(x,y,m))=alfam+beta*x
ppm(lansing, ~marks * x) #el símbolo "*" indica interacción en el sentido usual de los modelos lineales. El modelo ajustado es un MPP con log(lambda(x,y,m))=alfam+betam*x. La intensidad es log-lineal en x con una pendiente e intercepto diferente para cada marca.
#### 27. Métodos 13: modelos de Gibbs para patrones puntuales múltiples.
p185
### 28. Datos de segmento de línea
Comentarios
Publicar un comentario