基于ggplot的政经数据可视化-日期/时间对象

在可视化任务中，我们时常会遇到与日期/时间相关的信息。我们固然可以使用诸如"2019-07-09"或"2019年7月"这样的字符，来表示日期/时间信息，并用正则表达式来进行处理，但实际上，R为日期/时间对象提供了专门的对象类型和大量相关的函数。接下来，我们将依次介绍日期/时间表示法、时间序列对象、lubridate包中的若干函数，以及与日期/时间有关的缺失值填补。读者既可以现在学习这些内容，也可以等到读到后边的折线图的绘制方法时再学习这些内容。

一、日期

1. 生成日期

x=Sys.Date() # 查看当前日期

# [1] "2019-07-09"

class(x) # 查看对象类型

# [1] "Date" # 可见这里的日期对象的类型为Date（尽管看上去像一个character对象）

## 用as.Date将character对象转化为Date对象；用as.character将Date对象转回character对象

x=as.Date("1998-08-31")

y=as.character(x)

## 还可将Date对象转为数值，数值代表的是日期距离1970年1月1日的天数， 1970年1月1日以后的为正数，以前的为负数

x=as.Date(c("1970-01-03", "1969-12-30", "1998-08-31"))

as.numeric(x)

# [1] 2 -2 10469

## 注意：不能用as.Date转化不存在的日期

# 因此，as.Date("2019-2-29")或as.Date("2019-4-31")都将报错

## as.Date可识别的默认格式

# 年份用四位数表示；月和日可以写满两位也可以只写一位；间隔符可以是"-"也可以是"/"；但必须保证向量中的元素只包含一种格式

as.Date(c("1998-08-31", "1998-8-31"))

as.Date(c("1998/8/31", "1998/08/31"))

as.Date(c("1998-08-31", "1998-8-31", "1998/8/31", "1998/08/31")) # 结果中出现NA，因为向量包含两种间隔符

## 用format参数帮助函数理解字符中的日期信息

as.Date("19980831", format="%Y%m%d") # 若不指定format，则会报错

# 这里的"%Y"等都是代表日期元素的表达式，其中，"%Y"代表四位数年份，"%y"代表两位数年份，"%m"代表月份，"%d"代表日

as.Date("1998abc08xyz31", format="%Yabc%mxyz%d") # 只要指定format，即使输入字符包含无效内容也仍可被识别

as.Date("8/31 1998", format="%m/%d %Y")

as.Date("19/07 09", format="%y/%m %d") # 这里的"19"只会被理解成2019年，不会被理解成1919年

as.Date("07-09", "%m-%d") # 当不包含年份时，返回当前年份

# as.Date("2018-12", "%Y-%m") # 注意：必须加上日子，否则结果为NA

与日期相关并以"%"开头的表达式，在生成日期/时间对象或从日期/时间对象中提取信息时，常会被用到。它们的具体含义可通过?strptime查到。以下，我们对常用的表达式进行了总结，其中一些要到后文讲到时间对象时才会出现，到时候读者若遇到不明白的表达式可回到此处查看。

●%Y：四位数年份。

●%y：两位数年份。数字00-68将被识别为20xx，如：2000、2050、2068；数字69-99将被识别为19xx，如：1968、1990、1999。

●%m：月份，如：01、02、1、2、11、12。

●%B：月份全称，如：八月、August。

●%b：月份缩写，如：8月、AUG。

●%d：日，如：01、02、1、2、30。

●%A：星期全称，如：星期一、Monday。

●%a：星期缩写，如：周一、MON。

●%u：以数字形式出现的星期，使用1-7的数字，如：周一表示为1，周二表示为2。

●%w：以数字形式出现的星期，如用0-6的数字，如：周日表示为0，周一表示为1。

●%H：24小时制两位数小时，如：00、01、23，但形如"24:00:00"的表达方式也被接受。

●%I：12小时制两位数小时，如：01、02、12。

●%p：在12小时制下指示上午或下午，如：上午、下午、AM、PM。有的区域设置或操作系统无法使用此表达式。

