### SKATER application to spatil point process ### generating data example set.seed(1) x <- c(runif(50), runif(25,0,.3), runif(25,.7,1)) y <- c(runif(50), runif(25,0,.3), runif(25,.7,1)) ids <- list(which(x<.3 & y<.3), which(x>.7&y>.7)) ids[[3]] <- setdiff(1:length(x), unlist(ids)) plot(x,y, xlim=c(0,1), ylim=c(0,1), asp=1) ### Voronoi tesselation to make neighboorhod structure require(deldir) del <- deldir(x,y, rw=c(0,1,0,1)) plot(del, add=T, wl="tes") ### the area of Dirichlet (Voronoi) ### tessellation is proportional ### to inverse of incidence ratio ### creating data.frame with influence area of each point ### points with similar area have similar incidence dpad <- scale(data.frame(del$summary)$dir.area) require(spdep) ### neighboorhod list nb <- tapply(c(del$dirs[,5], del$dirs[,6]), c(del$dirs[,6], del$dirs[,5]), as.integer) class(nb) <- "nb" nb plot(nb, cbind(x,y)) plot(del, add=T, wl="tes") segments(c(0,1,1,0), c(0,0,1,1), c(1,1,0,0), c(0,1,1,0)) ### calculing costs (weights of edges) based on euclidean ### distance between areas of neighbour points lcosts <- nbcosts(nb, cbind(1,dpad)) ### making listw nb.w <- nb2listw(nb, lcosts, style="B") ### find a minimum spanning tree mst <- mstree(nb.w) ### the mstree plot ### This tree is reasonable with incidence par(mar=c(0,0,0,0)) plot(mst, cbind(x,y), col=2, cex.lab=.7, cex.circles=0.03, fg="blue") plot(del, add=T, wl="tes") segments(c(0,1,1,0), c(0,0,1,1), c(1,1,0,0), c(0,1,1,0)) ### twoo groups with no restriction res1 <- skater(mst[,1:2], data.frame(cbind(1,dpad)), 1) plot(res1, cbind(x,y), cex.circles=.03) plot(del, add=T, wl="tes") segments(c(0,1,1,0), c(0,0,1,1), c(1,1,0,0), c(0,1,1,0)) ### more one group, with no retriction res2 <- skater(res1, data.frame(cbind(1,dpad)), 1) plot(res2, cbind(x,y), cex.circles=.03) plot(del, add=T, wl="tes") segments(c(0,1,1,0), c(0,0,1,1), c(1,1,0,0), c(0,1,1,0)) ### concordance with original groups g.ori <- integer(length(x)) for (i in 1:length(ids)) g.ori[ids[[i]]] <- i table(res2$group, g.ori) ### labeling groups in adequate label reagr <- function(gr, ref) { tab <- table(gr, ref) id <- apply(tab, 2, function(x) which.max(x/sum(x))) tab <- rbind(tab[id, ], tab[-id,]) tab <- cbind(tab, matrix(0, nrow(tab), nrow(tab)-ncol(tab))) class(tab) <- "table" tab } tab.conc <- reagr(res2$group, g.ori) tab.conc ### compute the kappa statistic require(epibasix) epiKappa(tab.conc) ### 0.876 ### Discussão: ### 1) a distância euclideana entre as áreas de influência ### é uma boa medida para fazer a árvore ### 2) há um modelo de point pattern que pode ser utilizado ### para testar se os grupos gerados são estatisticamente ### diferentes? usar isso como critério de parada.