用于识别年龄频率数据相关性的滞后线性模型
Lagged linear model to identify correlation in age frequency data
我有这些数据，我正在尝试在 r 中进行滞后线性回归以确定 YOY 的数量是否与下一年的 1 岁儿童和后年的 2 岁儿童的数量显着相关...等等...
数据：
structure(list(Year = c("2008", "2009", "2010", "2011", "2012",
"2013", "2014", "2015", "2016", "2017", "2018", "2007", "2007",
"2007", "2007", "2008", "2008", "2008", "2009", "2009", "2009",
"2009", "2009", "2009", "2009", "2010", "2010", "2010", "2010",
"2010", "2011", "2011", "2011", "2011", "2011", "2011", "2011",
"2011", "2011", "2012", "2012", "2012", "2012", "2012", "2012",
"2012", "2012", "2013", "2013", "2013", "2013", "2013", "2013",
"2013", "2013", "2014", "2014", "2014", "2014", "2014", "2014",
"2014", "2014", "2014", "2015", "2015", "2015", "2015", "2015",
"2015", "2015", "2015", "2015", "2016", "2016", "2016", "2016",
"2016", "2016", "2016", "2017", "2017", "2017", "2017", "2017",
"2017", "2017", "2018", "2018", "2018", "2018", "2018", "2018",
"2018", "2018"), Age = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 3L, 6L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 2L, 3L, 4L, 5L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,
9L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L,
8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 10L, 2L, 3L, 4L, 5L, 6L, 7L, 10L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), .Label = c("0", "1",
"2", "3", "4", "5", "6", "7", "8", "9"), class = "factor"), n = c(166,
28, 34, 77, 170, 18, 3, 22, 43, 50, 151, 1, 8, 17, 1, 4, 19,
1, 1, 46, 37, 52, 5, 1, 1, 19, 41, 15, 16, 1, 1, 13, 4, 26, 12,
11, 1, 1, 1, 1, 87, 15, 13, 27, 13, 17, 1, 1, 32, 30, 3, 4, 1,
1, 1, 1, 24, 15, 23, 6, 2, 1, 2, 2, 4, 18, 13, 31, 28, 3, 3,
6, 1, 4, 6, 1, 5, 9, 1, 1, 1, 16, 16, 8, 1, 1, 4, 1, 12, 4, 7,
2, 1, 2, 1), id = c("YOY", "YOY", "YOY", "YOY", "YOY", "YOY",
"YOY", "YOY", "YOY", "YOY", "YOY", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult")), row.names = c(NA, -95L), class = "data.frame")
我制作了很棒的情节，表明这里确实看起来确实有东西。不完美，但有某种关系。
# Frequencey density plot of ages over year
ggplot(wi.age.count, aes(x=Year, y=Age)) +
geom_point(aes(cex = n, color = id)) +
#scale_fill_brewer(palette="Set1") +
labs(title = "Age frequency plot", subtitle = "Hogfish", y = "Age", x = "Year") +
scale_size(range = c(1,10), breaks=c(1,2, 5, 10, 20, 40, 60, 80, 110, 150)) +
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
作为奖励，如果有人知道如何将对角线连接到从年龄、年份到年龄+1、年份+1 等的数据...那就太好了。
我滞后的线性代码很糟糕，我花了数周时间阅读文献和堆栈问题。如果你愿意，我可以向你展示更多我失败的尝试。
这是一次尝试
# linear model
l.fit <- lm(wi.age.count$Year ~ wi.age.count$n + lag(wi.age.count$Year, +1)); par(mfrow=c(1,2))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
plot(wi.age.count$Year ~ wi.age.count$n, pch = 2, type="b", xlab = 'Year', ylab = 'Age Frequency', xlim=range(age.hog$Year), ylim=range(c(0,age.hog$n)), main="Hogfish")
abline(l.fit, lwd=3, lty=3); legend (0, 700, paste("AIC =", AIC.l.fit), bty = 'n')
hist(residuals(l.fit), xlab='Residuals', main='Quality check')
summary(l.fit)
我什至不确定哪个是最合适的，滞后线性模型或 ARIMA 或 acf() 或完全不同的东西。其中一个问题是我有 3 个维度……年、年龄和年龄。任何帮助将不胜感激。
除所有科学文献外，我还尝试效仿的来源。
R - predicting simple dyn model with one lag term
Iteratively forecasting dyn models
Issue when attempting to run a Distributed Lag model in R using dynlm
Comparing linear regressions with a factor and lagged predictors, using R
数据应该是这样的...没有前几年。
你能试试这个吗？
Years = sort(unique(as.numeric(wi.age.count$Year)))
ivs = data.frame(Year = NA)
for (i in 1:(length(Years)-1)) {
ivs$dum = NA
names(ivs)[ncol(ivs)] = paste("n.", i, sep="")
}
i = 0
while (i < length(Years)) {
i = i + 1
tmp = data.frame(Year = Years[i])
j = i
while (j < length(Years)) {
j = j + 1
tmp$dum = 0
if (nrow(wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]) > 0) {
tmp$dum = wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]$n
}
names(tmp)[ncol(tmp)] = paste("n.", j - i, sep="")
}
k = 0
while (k < i - 1) {
k = k + 1
tmp$dum = NA
names(tmp)[ncol(tmp)] = paste("n.", j - i + k, sep="")
}
ivs = rbind(ivs, tmp)
}
ivs = ivs[-1, ]
ivs = ivs[-(nrow(ivs)), ]
ivs[is.na(ivs)] = 0
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
让我解释一下我在这里做什么。你说你想要一个线性模型来确定 YOY 的数量是否与明年 1 岁儿童的数量相关，另一年与 2 岁儿童的数量相关，依此类推。因此，我首先使用以下行创建自变量：
Years = sort(unique(as.numeric(wi.age.count$Year)))
ivs = data.frame(Year = NA)
for (i in 1:(length(Years)-1)) {
ivs$dum = NA
names(ivs)[ncol(ivs)] = paste("n.", i, sep="")
}
i = 0
while (i < length(Years)) {
i = i + 1
tmp = data.frame(Year = Years[i])
j = i
while (j < length(Years)) {
j = j + 1
tmp$dum = 0
if (nrow(wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]) > 0) {
tmp$dum = wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]$n
}
names(tmp)[ncol(tmp)] = paste("n.", j - i, sep="")
}
k = 0
while (k < i - 1) {
k = k + 1
tmp$dum = NA
names(tmp)[ncol(tmp)] = paste("n.", j - i + k, sep="")
}
ivs = rbind(ivs, tmp)
}
ivs = ivs[-1, ]
ivs[is.na(ivs)] = 0
> ivs
Year n.1 n.2 n.3 n.4 n.5 n.6 n.7 n.8 n.9 n.10 n.11
2 2007 4 37 15 12 13 1 2 0 1 0 0
3 2008 46 41 26 27 1 1 6 0 0 0 0
4 2009 19 4 13 4 2 3 0 0 0 0 0
5 2010 13 15 3 6 3 1 0 0 0 0 0
6 2011 87 30 23 28 9 4 1 0 0 0 0
7 2012 32 15 31 5 1 2 0 0 0 0 0
8 2013 24 13 1 1 1 0 0 0 0 0 0
9 2014 18 6 8 2 0 0 0 0 0 0 0
10 2015 4 16 7 0 0 0 0 0 0 0 0
11 2016 16 4 0 0 0 0 0 0 0 0 0
12 2017 12 0 0 0 0 0 0 0 0 0 0
这里，n.1是明年1岁的人数，n.2是2岁的人数另一年等等。
我还创建了一个只有因变量的数据框：
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
> dv
Year n
1 2008 166
2 2009 28
3 2010 34
4 2011 77
5 2012 170
6 2013 18
7 2014 3
8 2015 22
9 2016 43
10 2017 50
11 2018 151
我在估计模型。这是您的模型：
但是，让我们先创建公式。我不知道你想要多少延迟。因此，我创建了一个以 3 个滞后作为自变量的公式。如果需要，您可以通过更改 for 循环中的 2:4 来更改公式。
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
> formula
[1] "n ~ n.1 + n.2 + n.3"
这是您的模型：
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
summary(l.fit)
Call:
lm(formula = formula, data = merge(dv, ivs))
Residuals:
Min 1Q Median 3Q Max
-40.389 -29.713 -0.262 25.390 44.063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.0023 19.8713 0.654 0.5372
n.1 -0.5888 0.7834 -0.752 0.4807
n.2 1.1125 1.5050 0.739 0.4877
n.3 4.2888 1.5825 2.710 0.0351 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 38.25 on 6 degrees of freedom
Multiple R-squared: 0.727, Adjusted R-squared: 0.5905
F-statistic: 5.326 on 3 and 6 DF, p-value: 0.03967
> AIC.l.fit
[1] 106
我将根据您在 2020 年 7 月 27 日发表的评论添加另一个答案。该图没有数字，但它给出了我应该在 IVS 矩阵中包含的数字的一些想法。请尝试以下代码，看看是否有意义。
tmp = wi.age.count[order(wi.age.count$Age), ]
ivs = reshape(tmp[which(tmp$Age != 0), -4], direction = "wide", idvar = "Year", timevar = "Age")
ivs[is.na(ivs)] = 0
> ivs
Year n.1 n.2 n.3 n.4 n.5 n.6 n.7 n.8 n.9
13 2007 8 17 0 0 1 0 0 0 0
16 2008 4 19 1 0 0 0 0 0 0
20 2009 46 37 52 5 1 1 0 0 0
26 2010 19 41 15 16 0 0 0 0 1
32 2011 13 4 26 12 11 1 1 1 0
41 2012 87 15 13 27 13 17 1 0 0
49 2013 32 30 3 4 1 1 1 0 0
57 2014 24 15 23 6 2 1 2 2 0
66 2015 18 13 31 28 3 3 6 0 1
74 2016 4 6 1 5 9 1 0 0 1
82 2017 16 16 8 1 1 4 0 0 0
89 2018 12 4 7 2 1 2 1 0 0
这是你的ivs矩阵。这看起来正确吗？
其他都一样。这是您的 dv 矩阵：
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
> dv
Year n
1 2008 166
2 2009 28
3 2010 34
4 2011 77
5 2012 170
6 2013 18
7 2014 3
8 2015 22
9 2016 43
10 2017 50
11 2018 151
你的公式有三个滞后。
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
> formula
[1] "n ~ n.1 + n.2 + n.3"
结果如下：
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
summary(l.fit)
Call:
lm(formula = formula, data = merge(dv, ivs))
Residuals:
Min 1Q Median 3Q Max
-60.367 -38.028 8.698 23.763 96.257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 99.8976 36.1761 2.761 0.028 *
n.1 1.1059 0.8388 1.318 0.229
n.2 -1.7339 1.5773 -1.099 0.308
n.3 -1.6346 1.2932 -1.264 0.247
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 59.48 on 7 degrees of freedom
Multiple R-squared: 0.3731, Adjusted R-squared: 0.1044
F-statistic: 1.389 on 3 and 7 DF, p-value: 0.3233
> AIC.l.fit
[1] 126
我有这些数据，我正在尝试在 r 中进行滞后线性回归以确定 YOY 的数量是否与下一年的 1 岁儿童和后年的 2 岁儿童的数量显着相关...等等...
数据：
structure(list(Year = c("2008", "2009", "2010", "2011", "2012",
"2013", "2014", "2015", "2016", "2017", "2018", "2007", "2007",
"2007", "2007", "2008", "2008", "2008", "2009", "2009", "2009",
"2009", "2009", "2009", "2009", "2010", "2010", "2010", "2010",
"2010", "2011", "2011", "2011", "2011", "2011", "2011", "2011",
"2011", "2011", "2012", "2012", "2012", "2012", "2012", "2012",
"2012", "2012", "2013", "2013", "2013", "2013", "2013", "2013",
"2013", "2013", "2014", "2014", "2014", "2014", "2014", "2014",
"2014", "2014", "2014", "2015", "2015", "2015", "2015", "2015",
"2015", "2015", "2015", "2015", "2016", "2016", "2016", "2016",
"2016", "2016", "2016", "2017", "2017", "2017", "2017", "2017",
"2017", "2017", "2018", "2018", "2018", "2018", "2018", "2018",
"2018", "2018"), Age = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 3L, 6L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 2L, 3L, 4L, 5L, 10L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,
9L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L,
8L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 10L, 2L, 3L, 4L, 5L, 6L, 7L, 10L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), .Label = c("0", "1",
"2", "3", "4", "5", "6", "7", "8", "9"), class = "factor"), n = c(166,
28, 34, 77, 170, 18, 3, 22, 43, 50, 151, 1, 8, 17, 1, 4, 19,
1, 1, 46, 37, 52, 5, 1, 1, 19, 41, 15, 16, 1, 1, 13, 4, 26, 12,
11, 1, 1, 1, 1, 87, 15, 13, 27, 13, 17, 1, 1, 32, 30, 3, 4, 1,
1, 1, 1, 24, 15, 23, 6, 2, 1, 2, 2, 4, 18, 13, 31, 28, 3, 3,
6, 1, 4, 6, 1, 5, 9, 1, 1, 1, 16, 16, 8, 1, 1, 4, 1, 12, 4, 7,
2, 1, 2, 1), id = c("YOY", "YOY", "YOY", "YOY", "YOY", "YOY",
"YOY", "YOY", "YOY", "YOY", "YOY", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult", "Adult", "Adult", "Adult",
"Adult", "Adult", "Adult", "Adult")), row.names = c(NA, -95L), class = "data.frame")
我制作了很棒的情节，表明这里确实看起来确实有东西。不完美，但有某种关系。
# Frequencey density plot of ages over year
ggplot(wi.age.count, aes(x=Year, y=Age)) +
geom_point(aes(cex = n, color = id)) +
#scale_fill_brewer(palette="Set1") +
labs(title = "Age frequency plot", subtitle = "Hogfish", y = "Age", x = "Year") +
scale_size(range = c(1,10), breaks=c(1,2, 5, 10, 20, 40, 60, 80, 110, 150)) +
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
作为奖励，如果有人知道如何将对角线连接到从年龄、年份到年龄+1、年份+1 等的数据...那就太好了。
我滞后的线性代码很糟糕，我花了数周时间阅读文献和堆栈问题。如果你愿意，我可以向你展示更多我失败的尝试。
这是一次尝试
# linear model
l.fit <- lm(wi.age.count$Year ~ wi.age.count$n + lag(wi.age.count$Year, +1)); par(mfrow=c(1,2))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
plot(wi.age.count$Year ~ wi.age.count$n, pch = 2, type="b", xlab = 'Year', ylab = 'Age Frequency', xlim=range(age.hog$Year), ylim=range(c(0,age.hog$n)), main="Hogfish")
abline(l.fit, lwd=3, lty=3); legend (0, 700, paste("AIC =", AIC.l.fit), bty = 'n')
hist(residuals(l.fit), xlab='Residuals', main='Quality check')
summary(l.fit)
我什至不确定哪个是最合适的，滞后线性模型或 ARIMA 或 acf() 或完全不同的东西。其中一个问题是我有 3 个维度……年、年龄和年龄。任何帮助将不胜感激。
除所有科学文献外，我还尝试效仿的来源。
R - predicting simple dyn model with one lag term
Iteratively forecasting dyn models
Issue when attempting to run a Distributed Lag model in R using dynlm
Comparing linear regressions with a factor and lagged predictors, using R
数据应该是这样的...没有前几年。
你能试试这个吗？
Years = sort(unique(as.numeric(wi.age.count$Year)))
ivs = data.frame(Year = NA)
for (i in 1:(length(Years)-1)) {
ivs$dum = NA
names(ivs)[ncol(ivs)] = paste("n.", i, sep="")
}
i = 0
while (i < length(Years)) {
i = i + 1
tmp = data.frame(Year = Years[i])
j = i
while (j < length(Years)) {
j = j + 1
tmp$dum = 0
if (nrow(wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]) > 0) {
tmp$dum = wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]$n
}
names(tmp)[ncol(tmp)] = paste("n.", j - i, sep="")
}
k = 0
while (k < i - 1) {
k = k + 1
tmp$dum = NA
names(tmp)[ncol(tmp)] = paste("n.", j - i + k, sep="")
}
ivs = rbind(ivs, tmp)
}
ivs = ivs[-1, ]
ivs = ivs[-(nrow(ivs)), ]
ivs[is.na(ivs)] = 0
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
让我解释一下我在这里做什么。你说你想要一个线性模型来确定 YOY 的数量是否与明年 1 岁儿童的数量相关，另一年与 2 岁儿童的数量相关，依此类推。因此，我首先使用以下行创建自变量：
Years = sort(unique(as.numeric(wi.age.count$Year)))
ivs = data.frame(Year = NA)
for (i in 1:(length(Years)-1)) {
ivs$dum = NA
names(ivs)[ncol(ivs)] = paste("n.", i, sep="")
}
i = 0
while (i < length(Years)) {
i = i + 1
tmp = data.frame(Year = Years[i])
j = i
while (j < length(Years)) {
j = j + 1
tmp$dum = 0
if (nrow(wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]) > 0) {
tmp$dum = wi.age.count[which(wi.age.count$Year == Years[j] & wi.age.count$Age == Years[j] - Years[i]), ]$n
}
names(tmp)[ncol(tmp)] = paste("n.", j - i, sep="")
}
k = 0
while (k < i - 1) {
k = k + 1
tmp$dum = NA
names(tmp)[ncol(tmp)] = paste("n.", j - i + k, sep="")
}
ivs = rbind(ivs, tmp)
}
ivs = ivs[-1, ]
ivs[is.na(ivs)] = 0
> ivs
Year n.1 n.2 n.3 n.4 n.5 n.6 n.7 n.8 n.9 n.10 n.11
2 2007 4 37 15 12 13 1 2 0 1 0 0
3 2008 46 41 26 27 1 1 6 0 0 0 0
4 2009 19 4 13 4 2 3 0 0 0 0 0
5 2010 13 15 3 6 3 1 0 0 0 0 0
6 2011 87 30 23 28 9 4 1 0 0 0 0
7 2012 32 15 31 5 1 2 0 0 0 0 0
8 2013 24 13 1 1 1 0 0 0 0 0 0
9 2014 18 6 8 2 0 0 0 0 0 0 0
10 2015 4 16 7 0 0 0 0 0 0 0 0
11 2016 16 4 0 0 0 0 0 0 0 0 0
12 2017 12 0 0 0 0 0 0 0 0 0 0
这里，n.1是明年1岁的人数，n.2是2岁的人数另一年等等。
我还创建了一个只有因变量的数据框：
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
> dv
Year n
1 2008 166
2 2009 28
3 2010 34
4 2011 77
5 2012 170
6 2013 18
7 2014 3
8 2015 22
9 2016 43
10 2017 50
11 2018 151
我在估计模型。这是您的模型：
但是，让我们先创建公式。我不知道你想要多少延迟。因此，我创建了一个以 3 个滞后作为自变量的公式。如果需要，您可以通过更改 for 循环中的 2:4 来更改公式。
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
> formula
[1] "n ~ n.1 + n.2 + n.3"
这是您的模型：
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
summary(l.fit)
Call:
lm(formula = formula, data = merge(dv, ivs))
Residuals:
Min 1Q Median 3Q Max
-40.389 -29.713 -0.262 25.390 44.063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.0023 19.8713 0.654 0.5372
n.1 -0.5888 0.7834 -0.752 0.4807
n.2 1.1125 1.5050 0.739 0.4877
n.3 4.2888 1.5825 2.710 0.0351 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 38.25 on 6 degrees of freedom
Multiple R-squared: 0.727, Adjusted R-squared: 0.5905
F-statistic: 5.326 on 3 and 6 DF, p-value: 0.03967
> AIC.l.fit
[1] 106
我将根据您在 2020 年 7 月 27 日发表的评论添加另一个答案。该图没有数字，但它给出了我应该在 IVS 矩阵中包含的数字的一些想法。请尝试以下代码，看看是否有意义。
tmp = wi.age.count[order(wi.age.count$Age), ]
ivs = reshape(tmp[which(tmp$Age != 0), -4], direction = "wide", idvar = "Year", timevar = "Age")
ivs[is.na(ivs)] = 0
> ivs
Year n.1 n.2 n.3 n.4 n.5 n.6 n.7 n.8 n.9
13 2007 8 17 0 0 1 0 0 0 0
16 2008 4 19 1 0 0 0 0 0 0
20 2009 46 37 52 5 1 1 0 0 0
26 2010 19 41 15 16 0 0 0 0 1
32 2011 13 4 26 12 11 1 1 1 0
41 2012 87 15 13 27 13 17 1 0 0
49 2013 32 30 3 4 1 1 1 0 0
57 2014 24 15 23 6 2 1 2 2 0
66 2015 18 13 31 28 3 3 6 0 1
74 2016 4 6 1 5 9 1 0 0 1
82 2017 16 16 8 1 1 4 0 0 0
89 2018 12 4 7 2 1 2 1 0 0
这是你的ivs矩阵。这看起来正确吗？
其他都一样。这是您的 dv 矩阵：
dv = wi.age.count[which(wi.age.count$id == "YOY"), c(1, 3)]
> dv
Year n
1 2008 166
2 2009 28
3 2010 34
4 2011 77
5 2012 170
6 2013 18
7 2014 3
8 2015 22
9 2016 43
10 2017 50
11 2018 151
你的公式有三个滞后。
formula = ""
for (i in 2:4) formula = paste(formula, "+", names(ivs)[i])
formula = paste("n ~", substr(formula, 4, nchar(formula)))
> formula
[1] "n ~ n.1 + n.2 + n.3"
结果如下：
l.fit = lm(formula, merge(dv, ivs))
AIC.l.fit <- signif(AIC(l.fit), digits = 3)
summary(l.fit)
Call:
lm(formula = formula, data = merge(dv, ivs))
Residuals:
Min 1Q Median 3Q Max
-60.367 -38.028 8.698 23.763 96.257
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 99.8976 36.1761 2.761 0.028 *
n.1 1.1059 0.8388 1.318 0.229
n.2 -1.7339 1.5773 -1.099 0.308
n.3 -1.6346 1.2932 -1.264 0.247
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 59.48 on 7 degrees of freedom
Multiple R-squared: 0.3731, Adjusted R-squared: 0.1044
F-statistic: 1.389 on 3 and 7 DF, p-value: 0.3233
> AIC.l.fit
[1] 126