r - 动物园包的零膨胀滚动贝塔回归？

Question

我正在尝试评估两个变量之间的关系随时间的变化。我使用 Zoo 包创建了 46 年的不规则时间序列对象。我的数据是取值为 0 和 1 的零膨胀比例。以下是数据:

edf
   Year     World        Ego
1  1760 1.0000000 0.00000000
2  1761 0.3055556 0.00000000
3  1762 0.3950617 0.11814413
4  1764 0.8677686 0.26984127
5  1766 0.0000000 0.00000000
6  1767 0.8580606 0.15407986
7  1769 0.7500000 0.00000000
8  1771 0.7416174 0.37698413
9  1772 0.6611570 0.53587372
10 1777 0.4375000 0.20000000
11 1778 0.9629630 0.36111111
12 1779 0.7229630 0.05291005
13 1781 0.0000000 0.00000000
14 1782 0.0000000 0.00000000
15 1783 0.7500000 0.00000000
16 1784 0.7966605 0.21893984
17 1785 0.8518519 0.12500000
18 1786 0.0000000 0.00000000
19 1787 0.2279036 0.00000000
20 1788 0.7425926 0.08585859
21 1789 0.4648760 0.17942337
22 1790 0.8888889 0.00000000
23 1791 0.7958546 0.35023819
24 1792 0.0000000 0.00000000
25 1794 0.8021333 0.65529337
26 1795 0.0000000 0.00000000
27 1800 0.9900000 0.10825397
28 1802 0.7866667 0.07500000
29 1803 0.0000000 0.00000000
30 1804 0.0000000 0.00000000
31 1805 0.7416026 0.34158521
32 1806 0.9420000 0.47337963
33 1810 0.7500000 0.00000000
34 1812 0.8397279 0.53089503
35 1818 0.4863946 0.31103450
36 1819 0.8636475 0.20591162
37 1820 0.8888889 0.00000000
38 1821 0.7197232 0.60557261
39 1822 0.7308806 0.27126586
40 1823 0.6113805 0.26487719
41 1824 0.6400000 0.00000000
42 1826 0.9086405 0.13932918
43 1827 0.7447051 0.16207173
44 1828 0.9183673 0.40000000
45 1830 0.9843750 0.50000000
46 1831 0.7053061 0.55736111

我正在使用 beta 回归，但使用手册中的建议转换因变量值:

y.transf.betareg <- function(y){
  n.obs <- sum(!is.na(y))
  (y * (n.obs - 1) + 0.5) / n.obs
}

然后使用 rollapply 计算移动回归。这是我的代码:

library(zoo)
library(betareg)
brol<-as.zoo(edf)
index1 <- rollapply(data = brol,  
                          width = 5,  
                          function(brr)  coef(betareg(y.transf.betareg(brr[3])~brr[2],
                                            data=as.data.frame(brr),
                                            na.action = na.omit
                                    ),
                      by.column = F,
                      align="right")) 

但我收到此错误:

