psysky
Reputation: 3195
Calculating Root Mean Square Percentage Error in R
I performed the forecast, using this dataset here the trainsample(part of)
train=structure(list(Store = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L,
10L, 11L, 12L, 13L, 14L, 15L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,
9L, 10L, 11L, 12L, 13L, 14L, 15L), DayOfWeek = c(5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), Date = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("30.01.2015",
"31.01.2015"), class = "factor"), Sales = c(5577L, 5919L, 6911L,
13307L, 5640L, 6555L, 11430L, 6401L, 8072L, 6350L, 10031L, 9156L,
7004L, 6491L, 8898L, 5577L, 5919L, 6911L, 13307L, 5640L, 6555L,
11430L, 6401L, 8072L, 6350L, 10031L, 9156L, 7004L, 6491L, 8898L
), Customers = c(616L, 624L, 678L, 1632L, 617L, 692L, 1077L,
747L, 643L, 602L, 1263L, 988L, 453L, 692L, 828L, 616L, 624L,
678L, 1632L, 617L, 692L, 1077L, 747L, 643L, 602L, 1263L, 988L,
453L, 692L, 828L), Open = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L), Promo = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L), StateHoliday = c(0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), SchoolHoliday = c(0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("Store",
"DayOfWeek", "Date", "Sales", "Customers", "Open", "Promo", "StateHoliday",
"SchoolHoliday"), class = "data.frame", row.names = c(NA, -30L
))
and here the test dataset
df2=structure(list(Id = 1:10, Store = 1:10, DayOfWeek = c(5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), Date = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "31.07.2015", class = "factor"),
Customers = c(555L, 625L, 821L, 1498L, 559L, 589L, 1414L,
833L, 687L, 681L), Open = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L), Promo = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), StateHoliday = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L
), SchoolHoliday = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L)), .Names = c("Id", "Store", "DayOfWeek", "Date", "Customers",
"Open", "Promo", "StateHoliday", "SchoolHoliday"), class = "data.frame", row.names = c(NA,
-10L))
Test sample has data only for 31/07/15, i have to predict the sales. So forecast was performed for each store separately.
Forecast performed via regression analysis
Forecast <- test
for(store in store1$Store) {
coeff <- lm(data = train[train$Store == store, ],
Sales ~ DateDiff)$coefficients
store1[store1$Store == store, 'reg_intercept']<- coeff[1]
store1[store1$Store == store, 'reg_slope'] <- coeff[2]
train[train$Store == store, 'LinearRegressionForecast'] <-
coeff[1] + coeff[2] * train[train$Store == store, 'DateDiff']
Forecast[Forecast$Store == store, 'LinearRegressionForecast'] <-
coeff[1] + coeff[2] * Forecast[Forecast$Store == store, 'DateDiff']
}
#set predictors
predictors <- c('Store', 'DayOfWeek', 'Promo')
modelForecast <- train %>%
group_by_(.dots=predictors) %>%
summarize(salesMinusForecast=mean(Sales - LinearRegressionForecast)) %>%
ungroup()
Forecast <- Forecast %>%
left_join(modelForecast, by=predictors) %>%
mutate(Sales=salesMinusForecast + LinearRegressionForecast) %>%
select(Id, Store, DayOfWeek, Date, Sales, Open, Promo, StateHoliday,
SchoolHoliday, WeekOfYear, DateDiff, LinearRegressionForecast)
#View(Forecast)
Forecast$Forecast <- 1
train$Forecast <- 0
x=train[, c(1:4, 6:13)]
y=Forecast[, 2:13]
testtrain=rbind(x,y)
testtrain[testtrain$Forecast == 1,
'Type'] <- "Forecast"
testtrain[testtrain$Forecast == 0,
'Type'] <- "Observed"
also i provide part of data for store dataset
store1=structure(list(Store = 1:13, StoreType = structure(c(2L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L), .Label = c("a",
"c", "d"), class = "factor"), Assortment = structure(c(1L, 1L,
1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L), .Label = c("a",
"c"), class = "factor")), .Names = c("Store", "StoreType", "Assortment"
), class = "data.frame", row.names = c(NA, -13L))
When forecast is performed, i just write data to csv file
Forecast <- data.frame(Id=Forecast$Id, Sales=Forecast$Sales)
So the main problem is that for each store i must calculate Root Mean Square Percentage Error (RMSPE).
The formula is
where y_i denotes the sales of a single store on a single day and ŷ_i denotes the corresponding prediction.
So output is simple number of store , sale, predicted sale and RMSPE, i.e. 4 column.
for example here the result for 173 and 174 store(without RMSPE)
testtrain=structure(list(Store = c(173L, 173L, 173L, 173L, 173L, 173L,
173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 174L, 174L, 174L,
174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L,
173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L,
173L, 173L, 173L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L,
174L, 174L, 174L, 174L, 174L, 174L, 173L, 173L, 173L, 173L, 173L,
173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 173L, 174L, 174L,
174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L, 174L,
174L), DayOfWeek = c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), Date = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L), .Label = c("15.07.2015", "16.07.2015", "17.07.2015"
), class = "factor"), Sales = structure(c(27L, 12L, 16L, 18L,
9L, 4L, 26L, 23L, 10L, 19L, 7L, 20L, 25L, 5L, 17L, 2L, 11L, 8L,
3L, 22L, 15L, 14L, 28L, 6L, 1L, 24L, 13L, 21L, 27L, 12L, 16L,
18L, 9L, 4L, 26L, 23L, 10L, 19L, 7L, 20L, 25L, 5L, 17L, 2L, 11L,
8L, 3L, 22L, 15L, 14L, 28L, 6L, 1L, 24L, 13L, 21L, 27L, 12L,
16L, 18L, 9L, 4L, 26L, 23L, 10L, 19L, 7L, 20L, 25L, 5L, 17L,
2L, 11L, 8L, 3L, 22L, 15L, 14L, 28L, 6L, 1L, 24L, 13L, 21L), .Label = c("10318.344",
"10725.268", "10765.647", "13546.236", "3418.328", "3939.406",
"4089.442", "4377.643", "5196.012", "5487.437", "5778.296", "6200.403",
"6216.929", "6331.589", "6404.693", "6472.833", "6693.678", "6751.922",
"6770.161", "7510.433", "7736.447", "7743.879", "8107.569", "8119.046",
"9087.104", "9326.839", "9718.452", "9855.327"), class = "factor"),
Promo = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), LinearRegressionForecast = structure(c(22L,
9L, 15L, 14L, 8L, 1L, 26L, 20L, 6L, 17L, 3L, 18L, 25L, 2L,
11L, 27L, 7L, 5L, 24L, 13L, 10L, 12L, 23L, 4L, 28L, 21L,
16L, 19L, 22L, 9L, 15L, 14L, 8L, 1L, 26L, 20L, 6L, 17L, 3L,
18L, 25L, 2L, 11L, 27L, 7L, 5L, 24L, 13L, 10L, 12L, 23L,
4L, 28L, 21L, 16L, 19L, 22L, 9L, 15L, 14L, 8L, 1L, 26L, 20L,
6L, 17L, 3L, 18L, 25L, 2L, 11L, 27L, 7L, 5L, 24L, 13L, 10L,
12L, 23L, 4L, 28L, 21L, 16L, 19L), .Label = c("10672.724",
"2286.724", "2940.339", "3038.273", "3265.624", "3387.729",
"3475.001", "3568.385", "4527.949", "5042.683", "5131.816",
"5196.835", "5204.855", "5239.113", "5572.545", "5605.564",
"5656.971", "6216.276", "6510.814", "6749.251", "6901.256",
"7248.194", "7310.538", "7549.539", "7585.489", "7842.506",
"8371.118", "8487.823"), class = "factor"), Type = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Forecast", "obser"
), class = "factor")), .Names = c("Store", "DayOfWeek", "Date",
"Sales", "Promo", "LinearRegressionForecast", "Type"), class = "data.frame", row.names = c(NA,
-84L))
LinearRegressionForecast column we don't touch. There is column Type. forecast is that we predicted, and obser is initial value.
How can i calculate this metric(RMSPE) for each store.
RMSPE must be calculted for each store and for each day
Upvotes: 0
Views: 3184
Answers (1)
Rui Barradas
Reputation: 76460
As posted, your data has two problems, column Date is not of class Date and columns Sales and LinearRegressionForecast are of class factor. So I will start by coercing those columns to appropriate classes.
testtrain$Date <- as.Date(testtrain$Date, "%d.%m.%Y")
testtrain$Sales <- as.numeric(levels(testtrain$Sales))[testtrain$Sales]
testtrain$LinearRegressionForecast <- as.numeric(levels(testtrain$LinearRegressionForecast))[testtrain$LinearRegressionForecast]
To compute the statistic of interest, define a function rmspe, split the dataset by Store and Date, then call the function passing each of these sub dataframes.
rmspe <- function(DF){
y <- DF[["Sales"]]
y_obs <- y[DF[["Type"]] == "obser"]
y_est <- y[DF[["Type"]] == "Forecast"]
sqrt(mean(((y_obs - y_est)/y_obs)^2))
}
sp <- split(testtrain, list(testtrain$Store, testtrain$Date))
r <- sapply(sp, rmspe)
result <- data.frame(Store = sub("\\..*$", "", names(r)),
Date = sub("^\\d+\\.([-[:digit:]]+).*$", "\\1", names(r)),
RMSPE = r)
row.names(result) <- NULL
rm(sp, r)
result
# Store Date RMSPE
#1 173 2015-07-15 0.7373282
#2 174 2015-07-15 0.3287756
#3 173 2015-07-16 0.7373282
#4 174 2015-07-16 0.3287756
#5 173 2015-07-17 0.7373282
#6 174 2015-07-17 0.3287756
Upvotes: 1
Related Questions
- calculating adjusted R-Squared in R
- How to obtain RMSE out of lm result?
- Root mean square error in R - mixed effect model
- RMSE (root mean square deviation) calculation in R
- How to go from root mean squared error to percentage accuracy in R?
- How do I find the root mean squared error in R?
- calculating mean squared prediction error based on simulated data in R
- R: Calculate predicted mean in a lm() function
- Equation for Standard Error
- Calculating the mean, standard error and % in R for a data frame