●%M：两位数分钟，即：00、01、59。

●%S：整数秒，取值为00-61的数字。

●%OS：可带小数的秒，"OS"后边的数字用于指定小数位数，如：对于13.348秒，"OS1"可提取到13.3，"OS2"可提取到13.34，"OS"后不加数字则只提取整数。

日期/时间对象的处理与区域设置有关。如果读者的操作系统的区域设置为简体中文的话，那么as.Date(x="1998/八月/31", format="%Y/%B/%d")将会得到正确的结果，而as.Date(x="1998/AUG/31", format="%Y/%B/%d")则无效。相反，如果读者的操作系统使用英文的话，后者有效，前者无效。注意：我们将使用"%B"和"%b"代表月份的完整名称和缩写。

## 那么，如何在不改变Windows系统区域设置的情况下在R内部作出更改

Sys.getlocale("LC_TIME") # 查询当前时间设置方式。在简体中文操作系统中,会得到"Chinese (Simplified)_China.936"或类似的结果

Sys.setlocale("LC_TIME", "English") # 我们现在将时间设置方式改为英语,以便使用英文名称

as.Date(c("1998/Aug/31", "1998/aug/31", "1998/AUG/31"), "%Y/%b/%d")

# 可忽略名称的大小写

as.Date(c("August3198", "august3198", "AUGUST3198"), " %B%d%y") # 可忽略名称的大小写

Sys.setlocale("LC_TIME", "Chinese (Simplified)_China.936") # 改回到原来的简体中文设置

## 当多种格式并存，并且我们能够罗列出所有可能的格式时，可通过以下方式指定多种格式

x=c("1998-08-31", "1998/08/31", "1998年8月31日", "19980831")

myformat=rep(NA, length(x))

for (i in 1: length(x)){

ii=x[i]

myformat[i]=if (grepl("\\-", ii)) "%Y-%m-%d"

else if (grepl("/", ii)) "%Y/%m/%d"

else if (grepl("年", ii)) "%Y年%m月%d日"

else "%Y%m%d"

}

as.Date(x, format=myformat)

## 通过指定天数和起始日期的时间确定日期

as.Date(2, origin="1998-08-31") # "1998-09-02"

as.Date(-2, origin="19980831", format="%Y%m%d") # "1998-08-29"

as.Date(1: 365, origin=as.Date("1998-08-31"))

## 通过seq生成等间隔时间

a=as.Date("1984-08-17"); b=as.Date("1998-08-31")

seq(a, b, length.out=10) # 生成等间隔的10个日期

seq(a, b, "2 month") # 以两个月为步长。注意："month"不需要使用复数

seq(a, b, "50 day") # "day"、"month"、"year"、"week"均可用来设置步长

seq(a, b, "2 week")

seq(a, by="2 year", length.out=8)

2. 比较和计算日期

## 日期比较

a=as.Date("1984-08-17"); b=as.Date("1998-08-31")

a < b # 靠后的时间较大

x=seq(a, b, by="1 year")

min(x); max(x); mean(x); median(x)

y=as.character(summary(x)) # 简单汇总

## 日期计算

c(a, b)+5 # 5天后

y=b-a # 相减

class(y) # 相减的结果是一个difftime对象

as.numeric(y) # 转化为普通数值

difftime(b, a) # 另一种相减的写法

difftime(b, a, units="hours") # 用units指定计算间隔时的单位，可选择"auto", "secs", "mins", "hours", "days", "weeks"

## 注意闰年带来的影响

as.Date("2020-7-9")-as.Date("2019-7-9") # 366

as.Date("2019-7-9")-as.Date("2018-7-9") # 365

## 顺序

x=as.Date(c("1970-01-03", "1969-12-30", "2019-07-09", "1900-0101"))

sort(x); order(x); rank(x)

二、时间

R自带的时间对象为POSIXct对象，它与Date对象的差异在于，前者除了年/月/日外还包括时钟时间。

1. 生成时间

x=Sys.time() # 获取当前时间

# "2019-07-09 16:08:31 CST"

