Подтвердить что ты не робот

Эффективный и точный расчет возраста (в годах, месяцах или неделях) в R, дата рождения и произвольная дата

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

После быстрого поиска в SO и Google я нашел 3 варианта:

Итак, вот мой игрушечный код:

# Some toy birthdates
birthdate <- as.Date(c("1978-12-30", "1978-12-31", "1979-01-01", 
                       "1962-12-30", "1962-12-31", "1963-01-01", 
                       "2000-06-16", "2000-06-17", "2000-06-18", 
                       "2007-03-18", "2007-03-19", "2007-03-20", 
                       "1968-02-29", "1968-02-29", "1968-02-29"))

# Given dates to calculate the age
givendate <- as.Date(c("2015-12-31", "2015-12-31", "2015-12-31", 
                       "2015-12-31", "2015-12-31", "2015-12-31", 
                       "2050-06-17", "2050-06-17", "2050-06-17",
                       "2008-03-19", "2008-03-19", "2008-03-19", 
                       "2015-02-28", "2015-03-01", "2015-03-02"))

# Using a common arithmetic procedure ("Time differences in days"/365.25)
(givendate-birthdate)/365.25

# Use the package lubridate
require(lubridate)
new_interval(start = birthdate, end = givendate) / 
                     duration(num = 1, units = "years")

# Use the package eeptools
library(eeptools)
age_calc(dob = birthdate, enddate = givendate, units = "years")

Давайте поговорим позже о точности и сосредоточимся в первую очередь на производительности. Вот код:

# Now let compare the performance of the alternatives using microbenchmark
library(microbenchmark)
mbm <- microbenchmark(
    arithmetic = (givendate - birthdate) / 365.25,
    lubridate = new_interval(start = birthdate, end = givendate) /
                                     duration(num = 1, units = "years"),
    eeptools = age_calc(dob = birthdate, enddate = givendate, 
                        units = "years"),
    times = 1000
)

# And examine the results
mbm
autoplot(mbm)

Вот результаты:

Microbenchmark results - table Microbenchmark results - plot

Итог: производительность функций lubridate и eeptools намного хуже, чем арифметический метод (/365.25 как минимум в 10 раз быстрее). К сожалению, арифметический метод недостаточно точен, и я не могу позволить себе несколько ошибок, которые этот метод допустит.

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

Как я читал в некоторых сообщениях, lubridate и eeptools не делают таких ошибок (хотя я не смотрел код/​​читал больше об этих функциях, чтобы узнать, какой метод они используют) и поэтому я хотел их использовать., но их производительность не работает для моего реального приложения.

Есть идеи по поводу эффективного и точного метода расчета возраста?

EDIT

Опс, кажется, lubridate также делает ошибки. И, очевидно, основываясь на этом игрушечном примере, он делает больше ошибок, чем арифметический метод (см. строки 3, 6, 9, 12). (я что-то не так делаю?)

toy_df <- data.frame(
    birthdate = birthdate,
    givendate = givendate,
    arithmetic = as.numeric((givendate - birthdate) / 365.25),
    lubridate = new_interval(start = birthdate, end = givendate) /
        duration(num = 1, units = "years"),
    eeptools = age_calc(dob = birthdate, enddate = givendate,
                        units = "years")
)
toy_df[, 3:5] <- floor(toy_df[, 3:5])
toy_df

    birthdate  givendate arithmetic lubridate eeptools
1  1978-12-30 2015-12-31         37        37       37
2  1978-12-31 2015-12-31         36        37       37
3  1979-01-01 2015-12-31         36        37       36
4  1962-12-30 2015-12-31         53        53       53
5  1962-12-31 2015-12-31         52        53       53
6  1963-01-01 2015-12-31         52        53       52
7  2000-06-16 2050-06-17         50        50       50
8  2000-06-17 2050-06-17         49        50       50
9  2000-06-18 2050-06-17         49        50       49
10 2007-03-18 2008-03-19          1         1        1
11 2007-03-19 2008-03-19          1         1        1
12 2007-03-20 2008-03-19          0         1        0
13 1968-02-29 2015-02-28         46        47       46
14 1968-02-29 2015-03-01         47        47       47
15 1968-02-29 2015-03-02         47        47       47
4b9b3361

Ответ 1

Хорошо, поэтому я нашел эту функцию в другом сообщении :

