resulting $R^{2}_{m_{0}}$ values are presented in Table 3. Observe that
$R_{m_{0}}^{2}$ for $m_{0}>1$ not necessarily is an increasing function of the
number of covariates included in the model. This is partly due to the fact
that all covariate vectors are weighted equally in the definition of
$R_{m_{0}}^{2}$, regardless of their number of observations.

Covariates | Link function | $g_{m}$ | $R^{2}_{1}$ | $R^{2}_{2}$ | $R^{2}_{3}$ | $R^{2}_{4}$
---|---|---|---|---|---|---
Age, sex, year |  | Legendre | 15.29 % | 56.99% | 56.94% | 58.13%
Identity | Hermite | 15.29% | 57.55% | 57.60% | 57.77%
| Laguerre | 15.29% | 50.22% | 55.63% | 56.52%
| Legendre | 15.30% | 56.90% | 56.84% | 57.94%
Logit | Hermite | 15.30% | 57.34% | 57.31% | 57.37%
| Laguerre | 15.30% | 49.71% | 55.27% | 56.15%
Year |  | Legendre | 2.89 % | 67.55% | 67.61% | 70.91%
Identity | Hermite | 2.89% | 70.88% | 71.01% | 71.05%
| Laguerre | 2.89% | 60.61% | 66.53% | 67.59%
| Legendre | 2.95% | 67.58% | 67.63% | 70.93%
Logit | Hermite | 2.95% | 70.85% | 70.98% | 71.02%
| Laguerre | 2.95% | 60.58% | 66.53% | 67.57%
Table 3: $R^{2}_{m_{0}}$-values for two choices of covariates, the identity
and logit link functions and the three polynomial systems of weight functions
defined in Examples 5-7. (A weighted version of) the classical coefficient of
determination, $R^{2}_{1}$, reflects that more of the variation in the
response is explained when age and sex are added as covariates, compared to
having only year as covariate. Notice though that $R^{2}_{m_{0}}$ for
$m_{0}>1$ is higher for a model with year as the only covariate compared to
using all three covariates. The reason is that the year only model has an
increased number of observations in each covariate specific subset, which
leads to a smoother $\hat{Q}(p|\mbox{\boldmath$x$})$ with smaller jumps at
each discontinuity point. This, in turn, leads to a closer fit of
$\hat{Q}_{\text{appr}}(p|\mbox{\boldmath$x$})$.

####  5.3.3 Covariate and $L$-functional estimates

In this section we will first illustrate how the $q=4$ parameter estimates of
(154) change as a function of the quantile $p$ (cf. Figure 6) over a grid of
100 values in $(0,1)$ for the case of an identity link function. Next we plot
the change in the location, scale, standardized skewness and standardized
kurtosis over time, for each subpopulation, i.e. each combination of age and
sex. These plots are shown in Figure 7. All estimates include approximate 95%
bootstrapped pointwise confidence intervals, represented by ribbons. The
negative estimates of the effect of year in Figure 6(d), and the negative
slopes over time of the expected value of conditional response distribution,
in Figure 7(a), both demonstrate that birds in recent years tend to arrive
earlier.

We also use the resampling results to present tables over the Mahalanobis
distance between different pairs of covariates’ $L$-functionals. For this we
need some additional notation. For any subset $I\subset\\{1,2,32,42\\}$ of
order numbers of location, scale, standardized skewness and standardized
kurtosis, we introduce the collection

$\mbox{\boldmath$\theta$}_{I}(\mbox{\boldmath$x$})=(T_{m}(F_{Y|\mbox{\scriptsize\boldmath$x$}});\,m\in
I)^{\top}$

of conditional (standardized) $L$-moments for covariate vector $x$. The
corresponding vector of estimated conditional $L$-moments is denoted

$\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$})=(T_{m}(\hat{F}_{Y|\mbox{\scriptsize\boldmath$x$}});\,m\in
I)^{\top}.$ (159)

Recall that each $x$ corresponds to a subpopulation (a combination of age
group and sex) at a specific time point. Assume that the data set is resampled
$B$ times, and let $F^{\ast b}_{Y|\mbox{\scriptsize\boldmath$x$}}$ be the
$b$:th resampled conditional response distribution ($b=1,\ldots,B$) for
covariate vector $x$. The corresponding resampled vector of conditional
$L$-moments is

$\mbox{\boldmath$\theta$}_{I}^{\ast b}(\mbox{\boldmath$x$})=(T_{m}(F^{\ast
b}_{Y|\mbox{\scriptsize\boldmath$x$}});\,m\in I)^{\top}.$ (160)

These $B$ vectors will be scattered around (159), with an estimated covariance
matrix

$\hat{\mbox{\boldmath$\Sigma$}}_{I}(\mbox{\boldmath$x$})=\frac{1}{B}\sum_{b=1}^{B}(\mbox{\boldmath$\theta$}_{I}^{\ast
b}(\mbox{\boldmath$x$})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))(\mbox{\boldmath$\theta$}_{I}^{\ast
b}(\mbox{\boldmath$x$})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))^{\top}.$

In particular, the Mahalanobis distance between (160) and the center point
(159) of the distribution is

$M_{I}^{b}(\mbox{\boldmath$x$})=\sqrt{(\mbox{\boldmath$\theta$}_{I}^{\ast
b}(\mbox{\boldmath$x$})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))^{T}\hat{\mbox{\boldmath$\Sigma$}}_{I}(\mbox{\boldmath$x$})^{-1}(\mbox{\boldmath$\theta$}_{I}^{\ast
b}(\mbox{\boldmath$x$})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))}.$
(161)

This distance can be interpreted as how many standard deviations away from the
center point $\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$})$ the
point $\mbox{\boldmath$\theta$}_{I}^{\ast b}(\mbox{\boldmath$x$})$ is. We
would however like to measure the distance between the center of the resampled
clouds for different pairs $\mbox{\boldmath$x$}_{1}$ and
$\mbox{\boldmath$x$}_{2}$ of covariate vectors, rather than each resampled
point’s distance to its center. The Mahalanobis distance between
$\mbox{\boldmath$x$}_{1}$ and $\mbox{\boldmath$x$}_{2}$, for the collection
$I$ of conditional $L$-moments, is

$M_{I}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2})=\\!\sqrt{(\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}_{1})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))^{T}[0.5(\hat{\mbox{\boldmath$\Sigma$}}_{I}(\mbox{\boldmath$x$}_{1})+\hat{\mbox{\boldmath$\Sigma$}}_{I}(\mbox{\boldmath$x$}_{2}))]^{-1}(\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}_{1})-\hat{\mbox{\boldmath$\theta$}}_{I}(\mbox{\boldmath$x$}))}.$
(162)

In Tables 4 and 5 we present values of the modified Mahalanobis distances
(162) for different combinations of $I$, $\mbox{\boldmath$x$}_{1}$, and
$\mbox{\boldmath$x$}_{2}$. From these tables it can be seen that the location
functional is most important for distinguishing the arrival distributions
between subpopulations (Table 4) and years (Table 5). Although the scale,
standardized skewness and standardized heavytailedness functionals are less
important, they still help to discriminate even more between the arrival time
distributions of these groups.

(a) The intercept.

(b) Effect of the bird being adult.

(c) Effect of the bird being male.

(d) Effect of year.

Figure 6: Plots of parameter estimates for the regression model with the best
fit, using the identity link. The black lines represent the $q=4$ estimates
$\hat{\beta}_{1}(p),\hat{\beta}_{2}(p),\hat{\beta}_{3}(p),\hat{\beta}_{4}(p)$,
where $p$ takes values on a grid of 100 equispaced points in $(0,1)$. The grey
ribbon of each subplot consists of approximate $95\%$ pointwise bootstrapped
confidence interval for $\\{\beta_{k}(p);\,0<p<1\\}$. Note that in (b), the
confidence interval is very broad for the lowest quantiles, which is due to
there being extremely few juvenile birds among the early arrivers. In (d), the
effect of year quickly approaches 0 in the upper quantiles. This might be due
to an increased chance of catching stray birds instead of migrating birds at
that time of the year. Lastly, smoothing $\hat{\boldsymbol{\beta}}(p)$ seems
like a good idea, given the wigglyness of the estimates.

(a) $T_{1}(\hat{F}_{Y|\mbox{\scriptsize\boldmath$x$}})$

(b) $T_{2}(\hat{F}_{Y|\mbox{\scriptsize\boldmath$x$}})$

(c) $T_{32}(\hat{F}_{Y|\mbox{\scriptsize\boldmath$x$}})$

(d) $T_{42}(\hat{F}_{Y|\mbox{\scriptsize\boldmath$x$}})$

