使用 rollmean 计算不包括 R 中第一个观察值的移动平均值
Using rollmean to calculate a moving average excluding the first observation in R
我目前正在将 Stata 中的一些时间序列数据命令翻译成 R。我正在使用 zoo 包来计算 R 中的移动平均值。这是我的数据的样子：
data <- cbind(c(1960:1970), c(95.5, 95.3, 95.3, 95.7, 95.7, 95.7,
95.1, 95.1, 95.1, 95, 95))
[,1] [,2]
[1,] 1960 95.5
[2,] 1961 95.3
[3,] 1962 95.3
[4,] 1963 95.7
[5,] 1964 95.7
[6,] 1965 95.7
[7,] 1966 95.1
[8,] 1967 95.1
[9,] 1968 95.1
[10,] 1969 95.0
[11,] 1970 95.0
我会把它变成 data.frame:
data <- as.data.frame(data)
现在，我可以使用 rollmean 函数用我的数据计算 turnout 的移动平均值：
data$turnout <- rollmean(data[,2], 1, fill = NA)
这就是我得到的：
V1 V2 turnout
1 1960 95.5 95.5
2 1961 95.3 95.3
3 1962 95.3 95.3
4 1963 95.7 95.7
5 1964 95.7 95.7
6 1965 95.7 95.7
7 1966 95.1 95.1
8 1967 95.1 95.1
9 1968 95.1 95.1
10 1969 95.0 95.0
11 1970 95.0 95.0
一切都很好，但我的问题是我希望我的列 turnout（移动平均线）从 1961 年而不是 1960 年开始。此代码不排除第一个观察值，即我正在努力。
作为参考，等效的 Stata 命令为：
tssmooth ma m1turnout = turnout, window (1 0)
我已经尝试过使用 align = "right" 函数，但这似乎并不能解决问题。有什么想法吗？
提前致谢！
编辑——为了澄清，我在不同的长度上做这件事。在 Stata 中，完整的代码是这样的，其中 since 是一个变量，它描述了自干预以来的年数。
foreach y of numlist 1(1)10{
tssmooth ma m`y'turnout = turnout, window (`y' 0)
}
gen dvturnout=.
foreach y of numlist 2(1)9{
replace dvturnout = l1.turnout if since==1
replace dvturnout = m`y'turnout if since==`y' & m`y'turnout!=.
replace dvturnout = m10turnout if (since==10 & m10turnout!=.) | (since==. & redist!=. & m10turnout!=.)
}
foreach y of numlist 1(1)10{
drop m`y'turnout
}
我的最终目标是这个 dvturnout 变量。
当我尝试我认为对应于 Stata 中代码的第一部分的内容时，即：
foreach y of numlist 1(1)10{
tssmooth ma m`y'turnout = turnout, window (`y' 0)
}
在 R 中，我这样做（其中 [,35] 是我开始向其添加变量的列）：
for (j in 1:10) {
data_countries[[i]][,35+j] <- rollmean(data_countries[[i]][,13], j, fill = NA, align = "right")
}
}
它为我吐出这个：
year since V36 V37 V38 V39 V40 V41 V42 V43 V44 V45
1 1960 NA 95.5 NA NA NA NA NA NA NA NA NA
2 1961 NA 95.3 95.40 NA NA NA NA NA NA NA NA
3 1962 NA 95.3 95.30 95.36667 NA NA NA NA NA NA NA
4 1963 NA 95.7 95.50 95.43333 95.450 NA NA NA NA NA NA
5 1964 NA 95.7 95.70 95.56667 95.500 95.50 NA NA NA NA NA
6 1965 NA 95.7 95.70 95.70000 95.600 95.54 95.53333 NA NA NA NA
7 1966 NA 95.1 95.40 95.50000 95.550 95.50 95.46667 95.47143 NA NA NA
8 1967 NA 95.1 95.10 95.30000 95.400 95.46 95.43333 95.41428 95.4250 NA NA
9 1968 NA 95.1 95.10 95.10000 95.250 95.34 95.40000 95.38571 95.3750 95.38889 NA
10 1969 NA 95.0 95.05 95.06667 95.075 95.20 95.28333 95.34286 95.3375 95.33333 95.35
11 1970 NA 95.0 95.00 95.03333 95.050 95.06 95.16667 95.24286 95.3000 95.30000 95.30
这些数字都很好，但比我希望的要低 "shifted"。这是相同操作在 Stata 中给我的结果：
year dvturnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout m10turnout
1960
1961 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5
1962 95.3 95.4 95.4 95.4 95.4 95.4 95.4 95.4 95.4 95.4
1963 95.3 95.3 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667
1964 95.7 95.5 95.43333 95.45 95.45 95.45 95.45 95.45 95.45 95.45
1965 95.7 95.7 95.56667 95.5 95.5 95.5 95.5 95.5 95.5 95.5
1966 95.7 95.7 95.7 95.6 95.54 95.53333 95.53333 95.53333 95.53333 95.53333
1967 95.1 95.39999 95.5 95.55 95.5 95.46667 95.47143 95.47143 95.47143 95.47143
1968 95.1 95.1 95.3 95.39999 95.46 95.43333 95.41428 95.425 95.425 95.425
1969 95.1 95.1 95.1 95.25 95.34 95.39999 95.38571 95.375 95.38889 95.38889
1970 95 95.05 95.06667 95.075 95.2 95.28333 95.34286 95.3375 95.33334 95.35
您需要的是一个不包括当前观测值的移动平均函数。谢天谢地，w_i_l_lwrote a function exactly like that。是什么让事情变得复杂：您论文的作者用上一列的结果填充了没有足够数据（例如，k = 4，但只有 3 个数据点）的移动平均线。我真的不建议这样做，因为如果没有非常明确地指出，这会（而且通常会）导致严重的混乱。
代码
# w_i_l_l's moving average function
mav <- function(x,n){filter(x,rep(1/n,n), sides=1)}
mavback <- function(x,n){
a<-mav(x,1)
b<-mav(x,(n+1))
c<-(1/n)*((n+1)*b - a)
return(c)
}
# Create 10 columns with moving averages of k = 1:10
result <- NULL
for(i in 1:10){
result <- cbind(result,mavback(test[,2], i))
}
# Give propers names to columns
colnames(result) <- paste0("m", 1:ncol(result)-1,"turnout")
# Combine result with base data
result <- cbind(test,data.frame(result))
# WONKY STATISTICS: If there is a NA (= not enough data for a
# moving average) fill it up with previous column's result
for(i in 4:ncol(result)){
# Nested loop starts from first row
for(j in 2:nrow(result)){
# Check for NA
if(is.na(result[j,i])){
result[j,i] <- result[j,i-1]
}
}
}
结果
> result
year turnout m0turnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout
1 1960 95.5 NA NA NA NA NA NA NA NA NA NA
2 1961 95.3 95.5 95.50 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000
3 1962 95.3 95.3 95.40 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000
4 1963 95.7 95.3 95.30 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667
5 1964 95.7 95.7 95.50 95.43333 95.45000 95.45000 95.45000 95.45000 95.45000 95.45000 95.45000
6 1965 95.7 95.7 95.70 95.56667 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000
7 1966 95.1 95.7 95.70 95.70000 95.60000 95.54000 95.53333 95.53333 95.53333 95.53333 95.53333
8 1967 95.1 95.1 95.40 95.50000 95.55000 95.50000 95.46667 95.47143 95.47143 95.47143 95.47143
9 1968 95.1 95.1 95.10 95.30000 95.40000 95.46000 95.43333 95.41429 95.42500 95.42500 95.42500
10 1969 95.0 95.1 95.10 95.10000 95.25000 95.34000 95.40000 95.38571 95.37500 95.38889 95.38889
11 1970 95.0 95.0 95.05 95.06667 95.07500 95.20000 95.28333 95.34286 95.33750 95.33333 95.35000
没有 "filling up"
的结果
> result
year turnout m0turnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout
1 1960 95.5 NA NA NA NA NA NA NA NA NA NA
2 1961 95.3 95.5 NA NA NA NA NA NA NA NA NA
3 1962 95.3 95.3 95.40 NA NA NA NA NA NA NA NA
4 1963 95.7 95.3 95.30 95.36667 NA NA NA NA NA NA NA
5 1964 95.7 95.7 95.50 95.43333 95.450 NA NA NA NA NA NA
6 1965 95.7 95.7 95.70 95.56667 95.500 95.50 NA NA NA NA NA
7 1966 95.1 95.7 95.70 95.70000 95.600 95.54 95.53333 NA NA NA NA
8 1967 95.1 95.1 95.40 95.50000 95.550 95.50 95.46667 95.47143 NA NA NA
9 1968 95.1 95.1 95.10 95.30000 95.400 95.46 95.43333 95.41429 95.4250 NA NA
10 1969 95.0 95.1 95.10 95.10000 95.250 95.34 95.40000 95.38571 95.3750 95.38889 NA
11 1970 95.0 95.0 95.05 95.06667 95.075 95.20 95.28333 95.34286 95.3375 95.33333 95.35
数据
test <- data.frame(cbind(year = c(1960:1970),
turnout = c(95.5, 95.3, 95.3, 95.7, 95.7,
95.7, 95.1, 95.1, 95.1, 95, 95)))
也许您正在寻找这样的东西：
library(zoo)
library(forecast)
data <- cbind(c(1960:1970), c(95.5, 95.3, 95.3, 95.7, 95.7, 95.7, 95.1, 95.1, 95.1, 95, 95))
x1 <- ts(data = data[, 2], start = 1960, end = 1970, frequency = 1)
x2 <- cbind(x1, turnout = zoo::rollmeanr(x1, k = 2))
打印时间序列对象：
x2
Time Series:
Start = 1960
End = 1970
Frequency = 1
x1 turnout
1960 95.5 NA
1961 95.3 95.40
1962 95.3 95.30
1963 95.7 95.50
1964 95.7 95.70
1965 95.7 95.70
1966 95.1 95.40
1967 95.1 95.10
1968 95.1 95.10
1969 95.0 95.05
1970 95.0 95.00
剧情：
forecast::autoplot(x2)
我发现最简单的方法是使用 lag 函数。
data$turnout <- lag(rollmean(data[,2], 1, fill = NA),1)
我目前正在将 Stata 中的一些时间序列数据命令翻译成 R。我正在使用 zoo 包来计算 R 中的移动平均值。这是我的数据的样子：
data <- cbind(c(1960:1970), c(95.5, 95.3, 95.3, 95.7, 95.7, 95.7,
95.1, 95.1, 95.1, 95, 95))
[,1] [,2]
[1,] 1960 95.5
[2,] 1961 95.3
[3,] 1962 95.3
[4,] 1963 95.7
[5,] 1964 95.7
[6,] 1965 95.7
[7,] 1966 95.1
[8,] 1967 95.1
[9,] 1968 95.1
[10,] 1969 95.0
[11,] 1970 95.0
我会把它变成 data.frame:
data <- as.data.frame(data)
现在，我可以使用 rollmean 函数用我的数据计算 turnout 的移动平均值：
data$turnout <- rollmean(data[,2], 1, fill = NA)
这就是我得到的：
V1 V2 turnout
1 1960 95.5 95.5
2 1961 95.3 95.3
3 1962 95.3 95.3
4 1963 95.7 95.7
5 1964 95.7 95.7
6 1965 95.7 95.7
7 1966 95.1 95.1
8 1967 95.1 95.1
9 1968 95.1 95.1
10 1969 95.0 95.0
11 1970 95.0 95.0
一切都很好，但我的问题是我希望我的列 turnout（移动平均线）从 1961 年而不是 1960 年开始。此代码不排除第一个观察值，即我正在努力。
作为参考，等效的 Stata 命令为：
tssmooth ma m1turnout = turnout, window (1 0)
我已经尝试过使用 align = "right" 函数，但这似乎并不能解决问题。有什么想法吗？
提前致谢！
编辑——为了澄清，我在不同的长度上做这件事。在 Stata 中，完整的代码是这样的，其中 since 是一个变量，它描述了自干预以来的年数。
foreach y of numlist 1(1)10{
tssmooth ma m`y'turnout = turnout, window (`y' 0)
}
gen dvturnout=.
foreach y of numlist 2(1)9{
replace dvturnout = l1.turnout if since==1
replace dvturnout = m`y'turnout if since==`y' & m`y'turnout!=.
replace dvturnout = m10turnout if (since==10 & m10turnout!=.) | (since==. & redist!=. & m10turnout!=.)
}
foreach y of numlist 1(1)10{
drop m`y'turnout
}
我的最终目标是这个 dvturnout 变量。
当我尝试我认为对应于 Stata 中代码的第一部分的内容时，即：
foreach y of numlist 1(1)10{
tssmooth ma m`y'turnout = turnout, window (`y' 0)
}
在 R 中，我这样做（其中 [,35] 是我开始向其添加变量的列）：
for (j in 1:10) {
data_countries[[i]][,35+j] <- rollmean(data_countries[[i]][,13], j, fill = NA, align = "right")
}
}
它为我吐出这个：
year since V36 V37 V38 V39 V40 V41 V42 V43 V44 V45
1 1960 NA 95.5 NA NA NA NA NA NA NA NA NA
2 1961 NA 95.3 95.40 NA NA NA NA NA NA NA NA
3 1962 NA 95.3 95.30 95.36667 NA NA NA NA NA NA NA
4 1963 NA 95.7 95.50 95.43333 95.450 NA NA NA NA NA NA
5 1964 NA 95.7 95.70 95.56667 95.500 95.50 NA NA NA NA NA
6 1965 NA 95.7 95.70 95.70000 95.600 95.54 95.53333 NA NA NA NA
7 1966 NA 95.1 95.40 95.50000 95.550 95.50 95.46667 95.47143 NA NA NA
8 1967 NA 95.1 95.10 95.30000 95.400 95.46 95.43333 95.41428 95.4250 NA NA
9 1968 NA 95.1 95.10 95.10000 95.250 95.34 95.40000 95.38571 95.3750 95.38889 NA
10 1969 NA 95.0 95.05 95.06667 95.075 95.20 95.28333 95.34286 95.3375 95.33333 95.35
11 1970 NA 95.0 95.00 95.03333 95.050 95.06 95.16667 95.24286 95.3000 95.30000 95.30
这些数字都很好，但比我希望的要低 "shifted"。这是相同操作在 Stata 中给我的结果：
year dvturnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout m10turnout
1960
1961 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5 95.5
1962 95.3 95.4 95.4 95.4 95.4 95.4 95.4 95.4 95.4 95.4
1963 95.3 95.3 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667
1964 95.7 95.5 95.43333 95.45 95.45 95.45 95.45 95.45 95.45 95.45
1965 95.7 95.7 95.56667 95.5 95.5 95.5 95.5 95.5 95.5 95.5
1966 95.7 95.7 95.7 95.6 95.54 95.53333 95.53333 95.53333 95.53333 95.53333
1967 95.1 95.39999 95.5 95.55 95.5 95.46667 95.47143 95.47143 95.47143 95.47143
1968 95.1 95.1 95.3 95.39999 95.46 95.43333 95.41428 95.425 95.425 95.425
1969 95.1 95.1 95.1 95.25 95.34 95.39999 95.38571 95.375 95.38889 95.38889
1970 95 95.05 95.06667 95.075 95.2 95.28333 95.34286 95.3375 95.33334 95.35
您需要的是一个不包括当前观测值的移动平均函数。谢天谢地，w_i_l_lwrote a function exactly like that。是什么让事情变得复杂：您论文的作者用上一列的结果填充了没有足够数据（例如，k = 4，但只有 3 个数据点）的移动平均线。我真的不建议这样做，因为如果没有非常明确地指出，这会（而且通常会）导致严重的混乱。
代码
# w_i_l_l's moving average function
mav <- function(x,n){filter(x,rep(1/n,n), sides=1)}
mavback <- function(x,n){
a<-mav(x,1)
b<-mav(x,(n+1))
c<-(1/n)*((n+1)*b - a)
return(c)
}
# Create 10 columns with moving averages of k = 1:10
result <- NULL
for(i in 1:10){
result <- cbind(result,mavback(test[,2], i))
}
# Give propers names to columns
colnames(result) <- paste0("m", 1:ncol(result)-1,"turnout")
# Combine result with base data
result <- cbind(test,data.frame(result))
# WONKY STATISTICS: If there is a NA (= not enough data for a
# moving average) fill it up with previous column's result
for(i in 4:ncol(result)){
# Nested loop starts from first row
for(j in 2:nrow(result)){
# Check for NA
if(is.na(result[j,i])){
result[j,i] <- result[j,i-1]
}
}
}
结果
> result
year turnout m0turnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout
1 1960 95.5 NA NA NA NA NA NA NA NA NA NA
2 1961 95.3 95.5 95.50 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000
3 1962 95.3 95.3 95.40 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000 95.40000
4 1963 95.7 95.3 95.30 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667 95.36667
5 1964 95.7 95.7 95.50 95.43333 95.45000 95.45000 95.45000 95.45000 95.45000 95.45000 95.45000
6 1965 95.7 95.7 95.70 95.56667 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000 95.50000
7 1966 95.1 95.7 95.70 95.70000 95.60000 95.54000 95.53333 95.53333 95.53333 95.53333 95.53333
8 1967 95.1 95.1 95.40 95.50000 95.55000 95.50000 95.46667 95.47143 95.47143 95.47143 95.47143
9 1968 95.1 95.1 95.10 95.30000 95.40000 95.46000 95.43333 95.41429 95.42500 95.42500 95.42500
10 1969 95.0 95.1 95.10 95.10000 95.25000 95.34000 95.40000 95.38571 95.37500 95.38889 95.38889
11 1970 95.0 95.0 95.05 95.06667 95.07500 95.20000 95.28333 95.34286 95.33750 95.33333 95.35000
没有 "filling up"
的结果> result
year turnout m0turnout m1turnout m2turnout m3turnout m4turnout m5turnout m6turnout m7turnout m8turnout m9turnout
1 1960 95.5 NA NA NA NA NA NA NA NA NA NA
2 1961 95.3 95.5 NA NA NA NA NA NA NA NA NA
3 1962 95.3 95.3 95.40 NA NA NA NA NA NA NA NA
4 1963 95.7 95.3 95.30 95.36667 NA NA NA NA NA NA NA
5 1964 95.7 95.7 95.50 95.43333 95.450 NA NA NA NA NA NA
6 1965 95.7 95.7 95.70 95.56667 95.500 95.50 NA NA NA NA NA
7 1966 95.1 95.7 95.70 95.70000 95.600 95.54 95.53333 NA NA NA NA
8 1967 95.1 95.1 95.40 95.50000 95.550 95.50 95.46667 95.47143 NA NA NA
9 1968 95.1 95.1 95.10 95.30000 95.400 95.46 95.43333 95.41429 95.4250 NA NA
10 1969 95.0 95.1 95.10 95.10000 95.250 95.34 95.40000 95.38571 95.3750 95.38889 NA
11 1970 95.0 95.0 95.05 95.06667 95.075 95.20 95.28333 95.34286 95.3375 95.33333 95.35
数据
test <- data.frame(cbind(year = c(1960:1970),
turnout = c(95.5, 95.3, 95.3, 95.7, 95.7,
95.7, 95.1, 95.1, 95.1, 95, 95)))
也许您正在寻找这样的东西：
library(zoo)
library(forecast)
data <- cbind(c(1960:1970), c(95.5, 95.3, 95.3, 95.7, 95.7, 95.7, 95.1, 95.1, 95.1, 95, 95))
x1 <- ts(data = data[, 2], start = 1960, end = 1970, frequency = 1)
x2 <- cbind(x1, turnout = zoo::rollmeanr(x1, k = 2))
打印时间序列对象：
x2
Time Series:
Start = 1960
End = 1970
Frequency = 1
x1 turnout
1960 95.5 NA
1961 95.3 95.40
1962 95.3 95.30
1963 95.7 95.50
1964 95.7 95.70
1965 95.7 95.70
1966 95.1 95.40
1967 95.1 95.10
1968 95.1 95.10
1969 95.0 95.05
1970 95.0 95.00
剧情：
forecast::autoplot(x2)
我发现最简单的方法是使用 lag 函数。
data$turnout <- lag(rollmean(data[,2], 1, fill = NA),1)