Forecasting Sleep Duration with Time-Series Models
Jackson Sbrissa
2026-10-09
Forecasting Sleep Duration with Time-Series Models
This tutorial explores daily sleep duration using generated demo data. It cleans the records, explores patterns, imputes missing nights, compares five forecast models on a 12-week holdout, checks rolling validation, tunes model variants, then makes a 12-week ETS forecast. The simulated data do not describe a real person. An optional section explains how to download Garmin data locally.
Question: Which forecasting approach gives useful predictions of daily sleep duration, and how does ETS forecast the next 12 weeks?
The generated values are for teaching only. They are not evidence about a real athlete, health advice, or a claim that these models perform well on personal data.
1. Get ready
Use RStudio. Create a new R Markdown file, replace
its contents with this tutorial, and save it as
sport-time-series-forecasting.Rmd. Run chunks in order
using each green triangle. Install these packages once in the
Console if prompted:
install.packages(c("tidyverse", "lubridate", "fpp3", "feasts", "ggtime", "timetk", "cowplot"))When all chunks run, click Knit to make an HTML page. The packages provide data cleaning, time-series plots, models, forecasts, and missing-value imputation.
library(tidyverse)
library(lubridate)
library(fpp3)
library(feasts)
library(fabletools)
library(ggtime)
library(timetk)2. Create a demo dataset and prepare it
For privacy and reproducibility, generate several years of made-up daily records. The data include a weekly pattern, gradual trend, yearly variation, random variation, and some missing nights marked as zero seconds to represent non-wear records.
set.seed(2026)
dates <- seq.Date(as.Date("2023-01-01"), as.Date("2026-09-28"), by = "day")
day_number <- seq_along(dates)
wake_day <- factor(wday(dates, label = TRUE, abbr = TRUE, week_start = 1),
levels = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))
weekday_effect <- c(Mon = 0.12, Tue = 0.08, Wed = 0.05, Thu = 0.02,
Fri = 0.00, Sat = -0.20, Sun = -0.14)
yearly_effect <- 0.12 * sin(2 * pi * yday(dates) / 365.25)
sleep_h_true <- 7.1 +
weekday_effect[as.character(wake_day)] +
0.00015 * day_number + yearly_effect +
rnorm(length(dates), mean = 0, sd = 0.45)
missing_night <- runif(length(dates)) < 0.04
sleep_raw <- tibble(
date = dates,
sleep_s = if_else(missing_night, 0, pmax(sleep_h_true, 3) * 3600)
)
sleep_raw |> slice_head(n = 10)## # A tibble: 10 × 2
## date sleep_s
## <date> <dbl>
## 1 2023-01-01 25907.
## 2 2023-01-02 24259.
## 3 2023-01-03 26097.
## 4 2023-01-04 25635.
## 5 2023-01-05 24592.
## 6 2023-01-06 21532.
## 7 2023-01-07 23705.
## 8 2023-01-08 23467.
## 9 2023-01-09 26247.
## 10 2023-01-10 25160.Clean and prepare the daily records. Durations are stored in seconds,
so convert them to hours. Treat zero duration as a missing night.
fill_gaps() makes the daily time index explicit, and
ts_impute_vec() fills missing values for modelling while
leaving the original sleep_h column available for
checking.
sleep_df <- sleep_raw |>
mutate(
date = as.Date(date),
sleep_s = na_if(sleep_s, 0),
sleep_h = sleep_s / 3600,
wake_day = wday(date, label = TRUE, week_start = 1),
weekend = wake_day %in% c("Sat", "Sun"),
month = month(date, label = TRUE),
short_night = sleep_h < 3
) |>
select(date, sleep_h, wake_day, weekend, month, short_night)
sleep_ts <- sleep_df |>
as_tsibble(index = date) |>
fill_gaps() |>
mutate(sleep_h_imp = timetk::ts_impute_vec(sleep_h, period = 7))
sleep_ts |> slice_head(n = 10)## # A tsibble: 10 x 7 [1D]
## date sleep_h wake_day weekend month short_night sleep_h_imp
## <date> <dbl> <ord> <lgl> <ord> <lgl> <dbl>
## 1 2023-01-01 7.20 Sun TRUE Jan FALSE 7.20
## 2 2023-01-02 6.74 Mon FALSE Jan FALSE 6.74
## 3 2023-01-03 7.25 Tue FALSE Jan FALSE 7.25
## 4 2023-01-04 7.12 Wed FALSE Jan FALSE 7.12
## 5 2023-01-05 6.83 Thu FALSE Jan FALSE 6.83
## 6 2023-01-06 5.98 Fri FALSE Jan FALSE 5.98
## 7 2023-01-07 6.58 Sat TRUE Jan FALSE 6.58
## 8 2023-01-08 6.52 Sun TRUE Jan FALSE 6.52
## 9 2023-01-09 7.29 Mon FALSE Jan FALSE 7.29
## 10 2023-01-10 6.99 Tue FALSE Jan FALSE 6.993. Explore the time series
Plot the imputed series first. Look for trend, repeating weekly pattern, gaps, and unusual values.
sleep_ts |> autoplot(sleep_h_imp) +
labs(title = "Simulated daily sleep duration", x = "Date", y = "Sleep (hours)")
The next plot uses STL decomposition to separate the long-term trend and weekly seasonality.
sleep_ts |>
model(STL(sleep_h_imp ~ trend(window = 365) + season(period = 7))) |>
components() |>
autoplot()
Check autocorrelation and weekly patterns. ACF and PACF show relationships with earlier observations; seasonal and subseries plots align observations by weekday.
acf_sleep <- sleep_ts |> ACF(sleep_h_imp, lag_max = 28) |> autoplot()
pacf_sleep <- sleep_ts |> PACF(sleep_h_imp, lag_max = 28) |> autoplot()
cowplot::plot_grid(acf_sleep, pacf_sleep)
sleep_ts |> gg_season(sleep_h_imp, period = "week")
sleep_ts |> gg_subseries(sleep_h_imp, period = "week")
The following summaries check stationarity, differencing, and whether a transformation may help. They are diagnostics, not automatic instructions to change the model.
sleep_ts |> features(sleep_h_imp, list(unitroot_kpss, unitroot_ndiffs, unitroot_nsdiffs))## # A tibble: 1 × 4
## kpss_stat kpss_pvalue ndiffs nsdiffs
## <dbl> <dbl> <int> <int>
## 1 0.694 0.0141 1 0sleep_ts |> features(sleep_h_imp, guerrero)## # A tibble: 1 × 1
## lambda_guerrero
## <dbl>
## 1 1.30Compare observed sleep duration by weekday, then show where missing records were imputed.
sleep_ts |>
as_tibble() |>
ggplot(aes(wake_day, sleep_h)) +
geom_boxplot() +
labs(x = "Wake day", y = "Sleep (hours)")## Warning: Removed 50 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
sleep_ts |>
ggplot(aes(date, sleep_h_imp, colour = is.na(sleep_h))) +
geom_point(size = 0.5) +
labs(x = "Date", y = "Sleep (hours)", colour = "Imputed")
4. Fit candidate models and hold out the newest 12 weeks
Compare a mean benchmark, seasonal-naive benchmark, time-series linear model (TSLM), exponential smoothing (ETS), and automatic ARIMA. Keep the newest 12 weeks out of model fitting so forecasts can be compared with observations the models have not seen.
train_ts <- sleep_ts |> filter_index(~ "2026-07-06")
test_ts <- sleep_ts |> filter_index("2026-07-07" ~ .)
nrow(train_ts)## [1] 1283nrow(test_ts)## [1] 84The training data end on 6 July 2026. The test data start the next day and cover the final 12 weeks in this demo.
sleep_fit <- train_ts |>
model(
mean = MEAN(sleep_h_imp),
snaive = SNAIVE(sleep_h_imp),
tslm = TSLM(sleep_h_imp ~ trend() + season()),
ets = ETS(sleep_h_imp),
arima = ARIMA(sleep_h_imp)
)
sleep_fit## # A mable: 1 x 5
## mean snaive tslm ets arima
## <model> <model> <model> <model> <model>
## 1 <MEAN> <SNAIVE> <TSLM> <ETS(M,N,M)> <ARIMA(1,1,2)(2,0,0)[7]>The models represent different assumptions: no pattern (mean), repeat the previous week (SNAIVE), fixed trend and weekday effects (TSLM), changing level/trend/seasonality (ETS), and an automatically selected ARIMA model.
5. Check model errors and holdout accuracy
Inspect ETS residuals for patterns the model has not explained. Residuals should look broadly like random noise; clear structure can signal a model has missed something.
sleep_fit |> select(ets) |> gg_tsresiduals()
Compare in-sample fit and out-of-sample forecast accuracy. RMSE, MAE, and MASE are error measures; lower is better for comparisons on the same data and scale.
sleep_fc <- sleep_fit |> forecast(new_data = test_ts)
model_tbl <- bind_rows(
sleep_fit |> accuracy(),
sleep_fc |> accuracy(sleep_ts)
) |>
select(.model, .type, RMSE, MAE, MASE) |>
arrange(.type, RMSE)
model_tbl## # A tibble: 10 × 5
## .model .type RMSE MAE MASE
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 ets Test 0.472 0.382 0.795
## 2 mean Test 0.489 0.393 0.818
## 3 tslm Test 0.497 0.402 0.837
## 4 arima Test 0.502 0.407 0.847
## 5 snaive Test 0.597 0.487 1.01
## 6 ets Training 0.450 0.359 0.748
## 7 tslm Training 0.454 0.364 0.757
## 8 arima Training 0.461 0.368 0.765
## 9 mean Training 0.473 0.378 0.787
## 10 snaive Training 0.607 0.481 1Plot forecasts over the recent part of the series. The test period is included so you can compare predicted values with observed demo values.
sleep_fc |>
autoplot(sleep_ts |> filter_index("2026-05-01" ~ .), level = NULL) +
labs(title = "Candidate model forecasts and held-out sleep data",
x = "Date", y = "Sleep duration (hours)", colour = "Model")
6. Test robustness with rolling 12-week windows
One holdout split can give a lucky or unlucky result. Repeat the comparison over six consecutive 12-week test windows. In each round, fit models only on dates before the test window, then score forecasts on that window.
origins <- as.Date("2026-07-06") - 84 * 0:5
cv_results <- purrr::map_dfr(origins, function(o) {
train <- sleep_ts |> filter(date <= o)
test <- sleep_ts |> filter(date > o, date <= o + 84)
train |>
model(
mean = MEAN(sleep_h_imp),
snaive = SNAIVE(sleep_h_imp),
tslm = TSLM(sleep_h_imp ~ trend() + season()),
ets = ETS(sleep_h_imp),
arima = ARIMA(sleep_h_imp)
) |>
forecast(new_data = test) |>
accuracy(sleep_ts) |>
mutate(origin = o)
})## Warning in sqrt(diag(best$var.coef)): NaNs producedcv_summary <- cv_results |>
group_by(.model) |>
summarise(RMSE = mean(RMSE), MAE = mean(MAE), .groups = "drop") |>
arrange(RMSE)
cv_summary## # A tibble: 5 × 3
## .model RMSE MAE
## <chr> <dbl> <dbl>
## 1 ets 0.457 0.368
## 2 tslm 0.459 0.368
## 3 mean 0.478 0.384
## 4 arima 0.479 0.386
## 5 snaive 0.623 0.502The average scores summarize performance across several periods. They can differ from the single holdout results because each window contains different observations.
augment(sleep_fit) |>
features(.innov, ljung_box, lag = 14)## # A tibble: 5 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arima 17.8 2.14e- 1
## 2 ets 13.9 4.60e- 1
## 3 mean 106. 4.44e-16
## 4 snaive 280. 0
## 5 tslm 52.7 2.13e- 67. Tune model variants
This chunk checks alternative ETS, TSLM, ARIMA, STL-ETS, and combined forecasts over the same rolling windows, plus models fit using only the most recent 52 weeks. It can take longer than earlier chunks.
tuning_results <- purrr::map_dfr(origins, function(o) {
train <- sleep_ts |> filter(date <= o)
train_52 <- train |> filter(date > o - 365)
test <- sleep_ts |> filter(date > o, date <= o + 84)
full_fit <- train |>
model(
ets = ETS(sleep_h_imp),
ets_ANA = ETS(sleep_h_imp ~ error("A") + trend("N") + season("A")),
tslm = TSLM(sleep_h_imp ~ trend() + season()),
tslm_notr = TSLM(sleep_h_imp ~ season()),
tslm_knot = TSLM(sleep_h_imp ~ trend(knots = as.Date(c("2023-07-01", "2025-07-01"))) + season()),
sarima = ARIMA(sleep_h_imp ~ PDQ(1, 0, 1)),
arima_four = ARIMA(sleep_h_imp ~ fourier(period = 7, K = 3) +
fourier(period = 365.25, K = 2) + PDQ(0, 0, 0)),
stl_ets = decomposition_model(
STL(sleep_h_imp ~ trend(window = 365) + season(period = 7)),
ETS(season_adjust ~ season("N")))
) |>
mutate(combo = (ets + tslm) / 2)
recent_fit <- train_52 |>
model(
ets_52w = ETS(sleep_h_imp),
tslm_52w = TSLM(sleep_h_imp ~ season()),
mean_52w = MEAN(sleep_h_imp)
)
bind_rows(forecast(full_fit, new_data = test),
forecast(recent_fit, new_data = test)) |>
accuracy(sleep_ts)
})## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `tslm_knot = (function (object, ...) ...`.
## Caused by warning:
## ! prediction from a rank-deficient fit may be misleadingtuning_summary <- tuning_results |>
group_by(.model) |>
summarise(RMSE = mean(RMSE), MAE = mean(MAE), .groups = "drop") |>
arrange(RMSE)
tuning_summary## # A tibble: 12 × 3
## .model RMSE MAE
## <chr> <dbl> <dbl>
## 1 combo 0.456 0.366
## 2 tslm_52w 0.456 0.363
## 3 ets_ANA 0.456 0.366
## 4 ets 0.457 0.368
## 5 tslm_notr 0.457 0.365
## 6 tslm 0.459 0.368
## 7 ets_52w 0.459 0.369
## 8 stl_ets 0.468 0.375
## 9 sarima 0.475 0.382
## 10 mean_52w 0.476 0.382
## 11 tslm_knot 0.488 0.391
## 12 arima_four 0.498 0.396This compares model variants on demo data. The best row is only the lowest average error for these simulated records and windows; it does not establish a best model for real sleep data.
8. Refit ETS and forecast 12 weeks
Refit ETS using all available observations, then forecast 12 weeks ahead.
final_fit <- sleep_ts |> model(ets = ETS(sleep_h_imp))
final_fc <- final_fit |> forecast(h = "12 weeks")
final_fc## # A fable: 84 x 4 [1D]
## # Key: .model [1]
## .model date
## <chr> <date>
## 1 ets 2026-09-29
## 2 ets 2026-09-30
## 3 ets 2026-10-01
## 4 ets 2026-10-02
## 5 ets 2026-10-03
## 6 ets 2026-10-04
## 7 ets 2026-10-05
## 8 ets 2026-10-06
## 9 ets 2026-10-07
## 10 ets 2026-10-08
## # ℹ 74 more rows
## # ℹ 2 more variables: sleep_h_imp <dist>, .mean <dbl>final_fc |>
autoplot(sleep_ts |> filter_index("2026-04-01" ~ .)) +
labs(title = "Demo 12-week sleep forecast from ETS",
x = "Date", y = "Sleep duration (hours)")
Summarise the 12-week forecast by weekday. These are averages of forecasts across future dates for each weekday, not targets or recommendations.
to_hm <- function(h) {
mins <- round(h * 60)
sprintf("%dh %02dm", mins %/% 60, mins %% 60)
}
n_days <- 84
weekday_tbl <- final_fit |>
forecast(new_data = new_data(sleep_ts, n = n_days)) |>
hilo(level = 80) |>
unpack_hilo(`80%`) |>
as_tibble() |>
mutate(wake_day = wday(date, label = TRUE, abbr = FALSE, week_start = 1)) |>
group_by(wake_day) |>
summarise(
nights = n(),
forecast_h = mean(.mean),
lower_80 = mean(`80%_lower`),
upper_80 = mean(`80%_upper`),
.groups = "drop"
) |>
mutate(
Forecast = to_hm(forecast_h),
`80% expected range` = paste(to_hm(lower_80), "to", to_hm(upper_80))
)
weekday_tbl |> select(`Wake day` = wake_day, Forecast, `80% expected range`) |> knitr::kable()| Wake day | Forecast | 80% expected range |
|---|---|---|
| Monday | 7h 12m | 6h 36m to 7h 48m |
| Tuesday | 7h 11m | 6h 35m to 7h 47m |
| Wednesday | 7h 10m | 6h 34m to 7h 46m |
| Thursday | 7h 03m | 6h 27m to 7h 40m |
| Friday | 7h 00m | 6h 24m to 7h 36m |
| Saturday | 6h 50m | 6h 14m to 7h 27m |
| Sunday | 6h 56m | 6h 20m to 7h 32m |
The shaded interval and table express model uncertainty. They are not guarantees about future sleep. Results describe generated values only.
Optional: download and analyse your own Garmin sleep data
How Python fits with this R tutorial: The data analysis and forecasts above use R. Garmin’s
garminconnectpackage uses Python, so this section is a separate Python script, not an R code chunk. Run it from Terminal or PowerShell on your own computer. It saves your sleep data to a local CSV, which you can then load into a private copy of this R Markdown file. The published tutorial continues to use simulated demo data only.
If you want to try the workflow with your own records, run the
following Python script on your own computer. It uses
the community-maintained garminconnect Python package to
request daily sleep summaries and save a CSV. This is an unofficial
wrapper for Garmin Connect; Garmin may change its service or the package
response format. Current package instructions require Python 3.12 or
later (package
guide).
Create a file named download_garmin_sleep.py, paste code
below into it, and run it in a terminal. Change the example start and
end dates to the period you want. Password prompt hides what you type.
Package stores sign-in tokens locally under
~/.garminconnect; protect that folder and never upload it
or your exported CSV. Do not put your Garmin password in the script or
in a published document.
from datetime import date, timedelta
from getpass import getpass
from pathlib import Path
import csv
import time
from garminconnect import Garmin
start_date = date(2025, 1, 1) # Change these dates on your computer
end_date = date(2025, 3, 31)
email = input("Garmin email: ").strip()
password = getpass("Garmin password (hidden): ")
client = Garmin(
email=email,
password=password,
prompt_mfa=lambda: getpass("Garmin verification code: "),
)
client.login(str(Path("~/.garminconnect").expanduser()))
password = None # Do not keep plaintext password in use longer than needed
columns = [
"date", "sleep_s", "deep_s", "light_s", "rem_s", "awake_s",
"sleep_score", "hrv_night_avg",
]
rows = []
number_of_days = (end_date - start_date).days + 1
for day_offset in range(number_of_days):
day = start_date + timedelta(days=day_offset)
response = client.get_sleep_data(day.isoformat()) or {}
sleep = response.get("dailySleepDTO") or {}
if sleep:
overall_score = (sleep.get("sleepScores") or {}).get("overall") or {}
rows.append({
"date": sleep.get("calendarDate", day.isoformat()),
"sleep_s": sleep.get("sleepTimeSeconds"),
"deep_s": sleep.get("deepSleepSeconds"),
"light_s": sleep.get("lightSleepSeconds"),
"rem_s": sleep.get("remSleepSeconds"),
"awake_s": sleep.get("awakeSleepSeconds"),
"sleep_score": overall_score.get("value"),
"hrv_night_avg": sleep.get("avgSleepHRV"),
})
time.sleep(0.5) # Small pause between requests
output_path = Path("garmin_sleep_local.csv")
with output_path.open("w", newline="", encoding="utf-8") as csv_file:
writer = csv.DictWriter(csv_file, fieldnames=columns)
writer.writeheader()
writer.writerows(rows)
print(f"Saved {len(rows)} nights to {output_path.resolve()}")Create a Python virtual environment in the folder containing the script. This keeps its packages separate from other Python projects. Open Terminal (macOS/Linux) or PowerShell (Windows), change to that folder, then run commands for your system:
# macOS / Linux
python3 -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade garminconnect curl_cffi
python3 download_garmin_sleep.py# Windows PowerShell
py -m venv .venv
.venv\Scripts\Activate.ps1
py -m pip install --upgrade garminconnect curl_cffi
py download_garmin_sleep.pyAt first login, type Garmin password and any verification code into
the local prompts. The package sends them to Garmin over HTTPS and saves
reusable sign-in tokens locally. Treat ~/.garminconnect
like a password; never publish or share it. Do not enter credentials in
RPubs.
To analyse your own Garmin data locally, make a separate private copy of this R Markdown file and replace the demo-data generation chunk with the import line below. Then run the cleaning and modelling chunks as written.
sleep_raw <- readr::read_csv("garmin_sleep_local.csv", show_col_types = FALSE)Keep that personal copy and its knitted HTML private. If you publish HTML after replacing demo data, plots or printed tables may reveal your sleep records. Publish the unchanged demo version only.
What to take away
- The workflow covers cleaning, exploration, model comparison, rolling validation, tuning, and forecasting.
- Generated demo data keep personal Garmin records out of this publication; numerical results depend on the data.
- Compare candidate models on the same held-out dates and across several rolling windows.
- Treat a personal sleep forecast as a planning aid, not a precise prediction or medical conclusion.
Reproduce this resource
- Install R and RStudio.
- Create an R Markdown document and paste this page into it.
- Install the packages listed above once using the Console.
- Run chunks from top to bottom, then click Knit. Model tuning may take several minutes.
- In the HTML preview, click Publish and choose RPubs.