In xliu-stat/NRRR: Multivariate Functional Regression via Nested Reduced-Rank Regularization
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)
Introduction: What is Nested Reduced-Rank Regression (NRRR)?
The nested reduced-rank regression (NRRR) is designed to fit a multivariate functional linear regression
$$
{\bf y}(t) = \int_{\mathcal{S}} {\bf C}0(s,t){\bf x}(s)ds + \epsilon(t),\ t\in\mathcal{T}
$$
where ${\bf y}(t)=(y_1(t),\ldots,y_d(t))^T \in \mathbb{R}^d$ with $t\in\mathcal{T}$ and ${\bf x}(s)=(x_1(s),\ldots,x_p(s))^T \in \mathbb{R}^p$ with $s\in\mathcal{S}$ are multivariate functional responses and predictors, respectively. $\epsilon(t)\in \mathbb{R}^d$ is a zero-mean random error function. The functional coefficient matrix ${\bf C}_0(s,t)=[c{k,l}(s,t)]{d\times p}$ consists of unknown bivariate functions $c{k,l}(s,t),\ k=1,\ldots,d, l=1,\ldots,p$, and NRRR aims to jointly estimate these functional surfaces by utilizing the potential associations among the functional variables.
First, a global reduced-rank structure assumes $$ {\bf C}0(s,t) = {\bf U}_0 {\bf C}^_0(s,t) {\bf V}_0^T,$$ where ${\bf U}_0 \in \mathbb{R}^{d\times r_y}$ with $r_y \leq d$ and ${\bf V}_0 \in \mathbb{R}^{p\times r_x}$ with $r_x \leq p$ provide the weights to form latent functional responses ${\bf y}^(t)={\bf U}_0^T {\bf y}(t)\in \mathbb{R}^{r_y}$ and latent functional predictors ${\bf x}^(s)={\bf V}_0^T {\bf x}(s)\in \mathbb{R}^{r_x}$, respectively. Therefore, the model is reduced to $${\bf y}^(t) = \int{\mathcal{S}} {\bf C}^_0(s,t){\bf x}^(s)ds + \epsilon^*(t),\ t\in\mathcal{T},$$ and a global dimension reduction is achieved once we have $r_y < d$ or $r_x < p$. In this way, ${\bf U}_0$ and ${\bf V}_0$ capture the association among the elements of ${\bf y}(t)$ and the association among the elements of ${\bf x}(s)$, respectively. This structure is particularly useful for simultaneously modeling a large number of functional responses and predictors that are highly correlated across $s$ or $t$.
Next, we deal with the latent regression surface ${\bf C}^_0(s,t)$. A basis expansion and truncation approach is first applied to ensure the smoothness and convert the infinite-dimensional problem to be finite-dimensional, i.e.,
$${\bf C}^0(s,t)\approx ({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf C}0^ ({\bf I}_{r_x}\otimes {\bf \Phi}(s)),\ {\bf C}_0^ \in \mathbb{R}^{J_yr_y\times J_xr_x},$$
where ${\bf \Psi}(t) = (\psi_1(t),\ldots,\psi{J_y}(t))$ consists $J_y$ basis functions and ${\bf \Phi}(s) = (\phi_1(s),\ldots,\phi_{J_x}(s))$ consists $J_x$ basis functions. Then, we impose the local reduced-rank structure on ${\bf C}_0^$ as $rank({\bf C}_0^)\leq r$ for $r \leq \min(J_xr_x, J_yr_y)$ and write it as a full-rank decomposition ${\bf C}_0^ = {\bf A}_0^{\bf B}_0^{T}$ for some ${\bf A}_0^ \in \mathbb{R}^{J_yr_y \times r}$ and ${\bf B}_0^* \in \mathbb{R}^{J_xr_x \times r}$. This structure induces the dependency between the latent responses and the latent predictors through their basis-expansion.
With these two layers of dimension reduction, the complexity of the model is greatly reduced. The basis functions used here are treated as given, such as spline, wavelet, and Fourier basis. In the package NRRR, the B-spline basis is mainly used. Also, we assume all the responses or the predictors share the same set of basis functions for simplicity.
Based on the selected basis functions, we obtain the integrated predictors and responses from
$${\bf x}=\int_{\mathcal{S}} ({\bf I}{p}\otimes {\bf \Phi}(s)){\bf x}(s)ds \in \mathbb{R}^{J_x p}, $$
and
$${\bf y}=({\bf I}_d \otimes {\bf J}{\psi\psi}^{-1/2})\int_{\mathcal{T}} ({\bf I}{d}\otimes {\bf \Psi}(t)){\bf y}(t)dt \in \mathbb{R}^{J_y d},$$
where ${\bf J}{\psi\psi}=\int_{\mathcal{T}}{\bf \Psi}(t){\bf \Psi}(t)^Tdt$ is a positive definite matrix. Thus, with all the aforementioned structures and quantities, the estimation criterion that minimizes the mean integrated squared error with respect to ${\bf C}(s,t)$, i.e.,
$$ \frac{1}{n} \sum_{i=1}^n \int_{\mathcal{T}} \|{\bf y}i(t) - \int{\mathcal{S}} {\bf C}(s,t){\bf x}i(s)ds \|^2 dt $$
can be written as
$$ \min{{\bf U, V, A^, B^}} { \frac{1}{n}\sum_{i=1}^n \| {\bf y}i - ({\bf I}_d\otimes {\bf J}{\psi\psi}^{1/2})({\bf U}\otimes {\bf I}{J_y}){\bf A^B^}^T ({\bf V}^T\otimes {\bf I}{J_x}){\bf x}i \|^2 }. $$ By properly rearranging the columns and rows of the data matrices and coefficient matrices, the problem finally boils down to
$$\min{{\bf C}} \| {\bf Y - XC} \|F^2, \ s.t.\ {\bf C}= ({\bf I}{J_x} \otimes {\bf V}) {\bf BA}^T ({\bf I}{J_y} \otimes {\bf U}^T),$$
where ${\bf Y} = ({\bf Y}{\cdot1}, \ldots, {\bf Y}{\cdot J_y})$ with ${\bf Y}{\cdot j} = (y_{ikj}){n\times d}$ for $j=1, \dots,J_y$ and ${\bf X} = ({\bf X}{\cdot1}, \ldots, {\bf X}{\cdot J_x})$ with ${\bf X}{\cdot j} = (x_{ilj}){n\times p}$ for $j=1, \dots,J_x$. That is, each ${\bf Y}{\cdot j}$ contains all the integrated values of the functional responses obtained based on $\psi_j(t)$ and ${\bf X}_{\cdot j}$ contains all the integrated values of the functional predictors obtained based on $\phi_j(s)$. This matrix approximation representation provides a clear illustration of the nested reduced-rank structure, i.e., ${\bf U}$ and ${\bf V}$ are designed to capture the shared column and row spaces among the blockwise sub-matrices of ${\bf C}$, which, as a whole, is also of low-rank. The applicability of this structure goes beyond the functional setup. For example, it can be applied in the high-dimensional vector autoregressive modeling in multivariate time series analysis, surveillance video processing and also the tensor-on-tensor regression.
This optimization problem is non-convex and has no explicit solution. We proposed a blockwise coordinate descent algorithm to obtain a local solution, and this algorithm is implemented in the package NRRR. For more details of NRRR, please read Liu, X., Ma, S., & Chen, K. (2020). Multivariate Functional Regression via Nested Reduced-Rank Regularization. arXiv: Methodology.
Next, we use the Adelaide electricity demand data as an example to show the usage of the package NRRR to investigate the functional association between weekly electricity demand trajectory and temperature trajectory. In the package NRRR, there is a function NRRR.func that takes functional observations as input and uses B-spline basis to conduct basis expansion and then fits a nested reduced-rank regression model. In the following, to provide more details about how to deal with discrete functional observations, we have not used NRRR.func. In practical application problems, the following introduction is helpful when other kinds of basis functions are preferred to conduct basis expansion.
Application: Adelaide Electricity Demand analysis
Adelaide is the capital city of the state of South Australia. In the summertime, the cooling in Adelaide mainly depends on air conditioning, which makes the electricity demand highly dependent on the weather conditions, and large volatility in temperature throughout the day could make stable electricity supply challenging. Therefore, for facilitating the supply management of electricity, it is important to understand the dependence and predictive association between the electricity demand and the temperature. Here we apply NRRR to perform a multivariate functional regression analysis between daily half-hour electricity demand profiles for the 7 days of a week and the corresponding temperature profiles for the 7 days of the same week.
Half-hourly temperature records at two locations, Adelaide Kent town and Adelaide airport, are available between 7/6/1997 and 3/31/2007. Also available are the half-hourly electricity demand records of Adelaide for the same period. The data is extracted from the R package fds.
library(fds)
# Electricity Demand Data
data(mondaydemand)
data(tuesdaydemand)
data(wednesdaydemand)
data(thursdaydemand)
data(fridaydemand)
data(saturdaydemand)
data(sundaydemand)
# Temperature at Kent town
data(mondaytempkent)
data(tuesdaytempkent)
data(wednesdaytempkent)
data(thursdaytempkent)
data(fridaytempkent)
data(saturdaytempkent)
data(sundaytempkent)
# Temperature at Adelaide airport
data(mondaytempairport)
data(tuesdaytempairport)
data(wednesdaytempairport)
data(thursdaytempairport)
data(fridaytempairport)
data(saturdaytempairport)
data(sundaytempairport)
As such, for each day during the period, there are three observed functional curves, each with 48 half-hourly observations. Here we plot the temperature and electricity demand profiles of all the Mondays from 7/6/1997 to 3/31/2007 as an illustration.
time_index <- c("00:00","00:30",
"1:00","1:30",
"2:00","2:30",
"3:00","3:30",
"4:00","4:30",
"5:00","5:30",
"6:00","6:30",
"7:00","7:30",
"8:00","8:30",
"9:00","9:30",
"10:00","10:30",
"11:00","11:30",
"12:00","12:30",
"13:00","13:30",
"14:00","14:30",
"15:00","15:30",
"16:00","16:30",
"17:00","17:30",
"18:00","18:30",
"19:00","19:30",
"20:00","20:30",
"21:00","21:30",
"22:00","22:30",
"23:00","23:30")
f1 <- rep(1:48,508)
prg1 <- vector()
prg2 <- vector()
prg3 <- vector()
f2 <- vector()
for(a in 1:508){
prg1 <- c(prg1,mondaydemand$y[,a])
prg2 <- c(prg2,mondaytempkent$y[,a])
prg3 <- c(prg3,mondaytempairport$y[,a])
f2 <- c(f2,rep(a, 48))
}
prgdata <- data.frame(b=f1,
day=f2,
PRG1=prg1, PRG2=prg2, PRG3=prg3)
prgdata$b <- factor(prgdata$b)
levels(prgdata$b) <- time_index
library(ggplot2)
# Electricity demand on Monday
ggplot(prgdata, aes(y = PRG1, x = b, group = day, colour = day)) +
geom_path(show.legend = FALSE) + #geom_tile() +
scale_x_discrete(breaks = time_index[seq(1,48,2)]) +
ggtitle("Electricity demand on Mondays") +
xlab("Time") +
ylab("Electricity demand") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(axis.text=element_text(size=8),axis.title=element_text(size=15),
plot.title = element_text(hjust = 0.5, size = 15))
# Temperature of Kent on Monday
ggplot(prgdata, aes(y = PRG2, x = b, group = day, colour = day)) +
geom_path(show.legend = FALSE) + #geom_tile() +
scale_x_discrete(breaks = time_index[seq(1,48,2)]) +
ggtitle("Temperature of Kent on Mondays") +
xlab("Time") +
ylab("Temperature") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(axis.text=element_text(size=8),axis.title=element_text(size=15),
plot.title = element_text(hjust = 0.5, size = 15))
# Temperature of Airport on Monday
ggplot(prgdata, aes(y = PRG3, x = b, group = day, colour = day)) +
geom_path(show.legend = FALSE) + #geom_tile() +
scale_x_discrete(breaks = time_index[seq(1,48,2)]) +
ggtitle("Temperature of Adelaide airport on Mondays") +
xlab("Time") +
ylab("Temperature") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(axis.text=element_text(size=8),axis.title=element_text(size=15),
plot.title = element_text(hjust = 0.5, size = 15))
Since our primary focus is on studying the general association between the within-day demand and temperature trajectories in a week, we center the 48 discrete observations of each daily curve to remove the between-day trend and seasonality of the data. Each week is then treated as a replication.
# response
mon <- as.matrix(mondaydemand$y)
tue <- as.matrix(tuesdaydemand$y)
wed <- as.matrix(wednesdaydemand$y)
thu <- as.matrix(thursdaydemand$y)
fri <- as.matrix(fridaydemand$y)
sat <- as.matrix(saturdaydemand$y)
sun <- as.matrix(sundaydemand$y)
# model dimension
ns <- 48
nt <- ns
n <- 508
p <- 7
d <- 7
WhoWeek <- array(dim = c(p, ns, n), NA)
WhoWeek[1, , ] <- apply(mon, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[2, , ] <- apply(tue, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[3, , ] <- apply(wed, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[4, , ] <- apply(thu, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[5, , ] <- apply(fri, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[6, , ] <- apply(sat, 2, scale, center = TRUE, scale = FALSE)
WhoWeek[7, , ] <- apply(sun, 2, scale, center = TRUE, scale = FALSE)
oriY <- array(dim = c(n, p, ns), NA)
for (i in 1:n) {
for (j in 1:p) {
oriY[i, j, ] <- WhoWeek[j, , i]
}
}
# predictor
WhoWeekX <- array(dim = c(d, nt, n), NA)
WhoWeekX[1, , ] <- apply(as.matrix(mondaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[2, , ] <- apply(as.matrix(tuesdaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[3, , ] <- apply(as.matrix(wednesdaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[4, , ] <- apply(as.matrix(thursdaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[5, , ] <- apply(as.matrix(fridaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[6, , ] <- apply(as.matrix(saturdaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
WhoWeekX[7, , ] <- apply(as.matrix(sundaytempkent$y), 2, scale,
center = TRUE, scale = FALSE)
oriX <- array(dim = c(n, d, nt), NA)
for (i in 1:n) {
for (j in 1:d) {
oriX[i, j, ] <- WhoWeekX[j, , i]
}
}
X <- oriX
Y <- oriY
Till now, we've obtained the discrete observations of the response and the predictor trajectories. Next, let's use B-spline basis with 30 degrees of freedom to convert the discrete observations into its integrated form as we introduced before.
library(splines)
jx <- 30
jy <- jx
sseq <- c(1:48)
phi <- bs(c(0, sseq), df = jx)[-1, ]
tseq <- c(1:48)
psi <- bs(c(0, tseq), df = jy)[-1, ]
# get integrated predictors
Xa <- array(dim = c(n, p, jx), NA)
sdiff <- (sseq - c(0, sseq[-ns]))
for (i in 1:n) {
for (l in 1:p) {
for (j in 1:jx) {
# integrate X over s
Xa[i, l, j] <- sum(phi[, j] * X[i, l, ] * sdiff)
}
}
}
# compute Jpsi and Jpsi^{-1/2}
Jpsi <- matrix(nrow = jy, ncol = jy, 0)
tdiff <- (tseq - c(0, tseq[-nt]))
for (t in 1:nt) Jpsi <- Jpsi + psi[t, ] %*% t(psi[t, ]) * tdiff[t]
eJpsi <- eigen(Jpsi)
Jpsihalf <- eJpsi$vectors %*% diag(sqrt(eJpsi$values)) %*% t(eJpsi$vectors)
Jpsihalfinv <- eJpsi$vectors %*% diag(1 / sqrt(eJpsi$values)) %*% t(eJpsi$vectors)
# get integrated responses
Ya <- array(dim = c(n, d, jy), NA)
tdiff <- (tseq - c(0, tseq[-nt]))
psistar <- psi %*% Jpsihalfinv
for (k in 1:d) {
for (i in 1:n) {
for (j in 1:jy) {
# integrate Y over t
Ya[i, k, j] <- sum(psistar[, j] * Y[i, k, ] * tdiff)
}
}
}
# obtain matrices Y and X in matrix approximation form
Yest <- Ya[, , 1]
for (j in 2:jy) Yest <- cbind(Yest, Ya[, , j])
Xest <- Xa[, , 1]
for (j in 2:jx) Xest <- cbind(Xest, Xa[, , j])
Here we use the daily half-hour electricity demand as the functional multivariate response with $d = 7$ (corresponding to 7 days in a week from Monday to Sunday). As for the predictors, we use the half-hour temperature data from Kent as the multivariate functional predictors, so that $p = 7$. The total sample size is $n = 508$, equaling the number of weeks in the study period.
Estimation
The package NRRR is available at Github, and we can install it with
#install.packages("devtools")
#devtools::install_github("xliu-stat/NRRR")
library(NRRR)
set.seed(6)
As fitting the NRRR model with cross validation could be time-consuming, we set eval=FALSE in this code chunk and display some previously saved fitting results here.
# set 'eval = TRUE' to reproduce the results
fit.NRRR <- NRRR.cv(Yest, Xest,
nfold = 10, norder = NULL, Ag0 = NULL, Bg0 = NULL,
jx = 30, jy = 30, p = 7, d = 7, n = 508,
maxiter = 300, conv = 1e-4, method = "RRR"
)
load("fit_nRRR.rda")
# the estimated r
fit.NRRR$rank
# the estimated rx
fit.NRRR$rx
# the estimated ry
fit.NRRR$ry
# the estimated U
fit.NRRR$Ag
# the estimated V
fit.NRRR$Bg
The estimated rank values are $\hat r = 4$, $\hat r_x = 1$, and $\hat r_y = 5$. The estimated loading matrix for the predictors is ${\bf V} = (-0.22, -0.39, -0.46, -0.52, -0.43, -0.28, -0.25)^T$. This shows that there is only one latent functional predictor that is driving the patterns of the electronic demands, and this factor can be roughly explained as the averaged daily temperature profile of the week. It appears that the days closer to the middle of the week load higher. On the response side, there is not much global reduction, as the estimated loading matrix ${\bf U}$ is of rank 5. To make sense of ${\bf U}$, it may be more convenient to examine the two basis vectors of its orthogonal complement, i.e., the first two singular vectors of ${\bf I - \hat U\hat U^T}$, which give the latent response factors that are not related to the temperatures at all. While the first loading vector $(-0.52, 0.36, 0.28, 0.25, −0.56, 0.34, −0.18)^T$ is hard to interpret, the second loading vector $(0.00, −0.68, 0.73, -0.01, −0.04, 0.04, −0.04))^T$ clearly indicates that the difference between the electronic demand profiles of Tuesday and Wednesday is mostly a noise process. In other words, the demand profiles of these two days are related to the temperature process in almost the same way.
Visualization of regression surface
The function NRRR.plot.reg is designed to plot the heatmap for each bivariate function in the functional coefficient matrix, e.g., $c_{k,l}(s,t),\ 1\leq k \leq d, 1 \leq l \leq p$ in ${\bf C}(s,t)$. Based on the fitting results from the nested reduced-rank regression, different kinds of regression surfaces can be visualized to give a clear illustration of the functional correlation between a user-specified predictor (or latent predictor) trajectory and response (or latent response) trajectory. With the nested reduced-rank structure, the functional regression becomes
$${\bf y}(t) \approx \int_{\mathcal{S}} {\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T {\bf x}(s)ds + {\bf \epsilon}(t),$$
based on which, we can visualize the bivariate function in:
-
${\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T \in \mathbb{R}^{d\times p}$ to see the effects of $x_i(s)$ on $y_j(t)$ with $1\leq i \leq p,\ 1\leq j \leq d$ by specifying method = "original", x_ind = i and y_ind = j.
-
${\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)) \in \mathbb{R}^{d\times r_x}$ to see the effects of the $i$-th element of ${\bf V}^T{\bf x}(s)$ on $y_j(t)$ with $1\leq i \leq r_x,\ 1\leq j \leq d$ by specifying method = "y_original", x_ind = i and y_ind = j.
-
$({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T \in \mathbb{R}^{r_y \times p}$ to see the effects of $x_i(s)$ on the $j$-th element of ${\bf U}^T{\bf y}(t)$ with $1\leq i \leq p,\ 1\leq j \leq r_y$ by specifying method = "x_original", x_ind = i and y_ind = j.
-
$({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)) \in \mathbb{R}^{r_y \times r_x}$ to see the effects of the $i$-th element of ${\bf V}^T{\bf x}(s)$ on the $j$-th element of ${\bf U}^T{\bf y}(t)$ with $1\leq i \leq r_x,\ 1\leq j \leq r_y$ by specifying method = "latent", x_ind = i and y_ind = j.
Based on the fitting results we obtained above, we've got only one latent predictor and five latent responses. Thus, next let's visualize the effect of the latent predictor on the original functional responses.
time_index <- c("00:00","00:30",
"1:00","1:30",
"2:00","2:30",
"3:00","3:30",
"4:00","4:30",
"5:00","5:30",
"6:00","6:30",
"7:00","7:30",
"8:00","8:30",
"9:00","9:30",
"10:00","10:30",
"11:00","11:30",
"12:00","12:30",
"13:00","13:30",
"14:00","14:30",
"15:00","15:30",
"16:00","16:30",
"17:00","17:30",
"18:00","18:30",
"19:00","19:30",
"20:00","20:30",
"21:00","21:30",
"22:00","22:30",
"23:00","23:30")
NRRR.plot.reg(-fit.NRRR$Ag, fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl,
fit.NRRR$rx, fit.NRRR$ry, sseq, phi, tseq, psi,
x_ind = 1, y_ind = 2, x_lab = "Temperature",
y_lab = "Electricity Demand",
tseq_index = time_index, sseq_index = time_index,
method = "y_original"
)
NRRR.plot.reg(-fit.NRRR$Ag, fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl,
fit.NRRR$rx, fit.NRRR$ry, sseq, phi, tseq, psi,
x_ind = 1, y_ind = 6, x_lab = "Temperature",
y_lab = "Electricity Demand",
tseq_index = time_index, sseq_index = time_index,
method = "y_original"
)
One thing that needs attention is the sign of the estimated matrices. Since all the elements in the estimated ${\bf V}$ here is negative and we're visualizing the relationship between ${\bf V^Tx}(s)$ and the components of ${\bf y}(t)$, we need to specify the input Ag in function NRRR.plot.reg as -fit.NRRR$Ag to avoid producing a heatmap of opposite signs.
Prediction
Here we use both the reduced-rank regression (RRR) and nested reduced-rank regression (NRRR) to conduct a comparison in terms of the prediction power. We use the first 400 samples as the fitting set and use the remaining 108 samples as the testing set.
TXindex <- 1:400
# training set, sample size 400
Xtrain <- Xest[TXindex, ]
Ytrain <- Yest[TXindex, ]
# testing set, sample size 108
PXindex <- (1:508)[-TXindex]
Xtest <- Xest[-TXindex, ]
Ytest <- Yest[-TXindex, ]
For the same reason as before, we set eval = FALSE for this code chunk and directly show some previously saved results.
# set 'eval = TRUE' to reproduce the results
##########################
## estimation with RRR
##########################
xr <- sum(svd(Xtrain)$d > 1e-2)
rankmax <- min(xr, d * jy, 20)
fit.RRR <- NRRR:::cv.rrr(Ytrain, Xtrain, maxrank = rankmax, nfold = 10)
norder <- fit.RRR$norder
###########################
## estimation with NRRR
###########################
fit.NRRR <- NRRR.cv(Ytrain, Xtrain,
nfold = 10, norder = norder,
Ag0 = NULL, Bg0 = NULL, jx, jy, p, d, n = 400,
maxiter = 300, conv = 1e-4,
method = "RRR"
)
load("pred_compare.rda")
# predict with RRR
if (sum(abs(fit.RRR$coef)) != 0) {
svdC <- svd(fit.RRR$coef, nu = fit.RRR$rank, nv = fit.RRR$rank)
Alrrr <- svdC$v
Blrrr <- svdC$u %*% diag(svdC$d[1:fit.RRR$rank], nrow = fit.RRR$rank,
ncol = fit.RRR$rank)
} else {
Alrrr <- matrix(nrow = d * jy, fit.RRR$rank, 0)
Blrrr <- matrix(nrow = p * jx, fit.RRR$rank, 0)
}
Ypred.RRR <- NRRR.pred(
tseq, X[PXindex, , ], sseq, diag(d),
diag(p), Alrrr, Blrrr, phi
)
# predict with NRRR
Ypred.NRRR <- NRRR.pred(
tseq, X[PXindex, , ], sseq, fit.NRRR$Ag,
fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl, phi
)
# relative prediction error of Y(t) for each sample in testing set
err1 <- vector()
err2 <- vector()
for (iy in 1:108) {
err1[iy] <- sum((Y[PXindex[iy], , ] -
Ypred.RRR$Ypred[iy, , ])^2) / sum((Y[PXindex[iy], , ])^2)
err2[iy] <- sum((Y[PXindex[iy], , ] -
Ypred.NRRR$Ypred[iy, , ])^2) / sum((Y[PXindex[iy], , ])^2)
}
c(mean(err1), mean(err2))
The averaged relative prediction error of RRR is 0.42, which is much larger than the one of NRRR. Finally, we use NRRR.plot.pred to plot the predicted response trajectory. In this plot, the black line is the truly observed response trajectory and the red line is the predicted trajectory.
NRRR.plot.pred(Ypred.NRRR$Ypred, Y = Y[PXindex, , ], i_ind = 1, yi_ind = 3,
tseq, t_index = time_index, x_lab = "Time",
y_lab = "Electricity demand")
xliu-stat/NRRR documentation built on Jan. 9, 2021, 3:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(rmarkdown.html_vignette.check_title = FALSE)
Introduction: What is Nested Reduced-Rank Regression (NRRR)?
The nested reduced-rank regression (NRRR) is designed to fit a multivariate functional linear regression
$$
{\bf y}(t) = \int_{\mathcal{S}} {\bf C}0(s,t){\bf x}(s)ds + \epsilon(t),\ t\in\mathcal{T}
$$
where ${\bf y}(t)=(y_1(t),\ldots,y_d(t))^T \in \mathbb{R}^d$ with $t\in\mathcal{T}$ and ${\bf x}(s)=(x_1(s),\ldots,x_p(s))^T \in \mathbb{R}^p$ with $s\in\mathcal{S}$ are multivariate functional responses and predictors, respectively. $\epsilon(t)\in \mathbb{R}^d$ is a zero-mean random error function. The functional coefficient matrix ${\bf C}_0(s,t)=[c{k,l}(s,t)]{d\times p}$ consists of unknown bivariate functions $c{k,l}(s,t),\ k=1,\ldots,d, l=1,\ldots,p$, and NRRR aims to jointly estimate these functional surfaces by utilizing the potential associations among the functional variables.
First, a global reduced-rank structure assumes $$ {\bf C}0(s,t) = {\bf U}_0 {\bf C}^_0(s,t) {\bf V}_0^T,$$ where ${\bf U}_0 \in \mathbb{R}^{d\times r_y}$ with $r_y \leq d$ and ${\bf V}_0 \in \mathbb{R}^{p\times r_x}$ with $r_x \leq p$ provide the weights to form latent functional responses ${\bf y}^(t)={\bf U}_0^T {\bf y}(t)\in \mathbb{R}^{r_y}$ and latent functional predictors ${\bf x}^(s)={\bf V}_0^T {\bf x}(s)\in \mathbb{R}^{r_x}$, respectively. Therefore, the model is reduced to $${\bf y}^(t) = \int{\mathcal{S}} {\bf C}^_0(s,t){\bf x}^(s)ds + \epsilon^*(t),\ t\in\mathcal{T},$$ and a global dimension reduction is achieved once we have $r_y < d$ or $r_x < p$. In this way, ${\bf U}_0$ and ${\bf V}_0$ capture the association among the elements of ${\bf y}(t)$ and the association among the elements of ${\bf x}(s)$, respectively. This structure is particularly useful for simultaneously modeling a large number of functional responses and predictors that are highly correlated across $s$ or $t$.
Next, we deal with the latent regression surface ${\bf C}^_0(s,t)$. A basis expansion and truncation approach is first applied to ensure the smoothness and convert the infinite-dimensional problem to be finite-dimensional, i.e., $${\bf C}^0(s,t)\approx ({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf C}0^ ({\bf I}_{r_x}\otimes {\bf \Phi}(s)),\ {\bf C}_0^ \in \mathbb{R}^{J_yr_y\times J_xr_x},$$ where ${\bf \Psi}(t) = (\psi_1(t),\ldots,\psi{J_y}(t))$ consists $J_y$ basis functions and ${\bf \Phi}(s) = (\phi_1(s),\ldots,\phi_{J_x}(s))$ consists $J_x$ basis functions. Then, we impose the local reduced-rank structure on ${\bf C}_0^$ as $rank({\bf C}_0^)\leq r$ for $r \leq \min(J_xr_x, J_yr_y)$ and write it as a full-rank decomposition ${\bf C}_0^ = {\bf A}_0^{\bf B}_0^{T}$ for some ${\bf A}_0^ \in \mathbb{R}^{J_yr_y \times r}$ and ${\bf B}_0^* \in \mathbb{R}^{J_xr_x \times r}$. This structure induces the dependency between the latent responses and the latent predictors through their basis-expansion.
With these two layers of dimension reduction, the complexity of the model is greatly reduced. The basis functions used here are treated as given, such as spline, wavelet, and Fourier basis. In the package NRRR, the B-spline basis is mainly used. Also, we assume all the responses or the predictors share the same set of basis functions for simplicity.
Based on the selected basis functions, we obtain the integrated predictors and responses from $${\bf x}=\int_{\mathcal{S}} ({\bf I}{p}\otimes {\bf \Phi}(s)){\bf x}(s)ds \in \mathbb{R}^{J_x p}, $$ and $${\bf y}=({\bf I}_d \otimes {\bf J}{\psi\psi}^{-1/2})\int_{\mathcal{T}} ({\bf I}{d}\otimes {\bf \Psi}(t)){\bf y}(t)dt \in \mathbb{R}^{J_y d},$$ where ${\bf J}{\psi\psi}=\int_{\mathcal{T}}{\bf \Psi}(t){\bf \Psi}(t)^Tdt$ is a positive definite matrix. Thus, with all the aforementioned structures and quantities, the estimation criterion that minimizes the mean integrated squared error with respect to ${\bf C}(s,t)$, i.e., $$ \frac{1}{n} \sum_{i=1}^n \int_{\mathcal{T}} \|{\bf y}i(t) - \int{\mathcal{S}} {\bf C}(s,t){\bf x}i(s)ds \|^2 dt $$ can be written as $$ \min{{\bf U, V, A^, B^}} { \frac{1}{n}\sum_{i=1}^n \| {\bf y}i - ({\bf I}_d\otimes {\bf J}{\psi\psi}^{1/2})({\bf U}\otimes {\bf I}{J_y}){\bf A^B^}^T ({\bf V}^T\otimes {\bf I}{J_x}){\bf x}i \|^2 }. $$ By properly rearranging the columns and rows of the data matrices and coefficient matrices, the problem finally boils down to $$\min{{\bf C}} \| {\bf Y - XC} \|F^2, \ s.t.\ {\bf C}= ({\bf I}{J_x} \otimes {\bf V}) {\bf BA}^T ({\bf I}{J_y} \otimes {\bf U}^T),$$ where ${\bf Y} = ({\bf Y}{\cdot1}, \ldots, {\bf Y}{\cdot J_y})$ with ${\bf Y}{\cdot j} = (y_{ikj}){n\times d}$ for $j=1, \dots,J_y$ and ${\bf X} = ({\bf X}{\cdot1}, \ldots, {\bf X}{\cdot J_x})$ with ${\bf X}{\cdot j} = (x_{ilj}){n\times p}$ for $j=1, \dots,J_x$. That is, each ${\bf Y}{\cdot j}$ contains all the integrated values of the functional responses obtained based on $\psi_j(t)$ and ${\bf X}_{\cdot j}$ contains all the integrated values of the functional predictors obtained based on $\phi_j(s)$. This matrix approximation representation provides a clear illustration of the nested reduced-rank structure, i.e., ${\bf U}$ and ${\bf V}$ are designed to capture the shared column and row spaces among the blockwise sub-matrices of ${\bf C}$, which, as a whole, is also of low-rank. The applicability of this structure goes beyond the functional setup. For example, it can be applied in the high-dimensional vector autoregressive modeling in multivariate time series analysis, surveillance video processing and also the tensor-on-tensor regression.
This optimization problem is non-convex and has no explicit solution. We proposed a blockwise coordinate descent algorithm to obtain a local solution, and this algorithm is implemented in the package NRRR. For more details of NRRR, please read Liu, X., Ma, S., & Chen, K. (2020). Multivariate Functional Regression via Nested Reduced-Rank Regularization. arXiv: Methodology.
Next, we use the Adelaide electricity demand data as an example to show the usage of the package NRRR to investigate the functional association between weekly electricity demand trajectory and temperature trajectory. In the package NRRR, there is a function NRRR.func that takes functional observations as input and uses B-spline basis to conduct basis expansion and then fits a nested reduced-rank regression model. In the following, to provide more details about how to deal with discrete functional observations, we have not used NRRR.func. In practical application problems, the following introduction is helpful when other kinds of basis functions are preferred to conduct basis expansion.
Application: Adelaide Electricity Demand analysis
Adelaide is the capital city of the state of South Australia. In the summertime, the cooling in Adelaide mainly depends on air conditioning, which makes the electricity demand highly dependent on the weather conditions, and large volatility in temperature throughout the day could make stable electricity supply challenging. Therefore, for facilitating the supply management of electricity, it is important to understand the dependence and predictive association between the electricity demand and the temperature. Here we apply NRRR to perform a multivariate functional regression analysis between daily half-hour electricity demand profiles for the 7 days of a week and the corresponding temperature profiles for the 7 days of the same week.
Half-hourly temperature records at two locations, Adelaide Kent town and Adelaide airport, are available between 7/6/1997 and 3/31/2007. Also available are the half-hourly electricity demand records of Adelaide for the same period. The data is extracted from the R package fds.
library(fds) # Electricity Demand Data data(mondaydemand) data(tuesdaydemand) data(wednesdaydemand) data(thursdaydemand) data(fridaydemand) data(saturdaydemand) data(sundaydemand) # Temperature at Kent town data(mondaytempkent) data(tuesdaytempkent) data(wednesdaytempkent) data(thursdaytempkent) data(fridaytempkent) data(saturdaytempkent) data(sundaytempkent) # Temperature at Adelaide airport data(mondaytempairport) data(tuesdaytempairport) data(wednesdaytempairport) data(thursdaytempairport) data(fridaytempairport) data(saturdaytempairport) data(sundaytempairport)
As such, for each day during the period, there are three observed functional curves, each with 48 half-hourly observations. Here we plot the temperature and electricity demand profiles of all the Mondays from 7/6/1997 to 3/31/2007 as an illustration.
time_index <- c("00:00","00:30", "1:00","1:30", "2:00","2:30", "3:00","3:30", "4:00","4:30", "5:00","5:30", "6:00","6:30", "7:00","7:30", "8:00","8:30", "9:00","9:30", "10:00","10:30", "11:00","11:30", "12:00","12:30", "13:00","13:30", "14:00","14:30", "15:00","15:30", "16:00","16:30", "17:00","17:30", "18:00","18:30", "19:00","19:30", "20:00","20:30", "21:00","21:30", "22:00","22:30", "23:00","23:30") f1 <- rep(1:48,508) prg1 <- vector() prg2 <- vector() prg3 <- vector() f2 <- vector() for(a in 1:508){ prg1 <- c(prg1,mondaydemand$y[,a]) prg2 <- c(prg2,mondaytempkent$y[,a]) prg3 <- c(prg3,mondaytempairport$y[,a]) f2 <- c(f2,rep(a, 48)) } prgdata <- data.frame(b=f1, day=f2, PRG1=prg1, PRG2=prg2, PRG3=prg3) prgdata$b <- factor(prgdata$b) levels(prgdata$b) <- time_index library(ggplot2) # Electricity demand on Monday ggplot(prgdata, aes(y = PRG1, x = b, group = day, colour = day)) + geom_path(show.legend = FALSE) + #geom_tile() + scale_x_discrete(breaks = time_index[seq(1,48,2)]) + ggtitle("Electricity demand on Mondays") + xlab("Time") + ylab("Electricity demand") + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + theme(axis.text=element_text(size=8),axis.title=element_text(size=15), plot.title = element_text(hjust = 0.5, size = 15)) # Temperature of Kent on Monday ggplot(prgdata, aes(y = PRG2, x = b, group = day, colour = day)) + geom_path(show.legend = FALSE) + #geom_tile() + scale_x_discrete(breaks = time_index[seq(1,48,2)]) + ggtitle("Temperature of Kent on Mondays") + xlab("Time") + ylab("Temperature") + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + theme(axis.text=element_text(size=8),axis.title=element_text(size=15), plot.title = element_text(hjust = 0.5, size = 15)) # Temperature of Airport on Monday ggplot(prgdata, aes(y = PRG3, x = b, group = day, colour = day)) + geom_path(show.legend = FALSE) + #geom_tile() + scale_x_discrete(breaks = time_index[seq(1,48,2)]) + ggtitle("Temperature of Adelaide airport on Mondays") + xlab("Time") + ylab("Temperature") + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + theme(axis.text=element_text(size=8),axis.title=element_text(size=15), plot.title = element_text(hjust = 0.5, size = 15))
Since our primary focus is on studying the general association between the within-day demand and temperature trajectories in a week, we center the 48 discrete observations of each daily curve to remove the between-day trend and seasonality of the data. Each week is then treated as a replication.
# response mon <- as.matrix(mondaydemand$y) tue <- as.matrix(tuesdaydemand$y) wed <- as.matrix(wednesdaydemand$y) thu <- as.matrix(thursdaydemand$y) fri <- as.matrix(fridaydemand$y) sat <- as.matrix(saturdaydemand$y) sun <- as.matrix(sundaydemand$y) # model dimension ns <- 48 nt <- ns n <- 508 p <- 7 d <- 7 WhoWeek <- array(dim = c(p, ns, n), NA) WhoWeek[1, , ] <- apply(mon, 2, scale, center = TRUE, scale = FALSE) WhoWeek[2, , ] <- apply(tue, 2, scale, center = TRUE, scale = FALSE) WhoWeek[3, , ] <- apply(wed, 2, scale, center = TRUE, scale = FALSE) WhoWeek[4, , ] <- apply(thu, 2, scale, center = TRUE, scale = FALSE) WhoWeek[5, , ] <- apply(fri, 2, scale, center = TRUE, scale = FALSE) WhoWeek[6, , ] <- apply(sat, 2, scale, center = TRUE, scale = FALSE) WhoWeek[7, , ] <- apply(sun, 2, scale, center = TRUE, scale = FALSE) oriY <- array(dim = c(n, p, ns), NA) for (i in 1:n) { for (j in 1:p) { oriY[i, j, ] <- WhoWeek[j, , i] } } # predictor WhoWeekX <- array(dim = c(d, nt, n), NA) WhoWeekX[1, , ] <- apply(as.matrix(mondaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[2, , ] <- apply(as.matrix(tuesdaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[3, , ] <- apply(as.matrix(wednesdaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[4, , ] <- apply(as.matrix(thursdaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[5, , ] <- apply(as.matrix(fridaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[6, , ] <- apply(as.matrix(saturdaytempkent$y), 2, scale, center = TRUE, scale = FALSE) WhoWeekX[7, , ] <- apply(as.matrix(sundaytempkent$y), 2, scale, center = TRUE, scale = FALSE) oriX <- array(dim = c(n, d, nt), NA) for (i in 1:n) { for (j in 1:d) { oriX[i, j, ] <- WhoWeekX[j, , i] } } X <- oriX Y <- oriY
Till now, we've obtained the discrete observations of the response and the predictor trajectories. Next, let's use B-spline basis with 30 degrees of freedom to convert the discrete observations into its integrated form as we introduced before.
library(splines) jx <- 30 jy <- jx sseq <- c(1:48) phi <- bs(c(0, sseq), df = jx)[-1, ] tseq <- c(1:48) psi <- bs(c(0, tseq), df = jy)[-1, ] # get integrated predictors Xa <- array(dim = c(n, p, jx), NA) sdiff <- (sseq - c(0, sseq[-ns])) for (i in 1:n) { for (l in 1:p) { for (j in 1:jx) { # integrate X over s Xa[i, l, j] <- sum(phi[, j] * X[i, l, ] * sdiff) } } } # compute Jpsi and Jpsi^{-1/2} Jpsi <- matrix(nrow = jy, ncol = jy, 0) tdiff <- (tseq - c(0, tseq[-nt])) for (t in 1:nt) Jpsi <- Jpsi + psi[t, ] %*% t(psi[t, ]) * tdiff[t] eJpsi <- eigen(Jpsi) Jpsihalf <- eJpsi$vectors %*% diag(sqrt(eJpsi$values)) %*% t(eJpsi$vectors) Jpsihalfinv <- eJpsi$vectors %*% diag(1 / sqrt(eJpsi$values)) %*% t(eJpsi$vectors) # get integrated responses Ya <- array(dim = c(n, d, jy), NA) tdiff <- (tseq - c(0, tseq[-nt])) psistar <- psi %*% Jpsihalfinv for (k in 1:d) { for (i in 1:n) { for (j in 1:jy) { # integrate Y over t Ya[i, k, j] <- sum(psistar[, j] * Y[i, k, ] * tdiff) } } } # obtain matrices Y and X in matrix approximation form Yest <- Ya[, , 1] for (j in 2:jy) Yest <- cbind(Yest, Ya[, , j]) Xest <- Xa[, , 1] for (j in 2:jx) Xest <- cbind(Xest, Xa[, , j])
Here we use the daily half-hour electricity demand as the functional multivariate response with $d = 7$ (corresponding to 7 days in a week from Monday to Sunday). As for the predictors, we use the half-hour temperature data from Kent as the multivariate functional predictors, so that $p = 7$. The total sample size is $n = 508$, equaling the number of weeks in the study period.
Estimation
The package NRRR is available at Github, and we can install it with
#install.packages("devtools") #devtools::install_github("xliu-stat/NRRR") library(NRRR) set.seed(6)
As fitting the NRRR model with cross validation could be time-consuming, we set eval=FALSE in this code chunk and display some previously saved fitting results here.
# set 'eval = TRUE' to reproduce the results fit.NRRR <- NRRR.cv(Yest, Xest, nfold = 10, norder = NULL, Ag0 = NULL, Bg0 = NULL, jx = 30, jy = 30, p = 7, d = 7, n = 508, maxiter = 300, conv = 1e-4, method = "RRR" )
load("fit_nRRR.rda") # the estimated r fit.NRRR$rank # the estimated rx fit.NRRR$rx # the estimated ry fit.NRRR$ry # the estimated U fit.NRRR$Ag # the estimated V fit.NRRR$Bg
The estimated rank values are $\hat r = 4$, $\hat r_x = 1$, and $\hat r_y = 5$. The estimated loading matrix for the predictors is ${\bf V} = (-0.22, -0.39, -0.46, -0.52, -0.43, -0.28, -0.25)^T$. This shows that there is only one latent functional predictor that is driving the patterns of the electronic demands, and this factor can be roughly explained as the averaged daily temperature profile of the week. It appears that the days closer to the middle of the week load higher. On the response side, there is not much global reduction, as the estimated loading matrix ${\bf U}$ is of rank 5. To make sense of ${\bf U}$, it may be more convenient to examine the two basis vectors of its orthogonal complement, i.e., the first two singular vectors of ${\bf I - \hat U\hat U^T}$, which give the latent response factors that are not related to the temperatures at all. While the first loading vector $(-0.52, 0.36, 0.28, 0.25, −0.56, 0.34, −0.18)^T$ is hard to interpret, the second loading vector $(0.00, −0.68, 0.73, -0.01, −0.04, 0.04, −0.04))^T$ clearly indicates that the difference between the electronic demand profiles of Tuesday and Wednesday is mostly a noise process. In other words, the demand profiles of these two days are related to the temperature process in almost the same way.
Visualization of regression surface
The function NRRR.plot.reg is designed to plot the heatmap for each bivariate function in the functional coefficient matrix, e.g., $c_{k,l}(s,t),\ 1\leq k \leq d, 1 \leq l \leq p$ in ${\bf C}(s,t)$. Based on the fitting results from the nested reduced-rank regression, different kinds of regression surfaces can be visualized to give a clear illustration of the functional correlation between a user-specified predictor (or latent predictor) trajectory and response (or latent response) trajectory. With the nested reduced-rank structure, the functional regression becomes
$${\bf y}(t) \approx \int_{\mathcal{S}} {\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T {\bf x}(s)ds + {\bf \epsilon}(t),$$
based on which, we can visualize the bivariate function in:
-
${\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T \in \mathbb{R}^{d\times p}$ to see the effects of $x_i(s)$ on $y_j(t)$ with $1\leq i \leq p,\ 1\leq j \leq d$ by specifying
method = "original",x_ind = iandy_ind = j. -
${\bf U}({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)) \in \mathbb{R}^{d\times r_x}$ to see the effects of the $i$-th element of ${\bf V}^T{\bf x}(s)$ on $y_j(t)$ with $1\leq i \leq r_x,\ 1\leq j \leq d$ by specifying
method = "y_original",x_ind = iandy_ind = j. -
$({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)){\bf V}^T \in \mathbb{R}^{r_y \times p}$ to see the effects of $x_i(s)$ on the $j$-th element of ${\bf U}^T{\bf y}(t)$ with $1\leq i \leq p,\ 1\leq j \leq r_y$ by specifying
method = "x_original",x_ind = iandy_ind = j. -
$({\bf I}{r_y}\otimes {\bf \Psi}^T(t)){\bf A^B^}^T({\bf I}{r_x}\otimes {\bf \Phi}(s)) \in \mathbb{R}^{r_y \times r_x}$ to see the effects of the $i$-th element of ${\bf V}^T{\bf x}(s)$ on the $j$-th element of ${\bf U}^T{\bf y}(t)$ with $1\leq i \leq r_x,\ 1\leq j \leq r_y$ by specifying
method = "latent",x_ind = iandy_ind = j.
Based on the fitting results we obtained above, we've got only one latent predictor and five latent responses. Thus, next let's visualize the effect of the latent predictor on the original functional responses.
time_index <- c("00:00","00:30", "1:00","1:30", "2:00","2:30", "3:00","3:30", "4:00","4:30", "5:00","5:30", "6:00","6:30", "7:00","7:30", "8:00","8:30", "9:00","9:30", "10:00","10:30", "11:00","11:30", "12:00","12:30", "13:00","13:30", "14:00","14:30", "15:00","15:30", "16:00","16:30", "17:00","17:30", "18:00","18:30", "19:00","19:30", "20:00","20:30", "21:00","21:30", "22:00","22:30", "23:00","23:30") NRRR.plot.reg(-fit.NRRR$Ag, fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl, fit.NRRR$rx, fit.NRRR$ry, sseq, phi, tseq, psi, x_ind = 1, y_ind = 2, x_lab = "Temperature", y_lab = "Electricity Demand", tseq_index = time_index, sseq_index = time_index, method = "y_original" ) NRRR.plot.reg(-fit.NRRR$Ag, fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl, fit.NRRR$rx, fit.NRRR$ry, sseq, phi, tseq, psi, x_ind = 1, y_ind = 6, x_lab = "Temperature", y_lab = "Electricity Demand", tseq_index = time_index, sseq_index = time_index, method = "y_original" )
One thing that needs attention is the sign of the estimated matrices. Since all the elements in the estimated ${\bf V}$ here is negative and we're visualizing the relationship between ${\bf V^Tx}(s)$ and the components of ${\bf y}(t)$, we need to specify the input Ag in function NRRR.plot.reg as -fit.NRRR$Ag to avoid producing a heatmap of opposite signs.
Prediction
Here we use both the reduced-rank regression (RRR) and nested reduced-rank regression (NRRR) to conduct a comparison in terms of the prediction power. We use the first 400 samples as the fitting set and use the remaining 108 samples as the testing set.
TXindex <- 1:400 # training set, sample size 400 Xtrain <- Xest[TXindex, ] Ytrain <- Yest[TXindex, ] # testing set, sample size 108 PXindex <- (1:508)[-TXindex] Xtest <- Xest[-TXindex, ] Ytest <- Yest[-TXindex, ]
For the same reason as before, we set eval = FALSE for this code chunk and directly show some previously saved results.
# set 'eval = TRUE' to reproduce the results ########################## ## estimation with RRR ########################## xr <- sum(svd(Xtrain)$d > 1e-2) rankmax <- min(xr, d * jy, 20) fit.RRR <- NRRR:::cv.rrr(Ytrain, Xtrain, maxrank = rankmax, nfold = 10) norder <- fit.RRR$norder ########################### ## estimation with NRRR ########################### fit.NRRR <- NRRR.cv(Ytrain, Xtrain, nfold = 10, norder = norder, Ag0 = NULL, Bg0 = NULL, jx, jy, p, d, n = 400, maxiter = 300, conv = 1e-4, method = "RRR" )
load("pred_compare.rda") # predict with RRR if (sum(abs(fit.RRR$coef)) != 0) { svdC <- svd(fit.RRR$coef, nu = fit.RRR$rank, nv = fit.RRR$rank) Alrrr <- svdC$v Blrrr <- svdC$u %*% diag(svdC$d[1:fit.RRR$rank], nrow = fit.RRR$rank, ncol = fit.RRR$rank) } else { Alrrr <- matrix(nrow = d * jy, fit.RRR$rank, 0) Blrrr <- matrix(nrow = p * jx, fit.RRR$rank, 0) } Ypred.RRR <- NRRR.pred( tseq, X[PXindex, , ], sseq, diag(d), diag(p), Alrrr, Blrrr, phi ) # predict with NRRR Ypred.NRRR <- NRRR.pred( tseq, X[PXindex, , ], sseq, fit.NRRR$Ag, fit.NRRR$Bg, fit.NRRR$Al, fit.NRRR$Bl, phi ) # relative prediction error of Y(t) for each sample in testing set err1 <- vector() err2 <- vector() for (iy in 1:108) { err1[iy] <- sum((Y[PXindex[iy], , ] - Ypred.RRR$Ypred[iy, , ])^2) / sum((Y[PXindex[iy], , ])^2) err2[iy] <- sum((Y[PXindex[iy], , ] - Ypred.NRRR$Ypred[iy, , ])^2) / sum((Y[PXindex[iy], , ])^2) } c(mean(err1), mean(err2))
The averaged relative prediction error of RRR is 0.42, which is much larger than the one of NRRR. Finally, we use NRRR.plot.pred to plot the predicted response trajectory. In this plot, the black line is the truly observed response trajectory and the red line is the predicted trajectory.
NRRR.plot.pred(Ypred.NRRR$Ypred, Y = Y[PXindex, , ], i_ind = 1, yi_ind = 3, tseq, t_index = time_index, x_lab = "Time", y_lab = "Electricity demand")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.