class(x) # "POSIXct" "POSIXt"

x=as.POSIXct("2019-07-09 16:08:31 CST") # 将字符转化为POSIXct对象

as.character(x) # 将POSIXct对象转为字符

## 字符包含的时区会被忽略，因为R会强制使用本地时区；如果字符不包含时区,会被自动补上本地时区

x=c("2019-07-09 12:24:16 UTC", "2019-07-09 12:24:16", "2019-07-09 12:24:16 EST")

y=as.POSIXct(x) # "2019-07-09 12:24:16 CST" "2019-07-09 12:24:16 CST" "2019-07-09 12:24:16 CST"

y=as.POSIXct(x, tz="UTC") # 要想指定时区，必须使用tz参数

# 注意：时区UTC和GMT是相同的

# 笔者的时区为CST；若要显示此时区，应设定tz="Asia/Shanghai"

y=as.POSIXct(x, tz="Asia/Shanghai")

## 修改format参数

## 以下我们将使用"%H"、"%M"、"%OS"来提取小时、分钟和秒的信息

## as.POSIXct函数会自动尝试以下格式："%Y-%m-%d %H:%M:%OS"、"%Y/%m/%d %H:%M:%OS"、"%Y-%m-%d %H:%M"、"%Y/%m/%d %H:%M"、"%Y-%m-%d"、"%Y/%m/%d"

as.POSIXct(x="2019-07-09 13:01") # 不加秒

as.POSIXct("7/9-2019", format="%m/%d-%Y") # 只有日期

as.POSIXct("8:20:01 2019-07-09", tz="HST", format="%H:%M:%OS %Y-%m%d") # 自定义format

## 注意：在as.POSIXct中指定origin参数时会因时区问题出错

as.POSIXct(3600, origin="2019-07-09 13:56:40")

# [1] "2019-07-09 22:56:40 CST" # 结果并不是预想的"2019-07-09 14:56:40 UTC"，这是因为这种写法相当于as.POSIXct(3600, origin=as.POSIXct("2019 07-09 13:56:40", tz="UTC"), tz="Asia/Shanghai")

## 解决办法一，在两个位置同时用tz参数指定同一个时区

as.POSIXct(3600, origin=as.POSIXct("2019-07-09 13:56:40", tz="UTC"), tz="UTC")

as.POSIXct(3600, origin=as.POSIXct("2019-07-09 13:56:40", tz="Asia/Shanghai"), tz="Asia/Shanghai")

## 解决办法二，用+，得到本地时区的结果

as.POSIXct("2019-07-09 13:56:40")+3600

## 12小时制与24小时制："%I"代表12小时制的时间格式，它必须与"%p"搭配使用

# 当Sys.getlocale("LC_TIME")的值为中文简体时，"%p"与"上午"或"下午"匹配

as.POSIXct("2019-07-09 下午 8:24:16", format="%Y-%m-%d %p %I:%M:%OS")

# 当Sys.getlocale("LC_TIME")的值为英文时，"%p"与"am"或"pm"匹配（忽略大小写）

Sys.setlocale("LC_TIME", "English")

as.POSIXct("2019-07-09 PM 8:24:16", format="%Y-%m-%d %p %I:%M:%OS")

# 完成上述操作后请将区域修改成原样：Sys.setlocale("LC_TIME", "Chinese(Simplified)_China.936")

## 用seq生成序列

a=as.POSIXct("2019-07-09 13:56:40"); b=as.POSIXct("2019-07-09 14:56:40")

seq(a, b, length.out=10)

seq(a, b, by="2 min")

## round.POSIXt，trunc.POSIXt，可选units有"secs", "mins", "hours", "days","months", "years"

x=c("2019-05-16", "2019-05-17", "2019-04-16", "2019-02-15", "202002-15")

round.POSIXt(as.Date(x), "months")

# [1] "2019-05-01 CST" "2019-06-01 CST" "2019-05-01 CST" "2019-0301 CST" "2020-02-01 CST" # 注意：结果跟所在月份的天数有关

