From e4979db2553a915448c7e8a622fc2badf763c8e0 Mon Sep 17 00:00:00 2001 From: Aritz Adin Date: Fri, 8 Sep 2023 10:17:13 +0200 Subject: [PATCH] Updated to version 0.5.3 --- DESCRIPTION | 6 ++-- NAMESPACE | 2 +- NEWS | 4 +++ R/CAR_INLA.R | 10 +++--- R/STCAR_INLA.R | 10 +++--- R/divide_carto.R | 94 ++++++++++++++++++++++++------------------------ README.md | 4 +++ 7 files changed, 70 insertions(+), 60 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 110780f..0b78940 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: bigDM Type: Package Title: Scalable Bayesian Disease Mapping Models for High-Dimensional Data -Version: 0.5.2 -Date: 2023-07-19 +Version: 0.5.3 +Date: 2023-09-08 Authors@R: c(person(given = "Aritz", family = "Adin", @@ -24,7 +24,7 @@ URL: https://github.com/spatialstatisticsupna/bigDM BugReports: https://github.com/spatialstatisticsupna/bigDM/issues Depends: R (>= 4.0.0) Additional_repositories: https://inla.r-inla-download.org/R/stable -Imports: crayon, doParallel, fastDummies, foreach, future, future.apply, MASS, Matrix, methods, parallel, RColorBrewer, Rdpack, sf, spatialreg, spdep, stats, utils, rlist +Imports: crayon, doParallel, fastDummies, foreach, future, future.apply, geos, MASS, Matrix, methods, parallel, RColorBrewer, Rdpack, sf, spatialreg, spdep, stats, utils, rlist Suggests: bookdown, INLA (>= 22.12.16), knitr, rmarkdown, testthat (>= 3.0.0), tmap RdMacros: Rdpack License: GPL-3 diff --git a/NAMESPACE b/NAMESPACE index 43c1b4c..273f744 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,7 @@ importFrom(fastDummies,dummy_cols) importFrom(foreach,"%dopar%") importFrom(foreach,foreach) importFrom(future.apply,future_mapply) +importFrom(geos,geos_prepared_intersects) importFrom(grDevices,boxplot.stats) importFrom(methods,as) importFrom(methods,is) @@ -35,7 +36,6 @@ importFrom(rlist,list.flatten) importFrom(sf,st_as_sf) importFrom(sf,st_bbox) importFrom(sf,st_centroid) -importFrom(sf,st_contains) importFrom(sf,st_distance) importFrom(sf,st_drop_geometry) importFrom(sf,st_geometry) diff --git a/NEWS b/NEWS index 413d415..40397c3 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,7 @@ +version 0.5.3 +- bug fixed +- faster implementation of divide_carto() function + version 0.5.2 - 'X' argument included to STCAR_INLA() function - changes in mergeINLA() function diff --git a/R/CAR_INLA.R b/R/CAR_INLA.R index 4e16951..1430b14 100644 --- a/R/CAR_INLA.R +++ b/R/CAR_INLA.R @@ -265,8 +265,8 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X= cat(sprintf("+ Model %d of %d",d,D),"\n") - Rs <- as(Rs,"TsparseMatrix") - Rs.Leroux <- as(Rs.Leroux,"TsparseMatrix") + Rs <- as(Rs,"Matrix") + Rs.Leroux <- as(Rs.Leroux,"Matrix") # Rs <- inla.as.sparse(Rs) # Rs.Leroux <- inla.as.sparse(Rs.Leroux) @@ -308,8 +308,8 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X= W <- aux$W n <- nrow(W) - Rs <- as(Diagonal(n,colSums(W))-W,"TsparseMatrix") - Rs.Leroux <- as(Diagonal(n)-Rs,"TsparseMatrix") + Rs <- as(Diagonal(n,colSums(W))-W,"Matrix") + Rs.Leroux <- as(Diagonal(n)-Rs,"Matrix") # Rs <- inla.as.sparse(Diagonal(n,colSums(W))-W) # Rs.Leroux <- inla.as.sparse(Diagonal(n)-Rs) @@ -471,5 +471,5 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X= } } -utils::globalVariables(c("inla.as.sparse")) +# utils::globalVariables(c("inla.as.sparse")) utils::globalVariables(c("inla.setOption")) diff --git a/R/STCAR_INLA.R b/R/STCAR_INLA.R index 766e263..008242b 100644 --- a/R/STCAR_INLA.R +++ b/R/STCAR_INLA.R @@ -215,7 +215,7 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro if(temporal=="rw1") dif <- 1 if(temporal=="rw2") dif <- 2 D <- diff(diag(T), differences=dif) - Rt <- as(t(D)%*%D, "TsparseMatrix") + Rt <- as(t(D)%*%D, "Matrix") # Rt <- inla.as.sparse(t(D)%*%D) ## Define hyperprior distributions ## @@ -352,8 +352,8 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro cat(sprintf("+ Model %d of %d",d,D),"\n") - Rs <- as(Rs,"TsparseMatrix") - Rs.Leroux <- as(Rs.Leroux,"TsparseMatrix") + Rs <- as(Rs,"Matrix") + Rs.Leroux <- as(Rs.Leroux,"Matrix") # Rs <- inla.as.sparse(Rs) # Rs.Leroux <- inla.as.sparse(Rs.Leroux) S <- nrow(Rs) @@ -435,8 +435,8 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro cat("STEP 2: Fitting global model with INLA (this may take a while...)\n") W <- aux$W - Rs <- as(Diagonal(S,colSums(W))-W, "TsparseMatrix") - Rs.Leroux <- as(Diagonal(S)-Rs, "TsparseMatrix") + Rs <- as(Diagonal(S,colSums(W))-W, "Matrix") + Rs.Leroux <- as(Diagonal(S)-Rs, "Matrix") # Rs <- inla.as.sparse(Diagonal(S,colSums(W))-W) # Rs.Leroux <- inla.as.sparse(Diagonal(S)-Rs) diff --git a/R/divide_carto.R b/R/divide_carto.R index 3aa10fe..e3bc844 100644 --- a/R/divide_carto.R +++ b/R/divide_carto.R @@ -11,8 +11,9 @@ #' #' @return List of \code{sf} objects with the spatial polygons of each subdomain. #' +#' @importFrom geos geos_prepared_intersects #' @importFrom RColorBrewer brewer.pal -#' @importFrom sf st_as_sf st_contains st_intersects st_set_geometry st_union +#' @importFrom sf st_as_sf st_set_geometry st_union #' @importFrom stats aggregate #' #' @examples @@ -45,62 +46,63 @@ #' @export divide_carto <- function(carto, ID.group=NULL, k=0, plot=FALSE){ - ## Transform 'SpatialPolygonsDataFrame' object to 'sf' class - carto <- sf::st_as_sf(carto) + ## Transform 'SpatialPolygonsDataFrame' object to 'sf' class + carto <- sf::st_as_sf(carto) - ## Construct the grouped 'sf' object by ID.group variable ## - carto$temp <- seq(1,nrow(carto)) - Data <- sf::st_set_geometry(carto, NULL) + ## Construct the grouped 'sf' object by ID.group variable ## + Data <- sf::st_set_geometry(carto, NULL) + carto.group <- stats::aggregate(carto[,"geometry"], list(ID.group=Data[,ID.group]), utils::head) + D <- nrow(carto.group) - carto.group <- stats::aggregate(carto[,"geometry"], list(ID.group=Data[,ID.group]), utils::head) - D <- nrow(carto.group) + ######################## + ## Disjoint partition ## + ######################## + group.names <- carto.group$ID.group + names(group.names) <- group.names - ######################## - ## Disjoint partition ## - ######################## - carto.k0 <- vector("list",D) - names(carto.k0) <- carto.group$ID.group + carto.k0 <- lapply(group.names, function(x) { + aux <- carto[st_set_geometry(carto[,ID.group],NULL)==x,] + rownames(aux) <- NULL + aux + }) - for(i in 1:D){ - loc <- sf::st_contains(carto.group$geometry[i], carto[,"geometry"]) - carto.k0[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc") - carto.k0[[i]]$temp <- NULL - } + if(k==0){ + if(plot) lapply(carto.k0, function(x) plot(x$geometry, main=unique(st_set_geometry(x,NULL)[,ID.group]))) + return(carto.k0) + } - if(k==0){ - if(plot) lapply(carto.k0, function(x) plot(x$geometry, main=unique(st_set_geometry(x,NULL)[,ID.group]))) - return(carto.k0) - } + ############################################ + ## Partition including k-order neighbours ## + ############################################ + if(k>0){ - ############################################ - ## Partition including k-order neighbours ## - ############################################ - if(k>0){ + carto.k <- vector("list",D) + names(carto.k) <- names(carto.k0) - carto.k <- vector("list",D) - names(carto.k) <- names(carto.k0) + if(plot) color <- RColorBrewer::brewer.pal(k+2,"Blues") - if(plot) color <- RColorBrewer::brewer.pal(k+2,"Blues") + for(i in 1:D){ + aux.carto <- carto.group$geometry[i] - for(i in 1:D){ - if(plot) plot(carto.k0[[i]]$geometry, col="lightgrey", main=sort(unique(Data[,ID.group]))[i], - xlim=sf::st_bbox(carto.group$geometry[i])[c(1,3)]*c(0.99,1.01), - ylim=sf::st_bbox(carto.group$geometry[i])[c(2,4)]*c(0.99,1.01)) + for(j in 1:k){ + # loc <- sf::st_intersects(aux.carto, carto[,"geometry"]) + # carto.k[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc") + loc <- geos::geos_prepared_intersects(aux.carto, carto[,"geometry"]) + carto.k[[i]] <- carto[loc,] + rownames(carto.k[[i]]) <- NULL - aux.temp <- carto.k0[[i]]$temp - aux.carto <- carto.group$geometry[i] + if(plot){ + plot(carto.k[[i]]$geometry, col=color[j+2], main=sort(unique(Data[,ID.group]))[i], + xlim=sf::st_bbox(carto.group$geometry[i])[c(1,3)]*c(0.99,1.01), + ylim=sf::st_bbox(carto.group$geometry[i])[c(2,4)]*c(0.99,1.01)) - for(j in 1:k){ - loc <- sf::st_intersects(aux.carto, carto[,"geometry"]) - carto.k[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc") + plot(carto.k0[[i]]$geometry, col=color[1], add=TRUE) + } - if(plot) plot(carto.k[[i]]$geometry[!(carto.k[[i]]$temp %in% aux.temp)], col=color[j+2], add=TRUE) - aux.temp <- carto.k[[i]]$temp - aux.carto <- sf::st_union(carto.k[[i]]) - } - carto.k[[i]]$temp <- NULL - } + aux.carto <- sf::st_union(carto.k[[i]]) + } + } - return(carto.k) - } + return(carto.k) + } } diff --git a/README.md b/README.md index a7c3b4b..e0e874d 100644 --- a/README.md +++ b/README.md @@ -91,6 +91,10 @@ When using this package, please cite the following papers: ``` news(package="bigDM") ``` +__Changes in version 0.5.3__ (2023 Sep 08) +* small bugs fixed +* faster implementation of `divide_carto()` function + __Changes in version 0.5.2__ (2023 Jun 14) * changes in `mergeINLA()` function * 'X' argument included to `STCAR_INLA()` function