使用 zoo 包零膨胀滚动 beta 回归？
Zero inflated rolling beta regression with zoo package?
我正在尝试评估两个变量之间的关系随时间的变化。我使用 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 link 的二项式 GLM 或 GAMLSS，但不适用于 betareg。
从一些研究来看，传递给函数的每条数据似乎都不是满秩的，但我不知道如何处理这个问题。谁能建议？非常非常感谢。
编辑：由我的朋友解决。作为记录，如果有人关心的话：这是一个 window 宽度的问题，我认为我已经尝试过了 - 但还不够。 beta 回归模型在估计 window 的每次迭代的系数时遇到问题，因此我们使用 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
然后我们设置window宽度为4，代码运行。
免责声明：我有几条评论 - 不仅仅是一个答案 - 但由于我想显示代码和输出并且需要更多 space，我以答案的形式进行。
首先，在您的 y.transf.betareg 中，您想使用 betareg 小插图中推荐的转换。由于 "number of observations" 您使用 46，即数据中的时间点数。但是，对于校正项，应该使用计算比例的观察次数（如果适用）。例如，在 1761 年，变量 World 是 0.3055556，这可能来自大约 11/36。如果是这样的话，那么观察次数应该是36。
其次，您正在拟合的 beta 回归模型具有三个参数（截距、斜率、精度），因此使用滚动 window 大小的四个观测值非常 - 比方说 - 乐观。我无法想象在您的应用程序中这不仅仅是随机噪声。
因此，我建议首先找到一个您希望大致适用于您的数据的模型。鉴于左删失为零，一个自然的候选者似乎是一个简单的托比特模型。数据的散点图以及托比特回归线是：
复制代码在最后。总体而言，该模型似乎与显着的斜率估计相当吻合：
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))
我正在尝试评估两个变量之间的关系随时间的变化。我使用 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 link 的二项式 GLM 或 GAMLSS，但不适用于 betareg。
从一些研究来看，传递给函数的每条数据似乎都不是满秩的，但我不知道如何处理这个问题。谁能建议？非常非常感谢。
编辑：由我的朋友解决。作为记录，如果有人关心的话：这是一个 window 宽度的问题，我认为我已经尝试过了 - 但还不够。 beta 回归模型在估计 window 的每次迭代的系数时遇到问题，因此我们使用 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
然后我们设置window宽度为4，代码运行。
免责声明：我有几条评论 - 不仅仅是一个答案 - 但由于我想显示代码和输出并且需要更多 space，我以答案的形式进行。
首先，在您的 y.transf.betareg 中，您想使用 betareg 小插图中推荐的转换。由于 "number of observations" 您使用 46，即数据中的时间点数。但是，对于校正项，应该使用计算比例的观察次数（如果适用）。例如，在 1761 年，变量 World 是 0.3055556，这可能来自大约 11/36。如果是这样的话，那么观察次数应该是36。
其次，您正在拟合的 beta 回归模型具有三个参数（截距、斜率、精度），因此使用滚动 window 大小的四个观测值非常 - 比方说 - 乐观。我无法想象在您的应用程序中这不仅仅是随机噪声。
因此，我建议首先找到一个您希望大致适用于您的数据的模型。鉴于左删失为零，一个自然的候选者似乎是一个简单的托比特模型。数据的散点图以及托比特回归线是：
复制代码在最后。总体而言，该模型似乎与显着的斜率估计相当吻合：
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))