Figure 7: Estimates of the four conditional $L$-functionals of location, scale, standardized skewness and standardized kurtosis for the bird migration data, for each combination $x$ of year and the binary covariates age and sex. These conditional $L$-functionals are estimated from a linear quantile regression model, cf. (95) and (109), using the Legendre weight functions $g_{m}(p)$. A subplot illustrates a specific conditional $L$-functional. It contains four curves, each one of which corresponds to a fixed combination of age and sex, whereas year varies along the horizontal axis. For the location plot (a), all polynomial series would give the same result (the conditional mean), since the first polynomial in each series is 1. For the scale plot (b) the lines correspond to Gini’s mean difference times a proportionality constant (cf. Example 5). The standardized skewness (c) and standardized kurtosis (d) are measured relative to a uniform distribution, since we use weight functions based on Legendre polynomials. The important information obtained from the figure is qualitative; how the various $L$-functionals change over time, rather than their actual values. From (c) we notice that the arrival distribution shifts from right-skewed towards symmetric over the study period, whereas from (d) we observe that the kurtosis increases. The conclusion would be that the left tail slowly becomes as thick as the right tail, meaning that early arrivers become more frequent over time. | Subpopulation comparison
---|---
Order | Juvenile female | Juvenile male | Juvenile female | Adult female
number set $I$ | v. Adult female | v. Adult male | v. Juvenile male | v. Adult male
$\\{1\\}$ | 12.41 | 10.46 | 19.27 | 32.44
$\\{2\\}$ | 7.14 | 6.21 | 2.97 | 4.36
$\\{32\\}$ | 0.03 | 0.03 | 0.33 | 0.53
$\\{42\\}$ | 1.71 | 1.77 | 1.13 | 1.44
$\\{1,2\\}$ | 14.83 | 11.81 | 19.43 | 33.22
$\\{1,32\\}$ | 12.59 | 10.47 | 19.53 | 32.44
$\\{1,42\\}$ | 12.49 | 10.47 | 19.36 | 33.19
$\\{2,32\\}$ | 7.29 | 6.21 | 2.97 | 4.37
$\\{2,42\\}$ | 7.39 | 6.31 | 3.25 | 4.60
$\\{32,42\\}$ | 1.75 | 1.77 | 1.14 | 1.54
$\\{1,32,42\\}$ | 15.36 | 11.82 | 19.65 | 33.23
$\\{1,2,42\\}$ | 14.89 | 11.81 | 19.51 | 33.98
$\\{1,32,42\\}$ | 12.75 | 10.48 | 19.69 | 33.19
$\\{2,32,42\\}$ | 7.47 | 6.31 | 3.26 | 4.61
$\\{1,2,32,42\\}$ | 15.56 | 11.82 | 19.79 | 33.98
Table 4: The Mahalanobis distance $M_{I}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2})$ between different pairs of subpopulations for the year 2000. These distances are computed for different collections $I$ of conditional $L$-moments. In accordance with Figure 7 the subpopulations mainly differ in location, although some of them have quite large values of $M_{I}$ for scale as well. It can be seen that skewness differs the least and kurtosis the second least between subpopulations.  | Subpopulation
---|---
$I$ | Juvenile female | Juvenile male | Adult female | Adult male
$\\{1\\}$ | 23.66 | 19.46 | 29.81 | 27.11
$\\{2\\}$ | 2.41 | 2.05 | 2.80 | 2.46
$\\{32\\}$ | 3.20 | 3.01 | 4.13 | 3.57
$\\{42\\}$ | 3.27 | 2.86 | 3.73 | 3.17
$\\{1,2\\}$ | 23.67 | 19.78 | 30.06 | 28.76
$\\{1,32\\}$ | 24.91 | 20.04 | 30.32 | 27.17
$\\{1,42\\}$ | 23.67 | 19.48 | 30.27 | 27.17
$\\{2,32\\}$ | 3.71 | 3.61 | 4.73 | 4.73
$\\{2,42\\}$ | 4.25 | 3.89 | 4.86 | 4.11
$\\{32,42\\}$ | 4.00 | 4.31 | 5.84 | 5.33
$\\{1,2,32\\}$ | 24.92 | 20.34 | 30.49 | 28.94
$\\{1,2,42\\}$ | 23.67 | 19.78 | 30.47 | 28.84
$\\{1,32,42\\}$ | 25.15 | 20.05 | 30.71 | 27.25
$\\{2,32,42\\}$ | 4.61 | 5.05 | 6.39 | 6.35
$\\{1,2,32,42\\}$ | 25.19 | 20.34 | 30.84 | 29.08
Table 5: The Mahalanobis distance
$M_{I}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2})$, for each
subpopulation, between the years 1982 and 2019. The distances are computed for
different collections $I$ of conditional $L$-functionals. As for Table 4 the
main differences are in location. However, standardized skewness and
standardized kurtosis both have values of $M_{I}$ slightly below 4. If we
regard 4 standard deviations as significant, these values are on the
borderline of demonstrating significant differences in skewness and kurtosis
between 1982 and 2019. Notice also that the Mahalanobis distances of
$I=\\{2,32,42\\}$ are larger than 4 for all four subpopulations, indicating
that there are other changes other than those in location in the arrival
distribution between 1982 and 2019.

##  6 Discussion

In this paper we developed a general theory of $L$-functionals of the response
variable distribution of regression models. Based on orthogonal series
expansions of the quantile functions of these distributions we generalized the
concept of $L$-moments and identified collections $\\{T_{m}(F)\\}_{m=1}^{4}$
of $L$-functionals that correspond to measures of location, scale,
unstandardized skewness and unstandardized heavytailedness of the response.
Different collections of $L$-functionals were introduced, depending on whether
the domain of the response variable is bounded, or unbounded in one or two
directions.

A number of generalizations of our work is possible. The first extension is to
study more systematically the standardized $L$-functionals of skewness,
heavytailedness, … $T_{m2}(F)=T_{m}(F)/T_{2}(F)$ when $m\geq 3$. In
particular, it is of interest to know which types of stochastic orderings
between response distributions (Oja, 1981) these functionals preserve. Second,
Takemura (1983) defined orthogonal series expansions of quantile functions for
arbitrary reference distributions $F_{0}$ for which $T_{m}(F_{0})$ vanish when
$m\geq 3$. In Examples 5-7 we considered expansions for uniform, normal, and
exponential distributions. It is also possible to define a system (68) of
$L$-functionals for Weibull, log-logistic, and other reference distributions
$F_{0}$ that are of interest in survival analysis. Third, a general way of
robustifying a collection $\\{T_{m}(F)\\}$ of $L$-functionals is to introduce
a weight function $w$ that downweights contributions from the lower and upper
tails of $F$. The weight density

$g_{m}(p)=w(p)P_{m-1}(Q_{0}(p))$

of order $m$ is a generalization of (68), using a variable rather than a
constant weight function $w(p)\equiv 1$. For instance, the trimmed $L$-moments
of Elamir and Seheult (2003) correspond to a weight function
$w(p)=(1-p)^{t}p^{t}$ for some positive integer $t$. Notice however that in
general the orthogonality property (69) is lost for such a system of weight
functions. Fourth, recursive estimation of quantiles (Stephanou et al., 2017)
could be extended to online estimation of $L$-functionals. Fifth, regression
quantiles have been used for time series (Cai and Zu, 2008, White et al.,
2008) in order to estimate quantiles and robust measures of skewness/kurtosis
of predictive distributions. It is of interest to analyze (ratios of)
$L$-functionals of such predictive distributions. Sixth, measures of location,
scale, skewness and heavytailedness of multivariate distributions can be
defined in terms of multivariate $L$-statistics (Liu, 1990, Liu et al. 1999,
Zuo et al., 2004, and Dang et al., 2009). In this context it is of interest to
study different polynomial systems of $L$-functionals and their order numbers.

## References

Bennett, C.A. (1952). Asymptotic properties of ideal linear estimators. Ph.D.
dissertation, University of Michigan.

Bennett, S. (1983). Analysis of survival data by the proportional odds model.
Statistics in Medicine 2, 273-277.

Bickel, P.J. (1965). On some robust estimators of location. Annals of
Mathematical Statistics 36, 847-858.

Bickel, P.J. (1973). On some analogues of linear combinations of order
statistics in the linear model. Annals of Statistics 1, 597-616.

Bickel, P.J. and Lehmann, E.L. (1975). Descriptive statistics for
nonparametric models II: Location. Annals of Statistics 3, 1045-1069.

Bickel, P.J. and Lehmann, E.L. (1976). Descriptive statistics for
nonparametric models III: Dispersion. Annals of Statistics 4, 1139-1158.

Billingsley, P. (1999). Convergence of Probability Measures, 2nd ed. Wiley,
New York.

Bottai, M., Cai, B. and McKeown, R.E. (2010). Logistic quantile regression for
bounded outcomes. Statistical Medicine 29, 309-317.

Bowley, A.L. (1920). Elements of Statistics, Schribner’s, New York.

Box, G.E.P. and Cox, D.R. (1964). An analysis of transformations. Journal of
the American Statistical Association 26, 211-252.

Cai, Z. and Xu, X. (2008). Nonparametric quantile estimations for dynamic
smooth coefficient models. Journal of the American Statistical Association
103, 1595-1608. This reference was removed since it deals time series data.

Chauduri, P. (1991). Global nonparametric estimation of conditional quantile
functions and their derivatives. Journal of Multivariate Analysis 39, 246-269.

Chauduri, P. and Loh, W.-L. (2002). Nonparametric estimation of conditional
quantiles using quantile regression trees. Bernoulli 8, 561-576.

