# Inladen van de benodigde libraries
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.2
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.5.1 ✔ purrr 0.3.4
## ✔ tibble 3.2.1 ✔ dplyr 1.1.4
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(tsibble) # Voor tijdreeksanalyse
## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr
##
## Attaching package: 'tsibble'
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(feasts) # Voor tijdreekscomponentanalyse
## Loading required package: fabletools
library(lubridate) # Voor datummanipulatie
##
## Attaching package: 'lubridate'
##
## The following object is masked from 'package:tsibble':
##
## interval
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(zoo) # Voor na.approx()
## Warning: package 'zoo' was built under R version 4.2.3
##
## Attaching package: 'zoo'
##
## The following object is masked from 'package:tsibble':
##
## index
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(dplyr)
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.2.3
library(fable)
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.2.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.2.3
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
# path_filedata <- "/Users/hannahbrands/Downloads/datasets/filedata.csv"
# path_weer <- "/Users/hannahbrands/Downloads/datasets/weer.csv"
path_filedata <- "C:/Users/oleja/Downloads/datasets/filedata.csv"
path_weer <- "C:/Users/oleja/Downloads/weer.csv"
# Inlezen en voorbereiden van de filedataset
# We beginnen met het inlezen van de dataset en het opschonen van de data.
# Ongeldige of negatieve waarden worden gefilterd om fouten in de analyse te voorkomen.
filedata <- read.csv(path_filedata, sep = ";") %>%
# Filter: Verwijder niet-numerieke en negatieve waarden
filter(
!grepl("[^0-9.]", duurMinuten), # Controleer op niet-numerieke tekens in duurMinuten
!grepl("[^0-9.]", zwaarteKmMin) # Controleer op niet-numerieke tekens in zwaarteKmMin
) %>%
# Formaatconversie: Zet kolommen om naar juiste typen
mutate(
BeginDatum = as_date(BeginDatum), # Zet BeginDatum om naar datumformaat
duurMinuten = abs(as.numeric(duurMinuten)), # Zet duurMinuten om naar absolute numerieke waarden
zwaarteKmMin = abs(as.numeric(zwaarteKmMin)) # Zet zwaarteKmMin om naar absolute numerieke waarden
) %>%
# Filter: Verwijder irrelevante of foutieve waarnemingen
filter(
zwaarteKmMin > 0, # Filezwaarte moet groter dan 0 zijn
duurMinuten > 0, # Fileduur moet groter dan 0 zijn
zwaarteKmMin / 2 >= duurMinuten # Filelengte moet minimaal 2 keer de duur zijn
)
# Berekening van dagelijkse filezwaarte
# De totale filezwaarte per dag wordt berekend door de waarden per datum te sommeren.
dag_zwaarte <- filedata %>%
group_by(BeginDatum) %>%
summarize(filezwaarte = sum(zwaarteKmMin), .groups = "drop") %>% # Bereken de som per dag
as_tsibble(index = BeginDatum) %>% # Zet om naar een tsibble (tijdreeks)
fill_gaps(filezwaarte = NA) %>% # Vul ontbrekende dagen in met NA
mutate(filezwaarte = na.approx(filezwaarte)) # Vervang NA's met lineaire interpolatie
# Visualisatie van de dagelijkse filezwaarte
# Een lijnplot wordt gemaakt om de dagelijkse trends in filezwaarte te visualiseren.
dag_zwaarte %>%
autoplot(filezwaarte) +
labs(
title = "Dagelijkse Filezwaarte",
x = "Datum",
y = "Filezwaarte (km∙min)"
)
# STL-decompositie uitvoeren
# De tijdreeks wordt gedecomposeerd in trend, seizoen en restcomponenten.
# Er is gekozen voor een seizoensperiode van 365 dagen om jaarlijkse patronen te analyseren.
decomp <- dag_zwaarte %>%
model(STL(filezwaarte ~ season(period = 365), robust = TRUE)) %>%
components()
# Identificatie van uitbijters
# Uitbijters worden geïdentificeerd in de restcomponent met behulp van een boxplot-methode (IQR-regel).
filezwaarte_outliers <- decomp %>%
filter(
remainder < quantile(remainder, 0.25) - 3 * IQR(remainder) | # Ondergrens
remainder > quantile(remainder, 0.75) + 3 * IQR(remainder) # Bovengrens
)
# Visualisatie van de uitbijters
# Een boxplot toont de restcomponent, waarbij uitbijters worden gemarkeerd.
decomp %>%
ggplot(aes(x = "", y = remainder)) + # Gebruik een lege x-as voor een enkele boxplot
geom_boxplot(outlier.color = "red", coef = 3, fill = "skyblue", color = "black") +
labs(
title = "Boxplot van de restcomponent",
x = NULL,
y = "Restcomponent"
) +
theme_minimal()
# Verwijderen van uitbijters
# Uitbijters worden verwijderd en de data wordt opnieuw geïnterpoleerd.
dag_zwaarte_clean <- dag_zwaarte %>%
anti_join(filezwaarte_outliers, by = "BeginDatum") %>% # Verwijder uitbijters
fill_gaps(filezwaarte = NA) %>% # Vul opnieuw hiaten in met NA
mutate(filezwaarte = na.approx(filezwaarte)) # Vervang NA's met lineaire interpolatie
# Nieuwe STL-decompositie
# Na het opschonen van de data wordt opnieuw een STL-decompositie uitgevoerd.
decomp_clean <- dag_zwaarte_clean %>%
model(STL(filezwaarte ~ season(period = 365), robust = TRUE)) %>%
components()
# Visualisatie van de nieuwe decompositie
# De decompositie wordt opnieuw gevisualiseerd, nu zonder de invloed van uitbijters.
decomp_clean %>%
autoplot() +
labs(
title = "STL-decompositie van dagelijkse filezwaarte (zonder uitbijters)",
subtitle = "Filezwaarte = trend + seizoen (365 dagen) + restcomponent",
x = "Datum",
y = "Filezwaarte"
) +
theme_minimal()
# Analyse van de restcomponent
# De restcomponent wordt verder geanalyseerd op statistische validiteit.
remainder <- decomp_clean$remainder
# ACF van de restcomponent
# Controleer of er autocorrelatie aanwezig is in de restcomponent.
acf(remainder, main = "ACF van de restcomponent")
# Ljung-Box-test
# Test voor significante autocorrelatie.
Box.test(remainder, lag = 20, type = "Ljung-Box")
##
## Box-Ljung test
##
## data: remainder
## X-squared = 16146, df = 20, p-value < 2.2e-16
# Q-Q-plot voor normaalverdeling
# Controleer visueel of de restcomponent normaal verdeeld is.
qq_plot <- ggplot(data.frame(sample = remainder), aes(sample = sample)) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "Q-Q-plot van de Restcomponent") +
theme_minimal()
grid.arrange(qq_plot, ncol = 1)
# Shapiro-Wilk-test voor normaalverdeling
# Neem een steekproef om de test uit te voeren op een grote dataset.
set.seed(123) # Zorg voor reproduceerbaarheid
sampled_remainder <- sample(remainder, size = 5000)
shapiro.test(sampled_remainder)
##
## Shapiro-Wilk normality test
##
## data: sampled_remainder
## W = 0.98297, p-value < 2.2e-16
# Breusch-Pagan-test voor homoscedasticiteit
# Controleer of de variantie van de restcomponent constant is.
bp_test <- bptest(remainder ~ seq_along(remainder))
print(bp_test)
##
## studentized Breusch-Pagan test
##
## data: remainder ~ seq_along(remainder)
## BP = 1.2641, df = 1, p-value = 0.2609
### boxplot maken om inzicht te krijgen in de algemene data ###
# Bereken de mediaan
median_value <- median(dag_zwaarte_clean$filezwaarte, na.rm = TRUE)
# Bereken de grenzen
lower_bound <- median_value - 20000
upper_bound <- median_value + 20000
# Bepaal welke waarden buiten het bereik liggen
outside_bounds <- dag_zwaarte_clean$filezwaarte < lower_bound |
dag_zwaarte_clean$filezwaarte > upper_bound
# Bereken het percentage buiten het bereik
percentage_outside <- sum(outside_bounds, na.rm = TRUE) / nrow(dag_zwaarte_clean) * 100
cat("Percentage buiten bereik:", percentage_outside, "%\n")
## Percentage buiten bereik: 53.53639 %
# Maak de boxplot met rode grenzen
ggplot(dag_zwaarte_clean, aes(y = filezwaarte)) +
geom_boxplot(fill = "skyblue", color = "darkblue", outlier.color = "red") +
geom_hline(yintercept = lower_bound, color = "red") + # Ondergrens
geom_hline(yintercept = upper_bound, color = "red") + # Bovengrens
labs(
title = "Boxplot van Dagelijkse Filezwaarte met Grenzen",
y = "Dagelijkse Filezwaarte (km·min)"
) +
theme_minimal()
### Bereken het percentage buiten de grens van 20000 ###
percentage_outside_remainder <- sum(decomp$remainder < -20000 | decomp$remainder > 20000) / nrow(decomp) * 100
cat("Percentage buiten bereik remainder:", percentage_outside_remainder, "%\n")
## Percentage buiten bereik remainder: 35.37531 %
### Maak de boxplot van de remainders (al eerder gedaan) ###
ggplot(decomp_clean, aes(y = remainder)) +
geom_boxplot(fill = "skyblue", color = "darkblue", outlier.color = "red") +
geom_hline(yintercept = -20000, color = "red") + # Ondergrens
geom_hline(yintercept = 20000, color = "red") + # Bovengrens
labs(
title = "Boxplot van Dagelijkse remainders met Grenzen",
y = "Dagelijkse remainders (km·min)"
) +
theme_minimal()
### algemene data opslpitsen in dagen aangezien daar waarschijnlijk een verschil tussen zit
# hier eerst een lijn diagram van maken ###
# Voeg kolommen toe voor weekdagen en jaar
dag_zwaarte_clean_dagen <- dag_zwaarte_clean %>%
mutate(
Weekdag = weekdays(BeginDatum),
Jaar = as.numeric(format(BeginDatum, "%Y"))
)
dag_zwaarte_clean_dagen$BeginDatum <- NULL
# Zet de weekdagen in de juiste volgorde
dag_zwaarte_clean_dagen$Weekdag <- factor(dag_zwaarte_clean_dagen$Weekdag,
levels = c("maandag", "dinsdag", "woensdag",
"donderdag", "vrijdag", "zaterdag", "zondag"))
# Bereken gemiddelde filezwaarte per jaar en weekdag
seasonal_summary <- dag_zwaarte_clean_dagen %>%
group_by(Jaar, Weekdag) %>%
summarise(mean_zwaarte = mean(filezwaarte), .groups = "drop")
# Maak de seizoensplot
ggplot(seasonal_summary, aes(x = Weekdag, y = mean_zwaarte, color = factor(Jaar), group = Jaar)) +
geom_line(linewidth = 1) +
scale_color_viridis_d() + # Mooie kleuren voor jaren
labs(
title = "Seizoensplot: Gemiddelde filezwaarte per weekdag",
x = "Weekdag",
y = "Gemiddelde filezwaarte (km-min)",
color = "Jaar"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
### Boxplot maken met grenzen van RWS per dag ###
# Bereken de stats per dag
stats_per_day <- dag_zwaarte_clean_dagen %>%
group_by(Weekdag) %>%
summarise(median_value = median(filezwaarte),
total_count = n(),
outside_count = sum(filezwaarte < (median(filezwaarte) - 20000) |
filezwaarte > (median(filezwaarte) + 20000))) %>%
mutate(lower_bound = median_value - 20000,
upper_bound = median_value + 20000,
percentage_outside = (outside_count / total_count) * 100)
# plot maken
ggplot(dag_zwaarte_clean_dagen, aes(x = Weekdag, y = filezwaarte)) +
geom_boxplot() +
geom_segment(data = stats_per_day,
aes(x = as.numeric(Weekdag) - 0.4, xend = as.numeric(Weekdag) + 0.4,
y = lower_bound, yend = lower_bound),
color = "red") +
geom_segment(data = stats_per_day,
aes(x = as.numeric(Weekdag) - 0.4, xend = as.numeric(Weekdag) + 0.4,
y = upper_bound, yend = upper_bound),
color = "red") +
labs(title = "Daily Traffic per Day of the Week",
x = "Day of the Week",
y = "Daily Traffic") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
### Percentage buiten de grens van RWS per dag
ggplot(stats_per_day, aes(x = Weekdag, y = percentage_outside)) +
geom_bar(stat = "identity", fill = "skyblue") +
labs(
title = "Percentage van data buiten mediaan ± 20.000 per dag",
x = "Dag van de week",
y = "Percentage buiten grenzen (%)"
) +
theme_minimal()
### Alles wat we hiervoor hebben gedaan voor de algemene data ook doen voor alleen de remainders ###
decomp_clean <- dag_zwaarte_clean %>%
model(STL(filezwaarte ~ season(period = 365), robust = TRUE)) %>%
components()
# Voeg kolommen toe voor weekdag en jaar
decomp_clean_dagen <- decomp_clean %>%
mutate(
Weekdag = weekdays(BeginDatum),
Jaar = as.numeric(format(BeginDatum, "%Y"))
)
decomp_clean_dagen <- decomp_clean_dagen %>%
select(-.model, -filezwaarte, -trend, -season_365, -season_adjust)
decomp_clean_dagen$BeginDatum <- NULL
# Zet de weekdagen in de juiste volgorde
decomp_clean_dagen$Weekdag <- factor(decomp_clean_dagen$Weekdag,
levels = c("maandag", "dinsdag", "woensdag",
"donderdag", "vrijdag", "zaterdag", "zondag"))
# Bereken gemiddelde filezwaarte per jaar en weekdag
seasonal_summary <- decomp_clean_dagen %>%
group_by(Jaar, Weekdag) %>%
summarise(mean_zwaarte = mean(remainder), .groups = "drop")
# Maak de seizoensplot
ggplot(seasonal_summary, aes(x = Weekdag, y = mean_zwaarte, color = factor(Jaar), group = Jaar)) +
geom_line(linewidth = 1) +
scale_color_viridis_d() + # Mooie kleuren voor jaren
labs(
title = "Seizoensplot: Gemiddelde remainders per weekdag",
x = "Weekdag",
y = "Gemiddelde remainders (km-min)",
color = "Jaar"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Bereken de stats per dag
stats_per_day <- decomp_clean_dagen %>%
group_by(Weekdag) %>%
summarise(total_count = n(),
outside_count = sum(remainder < -20000 | remainder > 20000)) %>%
mutate(percentage_outside = (outside_count / total_count) * 100)
ggplot(decomp_clean_dagen, aes(x = Weekdag, y = remainder)) +
geom_boxplot() +
geom_segment(data = stats_per_day,
aes(x = as.numeric(Weekdag) - 0.4, xend = as.numeric(Weekdag) + 0.4,
y = -20000, yend = -20000),
color = "red") +
geom_segment(data = stats_per_day,
aes(x = as.numeric(Weekdag) - 0.4, xend = as.numeric(Weekdag) + 0.4,
y = 20000, yend = 20000),
color = "red") +
labs(title = "Daily remainders per Day of the Week",
x = "Day of the Week",
y = "Daily remainders") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplot(stats_per_day, aes(x = Weekdag, y = percentage_outside)) +
geom_bar(stat = "identity", fill = "skyblue") +
labs(
title = "Percentage remainders ± 20.000 per dag",
x = "Dag van de week",
y = "Percentage buiten remainders (%)"
) +
theme_minimal()
# Maak trainings- en testset aan
train <- dag_zwaarte_clean %>% filter_index(. ~ "2022-12-31")
test <- dag_zwaarte_clean %>% filter_index("2023-01-01" ~ .)
# Maak meerdere trainingsets aan (sliding window van 3 jaar)
train_sets <- train %>% slide_tsibble(.size = 1095, .step = 365, .id = "id") # 3 jaar = 1095 dagen)
# fit de verschillende modellen en genereer forecasts; dit duurt ook wel even
forecasts <- train_sets %>%
model(
Drift = RW(filezwaarte ~ drift()),
Mean = MEAN(filezwaarte),
Naive = NAIVE(filezwaarte),
`Seasonal naive` = SNAIVE(filezwaarte),
ETS = ETS(filezwaarte),
ARIMA = ARIMA(filezwaarte)) %>%
forecast(h = 30)
# Voeg kolom toe met forecast horizon
forecasts <- forecasts %>%
group_by(.model, id) %>%
mutate(h = row_number()) %>%
ungroup()
# Voeg werkelijke waarden toe aan de voorspellingen
forecasts <- forecasts %>%
left_join(dag_zwaarte_clean, by = c("BeginDatum"))
# Bereken absolute fouten
forecasts <- forecasts %>%
mutate(error = .mean - filezwaarte.y) # Gebruik de werkelijke waarden (filezwaarte.y)
# rood als buiten de 20000 grens
forecasts <- forecasts %>%
group_by(.model, h) %>%
mutate(exceeded = ifelse(min(error, na.rm = TRUE) < -20000 | max(error, na.rm = TRUE) > 20000, TRUE, FALSE)) %>%
ungroup()
# tenslotte plotten we alles
forecasts %>%
ggplot(aes(color = exceeded, x = as.factor(h), y = error)) +
geom_boxplot() +
geom_hline(yintercept = 20000, linetype = "dashed", colour = "blue") +
geom_hline(yintercept = -20000, linetype = "dashed", colour = "blue") +
scale_color_manual(values = c("black", "red")) +
theme(legend.position = "none") +
labs(x = "Forecast horizon (in days)",
y = "Absolute error",
title = "Absolute errors depending on forecast horizon",
subtitle = "") +
facet_wrap(vars(.model))
# Bereken de wekelijkse filezwaarte, omdat we 1 week moeten voorspellen
# We transformeren de dagelijkse filezwaarte naar wekelijkse data door de totale filezwaarte per week te berekenen.
week_zwaarte <- filedata %>%
mutate(Week = floor_date(BeginDatum, unit = "week")) %>% # Rond de datum af naar het begin van de week
group_by(Week) %>%
summarize(filezwaarte = sum(zwaarteKmMin, na.rm = TRUE), .groups = "drop") %>% # Som filezwaarte per week
as_tsibble(index = Week) %>% # Zet de data om naar een tijdreeksformaat (tsibble)
fill_gaps(filezwaarte = NA) %>% # Vul ontbrekende weken in met NA
mutate(filezwaarte = na.approx(filezwaarte)) # Vervang de NA's met lineaire interpolatie
# Splits de wekelijkse data in training- en testsets
# We verdelen de data in een trainingsset (80%) en een testset (20%) om onze modellen te trainen en valideren.
n_total <- nrow(week_zwaarte) # Totaal aantal weken
n_train <- floor(0.8 * n_total) # Bereken de omvang van de trainingsset (80%)
file_train <- week_zwaarte[1:n_train, ] # Trainingsset
file_test <- week_zwaarte[(n_train + 1):n_total, ] # Testset
# Zet de data om naar tijdreeksformaat
# Dit is nodig voor het trainen van tijdreeksmodellen zoals ARIMA en ETS.
file_train_ts <- ts(file_train$filezwaarte, frequency = 52) # Tijdreeks met wekelijkse frequentie
file_test_ts <- ts(file_test$filezwaarte, frequency = 52) # Tijdreeks voor de testset
# Train een ARIMA-model
# We gebruiken een automatisch ARIMA-model omdat het flexibel is en goed presteert bij veel soorten tijdreeksen.
auto_model <- auto.arima(file_train_ts)
# Train een ETS-model
# We trainen een ETS-model met behulp van stlf() om seizoensvoorspellingen te genereren.
ets_stlf_model <- stlf(file_train_ts, method = "ets", h = length(file_test_ts))
# Voorspellingen maken voor de testperiode
# We maken voorspellingen voor zowel ARIMA als ETS en vergelijken de resultaten.
forecast_auto <- forecast(auto_model, h = length(file_test_ts)) # Voorspellingen met ARIMA
forecast_ets_stlf <- forecast(ets_stlf_model, h = length(file_test_ts)) # Voorspellingen met ETS
# RMSE berekenen
# RMSE (Root Mean Square Error) meet de nauwkeurigheid van de voorspellingen.
rmse_auto <- sqrt(mean((as.numeric(forecast_auto$mean) - as.numeric(file_test_ts))^2))
rmse_ets_stlf <- sqrt(mean((as.numeric(forecast_ets_stlf$mean) - as.numeric(file_test_ts))^2))
# Print RMSE-waarden en selecteer het beste model
cat("RMSE ARIMA:", rmse_auto, "\n") # RMSE van het ARIMA-model
## RMSE ARIMA: 143513.5
cat("RMSE ETS (stlf):", rmse_ets_stlf, "\n") # RMSE van het ETS-model
## RMSE ETS (stlf): 146666.5
# Kies het beste model op basis van de laagste RMSE
if (rmse_auto < rmse_ets_stlf) {
cat("ARIMA heeft een lagere RMSE dan ETS. Daarom kiezen we ARIMA voor de voorspellingen.\n")
best_model <- "ARIMA"
} else {
cat("ETS heeft een lagere RMSE dan ARIMA. Daarom kiezen we ETS voor de voorspellingen.\n")
best_model <- "ETS"
}
## ARIMA heeft een lagere RMSE dan ETS. Daarom kiezen we ARIMA voor de voorspellingen.
# Combineer voorspellingen en werkelijke waarden
# Om de prestaties van de modellen te visualiseren, combineren we de voorspellingen met de werkelijke waarden.
test_data <- data.frame(
Date = file_test$Week, # Datums van de testperiode
Actual = file_test$filezwaarte, # Werkelijke waarden
ARIMA = as.numeric(forecast_auto$mean), # Voorspellingen van ARIMA
ETS = as.numeric(forecast_ets_stlf$mean) # Voorspellingen van ETS
)
# Visualisatie van de voorspellingen versus de werkelijke waarden
# De lijnplot toont hoe goed de modellen de werkelijke waarden voorspellen.
ggplot(test_data, aes(x = Date)) +
geom_line(aes(y = Actual, color = "Actual"), linewidth = 1) + # Werkelijke waarden
geom_line(aes(y = ARIMA, color = "ARIMA"), linewidth = 1) + # Voorspellingen van ARIMA
geom_line(aes(y = ETS, color = "ETS"), linewidth = 1) + # Voorspellingen van ETS
labs(
title = "Voorspellingen versus werkelijke waarden",
y = "Filezwaarte",
x = "Datum",
color = "Legenda"
) +
theme_minimal()
# Analyse van residuen (alleen voor ARIMA)
# We analyseren de residuen van ARIMA om te controleren of er nog patronen of systematische afwijkingen zijn.
residuals_auto <- residuals(auto_model)
# ACF van de residuen
# Controleer op autocorrelatie in de residuen om de validiteit van het model te evalueren.
acf(residuals_auto, main = "ACF van residuen (ARIMA)")
# Bereken het gemiddelde van de residuen
# Een gemiddeld residu van ongeveer nul wijst op een goed fit model.
mean_residuals_auto <- mean(residuals_auto, na.rm = TRUE)
cat("Gemiddelde residuen (ARIMA):", mean_residuals_auto, "\n")
## Gemiddelde residuen (ARIMA): -741.5857
# Voorspelling maken voor 1 week vooruit met ARIMA
# Omdat ARIMA het beste model is, maken we een voorspelling voor 1 week vooruit.
forecast_auto_one_week <- forecast(auto_model, h = 1)
# Rond de voorspelling af op duizendtallen
# Dit maakt de voorspelling gemakkelijker te interpreteren.
forecasts_auto <- round(as.numeric(forecast_auto_one_week$mean) / 1000) * 1000
# Label de voorspelling met de datum van de volgende week
forecast_week_auto <- seq(from = max(file_test$Week) + 7, by = "week", length.out = 1)
names(forecasts_auto) <- format(forecast_week_auto, "%d %B %Y")
# Print de voorspelling
cat("Voorspelling voor 1 week vooruit (ARIMA):\n")
## Voorspelling voor 1 week vooruit (ARIMA):
print(forecasts_auto)
## 07 januari 2024
## 250000
# Weersdataset inlezen en voorbereiden
# De dataset met weergegevens wordt ingelezen en voorbereid om te combineren met de filezwaartegegevens.
weer <- read.csv(path_weer, sep = ";")
# Kolomnamen omvormen en data transformeren
# We zetten de datum om naar het juiste formaat en maken de dataset breder door metingen per locatie te pivoteren.
weersdata <- weer %>%
mutate(YYYYMMDD = as.Date(as.character(YYYYMMDD), format = "%Y%m%d")) %>% # Datum omzetten naar correct formaat
pivot_wider(names_from = meting, values_from = meetwaarde) %>% # Metingen per kolom groeperen
group_by(YYYYMMDD) %>%
filter(!(locatie == "De Bilt" & as.Date(YYYYMMDD) < as.Date("1906-01-01"))) %>% # Data vóór 1906 voor De Bilt verwijderen voor consistentie
ungroup()
# Verwijderen van kolommen met teveel missende waarden
# Kolommen met meer dan 40% missende waarden worden verwijderd om de dataset schoon en bruikbaar te houden.
na_percentage <- colMeans(is.na(weersdata)) * 100
kolommen_te_verwijderen <- names(weersdata)[na_percentage >= 40]
weersdata <- weersdata[, na_percentage < 40]
# Invullen van missende waarden
# Overgebleven missende waarden worden opgevuld met het gemiddelde van de betreffende kolom.
weersdata <- weersdata %>%
mutate(across(where(is.numeric), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .)))
# Bereken gemiddelde per dag (indien meerdere locaties per dag)
# Om consistente dagwaarden te krijgen, berekenen we het gemiddelde van numerieke kolommen per dag.
weersdata <- weersdata %>%
group_by(YYYYMMDD) %>%
summarize(across(where(is.numeric), mean, na.rm = TRUE)) %>%
ungroup()
## Warning: There was 1 warning in `summarize()`.
## ℹ In argument: `across(where(is.numeric), mean, na.rm = TRUE)`.
## ℹ In group 1: `YYYYMMDD = 1906-01-01`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
# Filter weersdata voor overlappende periode
# We beperken de data tot het bereik van 2000-01-01 tot 2023-12-31 om consistent te blijven met de filezwaartegegevens.
weersdata <- weersdata %>%
filter(YYYYMMDD >= as.Date("2000-01-01") & YYYYMMDD <= as.Date("2023-12-31"))
# Maak van de filezwaartegegevens een tijdreeks
# Dit is nodig om later te combineren met de weersgegevens.
dag_zwaarte <- dag_zwaarte %>%
as_tsibble(index = BeginDatum)
# Combineer de weersgegevens met de filezwaartegegevens
# Door de datasets te koppelen op de datumkolom, voegen we weersgegevens toe aan de filezwaartedata.
combined_data <- dag_zwaarte %>%
left_join(weersdata, by = c("BeginDatum" = "YYYYMMDD"))
# Bereken correlaties met filezwaarte
# We onderzoeken welke weergerelateerde variabelen een relatie hebben met de filezwaarte.
cor_data <- combined_data %>%
as_tibble() %>%
select(-BeginDatum, filezwaarte, FG, TG, SQ, DR, RH, NG, UG) # Selecteer relevante variabelen
# We berekenen de correlatie van elke variabele met filezwaarte.
cor_values <- cor(cor_data, use = "complete.obs")["filezwaarte", ]
# Selecteer variabelen met significante correlatie
# We kiezen variabelen met een absolute correlatie > 0.08, gebaseerd op onze drempel voor relevantie.
relevante_variabelen <- names(cor_values[abs(cor_values) > 0.08])
# Visualisatie van de relatie tussen filezwaarte en neerslagduur (DR)
# We gebruiken scatterplots om de verbanden tussen weerfactoren en filezwaarte te illustreren.
combined_data %>%
ggplot(aes(x = DR, y = filezwaarte)) +
geom_point(alpha = 0.5) + # Puntenplot met transparantie
geom_smooth(method = "lm", se = FALSE, color = "blue") + # Lineaire trendlijn toevoegen
labs(
title = "Relatie tussen DR en filezwaarte",
x = "Duur neerslag in 0,1 uur",
y = "Filezwaarte (km∙min)"
)
## `geom_smooth()` using formula = 'y ~ x'
# Visualisatie van de relatie tussen filezwaarte en neerslaghoeveelheid (RH)
combined_data %>%
ggplot(aes(x = RH, y = filezwaarte)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(
title = "Relatie tussen neerslag (RH) en filezwaarte",
x = "Neerslag (0.1 mm)",
y = "Filezwaarte (km∙min)"
)
## `geom_smooth()` using formula = 'y ~ x'
# Lineaire regressie met TSLM
# We voeren een regressieanalyse uit waarbij filezwaarte de afhankelijke variabele is en weerfactoren de voorspellers.
combined_data %>%
model(TSLM(filezwaarte ~ DR + RH + RHX + PG + PN + UG)) %>% # TSLM met geselecteerde variabelen
report() # Samenvatting van het model
## Series: filezwaarte
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -68600 -24396 -4986 18481 323153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92512.868 37134.063 2.491 0.01275 *
## DR 57.095 27.592 2.069 0.03855 *
## RH 142.458 30.449 4.679 2.93e-06 ***
## RHX -270.259 58.711 -4.603 4.22e-06 ***
## PG -8.833 14.570 -0.606 0.54438
## PN 1.682 13.675 0.123 0.90213
## UG 124.502 39.557 3.147 0.00165 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28530 on 8759 degrees of freedom
## Multiple R-squared: 0.02561, Adjusted R-squared: 0.02495
## F-statistic: 38.38 on 6 and 8759 DF, p-value: < 2.22e-16
# Combineer filedata met weersgegevens
model_data <- dag_zwaarte %>%
left_join(weersdata, by = c("BeginDatum" = "YYYYMMDD"))
# De dataset met filezwaarte wordt gecombineerd met weersdata op basis van de datum.
# Dit is nodig om weersinvloeden op de filezwaarte te kunnen analyseren.
# Splits de data in training- en testsets (80% voor training)
n_total <- nrow(model_data)
n_train <- floor(0.8 * n_total) # 80% van de data wordt gebruikt voor het trainen van het model.
train <- model_data[1:n_train, ]
test <- model_data[(n_train + 1):n_total, ]
# De training- en testdata worden gescheiden, zodat het model eerst kan leren en daarna wordt geëvalueerd.
# Zet training- en testsets om naar tijdreeksen (tsibble-formaat)
train <- train %>%
as_tsibble(index = BeginDatum)
test <- test %>%
as_tsibble(index = BeginDatum)
# Het omzetten naar tsibble-formaat is nodig voor tijdreeksanalyse in R.
# SARIMA-model met regressoren trainen
sarima_model <- train %>%
model(
SARIMA = ARIMA(filezwaarte ~ RH + RHX + UG) # ARIMA-model met relevante variabelen
)
# Hier gebruiken we het SARIMA-model om filezwaarte te voorspellen, rekening houdend met significante weersvariabelen
# (RH: dagelijkse regen, RHX: maximale regen in een uur, UG: gemiddelde windsnelheid).
# Dynamische regressie (TSLM) model trainen
tslm_model <- train %>%
model(
TSLM = TSLM(filezwaarte ~ RH + RHX + UG) # Lineair model met dezelfde variabelen als SARIMA
)
# TSLM-model wordt gebruikt als vergelijkingsmodel. Dit model werkt lineair en is minder complex dan SARIMA.
# Voorspellingen maken met SARIMA
sarima_forecast <- sarima_model %>%
forecast(new_data = test %>% select(BeginDatum, RH, RHX, UG, filezwaarte))
# SARIMA voorspelt de filezwaarte voor de testset, gebaseerd op de geselecteerde weersvariabelen.
# Bereken RMSE (Root Mean Square Error) voor SARIMA
sarima_accuracy <- accuracy(sarima_forecast, test)
cat("SARIMA RMSE:", sarima_accuracy$RMSE, "\n")
## SARIMA RMSE: 29860.17
# De RMSE voor SARIMA wordt berekend om de nauwkeurigheid van het model te evalueren. Hoe lager de RMSE, hoe beter het model.
# Voorspellingen maken met TSLM
tslm_forecast <- tslm_model %>%
forecast(new_data = test %>% select(BeginDatum, RH, RHX, UG, filezwaarte))
# TSLM voorspelt de filezwaarte voor dezelfde testset.
# Bereken RMSE voor TSLM
tslm_accuracy <- accuracy(tslm_forecast, test)
cat("Dynamic TSLM RMSE:", tslm_accuracy$RMSE, "\n")
## Dynamic TSLM RMSE: 32942.95
# De RMSE voor TSLM wordt berekend om te vergelijken met SARIMA.
# Vergelijk de RMSE-waarden
rmse_results <- tibble(
Model = c("SARIMA met regressoren", "TSLM met regressoren"),
RMSE = c(sarima_accuracy$RMSE, tslm_accuracy$RMSE)
)
print(rmse_results)
## # A tibble: 2 × 2
## Model RMSE
## <chr> <dbl>
## 1 SARIMA met regressoren 29860.
## 2 TSLM met regressoren 32943.
# De RMSE-waarden van beide modellen worden in een tabel gezet, zodat we eenvoudig kunnen zien welk model beter presteert.
# Selecteer het beste model op basis van de laagste RMSE
best_model <- rmse_results %>%
filter(RMSE == min(RMSE)) %>%
pull(Model)
cat("Het beste model is:", best_model, "\n")
## Het beste model is: SARIMA met regressoren
# Het model met de laagste RMSE wordt als het beste model geselecteerd. Dit wordt ook vermeld in de uitvoer.
# Voorspellingen maken met het beste model
if (best_model == "SARIMA met regressoren") {
best_forecast <- sarima_model %>%
forecast(new_data = test %>% select(BeginDatum, RH, RHX, UG, filezwaarte))
} else {
best_forecast <- tslm_model %>%
forecast(new_data = test %>% select(BeginDatum, RH, RHX, UG, filezwaarte))
}
# Het geselecteerde model (SARIMA of TSLM) wordt gebruikt om voorspellingen te maken op basis van de testdata.
# Beste model uit opgave 3 blijkt ARIMA zonder regressoren
# Splits de data in trainings- en testset
train <- dag_zwaarte %>%
filter(BeginDatum < as.Date("2020-03-01")) # Data vóór maart 2020
test <- dag_zwaarte %>%
filter(BeginDatum >= as.Date("2020-03-01") & BeginDatum <= as.Date("2020-03-31")) # Data van maart 2020
rest <- dag_zwaarte %>%
filter(BeginDatum >= as.Date("2020-04-01")) # Data na maart 2020
# Zet data om naar tsibble-formaat
train <- train %>% as_tsibble(index = BeginDatum)
test <- test %>% as_tsibble(index = BeginDatum)
# Pas ARIMA en ETS modellen toe
fit_ARIMA <- train %>%
model(ARIMA = ARIMA(filezwaarte))
# Genereer voorspellingen voor testset
forecasts <- fit_ARIMA %>% forecast(h = nrow(test))
# Bereken nauwkeurigheid voor beide modellen
accuracy_results <- forecasts %>%
accuracy(test) %>%
select(.model, RMSE, MAPE)
print(accuracy_results)
## # A tibble: 1 × 3
## .model RMSE MAPE
## <chr> <dbl> <dbl>
## 1 ARIMA 27978. 2020.
# Plot alleen de gefilterde data
forecasts %>%
autoplot(train, level = 95) +
scale_x_date(limits = c(as.Date("2019-10-01"), as.Date("2020-03-31"))) +
geom_line(data = test, aes(x = BeginDatum, y = filezwaarte)) + # Voeg de testlijn toe
#geom_line(data = rest, aes(x = BeginDatum, y = filezwaarte)) + # Voeg de testlijn toe
labs(
title = "Voorspelling van dagelijkse filezwaarte in maart 2020",
x = "Datum",
y = "Filezwaarte (km·min)",
caption = "Inclusief 95%-voorspellingsinterval"
) +
theme_minimal()
## Warning: Removed 7213 rows containing missing values or values outside the scale range
## (`geom_line()`).
# Voeg een kolom toe met de dag van de maand
train_march <- train %>%
filter(month(BeginDatum) == 3 & year(BeginDatum) < 2020 & year(BeginDatum) > 2015) %>%
mutate(Day = day(BeginDatum)) # Dag van de maand
forecast_march <- forecasts %>%
filter(lubridate::month(BeginDatum) == 3) %>%
mutate(Day = day(BeginDatum)) # Dag van de maand
forecast_march_95 <- forecast_march %>%
mutate(
.lower = .mean - 1.96 * sd(test$filezwaarte),
.upper = .mean + 1.96 * sd(test$filezwaarte)
)
# Plot de data
ggplot() +
geom_line(data = train_march, aes(x = Day, y = filezwaarte, color = as.factor(year(BeginDatum))), alpha = 1) + # Trainingsdata
geom_ribbon(data = forecast_march_95, aes(x = Day, ymin = .lower, ymax = .upper), fill = "blue", alpha = 0.1) +
geom_line(data = forecast_march_95, aes(x = Day, y = .mean)) + # Voorspellingen
labs(
title = "Vergelijking van dagelijkse filezwaarte in maart (2000-2012)",
x = "Dag van maart",
y = "Filezwaarte (km·min)",
color = "Jaar",
caption = "Voorspelde maart 2020 weergegeven in rood"
) +
theme_minimal() +
scale_x_continuous(breaks = 1:31) # Zorg dat de x-as alleen 1 t/m 31 toont
