Блог

Итак, хотя lag и lead в dplyr великолепны, я хочу смоделировать временные ряды чего-то вроде роста населения. Мой старый школьный код выглядел бы примерно так:

 tdf <- data.frame(time=1:5, pop=50)
for(i in 2:5){
  tdf$pop[i] = 1.1*tdf$pop[i-1]
}
  

что приводит к

   time    pop
1    1 50.000
2    2 55.000
3    3 60.500
4    4 66.550
5    5 73.205
  

Я чувствую dplyr tidyverse , что для этого должен быть способ or (как бы мне ни нравился мой цикл for).

Но что-то вроде

 tdf <- data.frame(time=1:5, pop=50) %>%
  mutate(pop = 1.1*lag(pop))
  

что было бы моим первым предположением, просто приводит

   time pop
1    1  NA
2    2  55
3    3  55
4    4  55
5    5  55
  

Я чувствую, что упускаю что-то очевидное…. что это?

Примечание — это тривиальный пример — в моих реальных примерах используется несколько параметров, многие из которых изменяются во времени (я моделирую прогнозы в разных сценариях GCM), поэтому tidyverse оказывается мощным инструментом для объединения моих симуляций.

Reduce (или его варианты purrr, если хотите) — это то, что вам нужно для кумулятивных функций, для которых еще не написана cum* версия:

 data.frame(time = 1:5, pop = 50) %>%
    mutate(pop = Reduce(function(x, y){x * 1.1}, pop, accumulate = TRUE))

##   time    pop
## 1    1 50.000
## 2    2 55.000
## 3    3 60.500
## 4    4 66.550
## 5    5 73.205
  

или с помощью purrr,

 data.frame(time = 1:5, pop = 50) %>%
    mutate(pop = accumulate(pop, ~.x * 1.1))

##   time    pop
## 1    1 50.000
## 2    2 55.000
## 3    3 60.500
## 4    4 66.550
## 5    5 73.205
  

Если начальное значение pop равно, скажем, 50, тогда pop = 50 * 1.1^(0:4) вы получите следующие четыре значения. С вашим кодом вы могли бы сделать:

 data.frame(time=1:5, pop=50) %>%
  mutate(pop = pop * 1.1^(1:n() - 1))
  

Или,

 base = 50

data.frame(time=1:5) %>%
  mutate(pop = base * 1.1^(1:n()-1))
  

Функция накопления Purrr может обрабатывать изменяющиеся во времени индексы, если вы передаете их в свою функцию моделирования в виде списка со всеми параметрами в нем. Тем не менее, требуется немного споров, чтобы заставить это работать правильно. Хитрость здесь в том, что accumulate() может работать как со списками, так и с векторными столбцами. Вы можете использовать tidyr функцию nest() для группировки столбцов в вектор списка, содержащий текущее состояние и параметры заполнения, а затем использовать accumulate() для результирующего столбца списка. Это немного сложно объяснить, поэтому я включил демонстрацию, имитирующую логистический рост либо с постоянной скоростью роста, либо с изменяющейся во времени скоростью стохастического роста. Я также включил пример того, как использовать это для имитации нескольких повторов для данной модели с использованием dpylr purrr tidyr.

 library(dplyr)
library(purrr)
library(ggplot2)
library(tidyr)

# Declare the population growth function. Note: the first two arguments
# have to be .x (the prior vector of populations and parameters) and .y,
# the current parameter value and population vector. 
# This example function is a Ricker population growth model. 
logistic_growth = function(.x, .y, growth, comp) {
  pop = .x$pop[1]
  growth = .y$growth[1]
  comp  = .y$comp[1]
  # Note: this uses the state from .x, and the parameter values from .y.
  # The first observation will use the first entry in the vector for .x and .y
  new_pop = pop*exp(growth - pop*comp)
  .y$pop[1] = new_pop
  return(.y)
}

# Starting parameters the number of time steps to simulate, initial population size,
# and ecological parameters (growth rate and intraspecific competition rate)
n_steps  = 100
pop_init = 1
growth = 0.5
comp = 0.05

#First test: fixed growth rates
test1 = data_frame(time = 1:n_steps,pop = pop_init, 
                   growth=growth,comp =comp)


# here, the combination of nest() and group_by() split the data into individual 
# time points and then groups all parameters into a new vector called state.
# ungroup() removes the grouping structure, then accumulate runs the function
#on the vector of states. Finally unnest transforms it all back to a
#data frame
out1 = test1 %>%
  group_by(time)%>%
  nest(pop, growth, comp,.key = state)%>%
  ungroup()%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

# This is the same example, except I drew the growth rates from a normal distribution
# with a mean equal to the mean growth rate and a std. dev. of 0.1
test2 = data_frame(time = 1:n_steps,pop = pop_init, 
                   growth=rnorm(n_steps, growth,0.1),comp=comp)

out2 = test2 %>%
  group_by(time)%>%
  nest(pop, growth, comp,.key = state)%>%
  ungroup()%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

# This demostrates how to use this approach to simulate replicates using dplyr
# Note the crossing function creates all combinations of its input values
test3 = crossing(rep = 1:10, time = 1:n_steps,pop = pop_init, comp=comp) %>%
  mutate(growth=rnorm(n_steps*10, growth,0.1))

out3 = test3 %>%
  group_by(rep)%>%
  group_by(rep,time)%>%
  nest(pop, growth, comp,.key = state)%>%
  group_by(rep)%>%
  mutate(
    state = accumulate(state,logistic_growth))%>%
  unnest()

print(qplot(time, pop, data=out1) 
        geom_line()  
        geom_point(data= out2, col="red") 
        geom_line(data=out2, col="red") 
        geom_point(data=out3, col="red", alpha=0.1) 
        geom_line(data=out3, col="red", alpha=0.1,aes(group=rep)))
  

Проблема здесь в том, что dplyr выполняется это как набор векторных операций, а не вычисление термина по одному за раз. Здесь 1.1*lag(pop) интерпретируется как «вычислите запаздывающие значения для всех pop, затем умножьте их все на 1.1». Поскольку вы set pop=50 запаздывали, значения для всех шагов были равны 50.

dplyr имеет некоторые вспомогательные функции для последовательной оценки; стандартная функция cumsum cumprod , и т.д. работают, и несколько новых (см. ?cummean ) все работают внутри dplyr . В вашем примере вы могли бы смоделировать модель с помощью:

 tdf <- data.frame(time=1:5, pop=50, growth_rate = c(1, rep(1.1,times=4)) %>%
    mutate(pop = pop*cumprod(growth_rate))


time    pop     growth_rate
1       50.000  1.0
2       55.000  1.1
3       60.500  1.1
4       66.550  1.1
5       73.205  1.1
  

Обратите внимание, что я добавил здесь скорость роста в качестве столбца, и я установил первую скорость роста равной 1. Вы также можете указать это следующим образом:

 tdf <- data.frame(time=1:5, pop=50, growth_rate = 1.1) %>%
    mutate(pop = pop*cumprod(lead(growth_rate,default=1))
  

Это делает явным, что столбец скорости роста относится к скорости роста на текущем временном шаге по сравнению с предыдущим.

Существуют ограничения на количество различных симуляций, которые вы можете выполнить таким образом, но должно быть возможно построить множество экологических моделей с дискретным временем, используя некоторую комбинацию кумулятивных функций и параметров, указанных в столбцах.

Как насчет функций отображения, т.Е.

 tdf <- data_frame(time=1:5)
tdf %>% mutate(pop = map_dbl(.x = tdf$time, .f = (function(x) 50*1.1^x)))