Chen, X., Linton, O. and van Kellegom, I. (2003). Estimation of semiparametric
functions when the criterion function is not smooth. Econometrica 71,
1591-1608.

Chernoff, H., Gastwirth, J.L. and Jones, M.V. (1967). Asymptotic distribution
of linear combinations of functions of order statistics with applications to
estimation. Annals of Mathematical Statistics 38, 52-72.

Chissom, B.S. (1970). Interpretation of the kurtosis statistics. American
Statistician 24(4), 19-23.

Cox, D.R. (1972). Regression models and life tables. Journal of the Royal
Statistical Society, Series B, 34, 187-220.

Dang, X., Serfling, R., and Zhou, W. (2009). Influence functions of some depth
functions, and application to depth-weighted L-statistics. Journal of
Nonparametric Statistics, 21(1), 49-66.

Doksum, K. and Gasko, M. (1990). On a correspondence between models in binary
regression and survival analysis. International Statistical Review 58,
243-252.

Elamir, E.A. and Seheult, A.H. (2003). Trimmed L-moments. Computational
Statistics and Data Analysis 43, 299-314.

Efron, B. (1991). Regression percentiles using asymmetric squared error loss.
Statistica Sinica 1(1), 93-125.

Efron, B. and Tibshirani, R. (1993). An introduction to the bootstrap, Chapman
and Hall, Boca-Raton, FL.

Fasiolo M., Goude Y., Nedellec R. and Wood S. N. (2017). Fast calibrated
additive quantile regression. URL: https://arxiv.org/abs/1707.03307

Fisher, R.A. and Cornish, E.A. (1960). The percentile points of distributions
having known cumulants. Technometrics 2, 209-225.

Forsythe, G. E. (1977). Computer methods for mathematical computations.
Prentice-Hall series in automatic computation, 259.

Frumento, P. and Bottai, M. (2016). Parametric modeling of quantile regression
coefficient models. Biometrics 72, 74-84.

Frumento, P. and Bottai, M. (2017). An estimating equation for censored and
truncated quantile regression. Computational Statistics and Data Analysis 113,
53-63.

Galton, F. (1883). Enquiry to human faculty and its development. London,
MacMillan.

Garcia, V.J., Martel-Escobar, M. and Vázquez-Polo, F.J. (2018). A note on
ordering probability distributions by skewness. Symmetry 10(7), 286.

García-Pareja, C. and Bottai, M. (2018). On mean decomposition for summarizing
conditional distributions. Statistics 2018:7:e208.

García-Pareja, C., Santacatterian, M., Ekström, A.M. and Bottai, M. (2019). In
Topics in mathematical statistics for medical applications: summary measures
and exact simulation of diffusions. PhD Thesis, Institute of Environmental
Medicine, Karolinska Institutet, Stockholm.

Gelfand, A.E., Ghosh, S.K., Christiansen, C., Somerai, S.B. and McLaughlin,
T.J. (2000). Proportional hazards models: a latent competing risk approach.
Applied Statistics 49, 385-397.

Gilchrist, W. (2000). Statistical Modelling with Quantile Functions. Chapman
and Hall.

Gilchrist, W. (2007). Modelling and fitting quantile distributions and
regressions. American Journal of Mathematical and Management Sciences 27,
401-439.

Greenwood, J.A., Landwehr, J.M., Matalas, N.C. and Wallis, J.R. (1979)
Probability weighted moments: definition and relation to parameters of several
distributions expressable in inverse form. Wat. Resour. Res. 15, 1049-1054.

Groeneveld, R.A. and Meeden, G. (1984). Measuring skewness and kurtosis. The
Statistician 33, 391-399.

Gutenbrunner, C. and Jurečková, J. (1992). Regression quantile and regression
rank score processes in the linear model and derived statistics. Annals of
Statistics 20, 305-330.

Haeusler, E. and Teugels, J.L. (1985). On asymptotic normality of Hill’s
estimator for the exponent of regular variation. Annals of Statistics 13(2),
743-756.

He, X. and Shi, P. (1994). Convergence rate of $B$-spline estimators of
nonparametric quantile functions. Journal of Nonparametric Statistics 33,
299-308.

Hill, B.M. (1975). A simple general approach to inference about the tail of a
distribution. Annals of Statistics 3, 1163-1174.

Hinkley, D.V. (1975). On power transformations to symmetry. Biometrika 62,
101-111.

Hogg, R.V. (1972). More light on the kurtosis and related statistics. Journal
of the American Statistical Association 67, 422-424.

Hogg, R.V. (1974). Adaptive robust procedures: A partial review of some
suggestions for future applications and theory. Journal of the American
Statistical Association 67, 422-424.

Hosking, J.R.M. (1990). L-moments: analysis and estimation of distributions
using linear combinations or order statistics. Journal of the Royal
Statistical Society Ser. B 52(1), 105-124

Hosking, J.R.M. (1992). Moments or $L$ moments? An example comparing two
measures of distributional shape. The American Statistician 46(3), 186-189.

Hosking, J.R.M. (2006). On the characterization of distributions by their
L-moments. Journal of Statistical Inference and Planning 136(1), 193-198.

H ssjer, O. (2008). On the coefficient of determination for mixed regression
models. Journal of Statistical Planning and Inference 138, 3022-3038.

Jung, J. (1955). On linear estimates defined by continuous weight function.
Arkiv für Matematik $d$Bd 3, 199-209.

Kalbfleish, J.D. and Prentice, R.L. (2002). Statistical Analysis of Failure
Time Data (2nd ed.), Wiley, Noboken, NJ.

Karian, Z.A. and Dudewicz, E.J. (2000). Fitting Statistical Distributions: The
Generalized Lambda Distribution and Generalized Bootstrap Methods. CRC Press,
Boca Raton, Florida.

Karvanen, J. (2006). Estimation of quantile mixtures via L-moments and trimmed
L-moments. Computational Statistics and Data Analysis 51(2), 947-959.

Karvanen, J. and Nuutinen, A. (2008). Characterizing the generalized lambda
distribution by L-moments. Computational Statistics and Data Analysis 52,
1971-1983.

Kim, M.-O. (2007). Quantile regression with varying coefficients. Annals of
Statistics 35(1), 92-108.

Kim, T.-H. and White, H. (2004). On more robust estimation of skewness and
kurtosis. Finance Research Letters 1, 56-73.

Koenker, R. (2005). Quantile Regression. Cambridge University Press,
Cambridge.

Koenker, R. (2020). quantreg: Quantile Regression. R package version 5.73.
https://CRAN.R-project.org/package=quantreg

Koenker, R. and Bassett Jr, G. (1978). Regression quantiles. Econometrica
46(1), 33-50.

Koenker, R. and Bassett Jr, G. (1982). Robust tests for heteroscedasticity
based on regression quantiles. Econometrica 50, 43-61.

Koenker, R. and Geling, O. (2001). Reappraising medfly longevity: A quantile
regression survival analysis. Journal of the American Statistical Association
96, 458-468.

Koenker, R. and Hallock, K. (2001). Quantile regression: An introduction.
Journal of Economic Perspectives 15, 143-156.

Koenker, R., Ng, N. and Portnoy, S. (1994). Quantile smoothing splines.
Biometrika 81, 673-680.

Koenker, R. and Portnoy, S. (1987). $L$-estimation for linear models. Journal
of the American Statistical Association 82, 85-1857.

Koenker, R. and Zhao, Q. (1994). $L$-estimation for linear heteroscedastic
models. Journal of Nonparametric Statistics 3, 223-235.

Lehikoinen, A., Lindén, A., Karlsson, M., Andersson, A., Crewe, T. L., Dunn,
E. H., Gregory, G., Karlsson, L., Kristiansen, V., Mackenzie, S. and others
(2019). Phenology of the avian spring migratory passage in Europe and North
America: Asymmetric advancement in time and increase in duration. Ecological
Indicators 101, 985-999.

Leng, C. and Tong, X. (2013). A quantile regression estimator for censored
data. Bernoulli 19(1), 344-361.

Lindén, A., Meller, K. and Knape, J (2017). An empirical comparison of models
for the phenology of bird migration. Journal of Avian Biology 48, 255-265.

Lindgren, A. (1997). Quantile regression with censored data using generalized
$L_{1}$ minimization. Computational Statistics and Data Analysis 23, 509-524.

Liu, R. Y. (1990). On a notion of data depth based on random simplices. The
Annals of Statistics, 405-414.

Liu, Y. and Bottai, M. (2009). Mixed-effects models for conditional quantiles
with longitudinal data. International Journal of Biostatistics 5(1), Article
28.

Liu, R. Y., Parelius, J. M., and Singh, K. (1999). Multivariate analysis by
data depth: descriptive statistics, graphics and inference (with discussion
and a rejoinder by liu and singh). The annals of statistics, 27(3), 783-858.

McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models, Chapman and
Hall/CRC, New York.

McGillivray, H.L. (1986). Skewness and asymmetry: Measures and orderings.
Annals of Statistics 14, 994-1011.

Moore, D.S. (1968). Annals of Mathematical Statistics 39(1), 263-265.

Moors, J.J.A. (1988). Q quantile alternative for kurtosis. The Statistician
37, 25-32.

