如何在面板数据上进行线性和趋势外推？

Question

如何在面板数据上进行线性和趋势外推？

3

我试图使用下面的代码来推断我数据中的缺失值 (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() 命令进行线性外推。

- Laurenz Hamel

请在您的问题中发布dput(df_complete)的结果。 - Cettt

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jay.sf · Accepted Answer

您需要指定您的xout。函数会处理NA。您可以查看approx函数，可以找到一些示例（用于插值，但类似）；键入?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

编辑： 如果要使用“趋势”来进行推断，您可以使用imputeTS包中的na_kalman。

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"))