零膨胀滚动β回归动物园包？

Question

随着时间的推移，我试图评估两个变量之间关系的变化。我用动物园的包装创造了一个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链接或GAMLSS的二项式GLM，而不是betareg。

从一些研究来看，似乎传递给函数的每一段数据都不是满级的，但我不知道如何处理。有人能告诉我吗？非常，非常感谢。

Answer 1

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

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

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

因此，我建议首先找到一个模型，您可以期望该模型对您的数据具有大致的支持作用。考虑到有左审查在零，一个自然的候选人似乎是一个简单的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测试的可视化(即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))

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，代码运行。