如何对面板数据执行线性和趋势外推？
How can I perform linear and trend extrapolation on panel data?
我正在尝试使用这行代码推断我的数据中的以下缺失值 (NA)，但它不起作用。
我的数据：
landkreis jahr deDomains
<chr> <dbl> <dbl>
1 Ahrweile… 2007 NA
2 Ahrweile… 2008 NA
3 Ahrweile… 2009 NA
4 Ahrweile… 2010 NA
5 Ahrweile… 2011 NA
6 Ahrweile… 2012 NA
7 Ahrweile… 2013 22224
8 Ahrweile… 2014 22460
9 Ahrweile… 2015 2379
10 Ahrweile… 2016 22769
11 Ahrweile… 2017 23268
12 Aichach-… 2007 NA
13 Aichach-… 2008 NA
14 Aichach-… 2009 NA
15 Aichach-… 2010 NA
16 Aichach-… 2011 NA
17 Aichach-… 2012 NA
18 Aichach-… 2013 21341
19 Aichach-… 2014 21393
20 Aichach-… 2015 21338
我正在尝试使用以下代码扩展 deDomains 变量上的 NA，但它不起作用
df_complete <- df_complete %>%
group_by(landkreis) %>%
mutate(`deDomains` = approxExtrap(which(!is.na(`deDomains`)),
`deDomains`[!is.na(`deDomains`)])$y)
我正在使用 Hmisc 包中的 approxExtrap() 命令进行线性外推。
您需要指定您的xout。 NAs 实际上是由函数处理的。您可能想查看 approx 函数，您可以在其中找到一些示例（尽管 interpolation，但它是相似的）；输入 ?approx.
library(dplyr)
library(Hmisc)
df_complete %>%
group_by(landkreis) %>%
mutate(`deDomains`=approxExtrap(x=jahr, y=deDomains, xout=jahr)$y)
# # A tibble: 20 x 3
# # Groups: landkreis [2]
# landkreis jahr deDomains
# <fct> <int> <dbl>
# 1 Ahrweile… 2007 22224
# 2 Ahrweile… 2008 22224
# 3 Ahrweile… 2009 22224
# 4 Ahrweile… 2010 22224
# 5 Ahrweile… 2011 22224
# 6 Ahrweile… 2012 22224
# 7 Ahrweile… 2013 22224
# 8 Ahrweile… 2014 22460
# 9 Ahrweile… 2015 2379
# 10 Ahrweile… 2016 22769
# 11 Ahrweile… 2017 23268
# 12 Aichach-… 2007 21341
# 13 Aichach-… 2008 21341
# 14 Aichach-… 2009 21341
# 15 Aichach-… 2010 21341
# 16 Aichach-… 2011 21341
# 17 Aichach-… 2012 21341
# 18 Aichach-… 2013 21341
# 19 Aichach-… 2014 21393
# 20 Aichach-… 2015 21338
或使用by:
library(Hmisc)
do.call(rbind, by(df_complete, df_complete$landkreis, function(x) {
transform(x,
deDomains=approxExtrap(x=x$jahr, y=x$deDomains, xout=x$jahr)$y
)
}))
# landkreis jahr deDomains
# Ahrweile….1 Ahrweile… 2007 22224
# Ahrweile….2 Ahrweile… 2008 22224
# Ahrweile….3 Ahrweile… 2009 22224
# Ahrweile….4 Ahrweile… 2010 22224
# Ahrweile….5 Ahrweile… 2011 22224
# Ahrweile….6 Ahrweile… 2012 22224
# Ahrweile….7 Ahrweile… 2013 22224
# Ahrweile….8 Ahrweile… 2014 22460
# Ahrweile….9 Ahrweile… 2015 2379
# Ahrweile….10 Ahrweile… 2016 22769
# Ahrweile….11 Ahrweile… 2017 23268
# Aichach-….12 Aichach-… 2007 21341
# Aichach-….13 Aichach-… 2008 21341
# Aichach-….14 Aichach-… 2009 21341
# Aichach-….15 Aichach-… 2010 21341
# Aichach-….16 Aichach-… 2011 21341
# Aichach-….17 Aichach-… 2012 21341
# Aichach-….18 Aichach-… 2013 21341
# Aichach-….19 Aichach-… 2014 21393
# Aichach-….20 Aichach-… 2015 21338
编辑： 要使用 "trend" 进行推断，您可以使用例如na_kalman 来自 imputeTS 包。
library(imputeTS)
res <- do.call(rbind, by(df_complete, df_complete$landkreis, function(x) {
transform(x,
deDomains.ex=na_kalman(x$deDomains, model = "StructTS", smooth = TRUE)
)
}))
# landkreis jahr deDomains deDomains.ex
# Ahrweile….1 Ahrweile… 2007 NA 21532.16
# Ahrweile….2 Ahrweile… 2008 NA 21186.24
# Ahrweile….3 Ahrweile… 2009 NA 20840.32
# Ahrweile….4 Ahrweile… 2010 NA 20494.40
# Ahrweile….5 Ahrweile… 2011 NA 20148.48
# Ahrweile….6 Ahrweile… 2012 NA 19802.56
# Ahrweile….7 Ahrweile… 2013 22224 22224.00
# Ahrweile….8 Ahrweile… 2014 22460 22460.00
# Ahrweile….9 Ahrweile… 2015 2379 2379.00
# Ahrweile….10 Ahrweile… 2016 22769 22769.00
# Ahrweile….11 Ahrweile… 2017 23268 23268.00
# Aichach-….12 Aichach-… 2007 NA 21344.52
# Aichach-….13 Aichach-… 2008 NA 21346.28
# Aichach-….14 Aichach-… 2009 NA 21348.04
# Aichach-….15 Aichach-… 2010 NA 21349.80
# Aichach-….16 Aichach-… 2011 NA 21351.55
# Aichach-….17 Aichach-… 2012 NA 21353.31
# Aichach-….18 Aichach-… 2013 21341 21341.00
# Aichach-….19 Aichach-… 2014 21393 21393.00
# Aichach-….20 Aichach-… 2015 21338 21338.00
可能有更好的数据来演示，不过还是看个图吧：
plot(deDomains ~ jahr, type="n", data=res)
sapply(seq(res$landkreis), function(x)
with(res[res$landkreis == unique(res$landkreis)[x], ],
{lines(jahr, deDomains.ex, col=x + 1)
points(jahr, deDomains, col=x + 1)}))
legend("bottomleft", legend=c(as.character(unique(res$landkreis)), "true points"),
col=c(2, 3, 1), lty=c(1, 1, NA), pch=c(NA, NA, 1))
您还可以查看 imputeTS::na_seadec 函数，其中 - 在卡尔曼中 - 可以选择其他算法，也可以检测频率。
数据：
df_complete <- structure(list(landkreis = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Ahrweile…",
"Aichach-…"), class = "factor"), jahr = c(2007L, 2008L, 2009L,
2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2007L,
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L), deDomains = c(NA,
NA, NA, NA, NA, NA, 22224L, 22460L, 2379L, 22769L, 23268L, NA,
NA, NA, NA, NA, NA, 21341L, 21393L, 21338L)), class = "data.frame", row.names = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20"))
我正在尝试使用这行代码推断我的数据中的以下缺失值 (NA)，但它不起作用。
我的数据：
landkreis jahr deDomains
<chr> <dbl> <dbl>
1 Ahrweile… 2007 NA
2 Ahrweile… 2008 NA
3 Ahrweile… 2009 NA
4 Ahrweile… 2010 NA
5 Ahrweile… 2011 NA
6 Ahrweile… 2012 NA
7 Ahrweile… 2013 22224
8 Ahrweile… 2014 22460
9 Ahrweile… 2015 2379
10 Ahrweile… 2016 22769
11 Ahrweile… 2017 23268
12 Aichach-… 2007 NA
13 Aichach-… 2008 NA
14 Aichach-… 2009 NA
15 Aichach-… 2010 NA
16 Aichach-… 2011 NA
17 Aichach-… 2012 NA
18 Aichach-… 2013 21341
19 Aichach-… 2014 21393
20 Aichach-… 2015 21338
我正在尝试使用以下代码扩展 deDomains 变量上的 NA，但它不起作用
df_complete <- df_complete %>%
group_by(landkreis) %>%
mutate(`deDomains` = approxExtrap(which(!is.na(`deDomains`)),
`deDomains`[!is.na(`deDomains`)])$y)
我正在使用 Hmisc 包中的 approxExtrap() 命令进行线性外推。
您需要指定您的xout。 NAs 实际上是由函数处理的。您可能想查看 approx 函数，您可以在其中找到一些示例（尽管 interpolation，但它是相似的）；输入 ?approx.
library(dplyr)
library(Hmisc)
df_complete %>%
group_by(landkreis) %>%
mutate(`deDomains`=approxExtrap(x=jahr, y=deDomains, xout=jahr)$y)
# # A tibble: 20 x 3
# # Groups: landkreis [2]
# landkreis jahr deDomains
# <fct> <int> <dbl>
# 1 Ahrweile… 2007 22224
# 2 Ahrweile… 2008 22224
# 3 Ahrweile… 2009 22224
# 4 Ahrweile… 2010 22224
# 5 Ahrweile… 2011 22224
# 6 Ahrweile… 2012 22224
# 7 Ahrweile… 2013 22224
# 8 Ahrweile… 2014 22460
# 9 Ahrweile… 2015 2379
# 10 Ahrweile… 2016 22769
# 11 Ahrweile… 2017 23268
# 12 Aichach-… 2007 21341
# 13 Aichach-… 2008 21341
# 14 Aichach-… 2009 21341
# 15 Aichach-… 2010 21341
# 16 Aichach-… 2011 21341
# 17 Aichach-… 2012 21341
# 18 Aichach-… 2013 21341
# 19 Aichach-… 2014 21393
# 20 Aichach-… 2015 21338
或使用by:
library(Hmisc)
do.call(rbind, by(df_complete, df_complete$landkreis, function(x) {
transform(x,
deDomains=approxExtrap(x=x$jahr, y=x$deDomains, xout=x$jahr)$y
)
}))
# landkreis jahr deDomains
# Ahrweile….1 Ahrweile… 2007 22224
# Ahrweile….2 Ahrweile… 2008 22224
# Ahrweile….3 Ahrweile… 2009 22224
# Ahrweile….4 Ahrweile… 2010 22224
# Ahrweile….5 Ahrweile… 2011 22224
# Ahrweile….6 Ahrweile… 2012 22224
# Ahrweile….7 Ahrweile… 2013 22224
# Ahrweile….8 Ahrweile… 2014 22460
# Ahrweile….9 Ahrweile… 2015 2379
# Ahrweile….10 Ahrweile… 2016 22769
# Ahrweile….11 Ahrweile… 2017 23268
# Aichach-….12 Aichach-… 2007 21341
# Aichach-….13 Aichach-… 2008 21341
# Aichach-….14 Aichach-… 2009 21341
# Aichach-….15 Aichach-… 2010 21341
# Aichach-….16 Aichach-… 2011 21341
# Aichach-….17 Aichach-… 2012 21341
# Aichach-….18 Aichach-… 2013 21341
# Aichach-….19 Aichach-… 2014 21393
# Aichach-….20 Aichach-… 2015 21338
编辑： 要使用 "trend" 进行推断，您可以使用例如na_kalman 来自 imputeTS 包。
library(imputeTS)
res <- do.call(rbind, by(df_complete, df_complete$landkreis, function(x) {
transform(x,
deDomains.ex=na_kalman(x$deDomains, model = "StructTS", smooth = TRUE)
)
}))
# landkreis jahr deDomains deDomains.ex
# Ahrweile….1 Ahrweile… 2007 NA 21532.16
# Ahrweile….2 Ahrweile… 2008 NA 21186.24
# Ahrweile….3 Ahrweile… 2009 NA 20840.32
# Ahrweile….4 Ahrweile… 2010 NA 20494.40
# Ahrweile….5 Ahrweile… 2011 NA 20148.48
# Ahrweile….6 Ahrweile… 2012 NA 19802.56
# Ahrweile….7 Ahrweile… 2013 22224 22224.00
# Ahrweile….8 Ahrweile… 2014 22460 22460.00
# Ahrweile….9 Ahrweile… 2015 2379 2379.00
# Ahrweile….10 Ahrweile… 2016 22769 22769.00
# Ahrweile….11 Ahrweile… 2017 23268 23268.00
# Aichach-….12 Aichach-… 2007 NA 21344.52
# Aichach-….13 Aichach-… 2008 NA 21346.28
# Aichach-….14 Aichach-… 2009 NA 21348.04
# Aichach-….15 Aichach-… 2010 NA 21349.80
# Aichach-….16 Aichach-… 2011 NA 21351.55
# Aichach-….17 Aichach-… 2012 NA 21353.31
# Aichach-….18 Aichach-… 2013 21341 21341.00
# Aichach-….19 Aichach-… 2014 21393 21393.00
# Aichach-….20 Aichach-… 2015 21338 21338.00
可能有更好的数据来演示，不过还是看个图吧：
plot(deDomains ~ jahr, type="n", data=res)
sapply(seq(res$landkreis), function(x)
with(res[res$landkreis == unique(res$landkreis)[x], ],
{lines(jahr, deDomains.ex, col=x + 1)
points(jahr, deDomains, col=x + 1)}))
legend("bottomleft", legend=c(as.character(unique(res$landkreis)), "true points"),
col=c(2, 3, 1), lty=c(1, 1, NA), pch=c(NA, NA, 1))
您还可以查看 imputeTS::na_seadec 函数，其中 - 在卡尔曼中 - 可以选择其他算法，也可以检测频率。
数据：
df_complete <- structure(list(landkreis = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Ahrweile…",
"Aichach-…"), class = "factor"), jahr = c(2007L, 2008L, 2009L,
2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2007L,
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L), deDomains = c(NA,
NA, NA, NA, NA, NA, 22224L, 22460L, 2379L, 22769L, 23268L, NA,
NA, NA, NA, NA, NA, 21341L, 21393L, 21338L)), class = "data.frame", row.names = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20"))