round.POSIXt(as.POSIXct("2019-07-09 15:30:18"), units="months") #注意，在保留月份时，结果后边仍然会附带上本月1日

trunc.POSIXt(as.POSIXct("2019-07-30 15:38:18"), units="months") #与round.POSIXct进行四舍五入不同，trunc.POSIXt只截取到所要求的位数，不会四舍五入

trunc.POSIXt(as.POSIXct("2019-05-16 15:38:18"), units="hours")

2. 比较和计算时间

对时间进行比较和计算的方式与对日期进行比较和计算的方式相仿。

a=as.POSIXct("2019-07-09 13:56:40"); b=as.POSIXct("2019-07-09 14:56:40")

x=seq(a, b, length.out=10)

b > a

mean(x); max(x); min(x); median(x)

y=as.character(summary(x)) # 简单汇总

a+c(3600, 7200)

difftime(b, a)

difftime(b, a, units="secs") # units默认为"auto"，即自动选择，可选择"secs"、"mins"、"hours"、"days"、"weeks"

3. 日期/时间对象的保存

char1=c("1998/8/31", "2019/7/9")

char2=c("1998-8-31 06:06:06", "2019-7-9 11:11:11")

date=as.Date(c("1998/8/31", "2019/7/9"))

time=as.POSIXct(c("1998-8-31 06:06:06", "2019-7-9 11:11:11"))

char3=paste(c("1998-8-31 06:06:06", "2019-7-9 11:11:11"), "abcde", sep="")

dat=data.frame(char1, char2, date, time, char3)

# write.csv(dat, "datetime.csv") # 打开csv文件后，会发现上述信息会以Excel默认的格式显示

# dat=read.csv("datetime.csv", row.names=1) # 但是读取文件后，发现格式正常

三、从日期/时间对象中提取信息

以上，我们用as.Date、as.POSIXct等函数生成日期/时间，但是，如果一个对象已经是日期/时间对象了，我们如何从中提取出年份、月份等信息并将其书写成我们需要的格式呢？

## 日期

x=as.Date("2019-07-09")

format(x) # "2019-07-09" # 相当于as.character(x)

format(x, "%m-%d") # "07-09" # 按照我们定义的格式输出

format(x, "这是月%m这是日%d")

format(x, "%m, %d, %A")

## 时间

x=as.POSIXct("2019-07-09 13:30:01")

format(x, "%Y-%m-%d")

format(x, "现在时间：%H时%M分")

format(x, "%p%I:%M") # 24小时和12小时制转换

## 时区的影响

x=as.POSIXct("2019-07-09 13:30:01") # "2019-07-09 13:30:01 CST" #时区为笔者所在的CST

format(x, "%H:%M:%OS") # "13:30:01" # 不修改时区(www.xing528.com)

format(x, "%H:%M:%OS", tz="UTC") # "05:30:01" # 将时区改为UTC，则会出现8个小时的变动

## 几个方便函数

x=as.POSIXct(c("1998-8-31 06:06:06", "2019-7-9 11:11:11", "2019-0709 15:24:16"))

weekdays(x, abbreviate=FALSE)

months(x, abbreviate=FALSE)

quarters(x) # 返回所在季度，取值为"Q1"、"Q2"、"Q3"、"Q4"

#==========

# 练习

#==========

## 对一组出生日期进行分组，将1970年1月1日以前出生的归为老年，将1990年1月1日以后出生的归为青年，将这两个日期之间的归为中年

birth=as.Date(c("1960-03-10", "1975-11-05", "1997-08-30", "1988-1230"))

limit1=as.Date("1990-01-01"); limit2=as.Date("1970-01-01")

g=rep(NA, length(birth))

for (i in 1: length(birth)){

ii=birth[i]

g[i]=if (ii < limit2) "老年" else if (ii >= limit2 & ii < limit1) "中年" else "青年"

}

四、时间序列对象

时间序列对象是时间与数值的结合。例如，为了储存5个月的失业率信息，我们可以使用数据框，数据框的第1列是采样时间，第2列是失业率数值；但是，我们在R中还可以使用时间序列对象。

