Interpreting unexpected values in kerneloverlapHR(method = "PHR") using adehabitat package in R

Question

I am trying to calculate home range overlap for multiple groups of animals using the kerneloverlaphr() function from the adehabitatHR package in R. I am hoping to use the "PHR" method; see package documentation:

The choice method="PHR" computes the volume under the UD of the animal j that is inside the home range of the animal i (i.e., the probability to find the animal j in the home range of i). That is:

As far as I can understand, this equation should produce values between 0 and 1, and should always approximate 1 when i = j. The output is a matrix of PHR values where rows are i and columns are j.

When calculating overlap for two groups, I have not faced any issues. However, with three or more groups, the diagonals often are far from 1 (I have seen both extremely low values, like in my simulated example below, and as high as 84), and the overlap values do not visibly correlate to actual overlap (as in the images provided). Often, but not always, I am met with the warning message:

Warning messages: 1: In vi * aj : longer object length is not a multiple of shorter object length

According to the answer at AdehabitatHR kerneloverlaphr, the argument same4all = TRUE should correct this but has made no difference in my experience.

The kernels seem to be properly constructed, as far as I can tell, and area estimates are reasonable. I am a bit lost as to where I've gone wrong in either preparing the data, calculating the home ranges, or calculating the overlap to receive such nonsense values. Specifically - why would I receive values other than 1 on the diagonal? and why do the estimates not visibly match the kernels?

In this first simulated example, I replicate the issue of values other than 0 for groups that do not visibly overlap:

#required packages
library(tidyverse)
library(hms)
library(adehabitatHR)
library(sp)
library(sf)

#simulate data

#set up data frame
date <- seq(from = as.Date("2023/01/01"), to = as.Date("2023/01/31"), by = "day")
time <- as_hms(seq(from = as_hms("06:00:00"), to = as_hms("17:00:00"), length.out = 500))
id <- seq(1, 3)

example <- data.frame(
  expand_grid(date, time, id)
) |> 
  arrange(id, date, time)

#simulate home ranges(lat and lon value distributions similar to those of my data)
set.seed(3)
lat_range <- list("1" = rnorm(n = 500, 
                            mean = 0.43, 
                            sd = 0.0038),
                  "2" = rnorm(n = 500,
                            mean = 0.40,
                            sd = 0.0043),
                  "3" = rnorm(n = 300,
                            mean = 0.39,
                            sd = 0.0040))
lon_range <- list("1" = rnorm(n = 500, 
                              mean = 23.01, 
                              sd = 0.0025),
                  "2" = rnorm(n = 500,
                              mean = 22.97,
                              sd = 0.0029),
                  "3" = rnorm(n = 500,
                              mean = 22.99,
                              sd = 0.0029))

example$lat <- numeric(46500)
example$lon <- numeric(46500)
for(id in seq(1, 3)){
  example$lat <- replace(example$lat, example$id == id, sample(lat_range[[id]], size = 15500, replace = TRUE))
  example$lon <- replace(example$lon, example$id == id, sample(lon_range[[id]], size = 15500, replace = TRUE))
}

#build kernel

example.sp <- data.frame(id = example$id, x = example$lat, y = example$lon)
coordinates(example.sp) <- ~x+y
proj4string(example.sp) <- CRS( "+init=epsg:4326" )
example.sp <- spTransform(example.sp, CRS("+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs"))
example.kernelref <- kernelUD(example.sp, h = "href")
example.kernel.poly <- getverticeshr(example.kernelref, percent = 95, unin = "m", unout = "km2") 

#area
example.kernel.poly$area

#calculate overlap
overlap_example <- kerneloverlaphr(example.kernelref, method = "PHR")

Here is a visualization of the three simulated homeranges (note they do not overlap visibly):

The matrix of overlap values for the three simulated groups is:

Interpreting this matrix as the documentation suggests, there is a 27% probability to find an animal from group 2 in the home range of group 1 despite the two ranges being farthest from each other.

As a second example, here is a visualization of some real home range data next to the overlap matrix:

In this instance, there is apparently a probability of 4.8 to find an animal from group 2 in the home range of group 3, approximately double the probability of finding an animal from group 2 in its own homerange (neither of these probabilities are valid in the first place, but the scaling doesn't match reality either).

Any insight is much appreciated!