Mosteller, F. (1946). On some useful “inefficient” statistics. Annals of
Mathematical Statistics 17, 377-408.

Mu, Y. and He, X. (2007). Power transformation toward a linear regression
quantile. Journal of the American Statistical Association 102, 269-279.

Mudholkar, G.S. and Hutson, A.D. (1998). $LQ$-moments: analogs of $L$-moments.
Journal of Statistical Planning and Inference 71, 191-208.

Neocleous, T., Vanden Branden, K. and Portnoy, S. (2006). Correction to
Portnoy (2003). Journal of the American Statistical Association 101, 860-861.

Newey, W.K. and McFadden, D.I. (1994). Large sample estimation and hypothesis
testing In: Engle, R.F., McFadden, D.I. (eds.) Handbook of Econometics 4,
North-Holland, Amsterdam, pp. 2211-2245.

Oberhofer, W. (1982). The consistency of nonlinear regression minimizing the
$L_{1}$ norm. Annals of Statistics 10, 316-319.

Oja, H. (1981). On location, scale, skewness and kurtosis of univariate
distributions. Scandinavian Journal of Statistics 8, 154-168.

Okagbue, H.I., Adamu, M.O. and Anake, T.A. (2019). Quantile mechanics: Issues
arising from critical review. International Journal of Advanced and Applied
Sciences 6(1), 9-23.

Olkin, I. and Pukelsheim, F. (1982). The distance between two random vectors
with given dispersion matrices. Linear Algebra Appl. 48, 257-263.

Parzen, E. (1979). Nonparametric statistical modelling (with commmets).
Journal of the American Statistical Association 74, 105-131.

Peng, L. and Huang, Y. (2008). Survival analysis with quantile regression
models. Journal of the American Statistical Association 103, 637-649.

Portnoy, S. (2003). Censored quantile regression. Journal of the American
Statistical Association 98, 1001-1012.

Powell, J.L. (1986). Censored regression quantiles. Journal of Econometrics
32, 143-155.

R Core Team (2021). R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/.

Rousseeuw, P.J. and Leroy, A. (1987). Robust regression and outlier detection,
Wiley, New York.

Royston, P. and Parmar, M.K.B. (2002). Flexible parametric proportional-
hazards an proportional-odds models for censored survival data, with
application to prognostic modelling and estimation of treatment effects.
Statistics in Medicine 21, 2175-2197.

Ruppert, D. and Carroll, R. (1980). Trimmed least squares estimation in the
linear model. Journal of the American Statistical Association 75, 828-838.

Serfling, R. (1980). Asymptotic Theory of Statistics. John Wiley, New York.

Sheather, S.J and Marron. J.S. (1990). Kernel quantile estimators. Journal of
the American Statistical Association 85, 410-416.

Sillito, G. (1969). Derivation of approximants to the inverse distribution
function of a continuous univariate population from the order statistics of a
sample. Biometrika 56(3), 641-650.

Stephanou, M., Varughese, M. and Mcdonald, I. (2017). Sequential quantiles via
Hermite series density estimation. arXiv:15017.05073v2.

Stigler, S. (1977). Do robust estimators work with real data? Annals of
Statistics 5, 1055-1098.

Takamura, A. (1983). Orthogonal expansion of quantile function and components
of the shapiro-statistic, Report No. TR-8, Stanford University, Department of
Statistics, California, USA.

Tukey, J.W. and McLaughlin, G. (1963). Less vulnerable confidence and
significance procedures for location based on a single sample. Sankya Ser. A
35, 331-352.

Wang, H.J. and Wang, L. (2009). Locally weighted censored quantile regression.
Journal of the American Statistical Association 104, 1117-1128.

van Zwet, W.R. (1964). Convex Transformations of Random Variables, Math.
Centrum, Amsterdam.

Welsh, A.H. (1987). The trimmed mean in the linear model. Annals of Statistics
15(1), 20-36.

Welsh, A.H. and Morrison, H.L. (1990). Robust L-estimation of scale with an
application to astronomy. Journal of the American Statistical Association 85,
729-743.

White, H., Tae-Hwan, K. and Manganelli, S. (2008). Modeling autoregressive
conditional skewness and kurtosis with multi-quantile CAViaR. European Central
Bank, Working paper 957, November 2008.

Yang, S. (1999). Censored median regression using weighted empirical survival
and hazard functions. Journal of the American Statistical Association 94,
137-145.

Ying, Z., Jung, S.H. and Wei, L.J. (1995). Survival analysis with median
regression models. Journal of the American Association 90, 178-184.

Younes, N. and Lachin, J. (1997). Link-based models for survival data with
interval and continuous time censoring. Biometrics 53, 1199-1211.

Zuo, Y., Cui, H., and Young, D. (2004). Influence function and maximum bias of
projection depth based estimators. The Annals of Statistics, 32(1), 189-218.


# Spike2Vec: An Efficient and Scalable Embedding Approach for COVID-19 Spike
Sequences

1st Sarwan Ali  Department of Computer Science
Georgia State University
Atlanta, GA, USA
<EMAIL_ADDRESS>2nd Murray Patterson  Department of Computer Science
Georgia State University
Atlanta, GA, USA
<EMAIL_ADDRESS>
###### Abstract

With the rapid global spread of COVID-19, more and more data related to this
virus is becoming available, including genomic sequence data. The total number
of genomic sequences that are publicly available on platforms such as GISAID
is currently several million, and is increasing with every day. The
availability of such Big Data creates a new opportunity for researchers to
study this virus in detail. This is particularly important with all of the
dynamics of the COVID-19 variants which emerge and circulate. This rich data
source will give us insights on the best ways to perform genomic surveillance
for this and future pandemic threats, with the ultimate goal of mitigating or
eliminating such threats. Analyzing and processing the several million genomic
sequences is a challenging task. Although traditional methods for sequence
classification are proven to be effective, they are not designed to deal with
these specific types of genomic sequences. Moreover, most of the existing
methods also face the issue of scalability. Previous studies which were
tailored to coronavirus genomic data proposed to use spike sequences
(corresponding to a subsequence of the genome), rather than using the complete
genomic sequence, to perform different machine learning (ML) tasks such as
classification and clustering. However, those methods suffer from scalability
issues.

In this paper, we propose an approach called Spike2Vec, an efficient and
scalable feature vector representation for each spike sequence that can be
used for downstream ML tasks. Through experiments, we show that Spike2Vec is
not only scalable on several million spike sequences, but also outperforms the
baseline models in terms of prediction accuracy, F1-score, etc. Since this
type of study on such huge numbers of spike sequences has not been done before
(to the best of our knowledge), we believe that it will open new doors for
researchers to use this data and perform different tasks to unfold new
information that was not available before. As an example of this, we use
information gain (IG) to compute the importance of each amino acid in the
spike sequence. The amino acids with higher IG values tend to be the same to
the ones reported by the USA based Centers for Disease Control and Prevention
(CDC) for different variants.

###### Index Terms:

COVID-19 Spike Sequences, Feature Vector Representation, k-mers,
Classification, Clustering

##  I Introduction

Very few fields of study remain untouched in the big data era, as massive
amounts of data are collected in every domain from finance [1, 2] to astronomy
[3, 4]. The field of biomedical and health informatics is no exception; one
which has had a recent and rather rapid growth spurt in the amount of
available data, due to the COVID-19 pandemic [5, 6]. One facet of this
increase is the amount of genomic data becoming available for COVID-19 in
databases such as GISAID111www.gisaid.org, where several million viral genome
(virome) sequences of COVID-19 — or more precisely, SARS-CoV-2222Severe Acute
Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), the virus which causes the
COVID-19 disease — are available.

Such data has a high _volume_ , as the SARS-CoV-2 virome has $\approx 30$K
nucleotide base-pairs, and there are more than 2.5 million such sequences
available in GISAID alone. While the number COVID-19 patients being sequenced
is a fraction of the actual number of cases, the sheer number of infections
(both now and in the past) means that the _velocity_ in which SARS-CoV-2
virome sequences are appearing is very high. For example, in March 2020, when
COVID-19 was declared a pandemic by the world health organization (WHO), there
were a few thousand sequences available. This grew to tens of thousands in the
late summer/fall, when the Alpha variant emerged in the UK [7]. By the end of
2020 it was hundreds of thousands, and in early 2021 it had reached 1 million;
today it is over 2.5 million. This will likely continue to increase
exponentially (see [8]) as many countries [9, 10] ramp up their sequencing
infrastructure for COVID-19 and future pandemics.

The available SARS-CoV-2 virome sequence data, in databases such as GISAID,
has a high _variety_ , since it comprises sequences from all over the world.
Since COVID-19 has spread all over the world for more than a year now, and
viruses continue to mutate over time, there are quite a number of variants of
the SARS-CoV-2 virome, and they continue to emerge. Because variants continue
to emerge and die off, some epidemiologists have even proposed a dynamic
nomenclature system similar to that used for the common cold or flu [11]. We
use this so-called “Pango Lineage” nomenclature to identify the variants we
study here, since only the very common variants of concern (VoCs) are named.
Examples of such named VoCs are the Alpha [7] (Pango Lineage B.1.1.7), Gamma
[12] (P.1) and Delta [13] (B.1.617.2) variants (see Table I for a more
complete list). The genomic variations (most of which happen in the spike
region, see Figure 1) that define these different variants have been
associated with increased transmissibility [14], and immune evasion [15].

