top clusters

[zaima1988]

Hi Ben, earlier, I've tried to get top clusters using this syntax
top10 <- top_clusters(res, zones, k = 10, overlapping = FALSE)
top10

but the result (top10), all clusters have gumble p value = 0, and altough I set overlapping = FALSE, the result is still overlapping. And then when I read your updates to top clusters and documentation, the result of top clusters are different than first syntax and all of MLC p value is 0.01. Beside that, when I use the syntax for show subregion in top10 cluster in Flexible Zones, there was error
Error: object of type 'closure' is not subsettable
What should I do? Thank you very much

Here the First syntax
knn_mat <- coords_to_knn(unique(data[,6:7]), 12)
zones <- knn_zones(knn_mat)

t<-length(unique(data$year))
m<-length(unique(data$subregion))
counts<-matrix(data$case,nrow=t, ncol=m)
View(counts)
population<-matrix(data$population,nrow=t, ncol=m)

res <- scan_pb_poisson(counts = counts,
zones = zones,
population = population,
n_mcsim = 99,
max_only = FALSE)

res$MLC

hotspot<-unique(data$id)[res$MLC$locations]
hotspot

#TOP Cluster
top10 <- top_clusters(res, zones, k = 10, overlapping = FALSE)
top10

#show subregion in top10 cluster
j=1
clustersubregion<-list()
for(i in top10$zone){
clustersubregion[[j]]<-unique(data$id)[zones[[i]]]
j<-j+1
}
clustersubregion

Second Syntax
knn_mat <- coords_to_knn(unique(data[,6:7]), 12)
zones <- knn_zones(knn_mat)

t<-length(unique(data$year))
m<-length(unique(data$subregion))
counts<-matrix(data$case,nrow=t, ncol=m)
#View(counts)#
population<-matrix(data$population,nrow=t, ncol=m)

res <- scan_pb_poisson(counts = counts,
zones = zones,
population = population,
n_mcsim = 99,
max_only = FALSE)

res$MLC

hotspot<-unique(data$id)[res$MLC$locations]
hotspot

#tOP CLUSTER P VALUE
mc_pvalue <- function(observed, replicates) {
if (length(replicates) == 0) {
return(NULL)
} else {
f <- Vectorize(
function(y) {
(1 + sum(replicates > y)) / (1 + length(replicates))
}
)

return(f(observed))

}
}

gumbel_pvalue <- function(observed, replicates, method = "ML", ...) {
if (length(replicates) < 2) {
stop("Need at least 2 observations to fit Gumbel distribution.")
}

Fit Gumbel distribution to Monte Carlo replicates

gumbel_mu <- NA
gumbel_sigma <- NA
if (method == "ML") {
gum_fit <- gum.fit(replicates, show = FALSE, ...)
gumbel_mu <- gum_fit$mle[1]
gumbel_sigma <- gum_fit$mle[2]
} else {
gumbel_sigma <- sqrt(6 * var(replicates) / pi^2)
gumbel_mu <- mean(replicates) + digamma(1) * gumbel_sigma
}

pvalue <- pgumbel(observed, gumbel_mu, gumbel_sigma, lower.tail = FALSE)

return(list(pvalue = pvalue,
gumbel_mu = gumbel_mu,
gumbel_sigma = gumbel_sigma))
}

mtop_clusters <- function(x, zones, k = 10, overlapping = FALSE, gumbel = FALSE,
alpha = NULL, ...) {
k <- min(k, nrow(x$observed))
if (overlapping) {
return(x$observed[seq_len(k), ])
} else {
row_idx <- c(1L, integer(k - 1))
seen_locations <- zones[[x$observed[1,]$zone]]
n_added <- 1L
i <- 2L
while (n_added < k && i <= nrow(x$observed)) {
zone <- x$observed[i, ]$zone
if (zone != x$observed[i-1, ]$zone &&
length(intersect(seen_locations, zones[[zone]])) == 0) {
seen_locations <- c(seen_locations, zones[[zone]])
n_added <- n_added + 1L
row_idx[n_added] <- i
}
i <- i + 1L
}
res <- x$observed[row_idx[row_idx > 0], ]

if (nrow(x$replicates) > 0) {
  res$MC_pvalue <- mc_pvalue(res$score, x$replicates$score)
  
  if (gumbel) {
    res$Gumbel_pvalue <- gumbel_pvalue(res$score, 
                                       x$replicates$score)$pvalue
  }
  if (!is.null(alpha) && alpha >= 0 && alpha <= 1) {
    res$critical_value <- quantile(x$replicates$score, 1 - alpha)
  }
}
return(res)

}
}

top10 <- mtop_clusters(res, zones, k = 10, overlapping = FALSE, gumbel=FALSE,alpha=0.05)
top10

#show subregion in top10 cluster
j=1
clustersubregion<-list()
for(i in top10$zone){
clustersubregion[[j]]<-unique(data$id)[zones[[i]]]
j<-j+1
}
clustersubregion