age <- function(from, to) {
    from_lt = as.POSIXlt(from)
    to_lt = as.POSIXlt(to)

    age = to_lt$year - from_lt$year

    ifelse(to_lt$mon < from_lt$mon |
               (to_lt$mon == from_lt$mon & to_lt$mday < from_lt$mday),
           age - 1, age)
}

Было опубликовано @Jim: "Следующая функция принимает векторы объектов Date и вычисляет возраст, правильно учитывая високосные годы. Кажется, это более простое решение, чем любой другой ответ".

Это действительно проще, и это трюк, который я искал. В среднем, это на самом деле быстрее, чем арифметический метод (примерно на 75% быстрее).

mbm <- microbenchmark(
    arithmetic = (givendate - birthdate) / 365.25,
    lubridate = interval(start = birthdate, end = givendate) /
        duration(num = 1, units = "years"),
    eeptools = age_calc(dob = birthdate, enddate = givendate, 
                        units = "years"),
    age = age(from = birthdate, to = givendate),
    times = 1000
)
mbm
autoplot(mbm)

enter image description hereenter image description here

И по крайней мере в моих примерах это не делает никакой ошибки (и это не должно быть ни в одном примере, это довольно простая функция, использующая ifelse s).

toy_df <- data.frame(
    birthdate = birthdate,
    givendate = givendate,
    arithmetic = as.numeric((givendate - birthdate) / 365.25),
    lubridate = interval(start = birthdate, end = givendate) /
        duration(num = 1, units = "years"),
    eeptools = age_calc(dob = birthdate, enddate = givendate,
                        units = "years"),
    age = age(from = birthdate, to = givendate)
)
toy_df[, 3:6] <- floor(toy_df[, 3:6])
toy_df

    birthdate  givendate arithmetic lubridate eeptools age
1  1978-12-30 2015-12-31         37        37       37  37
2  1978-12-31 2015-12-31         36        37       37  37
3  1979-01-01 2015-12-31         36        37       36  36
4  1962-12-30 2015-12-31         53        53       53  53
5  1962-12-31 2015-12-31         52        53       53  53
6  1963-01-01 2015-12-31         52        53       52  52
7  2000-06-16 2050-06-17         50        50       50  50
8  2000-06-17 2050-06-17         49        50       50  50
9  2000-06-18 2050-06-17         49        50       49  49
10 2007-03-18 2008-03-19          1         1        1   1
11 2007-03-19 2008-03-19          1         1        1   1
12 2007-03-20 2008-03-19          0         1        0   0
13 1968-02-29 2015-02-28         46        47       46  46
14 1968-02-29 2015-03-01         47        47       47  47
15 1968-02-29 2015-03-02         47        47       47  47

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

  • Я бы подождал, когда @Jim опубликует его как ответ.
  • Я буду ждать, если кто-нибудь придумает полное решение (эффективное, точное и продюсирование по годам, месяцам или неделям по желанию).

Ответ 2

Причина, по которой lubridate, по-видимому, допускает ошибки, заключается в том, что вы вычисляете продолжительность (точное время, которое происходит между двумя моментами, где 1 год = 31536000 с), а не периоды (изменение времени часов, которое происходит между двумя моментами).

Чтобы изменить время (в годах, месяцах, днях и т.д.), Вам нужно использовать

as.period(interval(start = birthdate, end = givendate))

который дает следующий вывод

 "37y 0m 1d 0H 0M 0S"   
 "37y 0m 0d 0H 0M 0S"   
 "36y 11m 30d 0H 0M 0S" 
 ...
 "46y 11m 30d 1H 0M 0S" 
 "47y 0m 0d 1H 0M 0S"   
 "47y 0m 1d 1H 0M 0S" 

Чтобы просто извлечь годы, вы можете использовать следующие

as.period(interval(start = birthdate, end = givendate))$year
 [1] 37 37 36 53 53 52 50 50 49  1  1  0 46 47 47

Примечание, к сожалению, появляется даже медленнее, чем описанные выше методы!

> mbm
Unit: microseconds
       expr       min        lq       mean    median         uq        max neval cld
 arithmetic   116.595   138.149   181.7547   184.335   196.8565   5556.306  1000  a 
  lubridate 16807.683 17406.255 20388.1410 18053.274 21378.8875 157965.935  1000   b

Ответ 3

Я собирался оставить это в комментариях, но я думаю, что это заслуживает отдельного ответа. Как указывает @Molx, ваш "арифметический" метод не так прост, как кажется - взгляните на код для -.Date, самое главное:

