【问题标题】：Efficient and accurate age calculation (in years, months, or weeks) in R given birth date and an arbitrary date在给定出生日期和任意日期的 R 中，有效且准确的年龄计算（以年、月或周为单位）
【发布时间】：2015-09-16 13:38:30
【问题描述】：

在给定出生日期和任意日期的情况下，我面临着计算年龄（以年、月或周为单位）的常见任务。问题是我经常需要对许多记录（>3 亿）执行此操作，因此性能是这里的关键问题。

在 SO 和 Google 中快速搜索后，我找到了 3 个替代方案：

一个常用的算术过程 (/365.25) (link)
使用包lubridate（link）中的函数new_interval()和duration()
函数age_calc()来自包eeptools（link，link，link）

所以，这是我的玩具代码：

# 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 倍）。不幸的是，算术方法不够准确，我无法承受这种方法会犯的一些错误。

“因为现代公历的方式构造，没有简单的算术产生一个人的年龄的方法，根据常用用法——常用用法意味着一个人的年龄应该始终是一个整数，恰好增加生日”。(link)

正如我在一些帖子中看到的，lubridate 和 eeptools 没有犯这样的错误（不过，我还没有查看代码/阅读更多关于这些函数的信息以了解它们使用哪种方法），这就是我想要的原因使用它们，但它们的性能不适用于我的实际应用程序。

关于计算年龄的有效且准确的方法有什么想法吗？

编辑

Ops，看来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

【问题讨论】：

如果有比lubridate 更快/更简单的方法，我会感到惊讶。如果您真的需要性能提升，我唯一的建议是先执行算术方法，然后使用 lubridate 方法重做所有“关闭调用”（例如，如果 abs(floor(age) - age) < 0.01) 然后使用 lubridate）
谢谢。你是lubridate 的经验丰富的用户吗？，...，正如我在编辑的问题中所说的那样，我发现它会出错（可能不仅仅是算术方法），但我在几篇文章中读到了 lubridate 是确实是能够准确计算年龄的 R 软件包之一。所以现在我想知道我是否做错了什么。（我认为不是，我基本上是按照示例进行的，这很简单，但只是仔细检查）
不是解决方案，但 difftime(givendate, birthdate) / 365.25 似乎比 (givendate - birthdate) / 365.25) 快 5% 左右。如果你最终得到算术，可能会很有用。
@Molx 好电话！ -.Date 只是更稳健地调用 difftime。 (unclass(givendate) - unclass(birthdate)) / 365.25 应该更快，因为这进一步跳过了 difftime 的开销。

标签： r lubridate

【解决方案1】：

好的，所以我在另一个post找到了这个功能：

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)

至少在我的示例中它没有犯任何错误（在任何示例中都不应该；这是一个使用ifelses 的非常简单的函数）。

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

遗憾的是，Note 似乎比上述方法还要慢！

> 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

【讨论】：

new_interval() 现在在 lubridate 包中已弃用，因此请改用 interval()。此外，“对因素没有意义”警告消息（由 as.period() 导致已解决，因此不再出现。

【解决方案3】：

我打算把它留在 cmets 中，但我认为它值得单独回答。正如@Molx 指出的那样，您的“算术”方法并不像看起来那么简单——看看-.Date 的代码，最重要的是：

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

因此，Date 类对象的“算术”方法实际上是 difftime 函数的包装器。 difftime 呢？如果您追求的是原始速度，这也会产生大量开销。

关键是 Date 对象存储为自/直到 1970 年 1 月 1 日的整数天数（尽管它们实际上并未存储为 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)，更不用说ifelses)

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】：

我一直在努力解决这个问题，最终得到了 a) 完美准确*（与迄今为止提出的所有其他选项相比）和 b ）相当快（请参阅我在另一个答案中的基准）。它依赖于我手工完成的大量算术运算以及来自 data.table 包的精彩 foverlaps 函数。

该方法的本质是从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 功能方面也取得了很大进展，但是a）它相当丑陋（需要0-1460 的非线性变换以避免嵌套@987654329 @ 语句等）和 b）最后一个 for 循环（在整个日期列表中以 apply 的形式）是不可避免的，所以我认为这会让事情变得太慢。（转换为x1 = (unclass(birthdays) - 59) %% 1461; x2 = x1 * (729 - x1) / 402232 + x1，为后代）

我已将此功能添加到my package。

*（对于不关心non-leap centuries 的日期范围；不过，我相信处理此类日期的扩展应该不会太繁琐）

【讨论】：

对我来说，似乎是唯一可行的解决方案。干得好！