Because databases such as GISAID collect sequences from all over the world,
they come from heterogeneous sources of sequencing technologies and centers,
leading to multiple levels of _veracity_. However, the largest source of
different veracity in the data is the widely varying degree to which different
populations are represented. For example, the UK sequences about 5% of its
population of $\approx 70$ million, the USA sequences about 1% of its
population of $\approx 300$ million, while India sequences only a fraction of
a percent of its population of $\approx 1.3$ billion [16, 17] (see Figure 3).
Because of this, for example, even though the Delta variant likely originated
in India, the majority of the available sequences of this variant are from the
UK and the USA, after it arrived in these countries [18].

Since the genomic sequence of a virus encodes all of its functions such as
virulence and transmissibility, the _value_ of such massive amounts of genomic
data is clear. It is variation in this genomic sequence itself which defines
the different variants of SARS-CoV-2 such as Alpha, Delta and Gamma. All of
these variants differ from each other in effect (due to their unique genomic
variations), yet they all descend from the original SARS-CoV-2 sequence [19].
It is only through a process of evolution and transmission to many parts of
the globe for over a year, has it diverged to this extent. The amount of
sequence data available today puts us in the age of _genomic surveillance_ :
tracking the spread of pathogens in terms of genomic content [20, 21].

Approaches for rapidly clustering and classifying sequences will be crucial in
these genomic surveillance efforts. A clustering method, when applied to the
data on a daily basis, for example, would identify a new and rapidly emerging
variant in terms of a cluster which grows abnormally quickly, allowing
scientists to focus on this cluster. Classification, on the other hand, would
allow us to track the spread of known variants in new municipalities, regions,
countries and continents. For example, the USA had a wave of the Alpha variant
from the UK in early 2021, and later, a wave of the Delta variant from India
and/or via other intermediaries, such as the UK (see Figure 5). Such patterns
of spread can reveal information about the underlying transmission networks
between different countries (the UK and USA, or the UK and India), or even
parts of different countries. This can help overcome some of the different
veracity in the data, such as the widely varying degree to which different
countries are represented in terms of sequencing data, due to sampling bias.
For example, even though India is very under-sampled compared to the UK, the
wave of the Delta variant in the UK, along with information about flights from
India to the UK near the beginning of this wave could give us insights on how
the Delta variant originated and spread in India.

The development of clustering and classification approaches needs several
important considerations, however. For one, the number of sequences is so huge
that any way of extracting useful features becomes even more critical. Since
the spike protein is the entry point of the virus to the host cell, it is an
important characterizing feature of a coronavirus [22, 23]. Most of the
variants of SARS-CoV-2 are characterized by mutations which happen
disproportionately in the spike region of the genome [7, 12, 13]. Even the
mRNA vaccines (_e.g._ , Pfizer and Moderna) for COVID-19 are designed to
encode/target only the SARS-CoV-2 spike protein [24] (unlike traditional
vaccines which comprise an entire virome). Since the spike region is
sufficient to characterize most of the important features of a viral sample,
yet much smaller in length, as depicted in Figure 1, we focus on an embedding
approach tailored to the spike region of the sequences.

Figure 1: The SARS-CoV-2 genome is composed of $\approx 30$Kb nucleotide base
pairs, which codes for several proteins, including the spike protein. The
region of the genome which corresponds to the spike protein is composed of
3821 ($25,384-21,563$) nucleotide base pairs, hence $3821$ (+ 1 stop codon *)
$/3=1274$ amino acids.

Previously, some efforts has been done to perform classification and
clustering of SARS-CoV-2 spike sequences [25, 26, 27]. However, those methods
are not scalable to the amount of data we use in this study. Although they
were successful in getting higher predictive accuracy, it is not clear if the
proposed methods are robust and will give same predictive performance on
larger datasets. In this paper, we propose Spike2Vec, an efficient and
scalable feature vector generation approach for SARS-CoV-2 spike sequences, to
which we can apply different machine learning tasks downstream, such as
classification and clustering. Our contributions in this paper are as follows:

  1. 1.

We propose an embedding approach, called Spike2Vec that outperforms the
baseline classification method in terms of predictive accuracy.

  2. 2.

We show that our method is scalable on larger datasets by using $\approx 2.5$
million spike sequences.

  3. 3.

We prove from the results that the machine learning models used in [25, 26,
27] are not scalable on these larger datasets. This robust checking help us to
analyze the machine learning models in detail in terms of their
appropriateness for SARS-CoV-2 spike sequences.

  4. 4.

We also show that in terms of clustering, our embedding approach is better
than the baseline model.

The rest of the paper is organized as follows: Section II contains a
discussion on the previous studies related to our research problem. Section
III contains a detailed description of our Spike2Vec approach. Section IV
contains the implementation details of the experimental evaluation of
Spike2Vec, along with the dataset statistics and discussion of the baseline
models. We present and discuss the results of this experimental evaluation in
Section V. Finally, we conclude our paper in Section VI.

##  II Literature Review

Because of the rapid spread of COVID-19 since December 2019, a lot of sequence
data is available for this virus. This new source of information has attracted
researchers from all the fields to perform analysis on this data to better
understand the diversity and dynamics of this virus. Authors in [27] propose a
one-hot encoding based approach to classify the different coronavirus hosts
using spike sequences alone. They shows that near-optimal prediction accuracy
can be achieved by considering only the spike portion of the genome sequence
rather than using the whole sequence. Ali  _et al._ , in [25] perform
classification of different variants of the human SARS-CoV-2. Although they
were successful in achieving higher accuracy than in [27], the kernel method
used in their approach, however, is not scalable to the size of the data we
use in this study. This drawback makes it difficult to use this approach in
real-world scenario, _i.e._ , the current scenario.

Supervised and unsupervised feature selection methods such as ridge regression
[28], lasso regression, and principal component analysis (PCA) [29], etc., are
very popular for not only reducing the runtime but also for improving the
predictive performance of the underlying machine learning algorithms. Authors
in [26] performs clustering on SARS-CoV-2 spike sequences and show that
clustering performance could be improved by simply using lasso and ridge
regression. Although they were also able to get significant improvement in
terms of clustering quality as compared to the baseline, using feature
selection methods like ridge regression and lasso regression scales very
poorly on the larger datasets, such as the one we use in this study. Melnyk
_et al._ , in [30] perform clustering of the entire SARS-CoV-2 genome (rather
than just the spike sequence) using CliqueSNV [31], a method originally
designed for identifying haplotypes in an intra-host viral population.
Although they obtained good overall F1-scores, our (clustering) approach tends
to obtain better overall F1-scores. It would be interesting to know whether
that is because of our feature vector representation, or because we leverage
more (and more up-to-date) data, or both.

Farhan  _et al._ , in [32] propose an efficient approach to compute a
similarity matrix (kernel matrix). The computed kernel matrix is proven to be
efficient for sequence classification. However, since their approach requires
to save an entire $n\times n$ dimensional kernel matrix (where $n$ is the
total number of sequences), this makes their method expensive in terms of
space. Authors in [33] use the random feature method to map the original input
into low dimensional feature space so that the inner product of the low
dimensional data is approximately equals to the inner product of original data
points. Peng  _et al._ , in [34] use the random feature attention model for
text classification. Their approach is linear in terms of runtime and space
and uses random feature methods to approximate the softmax function.

While dealing with _Big Data_ , it is also important to analyze the trade-off
between the prediction accuracy and the runtime [35]. Although Ali  _et al._ ,
in [25] use the kernel method for spike sequence classification, since the
kernel computation is, however, expensive in terms of time and space, their
approach is only a proof of concept, and not feasible in a real-world
scenario.

##  III Proposed Approach

In this section, we give a step by step description of the Spike2Vec approach.
Given the SARS-CoV-2 spike sequences, we first generate $k$-mers so that we
can preserve some ordering information of the sequences. An interesting side
note is that _de novo_ genome assembly (when no reference is present) involves
inferring sequence order by assembling $k$-mers obtained from short reads
[36]. After the $k$-mers are generated, to convert the alphabetical
information of $k$-mers into numerical representation (so that ML algorithms
can be applied), we generate frequency vectors, which count the number of
occurrences of each $k$-mer in the spike sequence. We then map the high
dimensional frequency vectors to low dimensional embedding using an
approximate kernel approach. Each step of Spike2Vec is explained in detail
below.

###  III-A k-mers Generation

The first step of Spike2Vec is to compute all $k$-mers for the spike
sequences. The main idea behind using k-mers is to allow some order
information of the sequence to be preserved. The total number of $k$-mers,
which we can be generated from a spike sequence of length $N$ is:

$N-k+1$ (1)

In our data, the value for $N$ is $1274$. In Equation (1), $k$ is a user-
defined parameter for the size of each mer. See Figure 2 for an example. In
our experiments we use $k=3$; this was decided using a standard validation set
approach [37].

Figure 2: Example of 4-mers of the amino acid sequence “MDPEGRKMLSV”.

