Please note that the description in this document applies only to effects of the double-nested type.
(Assuming the user has long-format independent variables and short-format dependent variables.)
The model supports the use of a time index in the following scenario: when the response variable has fewer observations than the explanatory variables, a time index column can be added to the input data to indicate the recording time of the response variable.
When using this feature, name the time index column
times and merge it with the input matrix. The model will
automatically extract time points where both independent and response
variables are present, using these to obtain the fitted parameters.
However, the objects of operations such as centralization and exp-smooth
remain all independent variables.
The following methods are provided for support:
Data Preparation:
# The code below is not working; it is for demonstration purposes only.
res <- build_flat_from_long(y, time_seq, X_ori, w_ori)Predict all response variables:
# The code below is not working; it is for demonstration purposes only.
pre_all <- gamFactory:::pre_all_times(object = out, X_ori = X_ori,w_ori = w_ori)Important Notes:
The data used for prediction and the data used for training should share the same initial data to meet the requirements of exp-smooth;
When calling gam_nl, both n_si and n_nexp must be
explicitly specified;
Consider the daily electricity demand of a given region over 365 days.
For each day, temperature observations are collected at three time points—morning, noon, and evening—across four subregions (East, South, West, and North), together with their corresponding sampling times.
The aim is to predict the total daily electricity demand.
In this case, the response variable (electricity demand) has one-third the length of the predictor variables (temperature observations). To match the model input format, the data are preprocessed as follows:
Inf (or
NA).time_diff (length = 365×3). The first
missing value is replaced by the mode of the column.time_diff. The expanded matrix thus has a shape of
365×2×3.times (length = 365). For example, if electricity demand is
recorded on the evening of day 1 and at noon on day 2, then the first
two elements of times are 3 and 5.Generate data:
# Before running the code here, please run the supporting code at the bottom first.
library(gamFactory)
sim <- simulate_elec_demand()Example of raw temperature data:
tail(sim$X_temp)Flattened temperature matrix:
tail(sim$X_flat,3)Final input matrix (including response and time features):
str(sim$dat)A nonlinear model of type si_nexp is employed for
estimation:
trans <- trans_linear_nexpsm(n_si = 4, n_nexp = 2)
out <- gam_nl( formula = list( y ~ s_nest(xw, trans = trans, k = 3, m = c(4, 2)),
~ 1),
data = sim$dat,
family = fam_gaussian(),
sp = 0.01,
control = list(trace = TRUE) )Inside the model, predictors are internally reshaped back into the long format, centered, and then processed through the exponential smoothing layer.
pre_all <- gamFactory:::pre_all_times(object = out, X_ori = sim$X_ori, w_ori = sim$w_ori)Plot:
# ---------------- plot：only first N_day ----------------
N_day <- 5 # only plot the first N days
n_day <- 365
n_time <- 3
N_day <- max(1L, min(N_day, n_day))
idxN <- 1:(N_day * n_time)
# prediction(black line + point)
plot(idxN, pre_all[idxN],
type = "o", pch = 16, col = "black", cex = 0.6, lwd = 1.2,
xlab = sprintf("Time index (first %d days)", N_day),
ylab = "Electricity demand",
main = sprintf("Predicted vs Observed — first %d days", N_day))
# true data (blue triangle)
mask_obsN <- sim$time_seq <= max(idxN)
points(sim$time_seq[mask_obsN], sim$dat$y[mask_obsN],
pch = 17, col = adjustcolor("blue", alpha.f = 0.5), cex = 1)
legend("topleft",
legend = c("Prediction (pre_all)", "Observed (dat$y)"),
col = c("black", "blue"), pch = c(16,17),
lty = c(1, NA), bty = "n")This effect can also converge without exponential smoothing when
using only linear or sinusoidal transformations. In such cases, the
intercept of alpha_nexp tends toward -Inf, causing the
smoothing weight ω → 0. For example:
With linear transformation only, the intercept of
alpha_nexp ≈ -77;
With linear + sinusoidal transformation, the intercept ≈ -22.
In both cases, ω ≈ 0, meaning the smoothing layer is effectively disabled.
# ================================================================
# Simulated electricity demand forecasting — testing the `times` parameter
# ================================================================
# Data structure:
# Explanatory variable matrix X_temp: 365×3 rows ("day × morning/noon/evening") × 4 columns (East/South/West/North);
# Response variable y_demand: length 365;
# Vector time_point: indicates the time of day for each daily observation, length 365, values in {1=morning, 2=noon, 3=evening};
# This script simulates electricity demand data and tests the handling of the `times` parameter (generation + processing + verification).
# ======================Example of use============================
# sim <- simulate_elec_demand()
#
# # useful outputs
# dat <- sim$dat
# X_temp <- sim$X_temp
# X_ori <- sim$X_ori
# X_flat <- sim$X_flat
# w <- sim$w
# w_ori <- sim$w_ori
# w_flat <- sim$w_flat
# L_temp <- sim$L_temp
# y_true <- sim$y_true
# y_demand <- sim$y_demand
# tail(dat,10)
# ==================================================================
simulate_elec_demand <- function(
n_day = 365,
n_time = 3, # 1=morning, 2=noon, 3=evening
regions = c("East","South","West","North"), # regions names
seed = 20251002,
# --- param for temperature generation ---
mu_r = c(East=13, South=15, West=12, North=10), # Regional mean (°C)
A_r = c(East=9, South=10, West=8, North=7), # Seasonal amplitude
phi_r = c(East=30, South=25, West=40, North=50), # Phase (days)
d_time = c(-1.5, 0, 1.0), # Diurnal offset: early/mid/late
sd_eps = 1.2, # Temperature noise standard deviation
# --- Morning/Afternoon/Evening Probability ---
prob_time = c(0.32, 0.38, 0.30),
# --- linear and expsmooth parameter ---
lin_coef = c(East = 1, South = 2, West = 3, North = 1),
beta = c(-1, 0.5), # expsmooth param: intercept, time_diff
sd_eta = 0.01, # noise for y_true
mark_inf_after_tlast = TRUE
){
stopifnot(length(d_time) == n_time)
set.seed(seed)
n_reg <- length(regions)
# 1) day-time index
idx <- expand.grid(time = 1:n_time, day = 1:n_day)
idx <- idx[order(idx$day, idx$time), ]
rownames(idx) <- NULL
# 2) tempreature matrix X_temp: (n_day*n_time) × n_reg
# T_{d,t,r} = mu_r + A_r * cos( 2π*(day - phi_r)/n_day ) + d_time[t] + ε
mu_r <- setNames(mu_r[regions], regions)
A_r <- setNames(A_r[regions], regions)
phi_r <- setNames(phi_r[regions], regions)
make_region_temp <- function(mu, A, phi) {
cos_season <- A * cos(2*pi*(idx$day - phi)/n_day)
mu + cos_season + d_time[idx$time] + rnorm(nrow(idx), 0, sd_eps)
}
X_temp <- do.call(
cbind,
lapply(regions, function(rr) make_region_temp(mu_r[[rr]], A_r[[rr]], phi_r[[rr]]))
)
colnames(X_temp) <- regions
stopifnot(nrow(X_temp) == n_day * n_time, ncol(X_temp) == n_reg)
rownames(X_temp) <- paste0("day", idx$day, "_t", idx$time)
X_ori <- X_temp # save for output
# 3) time point and time_seq
time_point <- sample(1:n_time, size = n_day, replace = TRUE, prob = prob_time)
time_seq <- (1:n_day - 1L) * n_time + time_point
# 4) put inf after t_last (optinal)
t_last <- time_seq[n_day]
if (mark_inf_after_tlast) {
bad <- (t_last + 1L):(n_day * n_time)
if (length(bad) > 0) X_temp[bad, ] <- Inf
}
# 5) linear transformation：L_temp = X_temp %*% lin_coef
lin_coef <- setNames(lin_coef[regions], regions)
L_temp <- X_temp %*% lin_coef
if (mark_inf_after_tlast) {
bad <- (t_last + 1L):(n_day * n_time)
if (length(bad) > 0) L_temp[bad, ] <- Inf
}
# 6) X_flat: n_day × (n_time*n_reg)
X_flat <- matrix(NA_real_, nrow = n_day, ncol = n_time * n_reg)
colnames(X_flat) <- as.vector(sapply(1:n_time, function(tt) paste0(regions, "_t", tt)))
rownames(X_flat) <- paste0("day", 1:n_day)
for (d in 1:n_day) {
rows_d <- which(idx$day == d)
rows_d <- rows_d[order(idx$time[rows_d])] # t1,t2,t3
block_3xR <- X_temp[rows_d, regions, drop = FALSE]
X_flat[d, ] <- as.vector(t(block_3xR)) # c(t1, t2, t3)
}
# 7) time_diff
time_diff <- rep(8, nrow(X_temp)) + rnorm(n_day * n_time, 0, 2)
# 8) w and w_flat
intercept <- rep(1, nrow(X_temp))
w <- cbind(intercept = intercept, time_diff = time_diff)
w_ori <- w
if (mark_inf_after_tlast) {
bad <- (t_last + 1L):(n_day * n_time)
if (length(bad) > 0) w[bad, ] <- Inf
}
w_flat <- matrix(NA_real_, nrow = n_day, ncol = n_time * ncol(w))
colnames(w_flat) <- as.vector(sapply(1:n_time, function(tt) paste0(colnames(w), "_t", tt)))
rownames(w_flat) <- paste0("day", 1:n_day)
for (d in 1:n_day) {
rows_d <- which(idx$day == d)
rows_d <- rows_d[order(idx$time[rows_d])]
block_3x2 <- w[rows_d, , drop = FALSE]
w_flat[d, ] <- as.vector(t(block_3x2))
}
# 9) y_true and y_demand
y_true <- gamFactory:::expsmooth(y = L_temp, Xi = w, beta = beta,
times = time_seq, deriv = 0)$d0
y_demand <- as.numeric(y_true + rnorm(n_day, 0, sd_eta))
# 10) combine dat = [y, xw]
xw <- cbind(X_flat, w_flat, times = time_seq)
dat <- data.frame(y = y_demand, xw = I(xw))
# return
list(
dat = dat, # data.frame(y, xw=I(matrix))
X_temp = X_temp, # may contain Inf
X_ori = X_ori, # without Inf
X_flat = X_flat,
w = w, # may contain Inf
w_ori = w_ori, # without Inf
w_flat = w_flat,
L_temp = L_temp, # may contain Inf
y_true = y_true,
y_demand = y_demand,
idx = idx,
time_point = time_point,
time_seq = time_seq,
t_last = t_last,
params = list(
seed = seed, mu_r = mu_r, A_r = A_r, phi_r = phi_r, d_time = d_time,
sd_eps = sd_eps, prob_time = prob_time, lin_coef = lin_coef,
beta = beta, sd_eta = sd_eta, mark_inf_after_tlast = mark_inf_after_tlast
)
)
}# # ===============================================================
# # # Test case:
# set.seed(1)
# n_day <- 365
# y <- rnorm(n_day, 100, 5)
# n_time <- 3
#
# # Step 1: Randomly assign each day a time_point ∈ {1, 2, 3}
# time_point <- sample(1:n_time, n_day, replace = TRUE)
#
# # Step 2: Compute time_seq
# time_seq <- (seq_len(n_day) - 1L) * n_time + time_point
#
# # Construct: long-format matrices with 1093 rows (2 rows fewer than 365*3 = 1095)
# n_row_long <- 1093
#
# X_ori <- matrix(rnorm(n_row_long * 4, 15, 5), n_row_long, 4)
# colnames(X_ori) <- c("East", "South", "West", "North")
#
# w_ori <- cbind(
# intercept = rep(1, n_row_long),
# time_diff = rnorm(n_row_long, 8, 2)
# )
#
# # Apply the function
# res <- build_flat_from_long(y, time_seq, X_ori, w_ori)
#
# # Check results
# res$n_rep # should be 3
# res$t_last # tail(time_seq, 1)
# dim(res$X_flat) # 365 x (3*4) = 365 x 12
# dim(res$w_flat) # 365 x (3*2) = 365 x 6
#
# tail(res$X_flat, 2)
# tail(res$w_flat, 2)
# # ===============================================================
# ===============================================================
# Build flat matrices (per-day slicing) from long-format X_ori / w_ori
# ===============================================================
# Inputs:
# y : numeric vector, length = n_day
# time_seq : integer vector, length = n_day, last element is t_last
# X_ori : numeric matrix/data.frame, long format ordered by day:
# day1: t1..t_nrep, day2: t1..t_nrep, ..., dayD: t1..t_nrep
# w_ori : same long-format ordering as X_ori
#
# Outputs:
# $dat : data.frame(y = y, xw = I(cbind(X_flat, w_flat, times = time_seq)))
# $X_flat : n_day x (n_rep * ncol(X_ori))
# $w_flat : n_day x (n_rep * ncol(w_ori))
# $t_last : last element of time_seq
# $n_rep : repetition count (ceiling)
# ===============================================================
build_flat_from_long <- function(y, time_seq, X_ori, w_ori) {
# ---------- 0) Basic checks ----------
stopifnot(is.numeric(y), length(y) >= 1L)
stopifnot(is.numeric(time_seq), length(time_seq) == length(y))
n_day <- length(y)
t_last <- tail(time_seq, 1L)
# Coerce to numeric matrices
X_ori <- as.matrix(X_ori); storage.mode(X_ori) <- "double"
w_ori <- as.matrix(w_ori); storage.mode(w_ori) <- "double"
nx <- nrow(X_ori); px <- ncol(X_ori)
nw <- nrow(w_ori); pw <- ncol(w_ori)
if (is.null(colnames(X_ori))) colnames(X_ori) <- paste0("X", seq_len(px))
if (is.null(colnames(w_ori))) colnames(w_ori) <- paste0("w", seq_len(pw))
# ---------- 1) Constraint on t_last ----------
max_rows_long <- max(nx, nw)
if (t_last > max_rows_long) {
stop(sprintf("t_last (%d) > max(nrow(X_ori), nrow(w_ori)) = %d. Please check time_seq or long-table lengths.",
t_last, max_rows_long))
}
# ---------- 2) Decide n_rep & pad the tail with Inf to full daily slices ----------
n_rep <- ceiling(max_rows_long / n_day)
total_needed <- n_day * n_rep
pad_matrix_to <- function(M, total_needed) {
n_now <- nrow(M); p_now <- ncol(M)
if (n_now < total_needed) {
n_pad <- total_needed - n_now
M <- rbind(M, matrix(Inf, nrow = n_pad, ncol = p_now,
dimnames = list(NULL, colnames(M))))
} else if (n_now > total_needed) {
# If longer, truncate the tail to n_day * n_rep
M <- M[seq_len(total_needed), , drop = FALSE]
warning("Input has more rows than needed; truncated to n_day * n_rep.")
}
M
}
X_pad <- pad_matrix_to(X_ori, total_needed)
w_pad <- pad_matrix_to(w_ori, total_needed)
# ---------- 3) Per-day slicing & flatten (t1..t_nrep in order) ----------
# Assumes row order = day1: t1..t_nrep, day2: t1..t_nrep, ...
make_flat_per_day <- function(M, n_day, n_rep, base_names) {
p <- ncol(M)
flat <- matrix(NA_real_, nrow = n_day, ncol = p * n_rep)
colnames(flat) <- as.vector(sapply(seq_len(n_rep), function(tt) paste0(base_names, "_t", tt)))
rownames(flat) <- paste0("day", seq_len(n_day))
for (d in seq_len(n_day)) {
rows_d <- ((d - 1L) * n_rep + 1L):(d * n_rep) # t1..t_nrep of current day
# Concatenate column blocks t1, t2, ...
flat[d, ] <- unlist(lapply(rows_d, function(rr) M[rr, , drop = FALSE]))
}
flat
}
X_flat <- make_flat_per_day(X_pad, n_day, n_rep, colnames(X_ori))
w_flat <- make_flat_per_day(w_pad, n_day, n_rep, colnames(w_ori))
# ---------- 4) Assemble xw & dat ----------
xw <- cbind(X_flat, w_flat, times = time_seq)
dat <- data.frame(y = as.numeric(y), xw = I(xw))
# ---------- 5) Return ----------
list(
dat = dat,
X_flat = X_flat,
w_flat = w_flat,
t_last = t_last,
n_rep = n_rep
)
}