Negative and larger than 1 p-values with local_gstar_perm

[jchanovas]

I computed the local Gi star for a sf multipolygon object with 1 feature. To my surprise, I was seeing negative and very large p_values. P_sim and p_folded_sim were between 0-1 as expected. I thought there was something wrong with my code so I run the exact same code as in the Getis-Ord-Gi section of the Emerging Hot Spot Analysis vignette and also found negative p values.The res values are different from the ones in the vignette. I tried going back to older versions of sfdep and spdep and still the same issue. I couldn't replicate results. Is anybody else having the same problem? I am using R 4.3.2. Thank you so much.

[JosiahParry]

Can you please provide a reproducible example? There should be a chance that a change in spdep broke sfdep column renaming and instead of looking at p_value you're actually looking at a z-value

[jchanovas]

Thank you so much for your prompt response! Here's the code that I am using. I took it from your website. What got me worried is that the values weren't the same as the ones posted in the vignette. When I plotted the regions, the hot spot map looked different.

library(sfdep)
library(dplyr)

guerry_nb <- guerry |> 
  mutate(
    nb = include_self(st_contiguity(geometry)),
         wt = st_weights(nb)
    ) 

donat_gistar <- guerry_nb |> 
  [transmute](https://dplyr.tidyverse.org/reference/mutate.html)(gi_star = [local_gstar_perm](https://sfdep.josiahparry.com/reference/local_gstar)(donations, nb, wt, nsim = 199)) |> 
  tidyr::[unnest](https://tidyr.tidyverse.org/reference/nest.html)(gi_star)

donat_gistar

[JosiahParry]

Yeah, I was correct. Theres a new field in the output.

library(sfdep)
library(spdep)
library(dplyr)

guerry_nb <- guerry |>
  mutate(
    nb = include_self(st_contiguity(geometry)),
    wt = st_weights(nb)
  )

donat_gistar <- guerry_nb |>
  transmute(gi_star = local_gstar_perm(donations, nb, wt, nsim = 199)) |>
  tidyr::unnest(gi_star)

listw <- recreate_listw(guerry_nb$nb, guerry_nb$wt)
res <- spdep::localG_perm(guerry_nb$donations, nsim = 199, listw)

attr(res, "internals") |> 
  as.data.frame() |> 
  head()
#>            Gi       E.Gi       Var.Gi  StdDev.Gi Pr(z != E(Gi))
#> 1 0.006257415 0.01167366 1.281468e-05 -1.5130195     0.13027470
#> 2 0.010280039 0.01215929 8.851174e-06 -0.6316603     0.52760887
#> 3 0.013774212 0.01277186 8.371282e-06  0.3464380     0.72901359
#> 4 0.005676119 0.01044402 1.125149e-05 -1.4214179     0.15519530
#> 5 0.007262258 0.01181977 1.156165e-05 -1.3403478     0.18013231
#> 6 0.004947662 0.01156478 8.731848e-06 -2.2393165     0.02513533
#>   Pr(z != E(Gi)) Sim Pr(folded) Sim  Skewness     Kurtosis
#> 1               0.06          0.030 0.6253990  0.003340486
#> 2               0.63          0.315 0.6453272 -0.212173645
#> 3               0.71          0.355 0.6044298  0.164908652
#> 4               0.07          0.035 0.8004458  0.198392901
#> 5               0.12          0.060 0.8551096  0.849952943
#> 6               0.02          0.010 0.6969997  0.901873526

The error occurs here: https://github.com/JosiahParry/sfdep/blob/649b90c7ee45fa1dd80bac79a8bf3ec622b16425/R/localG.R#L44-L51

Do you feel comfortable making the adjustment here? The fix is to extract the correct fields from the result. Right now i'm extracting based on position but since there's a new field thats causing a problem. One solution would be to just add a new clean name in the place where it should be.

[jchanovas]

Thank you so much! For now, I will extract the "internals" like you did and combine it with the sf object.

[jchanovas]

I figured out how to modify the source function. This is what I did:

1) Accessed function to modify it.

trace(local_g_perm, edit = TRUE)

2) Modify local_g_perm

localg_names <- c("gi", "e_gi", "var_gi", "std_gi", "p_value", 
    "p_sim", "p_folded_sim", "skewness", "kurtosis")
  gi <- as.numeric(res)
  stats::setNames(cbind(gi = gi, as.data.frame(attr(res, "internals")[, 
    2:9])), localg_names)

Thank you so much for your help! :-)

[JosiahParry]

Would you be kind enough to submit a PR with the modification? Due to conflict of interest im no longer able to maintain the package

[jchanovas]

My pleasure @JosiahParry!