###  III-B Frequency Vectors Generation

Since a k-mers is an alphabets-based representation of a spike sequences, we
need to convert the k-mers into a numerical representation. Therefore, we
design feature vectors that contain the frequency/counts of each $k$-mer in
its respective spike sequence. Each sequence $A$ is over an alphabet $\Sigma$.

###### Remark 1

Alphabets in our dataset represents amino acids of the spike sequence.

These fixed length frequency vector have length $|\Sigma|^{k}$, which
represents the number of possible $k$-mers of a spike sequence. Since the
total number of alphabets in our data are $21$ (the number of amino acids),
the length of each frequency vector becomes $21^{3}=9261$.

###  III-C Low dimensional Representation

In large scale machine learning (ML) tasks such as classification and
regression, typical supervised/unsupervised dimensionality reduction methods
such as principal component analysis, ridge regression, and lasso regression,
etc., are not suitable because they take a lot more time to execute.
Therefore, in a real world scenario where we can have a huge amount of data,
the scalability of any underlying algorithm could be one of the major issues.
One option is to use kernel based algorithms that compute a similarity matrix
which can later be used for the underlying ML tasks. To compute the kernel
matrix (gram matrix), the kernel trick is used.

######  Definition III.1 (Kernel Trick)

The Kernel Trick is used to generate features for an algorithm which depends
on the inner product between only the pairs of input data points. The main
idea is to avoid the need to map the input data (explicitly) to a high-
dimensional feature space.

The Kernel Trick relies on the following observation: Any positive definite
function f(a,b), where $a,b\in\mathcal{R}^{d}$, defines an inner product and a
lifting $\phi$ so that we can quickly compute the inner product between the
lifted data points [33]. More formally:

$\langle\phi(a),\phi(b)\rangle=f(a,b)$ (2)

The major drawback of kernel methods is that in case of large training data,
the kernel method suffers from large initial computational and storage costs.

To overcome these computational and storage problems, we use an approximate
kernel method called Random Fourier Features (RFF) [33], which maps the input
data to a randomized low dimensional feature space (euclidean inner product
space). More formally:

$z:\mathcal{R}^{d}\rightarrow\mathcal{R}^{D}$ (3)

In this way, we approximate the inner product between a pair of transformed
points. More formally:

$f(a,b)=\langle\phi(a),\phi(b)\rangle\approx z(a)^{\prime}z(b)$ (4)

In Equation (4), $z$ is low dimensional (unlike the lifting $\phi$). In this
way, we can transform the original input data with $z$, which acts as the
approximate low dimensional embedding for the original data. This low
dimensional representation is then used as an input for different ML tasks
like classification and regression.

##  IV Experimental Evaluation

We now detail the experiments we performed to evaluate Spike2Vec in terms of
both the downstream classification and clustering results obtained.

###  IV-A Experimental Setup

All experiments are conducted using an Intel(R) Xeon(R) CPU E7-4850 v4 @
$2.10$GHz having Ubuntu $64$ bit OS ($16.04.7$ LTS Xenial Xerus) with 3023 GB
memory. Implementation of Spike2Vec is done in Python and the code is
available online for
reproducibility333https://github.com/sarwanpasha/Spike2Vec. Our pre-processed
data is also available
online444https://drive.google.com/drive/folders/1-YmIM8ipFpj-
glr9hSF3t6VuofrpgWUa?usp=sharing, which can be used after agreeing to terms
and conditions of GISAID555https://www.gisaid.org/. For the classification
algorithms, we use $1\%$ data for training and $99\%$ for testing. The purpose
of using smaller training dataset is to show how much performance gain we can
achieve while using minimal training data.

###### Remark 2

Our data split and pre-processing follow those of [25].

###  IV-B Dataset Statistics

We used the (aligned) amino acid sequences corresponding to the spike protein
from the largest known database of SARS-CoV-2 sequences, GISAID. In our
dataset, we have 2,519,386 spike sequences along with the COVID-19 variant
information (in our data, we have $1327$ variants in total) for each spike
sequence. The information about some of the more well-represented variants is
given in Table I. Since most of the variants are new, we do not have all the
information available for all them. Therefore, we put “-” in the field in
Table I for which we do not have any information available online.

Pango Lineage | Region |  Labels |  Num. Mutations S-gene/Genome |  Num. of sequences
---|---|---|---|---
B.1.1.7 | UK [7] |  Alpha |  8/17 |  976077
B.1.351 | South Africa [7] |  Beta |  9/21 |  20829
B.1.617.2 | India [13] |  Delta |  8/17 |  242820
P.1 | Brazil [12] |  Gamma |  10/21 |  56948
B.1.427 | California [38] |  Epsilon |  3/5 |  17799
AY.4 | India [39] |  Delta |  - |  156038
B.1.2 | - |  - |  - |  96253
B.1 |  |  |  |  78741
B.1.177 | - |  - |  - |  72298
B.1.1 | - |  |  - |  44851
B.1.429 | - |  - |  - |  38117
AY.12 | India [39] |  Delta |  - |  28845
B.1.160 | - |  - |  - |  25579
B.1.526 | New York [40] |  Iota |  6/16 |  25142
B.1.1.519 | - |  - |  - |  22509
B.1.1.214 | - |  - |  - |  17880
B.1.221 | - |  - |  - |  13121
B.1.258 | - |  - |  - |  13027
B.1.177.21 | - |  - |  - |  13019
D.2 | - |  - |  - |  12758
B.1.243 | - |  - |  - |  12510
R.1 | - |  - |  - |  10034
TABLE I: The SARS-CoV-2 variants which were represented in more than 10,000
sequences (of the $\approx 2.5$ million sequences). The S/Gen. column
represents number of mutations on the Spike (S) gene / entire genome. Total
number of amino acid sequences in our dataset is 2,519,386. The variants
discussed in this Table comprise 1,995,195 sequences.

Figure 3 shows the total number of spike sequences for the top $10$ countries
worldwide. In our GISAID dataset, a total of $219$ countries are represented.
Since USA has the highest number of spike sequences, we use it as a case study
to analyze the spread patterns of different variants in Section IV-C1.

Figure 3: Country-wise distribution (for the top 10 countries) of spike
sequences.

###  IV-C Data Visualization

To see if there is any (hidden) clustering in the data, we mapped the data to
2D real vectors using the t-distributed stochastic neighbor embedding (t-SNE)
approach [41]. Since it was not possible to run the t-SNE algorithm on all
$\approx 2.5$ million spike sequences, we obtained a representative subset of
sequences containing 7000 randomly selected sequences such that the proportion
of each variant in this subset is equal to its proportion in the original
data. The t-SNE plot for Delta, Beta, Iota, Epsilon, and Gamma variants is
shown in Figure 4.

Figure 4: t-SNE embeddings of spike sequences

####  IV-C1 USA Case Study

Figure 5 shows the COVID-19 spread pattern for three variants in the USA from
March 2020 to July 2021. We can see in Figure 5 that after the coronavirus
spread hit the peak in April $2021$, the number of cases of the coronavirus
started decreasing. That was the point where a significant proportion of the
population of the USA was vaccinated (hence peak spread reduced).

Figure 5: Spread pattern of Alpha (blue line), Delta (orange line), and Gamma
(black line) variants in USA country from March 2020 to July 2021. Y-axis
shows the total number of COVID-19 infected patients.

###  IV-D Performance Evaluation

To evaluate Spike2vec, we perform classification and clustering on the low
dimensional feature vectors that it produces. For the classification task, we
use Naive Bayes, Logistic Regression, and Ridge Classifier. For the clustering
analysis, we use the $K$-means algorithm. To evaluate the clustering
algorithms, we report accuracy, precision, recall, weighted F1, macro F1, ROC-
AUC. We also show the runtime of different classification algorithms. To
evaluate the clustering method, we use weighted F1 score.

Approach |  ML Algo. | Acc. | Prec. | Recall | F1 (weighted) | F1 (Macro) | ROC-AUC | Training runtime (sec.)
---|---|---|---|---|---|---|---|---
OHE |  NB | 0.306370 | 0.583309 | 0.306370 | 0.383127 | 0.179708 | 0.594794 | 566.099
LR | 0.568250 | 0.498771 | 0.568250 | 0.495177 | 0.196033 | 0.576567 | 1309.060
RC | 0.563617 | 0.479080 | 0.563617 | 0.485127 | 0.174295 | 0.566908 | 110.766
Spike2Vec |  NB | 0.420112 | 0.797562 | 0.420112 | 0.521133 | 0.391974 | 0.685318 | 457.5417
LR | 0.688674 | 0.689724 | 0.688674 | 0.649236 | 0.490654 | 0.694196 | 830.6327
RC | 0.675313 | 0.680797 | 0.675313 | 0.629078 | 0.447041 | 0.674559 | 95.7315
TABLE II: Variants Classification Results (1% training set and 99% testing
set) for top $22$ variants ($1995195$ spike sequences) discussed in Table I.
Best values are shown in bold.

###  IV-E Baseline Algorithm

