This is the code for using the FQSE methodology. It includes two main files which can be loaded by calling:
source("R/fqse.R")
source("R/utils_fqse.R")The file fqse.R contains the core methodology, including the primary
iterative algorithm that alternately optimizes scores and factors. It
also includes a corresponding plot function. All functions are
documented using standard roxygen syntax.
The file utils_fqse.R provides supporting utilities, including basic
quantile operations, cross-validation routines, and an exploratory
plotting function used in the supplementary quantile error analysis.
The main user-facing function is fqse(), and results can be visualized
using the associated plot.fqse() method. An example demonstrating how
to use these functions and configure their parameters is provided below.
Before using the code, please ensure that all required package dependencies are installed.
#quantile-related dependencies, computational, plotting, and data manipulation tools
packages <- c("quantreg", "splines2", "parallel", "DoFuture",
"dplyr", "tidyr", "ggplot2", "colorspace")
#if needed call
#install.packages(packages)
#for all required dependencies check
#all_packages <- renv::dependencies()$PackageNow find a basic example on how to use FQSE and its main functions.
Let
$a_i, b_i\sim\chi^2(1)$ $U_i(t)\sim U[0,1]$ $\varepsilon_i\sim\chi^2(3)$ -
$\phi_{\epsilon}(t)\sim N(4,4)$ .
n=50
t=250
time.points = seq(0, 2*pi, length.out = t)
a = rchisq(n,1) #scaling factors
b = rchisq(n,1)
u1 = matrix(runif(n*t, min = 0, max=1), nrow = n) #uniform dt
e1 = apply(u1, MARGIN = 2, function(x){qchisq(x, 1)}) #error dt
l1 = sin(2*time.points)**2 #first pc
l2 = sin(time.points)**2 #second pc
lerr = dnorm(time.points, mean = 4, sd = 2) #error component
x = a*(3*u1**2)*matrix(rep(l1,n), byrow = TRUE, nrow = n) +
b*(2*u1)*matrix(rep(l2,n), byrow = TRUE, nrow = n)
x = x + sweep(e1, MARGIN = 2, lerr, '*')
matplot(t(x[1:20,]), type = "l")To apply FQSE on this dataset call
result = fqse(Y=x, npc=2, quantile.value =0.5, tau.grid = seq(0.1,0.9,0.1),
train.grid = c(0.05,seq(0.1, 0.9, 0.1),0.95))where we can adjust the parameters
- npc: the number of factors to be used.
- quantile.value: the base quantile to be used as reference quantile for the factors.
- tau.grid: vector of quantiles that should be included in the resulting surface.
- train.grid: vector of quantile levels to used for optimizing the score coefficients, i.e., the accumulated quantile loss over these levels will be minimized when calculating sample-specific score coefficients. It is recommended to include boundary levels.
Factors and scores of the result can be retrieved via
factors = result$factors
score.grid = result$score.gridscore.grid is a list of result$tau.grid, i.e., the
vector of quantiles that should be included in the surface.
Estimates can be accessed by calling
result_surface = collect_surface(score.grid, factors)
#transform into dataframe
df_fqse = collect_surface_df(result_surface,tau.grid = result$tau.grid, data=x)which performs score-factor matrix-multiplications for every entry in
score.grid and stores them in a list; this list is then transformed
into a dataframe. Plot the results:
gg = plot.fqse(df_fqse, tau.grid=result$tau.grid,samples=1:5)
gg
New surfaces
can be collected anytime using the predict_new_scoresfunction:
#the new grid needs to include EVERY desired quantile level
#(not just the ones you would like to add to the previous result).
new.grid = c(0.05, 0.45, 0.6, 0.8)
#result$score.expansion includes the score spline basis and coefficients,
#result$rotation.vector the identifiability information
score.grid.new = predict_new_scores(result$score.expansion, result$rotation.vector,
new.grid = new.grid)
new.surface = collect_surface(score.grid.new, factors)
new.surface.df = collect_surface_df(new.surface, tau.grid = new.grid, data=x)
gg_new = plot.fqse(new.surface.df, tau.grid=new.grid, samples=1:5)
gg_new