return(difftime(e1, e2, units = "days"))

Таким образом, "арифметический" метод для объектов класса Date действительно является оберткой для функции difftime. Как насчет difftime? Это также имеет много накладных расходов, если вы ищете грубую скорость.

Ключевым моментом является то, что объекты Date хранятся как целое число дней с/до 1 января 1970 года (хотя на самом деле они не сохраняются как integer, отсюда и возникновение класса IDate в data.table), поэтому мы можем просто вычесть их и покончить с этим, но чтобы избежать вызова метода -.Date, нам нужно unclass наши входные данные:

(unclass(birthdate) - unclass(givendate)) / 365.25

Что касается эффекта, то этот подход на несколько порядков быстрее, чем даже метод @Jim age.

Вот еще несколько расширенных тестовых данных:

set.seed(20349)
NN <- 1e6
birthdate <- as.Date(sprintf('%d-%02d-%02d',
                             sample(1901:2030, NN, TRUE),
                             sample(12, NN, TRUE),
                             sample(28, NN, TRUE)))

#average 30 years, most data between 20 and 40 years
givendate <- birthdate + as.integer(rnorm(NN, mean = 10950, sd = 1000))

(исключая eeptools, потому что это почти невозможно медленнее - взгляд на код для age_calc показывает, что код заходит так далеко, что создает последовательность дат для каждой пары дат (O(n^2) -ish), не говоря уже о ifelse s)

microbenchmark(
  arithmetic = (givendate - birthdate) / 365.25,
  lubridate = interval(start = birthdate, end = givendate) /
    duration(num = 1, units = "years"),
  age = age(from = birthdate, to = givendate),
  fastar = (unclass(givendate) - unclass(birthdate)) / 365.25,
  overlaps = get_age(birthdate, givendate),
  times = 50)
# Unit: milliseconds
#        expr        min         lq      mean     median         uq      max neval  cld
#  arithmetic  28.153465  30.384639  62.96118  31.492764  34.052991 180.9556    50  b  
#   lubridate  94.327968  97.233009 157.30420 102.751351 240.717065 265.0283    50   c 
#         age 338.347756 479.598513 483.84529 483.580981 488.090832 770.1149    50    d
#      fastar   7.740098   7.831528  11.02521   7.913146   8.090902 153.3645    50 a   
#    overlaps 316.408920 458.734073 459.58974 463.806255 470.320072 769.0929    50    d

Таким образом, мы также подчеркиваем глупость сравнительного анализа небольших данных.

Большая стоимость метода @Jim заключается в том, что as.POSIXlt становится все дороже с ростом ваших векторов.

Проблема неточности остается, но если эта точность не имеет первостепенного значения, кажется, что метод unclass не имеет аналогов.

Ответ 4

Я старался справиться с этим и, наконец, получил что-то, что а) совершенно точно * (в отличие от всех других представленных вариантов) и б) достаточно быстро (см. мои тесты в другом ответе). Он опирается на кучу арифметических действий, которые я делал вручную, и на замечательную функцию foverlaps из пакета data.table.

Суть этого подхода заключается в том, чтобы работать с целочисленным представлением Date s, а также признавать, что все даты рождения попадают в один из четырех циклов 1461 (= 365 * 4 + 1) -day, в зависимости от того, когда в следующем году ваш день рождения займет 366 дней.

Здесь функция:

library(data.table)
get_age <- function(birthdays, ref_dates){
  x <- data.table(bday <- unclass(birthdays),
                  #rem: how many days has it been since the lapse of the
                  #  most recent quadrennium since your birth?
                  rem = ((ref <- unclass(ref_dates)) - bday) %% 1461)
  #cycle_type: which of the four years following your birthday
  #  was the one that had 366 days? 
  x[ , cycle_type := 
       foverlaps(data.table(start = bdr <- bday %% 1461L, end = bdr),
                 #these intervals were calculated by hand;
                 #  e.g., 59 is Feb. 28, 1970. I made the judgment
                 #  call to say that those born on Feb. 29 don't
                 #  have their "birthday" until the following March 1st.
                 data.table(start = c(0L, 59L, 424L, 790L, 1155L), 
                            end = c(58L, 423L, 789L, 1154L, 1460L), 
                            val = c(3L, 2L, 1L, 4L, 3L),
                            key = "start,end"))$val]
  I4 <- diag(4L)[ , -4L] #for conciseness below
  #The 'by' approach might seem a little abstruse for those
  #  not familiar with 'data.table'; see the edit history
  #  for a more palatable version (which is also slightly slower)
  x[ , extra := 
       foverlaps(data.table(start = rem, end = rem),
                 data.table(start = st <- cumsum(c(0L, rep(365L, 3L) +
                                                     I4[.BY[[1L]],])),
                            end = c(st[-1L] - 1L, 1461L),
                            int_yrs = 0:3, key = "start,end")
       )[ , int_yrs + (i.start - start) / (end + 1L - start)], by = cycle_type]
  #grand finale -- 4 years for every quadrennium, plus the fraction:
  4L * ((ref - bday) %/% 1461L) + x$extra
}