For the baseline approach, we use the one-hot encoding based approach proposed
in [27] (for reference, we call this method as One Hot Embedding “OHE” in rest
of the paper). Authors in [27] use a typical one-hot encoding approach to
convert the spike sequences into numerical representations. In spike
sequences, we have $21$ unique amino acids (unique alphabets forming $\Sigma$)
namely “ACDEFGHIKLMNPQRSTVWXY”. Also, the length of each spike sequence is
$1273$ plus an ending character $*$ at the $1274^{th}$ location. After getting
the one-hot encoding based numerical representation for each spike sequence,
we will get a feature vector of length $26,733$ corresponding to each spike
sequence ($21\times 1273=26,733$). After getting this numerical
representation, authors in [27] use the typical principal component analysis
(PCA) approach to reduce the dimensionality of the data. In our case, since
the size of our data is too huge, simply using PCA would take a lot of
computational time. Therefore, in OHE, we use RFF on the one-hot embeddings to
get the low dimensional feature vector representations for our baseline model.

###### Remark 3

If we directly apply the classification algorithms on the one-hot embedding,
the underlying classification algorithms simply do not run and exceed the
amount of available memory ($\approx$3TB). This shows that scalability is the
major problem with these typical ML algorithms.

###  IV-F Machine Learning Models

For the classification task, we use Naive Bayes (NB), Logistic Regression
(LR), and Ridge Classifier (RC) [42] with default parameters. For clustering,
we use simple $k$-means algorithm with default parameters.

###### Remark 4

For $k$-means, we use $22$ as the number of clusters. We selected the number
of clusters using the Elbow method [26]. To do this, we perform clustering
with different number of clusters ranging from $2$ to $100$ and then see the
trade-off between the sum of squared error (distortion score) and the runtime.
After analyzing the trade-off, we use “knee point detection algorithm (KPDA)”
[43] to find the optimal value of $k$ for the k-means algorithm.

##  V Results and Discussion

In this section, we first present the classification results for Spike2Vec and
show that it significantly outperforms the baseline in terms of prediction
accuracy. We then show the results for the clustering using weighted F1. In
the end, we report the importance of each amino acid by computing information
gain.

###  V-A Classification Results

Results for different classification algorithms are shown in Table II. We can
observe that overall Logistic regression is almost a clear winner in case of
Spike2vec. All the classifiers in case of Spike2Vec clearly outperform the
corresponding classifiers with OHE. This performance for different evaluation
metrics shows the effectiveness of using $k$-mers instead of one hot encoding
for feature vector representation of the spike sequences. Also, we can observe
that although the performance of RC in case of Spike2Vec is not better than
LR, it is significantly better than LR and NB however in terms of training
runtime. Therefore, we can conclude that overall LR is better in terms of
prediction performance with Spike2Vec while RC is better in terms of runtime
along with comparable performance to LR.

###  V-B Clustering Results

We also test the performance of Spike2Vec using the $k$-means clustering
method. Results for the clustering methods are shown in Table III. We can
observe that Spike2Vec clearly outperforms the baseline models in case of all
but $1$ variant. The reason for the bad performance in case of Beta and
Epsilon variants is due to the fact that they are in comparatively less
proportion in the dataset (see Table I). Because of the less information,
Spike2Vec is not able to design rich feature vector representation for these
variants.

| F1 Score (Weighted) for Different Variants
---|---
Methods | Alpha | Beta | Delta | Gamma | Epsilon
OHE | 0.0410 | 0.0479 | 0.5942 | 0.6432 | 0.0571
Spike2Vec | 0.9997 | 0.0300 | 0.8531 | 0.9680 | 0.2246
TABLE III: F1 score by applying the $k$-means clustering algorithm on all
$1327$ variants (2519386 spike sequences) in the GISAID dataset. Best values
are shown in bold.

The contingency table for some of the more well represented variants computed
using Spike2Vec is given in Table IV. Since we have a total of $1327$
variants, it is not possible to show the contingency table for all of the
variants. Therefore, we only show relationship between variants and clusters
for some of the popular variants in Table IV.

| K-means (Cluster IDs)
---|---
Variant |  0 |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 |  10 |  11 |  12 |  13 |  14 |  15 |  16 |  17 |  18 |  19 |  20 |  21
Epsilon |  109 |  0 |  3186 |  432 |  113 |  219 |  67 |  0 |  160 |  9 |  134 |  0 |  18 |  14 |  0 |  78 |  9 |  3792 |  113 |  0 |  48 |  41
Alpha |  6061 |  1 |  175923 |  23353 |  5846 |  11754 |  3376 |  0 |  9466 |  1041 |  6889 |  0 |  734 |  1281 |  0 |  4160 |  329 |  205848 |  5730 |  0 |  3193 |  2136
Gamma |  344 |  0 |  10403 |  1312 |  327 |  686 |  205 |  0 |  534 |  63 |  390 |  0 |  52 |  77 |  0 |  254 |  14 |  11977 |  324 |  0 |  182 |  137
Beta |  144 |  0 |  3853 |  436 |  115 |  237 |  64 |  0 |  191 |  8 |  148 |  0 |  19 |  25 |  0 |  81 |  7 |  4435 |  119 |  0 |  71 |  44
Delta |  1432 |  1 |  43691 |  5732 |  1391 |  2832 |  831 |  0 |  2342 |  241 |  1777 |  0 |  172 |  315 |  0 |  1016 |  77 |  51596 |  1400 |  0 |  836 |  541
TABLE IV: Contingency tables of variants vs clusters.

###  V-C Importance of Each Amino Acid

Authors in [25] compute Information Gain (IG) between the variants and each
amino acid (attribute) separately. The main goal for computing IG was to see
the importance of each amino acid. Since they are using very small dataset
($7000$ spike sequences), it is not clear if the information gain values they
computed will be the same for our $\approx 2.5$ million sequences. Therefore,
a “proof of concept” is needed to verify their results (on a larger dataset).
The IG is defined as follows:

$IG(Class,position)=H(Class)-H(Class|position)$ (5) $H=\sum_{i\in
Class}-p_{i}\log p_{i}$ (6)

where $H$ is the entropy, and $p_{i}$ is the probability of the class $i$. We
extracted $10$ sample datasets of $20000$ spike sequences each and computed IG
for each of the dataset separately. For each of the dataset, we got same amino
acids with the maximum IG values. The values for different amino acids for one
of the $10$ datasets is given in Figure 6. The Centers for Disease Control and
Prevention (CDC) located in USA pointed out mutations at a few positions that
take place in different variants [39]. We compare the mutation information
from the CDC with the (high) IG values that we got for different amino acids.
According to the CDC, R452L is present in Epsilon and Delta lineages and sub-
lineages while K417N, E484K, and N501Y substitutions are present in the Beta
variant.

###### Remark 5

Note that R452L means Amino acid at position 452 was ‘R’ before and after
mutation, it changed to ‘L’.

Similarly, K417T, E484K, and N501Y substitutions are present in the Gamma
variant [39]. We can see in Figure 6 that we obtained the maximum IG values
for the same amino acids positions mentioned by CDC. We also made the IG
values for all amino acids available
online666https://github.com/sarwanpasha/Spike2Vec/blob/main/correlation˙data.csv.

Figure 6: Information gain for each amino acid position with respect to
variants. The $x$-axis corresponds to amino acid positions in the spike
sequences.

##  VI Conclusion

We propose an efficient and scalable embedding approach in this paper that can
be used to perform different machine learning tasks on the SARS-CoV-2 spike
sequences. We show that our model can scale to several million sequences, and
it also outperforms the baseline models significantly. Since the COVID-19
disease is relatively new, we do not have enough information available for
different coronavirus variants so far. We will explore the new (and existing)
variants in more detail in the future. We will also use deep learning models
to enhance the prediction performance of Spike2Vec. Using adversarial examples
to test the robustness of the ML models in case of spike sequences is another
potential future extension.

##  VII Acknowledgments

This research was supported by a Georgia State University startup grant.

