¿Visualiza la dirección más probable de migración dados puntos con edades usando R?

Tengo una pregunta sobre la mejor manera de visualizar la dirección prevista de la migración dado un conjunto de puntos con edades estimadas.

Imagine que tengo un conjunto de coordenadas geográficas que representan comunidades. Para cada punto, tengo una estimación de la edad de esa comunidad. Quiero dibujar una flecha en el mapa que indique la dirección más probable de la migración general a lo largo del tiempo. También quiero que la magnitud de la flecha represente cuán consistente es el gradiente / qué tan seguros estamos de que hay una alineación significativa entre el tiempo y el espacio. No estoy preguntando sobre la teoría del proceso de migración, sino sobre cómo decidirías sobre el ángulo, la posición y la magnitud de la flecha.

Por ejemplo, debajo de los puntos rojos son más antiguos que los puntos verdes. En el primer recuadro, parece una conclusión sensata que la migración fue hacia el sudeste. En el segundo cuadro, no hay un patrón claro, por lo que la flecha es más pequeña.

¿Existe una forma estándar y de principios de hacer esto? Pensé en una forma de dibujar un vector entre el punto más antiguo y el segundo más antiguo, luego entre el segundo más antiguo y el tercero más antiguo, y así sucesivamente. Luego suma los vectores (¿pero cómo decidir la ubicación de inicio?). ¿O tal vez necesita tomar los vectores entre todos los pares de puntos y ponderarlos según las edades relativas? ¿O tal vez necesito resolver los contornos, luego dibujar un camino que comience en el punto más alto y vaya cuesta abajo?

Una solución basada en R sería ideal.

Aquí está mi solución actual, basada en la sugerencia anterior para ajustar una superficie. Esta solución no aborda la escala correctamente, por lo que las coordenadas y los valores z deben ser pequeños (por ejemplo, entre -1 y 1).

drawDirection = function(x,y,z){
  # use linear model to fit surface
  m = lm(z ~ x+y)

  # Treat coefficeints like vectors and 
  # define arrow start and end points
  arrow = c(0,0,
            m$coefficients["x"],
            m$coefficients['y'])
  # Move arrow to centroid
  adj.x = mean(x)- (m$coefficients["x"]/2)
  adj.y = mean(y)- (m$coefficients["y"]/2)
  arrow = arrow + c(adj.x,adj.y,adj.x,adj.y)

  # colours for points
  # Yellow = higher = more recent
  z.col = heat.colors(10)[as.numeric(cut(z,breaks=10))]
  # Plot points with some extra space around
  plot(x,y, col=z.col, pch=16,
       xlim=c(min(x)-sd(x)*2,max(x)+sd(x)*2),
       ylim=c(min(y)-sd(y)*2,max(y)+sd(y)*2))
  # Plot arrow
  # (scale width by 10 x the R squared)
  arrows(arrow[1],arrow[2],arrow[3],arrow[4],
         lwd = 1+(10*summary(m)$adj.r.squared))
}

par(mfrow=c(1,2))

# Simulate some data
n = 30
# Correlated North-East
x = runif(n,0,1)
y = jitter(x,amount=1)
z = jitter(x+y,amount=1)
x = x
y = y
drawDirection(x,y,z)

# Uncorrelated
x2 = runif(n,0,1)
y2 = runif(n,0,1)
z2 = runif(n,0,1)

drawDirection(x2,y2,z2)

Sin embargo, creo que estoy haciendo el escalado incorrectamente: si aumento el rango de x e y, la flecha no se escala correctamente:

# Correlated North-East, 
#  with x and y variables between 0 and 360
x = runif(n,0,1)
y = jitter(x,amount=1)
z = jitter(x+y,amount=1)
x = x*360
y = y*360
drawDirection(x,y,z)

Entonces, por ahora, podría escalar las coordenadas geográficas y el valor z.