spatialRF_talk/functions.R at main · BlasBenito/spatialRF_talk · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
moran_custom <- function(
  x = NULL,
  distance.matrix = NULL,
  distance.threshold = NULL,
  verbose = TRUE
){

  #extracting weights from distance matrix
  x.distance.weights <- weights_from_distance_matrix(
    distance.matrix = distance.matrix,
    distance.threshold = distance.threshold
  )

  #length of the input vector
  x.length <- length(x)

  #computing expected Moran I
  moran.i.expected <- -1/(x.length - 1)

  #sum of weights
  x.distance.weights.sum <- sum(x.distance.weights)

  #centering x
  x.mean <- mean(x)
  x.centered <- x - x.mean

  #upper term of the Moran's I equation
  #sum of cross-products of the lags
  #x.centered %o% x.centered equals (xi - x.mean) * (xj - x.mean)
  cross.product.sum <- sum(x.distance.weights * x.centered %o% x.centered)

  #lower term of the Moran's I equation
  #variance of the centered x
  x.centered.variance <- sum(x.centered^2)

  #observed Moran's I
  moran.i.observed <- (x.length / x.distance.weights.sum) * (cross.product.sum/x.centered.variance)

  #preparing moran plot
  #computing weighted mean of the lag
  x.lag <- apply(
    X = x.distance.weights,
    MARGIN = 1,
    FUN = function(y){y %*% x}
  )

  #preparing title
  plot.title <- paste0(
    "Moran's I = ",
    round(moran.i.observed, 3)
  )

  #preparing data frame
  plot.df <- data.frame(
    x = spatialRF::rescale_vector(x),
    x.lag = spatialRF::rescale_vector(x.lag)
  )

  #moran plot
  #plot
  moran.plot <- ggplot2::ggplot(data = plot.df) +
    ggplot2::aes(
      x = x,
      y = x.lag
    ) +
    ggplot2::geom_point(alpha = 0.75, color = "#f1605d") +
    ggplot2::geom_smooth(
      method = "lm",
      color = "#0d0887",
      fill = "#0d0887",
      se = TRUE,
      alpha = 0.2,
      formula = y ~ x,
      size = 0.6,
      level = 0.99
    ) +
    ggplot2::xlab("Value") +
    ggplot2::ylab("Weighted value of neighbors") +
    ggplot2::theme(
      legend.position = "none",
      plot.margin = ggplot2::margin(t = 0, r = 0, b = 0, l = 0),
      plot.title = element_text(hjust = 0.5)
      ) +
    ggplot2::ggtitle(plot.title) +
    ggplot2::theme_bw(base_size = 10) +
    ggplot2::coord_fixed(xlim = c(0, 1), ylim = c(0, 1)) +
    theme_bw(base_size = 25)

  moran.plot

}