dat=ts(c(0.06, 0.07, 0.05, 0.04, 0.06), frequency=12, start=1)

# Jan Feb Mar Apr May

# 1 0.06 0.07 0.05 0.04 0.06

class(dat) # "ts"

## 使用ts函数可以自动将采样数值与月份和季度匹配起来

ts(1: 15, frequency=12, start=c(1998, 2)) # 这里有15个数值需要分配给15个时间点，ts函数对这里的参数值的理解是：每年采样12次，也就是按月采样,从1998年2月开始采样

ts(1: 15, frequency=4, start=c(1998, 2)) # 与上例不同，这里的frequency= 4被ts函数视为一年采样4次，也就是按季度采样，而c(1998, 2)则代表从1998年第2季度采样

## 使用ts函数，但不让数值与月和季度匹配

ts(1: 14, frequency=6, start=3) # 从第3个时间位置开始采样，每个时间位置采样6次

# Start = c(3, 1)

# End = c(5, 2)

# Frequency = 6

# 输出结果的含义是：第1个数值来自第3个时间位置的第1次采样，最后一个数值是第5个时间位置的第2次采样，因此共有6+6+2=14个数值

ts(1: 14, frequency=6, start=c(3, 3)) # 现在改成，第1个数值是第3个时间位置的第3次采样

# Start = c(3, 3)

# End = c(5, 4)

# Frequency = 6

ts(1: 40, frequency=1, start=1978) # 一年一个数值

## 如果输入的是数据框，ts函数会将每一列当成一套独立的数据

x=data.frame(a=1: 12, b=101: 112)

dat=ts(x, frequency=4, start=c(2008, 1))

class(dat) # "mts" "ts" "matrix"

## 截取时间序列对象的一部分

x=ts(1: 31, frequency=12, start=c(1998, 1))

window(x, start=c(1998, 3), end=c(2000, 5)) # 截取1998年3月至2000年5月的数据

x=ts(1: 23, frequency=6, start=3)

window(x, start=c(3, 2), end=c(6, 5)) # 截取第3个时间位置第2次采样至第6个时间位置第5次采样之间的数据

五、lubridate包

lubridate包提供了很多用于处理日期/时间对象的方便函数。

# install.packages("lubridate")

library(lubridate)

## lubridate包提供了一些用于把字符转化成日期/时间对象的函数，这些函数比as.Date和as.POSIXct灵活得多

ymd(c("1998-08-31", "19980831", "1998/8/31", "1998-8/31", "1998, AUG 31st"))

# ymd(c("1998年8月31日", "1998年八月31日")) # 对中文的支持存在不确定性，所以请尽量用英文

mdy(c("8-31-1998", "08311998"))

dmy(c("31/8/1998", "31st/Aug/1998"))

ydm("1998 31st, 8")

ymd_hms(c("2019-07-09 20:08:03", "2019-07-09 8:08:03 pm", "2019-0709 8:08:03")) # 注意：lubridate包默认使用UTC时区

ymd_hms("2019-07-09 20:08:00", tz = "Asia/Shanghai") # 要使用CST时区，需手动设置tz函数

ymd_hm("2019-07-09 20:08")

ymd_h("2019-07-09 20")

# 同类函数还有dmy_hms、dmy_hm、dmy_h、mdy_hms、mdy_hm、mdy_h、ydm_hms、ydm_hm、ydm_h，我们根据名称就可猜到它们的用途

floor_date(ymd("2019-07-09"), "month") # "2019-07-01" # 向下取整

ceiling_date(ymd("2019-07-09"), "year") # "2020-01-01" # 向上取整

## 提取信息

x=ymd_hms("2019-07-09 20:08:03")

second(x) # 提取秒。返回数值而不是字符

## 同类函数还有minute、hour、day、year、tz

wday(x) # 返回数值形式的星期，1为周日

wday(x, label=TRUE) # 返回定序变量周日、周一、周二……

yday(x) # 一年中的第几天

week(x) # 一年中的第几周

month(x); month(x, label=TRUE) # 月份