Сравнение с вашим основным примером:

toy_df <- data.frame(
  birthdate = birthdate,
  givendate = givendate,
  arithmetic = as.numeric((givendate - birthdate) / 365.25),
  lubridate = interval(start = birthdate, end = givendate) /
    duration(num = 1, units = "years"),
  eeptools = age_calc(dob = birthdate, enddate = givendate,
                      units = "years"),
  mine = get_age(birthdate, givendate)
)

toy_df
#     birthdate  givendate arithmetic lubridate   eeptools       mine
# 1  1978-12-30 2015-12-31 37.0020534 37.027397 37.0027397 37.0027322 #eeptools wrong: will be 366 days until 12/31/16, so fraction is 1/366
# 2  1978-12-31 2015-12-31 36.9993155 37.024658 37.0000000 37.0000000
# 3  1979-01-01 2015-12-31 36.9965777 37.021918 36.9972603 36.9972603
# 4  1962-12-30 2015-12-31 53.0020534 53.038356 53.0027397 53.0027322 #same problem
# 5  1962-12-31 2015-12-31 52.9993155 53.035616 53.0000000 53.0000000
# 6  1963-01-01 2015-12-31 52.9965777 53.032877 52.9972603 52.9972603
# 7  2000-06-16 2050-06-17 50.0013689 50.035616 50.0000000 50.0027397 #eeptools wrong: not exactly the birthday
# 8  2000-06-17 2050-06-17 49.9986311 50.032877 50.9972603 50.0000000 #eeptools wrong: _is_ exactly the birthday
# 9  2000-06-18 2050-06-17 49.9958932 50.030137 49.9945205 49.9972603 #eeptools wrong: fraction should be 364/365
# 10 2007-03-18 2008-03-19  1.0047912  1.005479  1.0027322  1.0027397 #eeptools wrong: 2/29 already passed, only 365 days until 3/19/2009
# 11 2007-03-19 2008-03-19  1.0020534  1.002740  1.0000000  1.0000000
# 12 2007-03-20 2008-03-19  0.9993155  1.000000  0.9966839  0.9972678 #eeptools wrong: we passed 2/29, so should be 365/366
# 13 1968-02-29 2015-02-28 46.9979466 47.030137 46.9977019 46.9972603 #my judgment: birthday occurs on 3/1 for 2/29 babies, so 364/365 the way there
# 14 1968-02-29 2015-03-01 47.0006845 47.032877 47.0000000 47.0000000
# 15 1968-02-29 2015-03-02 47.0034223 47.035616 47.0027397 47.0027322

Этот подход может быть расширен, чтобы довольно легко обрабатывать месяцы/недели. Месяцы будут немного скучными (нужно указать длину месяца в 4 года), поэтому я не стал беспокоиться; недели - это просто (недели не зависят от соображений високосного года, поэтому мы можем просто разделить их на 7).

Я также добился большого прогресса в выполнении этого с функциями base, но а) это было довольно уродливо (необходимо нелинейное преобразование 0-1460, чтобы избежать выполнения вложенных операторов ifelse и т.д.) И б) в конец цикла for (в форме apply по всему списку дат) был неизбежен, поэтому я решил, что это слишком сильно замедлит ход событий. (трансформация x1 = (unclass(birthdays) - 59) %% 1461; x2 = x1 * (729 - x1) / 402232 + x1 для потомков)

Я добавил эту функцию в свой пакет.

* (для диапазонов дат, когда не високосные столетия не являются проблемой; однако я считаю, что расширение для обработки таких дат не должно быть слишком обременительным)