library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(spdep) # for spatial autocorrelation tests
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(sp)
load('DB_RSdata.Rda')
data_M <- DB_RSdata %>% filter(sex=="M")
territorial definition: ifelse(recap>3 &
call>1, 1, 0)home.X and home.Y: average of all capture
locationscall.X and call.Y: average of all calling
locationscolor: B, IB, I, IR, Rcolor2: B, I, R ([IB, I, IR] -> I)color_num: 0,1,2,3,4 based on colorcolor2_num: 0, 1, 2 based on color2data_sp <- data_M %>% select(home.X, home.Y, color_num) %>%
na.omit() %>% # remove those with NA coordinates
mutate(home.X = jitter(home.X), home.Y = jitter(home.Y)) #add some noise to avoid overlapping points
ggplot(data_sp, aes(x = home.X, y = home.Y, color = color_num)) +
geom_point() +
scale_color_gradient(low = "cornflowerblue", high = "indianred") +
labs(title = "Male locations", x = "X Coordinate", y = "Y Coordinate", color = "Color score") +
theme_bw(base_size = 15) +
coord_fixed()
data_sp <- data_M %>% filter(territorial == 1) %>%
select(call.X, call.Y, color_num) %>%
na.omit() %>% # remove those with NA coordinates
mutate(call.X = jitter(call.X), call.Y = jitter(call.Y)) #add some noise to avoid overlapping points
ggplot(data_sp,
aes(x = call.X, y = call.Y, color = color_num)) +
geom_point() +
scale_color_gradient(low = "cornflowerblue", high = "indianred") +
labs(title = "Territorial male locations", x = "X Coordinate", y = "Y Coordinate", color = "Color score") +
theme_bw(base_size = 15) +
coord_fixed()
# Compute pairwise distances using the Euclidean distance metric
distances <- data_M %>% select(home.X, home.Y, color2) %>%
dist(method = "euclidean") %>% as.matrix
## Warning in dist(., method = "euclidean"): NAs introduced by coercion
# Set the diagonal of the distance matrix to a large value to exclude self-distance
diag(distances) <- Inf
# Find the index of the nearest neighbor for each data point
nearest_neighbor_indices <- apply(distances, 1, which.min)
# Extract the colors of the nearest neighbors
nearest_neighbor_colors <- data_M$color2[nearest_neighbor_indices]
# Create a data frame with the original data, the nearest neighbor coordinates, and colors
data_NN <- data_M %>%
mutate(NN.X = data_M$home.X[nearest_neighbor_indices],
NN.Y = data_M$home.Y[nearest_neighbor_indices],
NN.color2 = nearest_neighbor_colors)
not really different in terms of proportion
data_NN %>% select(color2) %>% table()%>% prop.table
## color2
## B I R
## 0.1184211 0.5394737 0.3421053
data_NN %>% select(color2, NN.color2) %>% table() %>% prop.table(1)
## NN.color2
## color2 B I R
## B 0.16666667 0.55555556 0.27777778
## I 0.09756098 0.56097561 0.34146341
## R 0.13461538 0.57692308 0.28846154
# Compute pairwise distances using the Euclidean distance metric
distances <- data_M %>% filter(territorial==1) %>%
select(call.X, call.Y, color2) %>%
dist(method = "euclidean") %>% as.matrix
## Warning in dist(., method = "euclidean"): NAs introduced by coercion
# Set the diagonal of the distance matrix to a large value to exclude self-distance
diag(distances) <- Inf
# Find the index of the nearest neighbor for each data point
nearest_neighbor_indices <- apply(distances, 1, which.min)
# Extract the colors of the nearest neighbors
nearest_neighbor_colors <- data_M$color2[nearest_neighbor_indices]
# Create a data frame with the original data, the nearest neighbor coordinates, and colors
data_NN <- data_M %>% filter(territorial==1)%>%
mutate(NN.X = data_M$call.X[nearest_neighbor_indices],
NN.Y = data_M$call.Y[nearest_neighbor_indices],
NN.color2 = nearest_neighbor_colors)
sample size a bit too small
data_NN %>% select(color2) %>% table()%>% prop.table
## color2
## B I R
## 0.15 0.55 0.30
data_NN %>% select(color2, NN.color2) %>% table() %>% prop.table(1)
## NN.color2
## color2 B I R
## B 0.0000000 0.5000000 0.5000000
## I 0.1818182 0.5909091 0.2272727
## R 0.1666667 0.4166667 0.4166667
not different from random distribution
data_sp <- data_M %>% select(home.X, home.Y, color_num) %>%
na.omit() %>% # remove those with NA coordinates
mutate(home.X = jitter(home.X), home.Y = jitter(home.Y)) #add some noise to avoid overlapping points
# Create a spatial object with X and Y coordinates
spatial_data <- SpatialPoints(data_sp[, c("home.X", "home.Y")])
# Create a distance-based spatial weights matrix (use knn, 3 nearest neighbors)
W <- knn2nb(knearneigh(spatial_data, k = 3))
# Convert the neighbor list (W) to a listw object
W_listw <- nb2listw(W)
# Calculate Moran's I for the "color" variable
moran.test(data_sp$color_num, W_listw)
##
## Moran I test under randomisation
##
## data: data_sp$color_num
## weights: W_listw
##
## Moran I statistic standard deviate = 0.19218, p-value = 0.4238
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.005071310 -0.006711409 0.003758949
color as numeric scores (B = 0, R = 4, and everything in between) not different from random distribution
data_sp <- data_M %>% filter(territorial == 1) %>%
select(call.X, call.Y, color_num) %>%
na.omit() %>% # remove those with NA coordinates
mutate(call.X = jitter(call.X), call.Y = jitter(call.Y)) #add some noise to avoid overlapping points
# Create a spatial object with X and Y coordinates
spatial_data <- SpatialPoints(data_sp[, c("call.X", "call.Y")])
# Create a distance-based spatial weights matrix (use knn, 3 nearest neighbors)
W <- knn2nb(knearneigh(spatial_data, k = 3))
# Convert the neighbor list (W) to a listw object
W_listw <- nb2listw(W)
# Calculate Moran's I for the "color" variable
moran.test(data_sp$color_num, W_listw)
##
## Moran I test under randomisation
##
## data: data_sp$color_num
## weights: W_listw
##
## Moran I statistic standard deviate = -0.15084, p-value = 0.5599
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04275336 -0.02564103 0.01287056
not different from random distribution
data_sp <- data_M %>% filter(territorial == 1) %>%
select(home.X, home.Y, color_num) %>%
na.omit() %>% # remove those with NA coordinates
mutate(home.X = jitter(home.X), home.Y = jitter(home.Y)) #add some noise to avoid overlapping points
# Create a spatial object with X and Y coordinates
spatial_data <- SpatialPoints(data_sp[, c("home.X", "home.Y")])
# Create a distance-based spatial weights matrix (use knn, 3 nearest neighbors)
W <- knn2nb(knearneigh(spatial_data, k = 3))
# Convert the neighbor list (W) to a listw object
W_listw <- nb2listw(W)
# Calculate Moran's I for the "color" variable
moran.test(data_sp$color_num, W_listw)
##
## Moran I test under randomisation
##
## data: data_sp$color_num
## weights: W_listw
##
## Moran I statistic standard deviate = -0.64614, p-value = 0.7409
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.10015482 -0.02564103 0.01329916