leap_year(2020) # 判断是否是闰年

## 时间段

period_a=ymd("20190630") %--% ymd("20190731") # 生成第1个类型为Interval的时间段

period_b=ymd("20190715") %--% ymd("20190804") # 生成第2个时间段

int_overlaps(period_a, period_b) # 两个时间段是否有重合

intersect(period_a, period_b) # 如果有重合，重合的部分是什么

union(period_a, period_b) # 合并两个时间段

## Duration对象和Period对象

## 按照lubridate包的文档的解释，Duration对象是以秒计算的时间

dseconds(x=1); dminutes(x=1); dhours(x=1)

ddays(x=1); dweeks(x=1); dyears(x=1) #有的版本需要使用lubridate::: dmonths

## Period对象的表示方法有所不同

seconds(x=1); minutes(x=1); hours(x=1); days(x=1)

weeks(x=1); months(x=1); years(x=1)

ymd("2012-01-01")+dyears(1)

# "2012-12-31 06:00:00 UTC" # Duration对象不随闰年改变，所以这里相当于加上了365.25天

ymd("2012-01-01")+years(1)

# "2013-01-01" # Period对象随闰年改变，所以这里正确地加上了366天而不是365天

## 相加时出现日期不存在的情况

ymd("2012-01-31")+months(1) # NA # 加上1个月被认为是加上31天，而2月31日不存在

ymd("2012-01-31")+months(2) # "2012-03-31"

ymd("2012-01-31")+months(3) # NA # 4月31日也不存在

## 解决办法：用%m+%，确保输出确实存在的日期

ymd("2012-01-31") %m+% months(1) # "2012-02-29"

ymd("2012-01-31") %m+% months(3) # "2012-04-30"

## 修改单个元素

x=ymd_hms("2019-07-09 20:08:03")

second(x)=33

year(x)=2008

tz(x)="Asia/Shanghai" # 等同于force_tz(x, "Asia/Shanghai")

with_tz(x, "America/Chicago") # "2019-07-09 07:08:33 CDT" # 转化成另一个时区的时间

六、填充缺失值

1. 固定值填充和线性填充

尽管我们可以去掉包含缺失值的个案，但有时，为了得到完整或美观的图表，我们还需保留这些个案并把缺失值填充上。在处理跟日期和时间相关的缺失值时，我们可使用zoo包。

# install.packages("zoo")

library(zoo)

## 我们使用na.approx函数进行线性填充，但在这之前，必须把待填充的向量转化成zoo对象

x=zoo(c(NA, 1, NA, NA, NA, 3, NA, 7, 4, NA))

class(x) # "zoo"

y=na.approx(x, na.rm=FALSE, method="linear")

# 注意：有时第一个值或最后一个值是缺失值，默认情况下，这些值不但不会被填充而且会被删除；为保持结果的项数与输入的项数一致，务必设置na.rm=FALSE

as.numeric(y) # 将填充结果转为数值向量

## 使用固定值填充。默认设置是，将数值按顺序排列，用缺失值左边的非缺失值填充

na.approx(x, na.rm=FALSE, method="constant")

# 当把默认的f=0改为f=1时，则用缺失值右边的非缺失值填充

na.approx(x, na.rm=FALSE, method="constant", f=1)

# 事实上，这里的f可以改成0与1之间的任意数

na.approx(x, na.rm=FALSE, method="constant", f=0.25) # 以1与3之间的缺失值为例，(3-1)*0.25=0.5，1+0.5=1.5，所以这些缺失值被替换成了1.5

## zoo函数的索引号参数order.by

# 默认情况下，索引号是连续的

x=zoo(c(10, NA, 30, 40, 50)) # 相当于zoo(c(10, NA, 30, 40, 50), order.

by=1: 5)

na.approx(x)

# 但索引号也可以是非连续的

x=zoo(c(10, NA, 50, 30, 40), order.by=c(1, 2, 5, 3, 4))

na.approx(x)

# 索引号不连续的情境是，例如，尽管我们按照时间点1、2、3、4、5采样，但是当把数值录入到表格里时，并没有按顺序录入

