使用 zoo 包进行零膨胀滚动 beta 回归？

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

免责声明：我有几个 cmets - 不仅仅是一个答案 - 但因为我想显示代码和输出并且需要更多空间，所以我以答案的形式来做。

首先，在您的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（移动和）测试。下面我展示了带宽为 15%（即 7 个​​观察值）及其 5% 临界值的 MOSUM 测试的可视化。相应的 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))

Answer 2

编辑：由我的朋友解决。为了记录，如果有人关心：这是一个窗口宽度的问题，我以为我已经玩过了——但还不够。 Beta 回归模型在估计窗口的每次迭代的系数时遇到问题，因此我们使用 Beta 回归创建了一个循环，跟踪其进度，并查看错误何时出现：

brr.function <- function(brr) {
  coef(betareg(y.transf.betareg(brr[,3])~brr[,2],
              data=as.data.frame(brr),
              na.action = na.omit))
}

a <- NULL 
total.obs <- nrow(brol)
mw <- 5  # window length

for (i in 1:c(total.obs-mw)){
  a<-c(a,brr.function(brol[i:c(i+mw),]))
  cat("i=",i,"\n") # this code tracks the progress of the loop
}

我们看到它停在 12 点，所以我们检查了那段数据：

i <- 12
brol[i:c(i+mw),]

Year    World        Ego
12 1779 0.722963 0.05291005
13 1781 0.000000 0.00000000
14 1782 0.000000 0.00000000
15 1783 0.750000 0.00000000

然后我们设置窗口宽度为4，代码运行。