Error in optim(par = start, fn = loglikfun, gr = gradfun, method = method,  : 
  non-finite value supplied by optim

当我尝试将线性样条回归与 betareg 结合使用时，我遇到了相同的错误。

我编写的代码适用于我尝试过的其他模型，例如带有 logit 链接的二项式 GLM 或 GAMLSS，但不适用于 betareg。

经过一番研究，似乎传递给函数的每条数据可能都不是满级的，但我不知道如何处理这个问题。有人可以建议吗？非常非常感谢。

Answer 1

免责声明:我有几条评论 - 不仅仅是一个答案 - 但由于我想显示代码和输出并且需要更多空间，所以我以答案的形式进行。

首先，在您的 y.transf.betareg 中，您希望使用 betareg 小插图中推荐的转换。您使用 46 作为“观察数”，即数据中的时间点数量。然而，对于校正项，应使用计算比例的观测值数量(如果适用)。例如，在 1761 年，变量 World 为 0.3055556，可能来自大约 11/36。如果是这样的话，那么观测数应该是 36。

其次，您正在拟合的 beta 回归模型具有三个参数(截距、斜率、精度)，因此使用四个观测值的滚动窗口大小是非常乐观的。我无法想象在您的应用程序中这只是随机噪声。

因此，我建议首先找到一个您期望大致适用于您的数据的模型。鉴于左审查为零，自然的候选者似乎是简单的 tobit 模型。数据散点图以及 tobit 回归线为:

复制代码包含在末尾。总体而言，该模型似乎与显着的斜率估计相当吻合:

           Estimate Std. Error  z value  Pr(>|z|)    
(Intercept) -0.27869    0.12527  -2.2247 0.0261034 *  
World        0.61259    0.16398   3.7358 0.0001871 ***
Log(scale)  -1.40915    0.14001 -10.0645 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

然后，您似乎有兴趣评估参数在采样期间是否稳定。然而，由于样本量有限，对每个子样本重新估计模型会导致大量的随机波动。相反，我们可以采用模型分数的滚动总和(如果您愿意的话，可以是一种残差)作为对整个样本的估计。如果任何参数发生系统性变化，您将在滚动总和的系统性变化中看到它们。这种类型的测试在结构变化文献中也称为基于分数的 MOSUM(移动和)测试。下面我展示了 MOSUM 测试的可视化，带宽为 15%(即 7 个​​观测值)及其 5% 临界值。相应的 p 值为 49.6%，因此显然不显着。该图显示没有系统偏离零。

因此，在样本大小适中的情况下，我们无法检测到与上面给出的参数拟合的单个模型的显着偏差。 (增加或减少带宽的 MOSUM 测试会得到相同的结果。)

复制代码:

数据

library("zoo")
edf <- read.zoo(textConnection("   Year     World        Ego
1  1760 1.0000000 0.00000000
2  1761 0.3055556 0.00000000
3  1762 0.3950617 0.11814413
4  1764 0.8677686 0.26984127
5  1766 0.0000000 0.00000000
6  1767 0.8580606 0.15407986
7  1769 0.7500000 0.00000000
8  1771 0.7416174 0.37698413
9  1772 0.6611570 0.53587372
10 1777 0.4375000 0.20000000
11 1778 0.9629630 0.36111111
12 1779 0.7229630 0.05291005
13 1781 0.0000000 0.00000000
14 1782 0.0000000 0.00000000
15 1783 0.7500000 0.00000000
16 1784 0.7966605 0.21893984
17 1785 0.8518519 0.12500000
18 1786 0.0000000 0.00000000
19 1787 0.2279036 0.00000000
20 1788 0.7425926 0.08585859
21 1789 0.4648760 0.17942337
22 1790 0.8888889 0.00000000
23 1791 0.7958546 0.35023819
24 1792 0.0000000 0.00000000
25 1794 0.8021333 0.65529337
26 1795 0.0000000 0.00000000
27 1800 0.9900000 0.10825397
28 1802 0.7866667 0.07500000
29 1803 0.0000000 0.00000000
30 1804 0.0000000 0.00000000
31 1805 0.7416026 0.34158521
32 1806 0.9420000 0.47337963
33 1810 0.7500000 0.00000000
34 1812 0.8397279 0.53089503
35 1818 0.4863946 0.31103450
36 1819 0.8636475 0.20591162
37 1820 0.8888889 0.00000000
38 1821 0.7197232 0.60557261
39 1822 0.7308806 0.27126586
40 1823 0.6113805 0.26487719
41 1824 0.6400000 0.00000000
42 1826 0.9086405 0.13932918
43 1827 0.7447051 0.16207173
44 1828 0.9183673 0.40000000
45 1830 0.9843750 0.50000000
46 1831 0.7053061 0.55736111"), header = TRUE)

完整示例模型

library("AER")
m <- tobit(Ego ~ World, data = edf)
coeftest(m)

散点图

plot(jitter(Ego, 10) ~ jitter(World, 10), data = edf,
  xlab = "World (jittered)", ylab = "Ego (jittered)")
abline(m)
legend("topleft", "Tobit model", lwd = 1, bty = "n")

MOSUM 测试

library("strucchange")
sctest(m, order.by = time(edf), functional = maxMOSUM(0.15), 
  plot = TRUE, aggregate = FALSE, ylim = c(-1.5, 1.5))