dat=data.frame(v=c(10, NA, 50, 30, 40), time=c(1, 2, 5, 3, 4))

x=zoo(dat$v, order.by=dat$time)

na.approx(x)

# 注意：索引号不同时，结果也不同。所以，为防止出错，尽量使用连续的索引号

x=zoo(c(0, NA, NA, 1), order.by=c(1, 2, 3, 4)); na.approx(x)

x=zoo(c(0, NA, NA, 1), order.by=c(1, 2, 3, 5)); na.approx(x)

## 使用na.locf填充头部和尾部的缺失值

x=zoo(c(NA, NA, 1, NA, 2, 3, NA, NA))

y=na.approx(x, na.rm=FALSE) # 此时向量头部和尾部仍是缺失值

y=na.locf(y, na.rm=FALSE) # 填充尾部缺失值

y=na.locf(y, na.rm=FALSE, from Last=TRUE) # 填充头部缺失值

# 可见，要填充全部缺失值，我们需同时使用na.approx和na.locf

## 现在假设我们要对1至6月的数据作图，但实际上我们只有2、3、5月的数据,这意味着我们需要填充1、4、6月的缺失值

d1=as.Date(c("2019-02-01", "2019-03-01", "2019-05-01")) # 3个有数据的月份

x1=c(22, 33, 55) # 与上述3个月份相对应的数值

# 接下来开始填充：第1步，生成完整的日期

d2=seq(as.Date("2019-01-01"), as.Date("2019-06-01"), by="month")

# 第2步，生成完整的带缺失值的数据，并将非缺失值填入

x2=rep(NA, length(d2))

pos=match(d1, d2)

x2[pos]=x1

# 第3步，填充

# 方法1：以日期为索引号

z=zoo(x2, order.by=d2)

result=na.approx(z, na.rm=FALSE)

# 方法2：以1、2、3……为索引号，结果与方法1不同

z=zoo(x2, order.by=1: length(x2))

result=na.approx(z, na.rm=FALSE)

# 第4步，填充两边的缺失值

result=na.locf(result, na.rm=FALSE)

result=na.locf(result, na.rm=FALSE, from Last=TRUE)

## 上述方法1和方法2的区别在以下例子中更为明显

z=zoo(c(1, NA, 3, NA, 5), order.by=as.Date(c("2019-01-01", "2019-0102", "2019-01-03", "2019-01-05", "2019-01-31")))

na.approx(z, na.rm=FALSE) # 1月2日与1月1日、与1月3日是等间隔的,所以缺失值被替换成了1日与3日的中间值2；相反，1月5日与1月3日、与1月31日不是等间隔的，所以缺失值没有被替换成中间值4

z=zoo(c(1, NA, 3, NA, 5), order.by=1: 5)

na.approx(z, na.rm=FALSE) # 如果只以1、2、3……为索引号，则日期之间的真实间隔被忽略，所以缺失值被替换成了2和4

2. 其他填充方法

用zoo包中的na.spline函数可实现样条填充，使用方法与na.approx相同。impute TS包中的na.interpolation函数也可实现以线性方式填充。

# install.packages("impute TS")

library(impute TS)

x=c(NA, 1, 3, NA, 5, 6, NA, NA, NA, 9, NA)

na.interpolation(x, option="linear") # option还可选"spline"、"stine"

此外，impute TS包还包含其他填充方法。

na.kalman(x, model="Struct TS") # 用Struct TS函数拟合模型并填充

na.kalman(x, model="auto.arima") # 用auto.arima函数拟合模型并填充

na.ma(x, k=2, weighting="exponential") # 加权移动平均法。参数k指定使用多少个非缺失值，例如k=2代表使用左边2个+右边2个=4个非缺失值；weighting为加权方法，默认为"exponential"，还可选择"simple"或"linear"

na.mean(x, option="mean") # option还可以是"mean"、"median"、"mode",代表以均值、中位数或众数填充

na.replace(x, fill=9999) # 填充一个特定数