Эффективный и точный расчет возраста (в годах, месяцах или неделях) в R с учетом даты рождения и произвольной даты
Передо мной стоит общая задача расчета возраста (в годах, месяцах или неделях) с учетом даты рождения и произвольной даты. Дело в том, что довольно часто мне приходится делать это на множестве записей (>300 миллионов), поэтому производительность является ключевым вопросом здесь.
После быстрого поиска в SO и Google я нашел 3 варианта:
- Общая арифметическая процедура (/365.25) ( ссылка)
- Использование функций
new_interval()а такжеduration()из пакетаlubridate( ссылка) - функция
age_calc()из пакетаeeptools( ссылка, ссылка, ссылка)
Итак, вот мой игрушечный код:
# 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's 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)
Вот результаты:
Итог: производительность lubridate а также eeptools функции намного хуже, чем арифметический метод (/365.25 как минимум в 10 раз быстрее). К сожалению, арифметический метод недостаточно точен, и я не могу позволить себе несколько ошибок, которые этот метод сделает.
"из-за того, как построен современный григорианский календарь, не существует простого арифметического метода, который бы определял возраст человека, определяемый в соответствии с обычным использованием - общее использование означает, что возраст человека всегда должен быть целым числом, которое увеличивается ровно в день рождения". ( ссылка)
Как я читаю на некоторых постах, lubridate а также eeptools не делайте таких ошибок (хотя я не смотрел код / читал больше об этих функциях, чтобы узнать, какой метод они используют), и именно поэтому я хотел их использовать, но их производительность не работает для моего реального приложения.
Есть идеи по поводу эффективного и точного метода расчета возраста?
РЕДАКТИРОВАТЬ
Опс, похоже 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
4 ответа
Итак, я нашел эту функцию в другом посте:
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)
И по крайней мере в моих примерах это не делает никакой ошибки (и это не должно происходить ни в одном примере; это довольно простая функция, использующая ifelseс).
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 опубликует это как ответ.
- Я буду ждать, чтобы увидеть, придет ли кто-нибудь еще с полным решением (эффективный, точный и производящий возраст в годах, месяцах или неделях по желанию).
Причина, по которой lubridate, по-видимому, допускает ошибки, состоит в том, что вы вычисляете продолжительность (точное время, которое происходит между двумя моментами, где 1 год = 31536000 с), а не периоды (изменение времени часов, которое происходит между двумя моментами).
Для изменения времени (в годах, месяцах, днях и т. Д.) Необходимо использовать
as.period(new_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(new_interval(start = birthdate, end = givendate))$year
[1] 37 37 36 53 53 52 50 50 49 1 1 0 46 47 47
Обратите внимание, что это вызовет следующее предупреждающее сообщение (не знаю почему):
Warning message:
In Ops.factor(left, right) : ‘-’ not meaningful for factors
и, к сожалению, появляется даже медленнее, чем методы выше!
> 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
Я собирался оставить это в комментариях, но я думаю, что это заслуживает отдельного ответа. Как указывает @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's 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)иш), не говоря уже о том, чтобы ifelse с)
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
Таким образом, мы также подчеркиваем глупость сравнительного анализа небольших данных.
Большая стоимость метода @ Джима заключается в том, что as.POSIXlt становится все дороже, поскольку ваши векторы растут.
Проблема неточности остается, но если эта точность не имеет первостепенного значения, кажется, unclass метод не имеет аналогов.
Я старался справиться с этим и, наконец, получил что-то, что а) совершенно точно * (в отличие от всех других представленных вариантов) и б) достаточно быстро (см. Мои тесты в другом ответе). Он опирается на кучу арифметики, которую я сделал вручную, и замечательный foverlaps функция от data.table пакет.
Суть подхода заключается в том, чтобы работать от целочисленного представления Dates, а также признать, что все даты рождения попадают в один из четырех циклов 1461 (= 365 * 4 + 1) дней, в зависимости от того, когда наступит следующий год, когда наступит ваш день рождения в 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для потомков)
Я добавил эту функцию в свой пакет.
* (для диапазонов дат, когда не високосные столетия не являются проблемой; однако я считаю, что расширение для обработки таких дат не должно быть слишком обременительным)