Использование rollmean для вычисления скользящего среднего, исключая первое наблюдение в R - PullRequest
Новая
       53

Использование rollmean для вычисления скользящего среднего, исключая первое наблюдение в R

0 голосов
/ 05 октября 2018

В настоящее время я работаю над переводом некоторых команд для данных временных рядов в 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

Все эти числа в порядке, но они "сдвинуты" вниз от того места, где я хочу их видетьбыть.Вот что дает мне та же операция в 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

Ответы [ 3 ]

0 голосов
/ 06 октября 2018

Может быть, вы ищете что-то вроде этого: 

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)

enter image description here

0 голосов
/ 08 октября 2018

Я нашел самый простой способ сделать это с помощью функции lag. 

data$turnout <- lag(rollmean(data[,2], 1, fill = NA),1)
0 голосов
/ 05 октября 2018

Вам нужна функция скользящего среднего , которая не включает текущее наблюдение .К счастью, w_i_l_l написал функцию, подобную этой .Что усложняло ситуацию: автор вашей статьи заполнил скользящее среднее, у которого недостаточно данных (например, 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

Результат без "заполнения" 

> 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)))
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...