## References

  * [1] R. Sharma, A. Mateush, and J. Übi, “Tale of three states: analysis of large person-to-person online financial transactions in three baltic countries,” in _IEEE BigData_ , 2019, pp. 1497–1505.
  * [2] Y. Dong, D. Yan, A. I. Almudaifer, S. Yan, Z. Jiang, and Y. Zhou, “BELT: A pipeline for stock price prediction using news,” in _IEEE BigData_ , 2020, pp. 379–410.
  * [3] J. Kremer, K. Stensbo-Smidt, F. Gieseke, K. S. Pedersen, and C. Igel, “Big universe, big data: Machine learning and image analysis for astronomy,” _IEEE Intelligent Systems_ , vol. 32, no. 2, pp. 16–22, 2017.
  * [4] M. Mace, “Comet NEOWISE sizzles as it slides by the sun, providing a treat for observers,” _Infrared Processing and Analysis Center_ , 2020.
  * [5] C. K. Leung, Y. Chen, C. S. H. Hoi, S. Shang, and A. Cuzzocrea, “Machine learning and olap on big covid-19 data,” in _IEEE BigData_ , 2020.
  * [6] C. K. Leung, Y. Chen, S. Shang, and D. Deng, “Big data science on covid-19 data,” in _IEEE BigData_ , 2020.
  * [7] S. Galloway _et al._ , “Emergence of sars-cov-2 b. 1.1. 7 lineage,” _Morbidity and Mortality Weekly Report_ , vol. 70, no. 3, p. 95, 2021.
  * [8] S. ZD, L. SY, F. F, C. RH, Z. C, E. MJ _et al._ , “Big data: Astronomical or genomical?” _PLoS Biology_ , vol. 13, no. 7, p. e1002195, 2015.
  * [9] Genome Web, https://www.genomeweb.com/infectious-disease/cdc-commits-90m-create-public-health-pathogen-genomics-research-centers, [Online; accessed 5-September-2021].
  * [10] Reuters, https://www.reuters.com/article/us-china-genomics-state/chinese-state-fund-invests-in-gene-firm-bgi-idUSKBN2AM0AT, [Online; accessed 5-September-2021].
  * [11] A. Rambaut, E. C. Holmes, A. O’Toole, V. Hill, J. T. McCrone, C. Ruis, L. du Plessis, and O. G. Pybus, “A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology,” _Nature Microbiology_ , vol. 5, no. 11, pp. 1403–1407, 2020.
  * [12] F. Naveca _et al._ , “Phylogenetic relationship of sars-cov-2 sequences from amazonas with emerging brazilian variants harboring mutations e484k and n501y in the spike protein,” _Virological. org_ , vol. 1, 2021.
  * [13] P. Yadav _et al._ , “Neutralization potential of covishield vaccinated individuals sera against b. 1.617. 1,” _bioRxiv_ , vol. 1, 2021.
  * [14] E. Volz, S. Mishra, M. Chand, J. C. Barrett, R. Johnson, L. Geidelberg, W. R. Hinsley, D. J. Laydon, G. Dabrera, Á. O’Toole _et al._ , “Assessing transmissibility of sars-cov-2 lineage b. 1.1. 7 in england,” _Nature_ , vol. 593, pp. 266–269, 2021.
  * [15] M. McCallum, J. Bassi, A. Marco, A. Chen, A. Walls, J. Iulio, M. Tortorici, M. Navarro, C. Silacci-Fregni, C. Saliba, M. Agostini, D. Pinto, K. Culap, S. Bianchi, S. Jaconi, E. Cameroni, J. Bowen, S. Tilles, M. Pizzuto, S. Guastalla, G. Bona, A. Pellanda, C. Garzoni, W. Van Voorhis, L. Rosen, G. Snell, A. Telenti, H. Virgin, L. Piccoli, D. Corti, and D. Veesler, “Sars-cov-2 immune evasion by variant b.1.427/b.1.429,” vol. 1, 2021, doi:10.1101/2021.03.31.437925.
  * [16] Y. Furuse, “Genomic sequencing effort for sars-cov-2 by country during the pandemic,” _International Journal of Infectious Diseases_ , vol. 103, pp. 305–307, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1201971220325571
  * [17] A. Maxmen, “Why US coronavirus tracking can’t keep up with variants,” _Nature_ , vol. 592, 2021.
  * [18] GISAID Website, https://www.gisaid.org/, [Online; accessed 5-September-2021].
  * [19] F. Wu, S. Zhao, B. Yu, Y.-M. Chen, W. Wang, Z.-G. Song, Y. Hu, Z.-W. Tao, J.-H. Tian, Y.-Y. Pei _et al._ , “A new coronavirus associated with human respiratory disease in china,” _Nature_ , vol. 579, pp. 265–269, 2020.
  * [20] F. Wu, A. Xiao, J. Zhang, K. Moniz, N. Endo, F. Armas, M. Bushman, P. R. Chai, C. Duvallet, T. B. Erickson, K. Foppe, N. Ghaeli, X. Gu, W. P. Hanage, K. H. Huang, W. L. Lee, K. A. McElroy, S. F. Rhode, M. Matus, S. Wuertz, J. Thompson, and E. J. Alm, “Wastewater surveillance of sars-cov-2 across 40 u.s. states from february to june 2020,” _Water Research_ , vol. 202, p. 117400, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0043135421005984
  * [21] J. Gardy and N. Loman, “Towards a genomics-informed, real-time, global pathogen surveillance system,” _Nature Reviews Genetics_ , vol. 19, pp. 9–20, 2018.
  * [22] F. Li, “Structure, function, and evolution of coronavirus spike proteins,” _Annual Review Virology_ , vol. 3, no. 1, pp. 237–261, 2016.
  * [23] A. Walls, Y. Park, and M. Tortorici, “Structure, function and antigenicity of the sars-cov-2 spike glycoprotein,” _Cell_ , vol. 181, no. 2, pp. 281–292, 2020.
  * [24] K. S. Corbett, D. K. Edwards, S. R. Leist _et al._ , _Nature_ , vol. 586, pp. 567–571, 2020.
  * [25] S. Ali, B. Sahoo, N. Ullah, A. Zelikovskiy, M. Patterson, and I. Khan, “A k-mer based approach for sars-cov-2 variant identification,” _To Appear at: International Symposium on Bioinformatics Research and Applications (ISBRA)_ , 2021.
  * [26] S. Ali, T. e Ali, M. A. Khan, I. Khan, and M. Patterson, “Effective and scalable clustering of sars-cov-2 sequences,” _To appear at: International Conference on Big Data Research (ICBDR)_ , 2021.
  * [27] K. Kuzmin _et al._ , “Machine learning methods accurately predict host specificity of coronaviruses based on spike sequences alone,” _Biochemical and Biophysical Research Communications_ , vol. 533, pp. 553–558, 2020.
  * [28] S. Ali, S. Ciccolella, L. Lucarella, G. Della Vedova, and M. D. Patterson, “Simpler and faster development of tumor phylogeny pipelines,” _To appear at: Journal of Computational Biology (JCB)_ , 2021.
  * [29] S. Ali, H. Mansoor, N. Arshad, and I. Khan, “Short term load forecasting using smart meter data,” in _International Conference on Future Energy Systems (e-Energy)_ , 2019, pp. 419–421.
  * [30] A. Melnyk, F. Mohebbi, S. Knyazev, B. Sahoo, R. Hosseini, P. Skums, A. Zelikovskiy, and M. D. Patterson, “From alpha to zeta: Identifying variants and subtypes of sars-cov-2 via clustering,” _To appear at: Journal of Computational Biology (JCB)_ , 2021.
  * [31] S. Knyazev, V. Tsyvina, A. Shankar, A. Melnyk, A. Artyomenko, T. Malygina, Y. B. Porozov, E. M. Campbell, S. Mangul, W. M. Switzer _et al._ , “Accurate assembly of minority viral haplotypes from next-generation sequencing through efficient noise reduction,” _Nucleic Acids Research_ , p. 264242, 2020.
  * [32] M. Farhan, J. Tariq, A. Zaman, M. Shabbir, and I. Khan, “Efficient approximation algorithms for strings kernel based sequence classification,” in _Advances in neural information processing systems (NeurIPS)_. ., 2017, pp. 6935–6945.
  * [33] A. Rahimi, B. Recht _et al._ , “Random features for large-scale kernel machines.” in _NIPS_ , vol. 3, no. 4, 2007, p. 5.
  * [34] H. Peng, N. Pappas, D. Yogatama, R. Schwartz, N. A. Smith, and L. Kong, “Random feature attention,” in _International Conference on Learning Representations (ICLR)_ , 2021.
  * [35] S. Ali, “Cache replacement algorithm,” _arXiv preprint arXiv:2107.14646_ , 2021.
  * [36] R. Rizzi, S. Beretta, M. Patterson, Y. Pirola, M. Previtali, G. D. Vedova, and P. Bonizzoni, “Overlap graphs and de Bruijn graphs: data structures for de novo genome assembly in the big data era,” _Quantitative Biology_ , vol. 7, no. 4, pp. 278–292, 2019.
  * [37] P. Devijver and J. Kittler, “Pattern recognition: A statistical approach,” in _London, GB: Prentice-Hall_ , 1982, pp. 1–448.
  * [38] W. Zhang _et al._ , “Emergence of a novel sars-cov-2 variant in southern california,” _Jama_ , vol. 325, no. 13, pp. 1324–1326, 2021.
  * [39] SARS-CoV-2 Variant Classifications and Definitions, https://www.cdc.gov/coronavirus/2019-ncov/variants/variant-info.html, [Online; accessed 5-September-2021].
  * [40] A. West Jr _et al._ , “Detection and characterization of the sars-cov-2 lineage b. 1.526 in new york,” _bioRxiv_ , 2021.
  * [41] L. Van der M. and G. Hinton, “Visualizing data using t-sne.” _Journal of Machine Learning Research (JMLR)_ , vol. 9, no. 11, 2008.
  * [42] A. Singh, B. S. Prakash, and K. Chandrasekaran, “A comparison of linear discriminant analysis and ridge classifier on twitter data,” in _2016 International Conference on Computing, Communication and Automation (ICCCA)_ , 2016, pp. 133–138.
  * [43] V. Satopaa, J. Albrecht, D. Irwin, and B. Raghavan, “Finding a” kneedle” in a haystack: Detecting knee points in system behavior,” in _International conference on distributed computing systems workshops_. IEEE, 2011, pp. 166–171.