# What factors influence the popularity of user-generated text in the creative
domain? A case study of book reviews

Salim Sazzed
Department of Computer Science,Old Dominion University, Norfolk, VA USA
Department of Computer Science,University of Memphis, Memphis, TN, USA
<EMAIL_ADDRESS>
###### Abstract

This study investigates a range of psychological, lexical, semantic, and
readability features of book reviews to elucidate the factors underlying their
perceived popularity. To this end, we conduct statistical analyses of various
features, including the types and frequency of opinion and emotion-conveying
terms, connectives, character mentions, word uniqueness, commonness, and
sentence structure, among others. Additionally, we utilize two readability
tests to explore whether reading ease is positively associated with review
popularity. Finally, we employ traditional machine learning classifiers and
transformer-based fine-tuned language models with n-gram features to
automatically determine review popularity. Our findings indicate that, with
the exception of a few features (e.g., review length, emotions, and word
uniqueness), most attributes do not exhibit significant differences between
popular and non-popular review groups. Furthermore, the poor performance of
machine learning classifiers using the word n-gram feature highlights the
challenges associated with determining popularity in creative domains.
Overall, our study provides insights into the factors underlying review
popularity and highlights the need for further research in this area,
particularly in the creative realm.

##  I Introduction

User-generated content has an indispensable role in circulating information
regarding products and various other types of entities on social networks and
e-commerce sites [1, 2]. As predicting review usefulness has practical
implications, in this study, we scrutinize how various kinds of features of a
review are connected to its usefulness. However, instead of exploring review
usefulness, or, more precisely helpfulness for the e-commerce domain (e.g.,
product review) that has been studied by many earlier studies [2, 3, 4, 5],
here, we analyze book reviews that represent the creative domain. Note that
exist broad distinctions between helpful product reviews in the e-commerce
domain and the popularity of the book review. Although some elements may
overlap in both types, the core components of book reviews are likely to be
different than those of e-commerce product reviews. In the former, a reviewer
narrates and describes his feelings and opinions toward various aspects of the
book, such as the story, plot, and characters influenced by the reviewer’s
writing style (i.e., elements of creativity). In contrast, a product review
depicts the features and functionality of the product and may even provide a
comparison with the other products, which is highly unlikely in book reviews.
On the feedback level, in the e-commerce domain, the potential customer may
label a review as helpful if the product pros or cons are available in the
review, which is not the point for the creative domain; rather, the writing
style may be concomitant to the popularity.

As exploring the factors underlying the perceived popularity of user-generated
text in the creative domain can provide insights for content creators and
platforms, contribute to knowledge in the creative realm, and identify
challenges that warrant further research, here we explore a set of book
reviews. By analyzing a set of 2259 popular and 2555 non-popular book reviews,
we aim to uncover any association between the popularity of a review and its
content features. Note that the popularity of a review specifies how a user
review is accepted by fellow readers, whether they like it or not. Various
factors, including review helpfulness, clarity, readability, and information
and sentiment conveyed, among others, may help determine whether an individual
likes a particular review or not. We examine various affective and linguistic
features and perform readability tests. Various features, such as the presence
of sentiment and emotion, usage of common and unique words, connectives, and
mentions of characters, are studied. Besides, we analyze the review length and
readability of the reviews using two well-known readability tests, FRE (Flesch
Reading Ease) and FKR (Flesch–Kincaid readability). Finally, we employ several
traditional machine learning classifiers and transformer-based fine-tuned
language models for classifying reviews into popular and non-popular groups.
We observe that although some of the features, such as emotion or review
length, demonstrate some distinguishing capabilities, most of them do not
exhibit any noticeable difference between the two groups in seclusion. The FRE
and FKG readability tests, which are based on statistical syntactic
attributes, do not reveal any discernible variations between the two groups.
Classifiers show poor performances utilizing unigram and bigram-based word
features for this binary-level popularity classification task (macro F1 scores
between 0.60 - 0.65), which suggests determining the popularity of book
reviews is a highly complex task.

Contributions:

The main contributions of this study can be summarized as follows-

  * •

We investigate diverse types of textual attributes in book reviews to find
their relationships with attained popularity.

  * •

We employ two readability tests, FRE and FKG, to find whether reading ease is
connected to the popularity of the book review.

  * •

We study the performance of various machine learning classifiers and
transformer-based language models for the automatic classification of popular
and non-popular book reviews.

##  II Related Work

Various aspects of user-generated reviews, including helpfulness, sentiment,
and linguistic traits, have been investigated in many earlier studies [2, 6,
7, 4, 8, 9]. Similar to review helpfulness, review popularity is crucial for
understanding the significance of reviews. Although the two terms, review
helpfulness and review popularity, are closely related concepts, differences
exist between them [10]. The various attributes of reviews may have
contrasting relationships with review popularity and helpfulness. As mentioned
earlier, existing studies related to review content analysis mainly focused on
predicting review helpfulness in the e-commerce domain [2, 3, 4, 5]. The
findings of these works are often inconsistent, as reported by Hong et al.
[2]. Lutz et al. [3] found that longer reviews (e.g., Amazon customer reviews)
with frequent shifts between positive and negative arguments are perceived as
less helpful. Malik and Hussain [5] utilized six machine learning classifiers
and three real-life Amazon review datasets in their study. They discovered
that the number of syllables, spaces, auxiliary verbs, and driving words in
review text and reviewer’s productivity score are effective predictors of a
review helpfulness. Huang et al. [4] discovered a word count threshold beyond
which the word count negatively affects the helpfulness of a review.
Krishnamoorthy [11] used linguistic features from reviews, alongside metadata,
subjectivity, and readability features, to predict helpfulness. Experiments on
two real-world datasets showed the efficacy of proposed approach. Eslami et
al. [12] noticed that most helpful online reviews are often associated with
medium length, lower review scores, and negative or neutral argument frames.
Moreover, their research revealed that review length is one of the most
important factors in predicting the helpfulness of an online consumer review.
Hong et al. [2] performed a meta-analysis on the findings of the influence of
determinants on review helpfulness. The authors found most of the review-
related determinants (e.g., review depth, review readability, linear review
rating, review age) and two reviewer-related determinants (i.e., reviewer
information disclosure and reviewer expertise) have inconsistent conclusions(
i.e., perceived review helpfulness) based on existing literature. They found
several factors such as the criteria for helpfulness determination, the source
of the reviews, and varied product domains are responsible for most of the
mixed findings. Chua and Banerjee [13] investigated the relationship between
the review helpfulness with respect to reviewer reputation, review rating, and
review depth. Utilizing Amazon review data, the authors found a positive
correlation between the review helpfulness and the reviewer profile and review
depth. However, when rating is concerned, they observed negative relations.
Besides, they also noticed reviews having harsh criticism and adequate depth
are usually received positively by the user.

##  III Dataset

The review corpus used here consists of a collection of Goodreads book reviews
obtained from the Kaggle website
111https://www.kaggle.com/datasets/notkrishna/top-100-science-fiction-books-
and-their-reviews. Originally, it consists of around 10,000 reviews written
towards top 100 science fiction books (ranked based on Goodreads ratings).
However, upon examination, we find that a significant number of reviews
comprise various types of issues, including missing values for the review text
and like count, as well as duplicate sentences. We remove duplicate sentences
from the reviews and exclude problematic reviews with missing fields. After
clean-up, the corpus consists of 8845 reviews
222https://github.com/sazzadcsedu/ReviewPopularity.

As each book has different distributions of the like count, it is more
reasonable to decide the popularity of a review based on the like count
distribution of the individual book than considering all the books in the
corpus. We split reviews of each book into four quartiles (Q1-Q4) based on the
like count. The reviews of Q4 (i.e., like count above 75% percentiles) and Q1
(i.e., like count below 25% percentiles) are considered popular and non-
popular reviews, respectively. Our initial analysis (not reported here)
incorporating reviews of Q3 and Q4 quartiles in the popular group and Q1 and
Q2 in the non-popular group, in general, do not yield any distinguishing
patterns. Therefore, in this study, we exclude reviews from the Q2 and Q3
quartiles to investigate whether a large margin between popular and non-
popular categories indeed yields better distinguishing signals. The final
dataset contains 2259 popular (all reviews from Q4 ) and 2555 non-popular
reviews (all reviews from Q1). A few examples are shown in Fig 1.

Figure 1: Examples of books reviews from popular (Q4) and non-popular groups
(Q1)

##  IV Feature Analysis

We consider various types of features in our analysis, including lexical,
affective, discourse, and readability metrics that have been explored in prior
studies for related research [14, 15, 16].

###  IV-A Length Feature

The length-related statistics provide insight into various textual aspects,
such as the level of detail, the level of engagement, and the complexity of a
review [4]. We compute the distribution of the following length-related
features: i) review length (# sentences), ii) review length (# words), and
iii) sentence length (# words) in popular and non-popular groups to determine
the similarities and differences.

###  IV-B Affective Feature

Psychological features, specifically affective features in a text, encompass
attributes related to emotions and sentiments. We examine the distributions of
two types of affective features—sentiment and emotion—in two groups, as
earlier studies have indicated associations with product helpfulness [17].

####  Presence of Sentiment

We scrutinize the sentiment in the reviews of two groups, such as the presence
of sentiment and opinion words based on two popular English sentiment
lexicons: Opinion Lexicon [18] and VADER [19]. This assessment unveils how
sentiments and opinions influence the popularity of user-provided feedback.

#### Presence of Emotion

To recognize the presence of emotion in the text, we utilize the NRC Emotion
Lexicon [20]. The following six types of emotions (i.e., anger, anticipation,
surprise, sadness, joy, and disgust) are investigated in the reviews of two
groups.

###  IV-C Lexical Feature

#### Presence of Capitalization

Capitalization is a common way to strengthen an expression or incorporate a
louder tone within the text [21] and labels certain parts of an utterance as
more salient than others [22].

#### Presence of Common Words

The coverage of two common wordlists, one with 10,000 most common English
words (large list) and the other with 5,000 most common English words (small
list) are computed for the popular and non-popular review groups.
333https://github.com/sazzadcsedu/ReviewPopularity. The purpose is to find
whether the presence of rarely used words inflicts any negative consequence on
the popularity of the reviews.

#### Lexical Diversity

Lexical diversity investigates the presence of unique words in a text corpus
[23]. We compute the ratios of distinct words in each review for both groups
and then report the mean, median, and standard deviations. The goal is to find
whether lexical diversity has any positive or negative influence on the
popularity of the reviews.

###  IV-D Discourse Feature

We further examine various discourse features, such as connectives and
conjunctives, that are related to the coherence and cohesion of the text.

#### Connectives

The use of connecting words creates cohesive links between ideas and clauses
and provides clues about text organization [24]. We consider the following
five types of connectives-

  * •

Adding (ADD): and, also, as well as, moreover, furthermore, besides, in
addition.

  * •

Cause and Effect (CE): because, so, therefore, thus, consequently, as a result
of.

  * •

Emphasizing (EMP): above all, in particular, especially, significantly,
indeed, notably, most of all.

  * •

Illustrating (ILL): for example, such as, for instance, as revealed by, in the
case of, as shown by.

  * •

Sequence (SEQ): next, then, first, second, finally, meanwhile, after.

#### Subordinating Conjunctions

In addition, we investigate the presence of subordinating conjunctions that
are related to the complex sentence. Complex sentences are usually more
difficult to process than simple sentences; nevertheless, they convey a clear
and informative message. A list of 50 commonly used subordinating conjunctions
is considered in this study
444https://github.com/sazzadcsedu/ReviewPopularity.

###  IV-E Semantic Feature

#### Presence of Named Entity

A named entity (NE) represents a real-world object, such as the name of a
person or organization. The three most common types of NE are person,
location, and organization. We analyze the presence of these three types of NE
in the reviews, mentions of the character names (PER), geographical entities
(GE) such as the mention of countries, cities, or similar things, and
organizations (ORG). The spaCy [25] library is used to identify the presence
of NE in the text.

##  V Readability

Readability refers to the reading easiness and understandability of a piece of
text. Readability is related to the cognitive effort that is needed by an
individual to understand and comprehend a piece of text [26]. A text with high
readability improves the likelihood that the reader will clearly understand
the thoughts and ideas conveyed. Previous studies yielded mixed findings
regarding the relationship between readability and the effectiveness of
product reviews. For instance, Liu et al. [27] found positive effects of
readability on the usefulness of product reviews, whereas Hong et al. [2]
reported no significant association between these two factors. Thus, in this
study, we aim to investigate whether readability influences the popularity of
reviews in the creative domain. Two well-known readability tests, FRE and FKG,
are applied to assess the reading easiness of the reviews and their connection
with popularity. FRE provides readability scores on a scale from 1 to 100,
where a higher score indicates that a chunk of English text is comparatively
easier to read. The FKG level, on the other hand, is equivalent to the US
grade level of education.

##  VI  Feature Significance Test and Classification

###  VI-A Mann-Whitney U Test for Significance

In addition to providing quantitative statistics of various attributes, we
perform the significance test to see whether the differences are meaningful
for the likely candidates. We apply the Mann-Whitney U Test to determine if
there is a significant difference between the two groups, popular and non-
popular, based on individual attributes. The Mann-Whitney U Test is a
nonparametric test that makes no assumptions about the distribution of data
(e.g., normal distribution). If the probability value (p-value) is less than
the significance level (we use 0.05), then we reject the null hypothesis that
there is no difference between the two groups and conclude that there is a
significant difference between the two groups.

###  VI-B N-gram Classification

In addition to analyzing various linguistic features individually, we further
investigate whether word n-grams show distinguishing capabilities when used as
inputs of the classifiers. Four popular classical machine learning (CML)
classifiers: Logistic Regression (LR), Support Vector Machine (SVM), K-NN
(K-nearest Neighbor), and Gradient Boosting Tree (XGBoost), are employed. The
default parameter settings of the scikit-learn library for various classifiers
are used. For the K-NN classifier, the value of K is set to 5. In addition, we
fine-tune two variants of transformer-based pre-trained language models: BERT
[28] and RoBERTa [29]. A mini-batch size of 16 and a learning rate of 0.00002
are used. During the training, 20% samples are utilized as a validation set.
The Adam optimizer [30] and cross-entropy loss are used as optimizer and loss
functions, respectively. The training process runs for 3 epochs, and an early
stopping criterion is employed. The performances of various classifiers are
evaluated based on 5-fold cross-validation and utilizing precision, recall,
macro F1 score, and accuracy metrics.

##  VII Results and Discussion

TABLE I: Various length related features in the reviews of two groups Attribute | Popular | Non-popular
---|---|---
| median/mean/std. | median/mean/std.
#Sentences / review | 18/24.60/23.08 | 11/16.13/18.21
#Words / review | 332/441.64/398.18 | 194/290.59/316.83
#Words / sentence | 17.82/17.95/7.45 | 16.96/18.02/8.22

Table I shows the distributions of various length-related features in two
groups. We notice some differences in the review length of popular and non-
popular groups. The mean and median word-level review length for the popular
group is much higher (almost 1.6 times) than for the non-popular groups, with
a median value of 332 versus 194. This observation suggests that longer
reviews are likely to be more popular. However, we find that the sentence
length is similar in both groups, which indicates the shorter or longer
sentence has no noticeable impact on the popularity of the review.

TABLE II: Percentage of sentiment words in the reviews of popular and non-popular groups Sentiment | Sentiment | Popular | Non-Popular
---|---|---|---
Lexicon | Orientation | Median/Mean/Std. | Median/Mean/Std.
Opinion | Negative | 3.30/3.45/2.06 | 3.0/3.29/2.46
Lexicon | Positive | 4/4.30/2.40 | 4.2/4.91/3.55
VADER | Negative | 2.8/ 2.97/1.95 | 2.6/2.77/2.21
| Positive | 5.0/5.37/2.60 | 5.2/5.83/3.55
TABLE III: Percentage of emotion words in the reviews of popular and non-popular groups Emotion | Popular | Non-popular
---|---|---
| Median/Mean/Std. | Median/Mean/Std.
Anger | 1.185/1.341/1.257 | 1.03/1.24/1.42
Anticipaton | 2.286/2.473/1.610 | 2.22/2.54/2.12
Disgust | 0.869/1.041/1.13 | 0.75/1.0/1.37
Joy | 1.8210/2.0610/1.6204 | 1.82/2.18/2.21
Sadness | 1.333/1.465/1.151 | 1.15/1.35/1.43
Surprise | 0.984/1.126/1.088 | 0.94/1.17/ 1.38
TABLE IV: The presence of two sets of commonly-used English words (10,000 and 5,000) in the reviews of two groups #Common Words | Group | Median | Mean | Std.
---|---|---|---|---
10000 | Popular | 86.76% | 86.11% | 5.688%
| Non-popular | 87.86% | 87.34% | 5.7%
5000 | Popular | 73.03% | 72.49% | 6.60%
| Non-popular | 74.00% | 73.69% | 6.80%
TABLE V: The percentage of unique and capitalized words in the reviews of two groups | Unique words
---|---
Group | Median | Mean | Std.
Popular | 62.80% | 64.33% | 11.86%
Non-popular | 67.86% | 69.52% | 13.29%
| Capitalized words
|  (%) of reviews | (%) in words (mean) |
Popular | 787/2259 (34.83%) | 1.63% |
Non-popular | 559/2555 (21.87%) | 0.67% |

As we can see from Table II, the presence of the sentiment words is similar to
both review groups based on mean and median values. We notice that this
observation is true for both types of semantic orientations (positive and
negative) and both lexicons (Opinion Lexicon and VADER). Table III shows that,
overall, the percentage of the emotion word is a bit higher in popular groups.
Especially, the emotions related to Anger, Disgust, and Sadness are at least
10% higher in the popular group than the non-popular group.

TABLE VI: Presence of various types of connectives in the reviews of popular and non-popular groups (shown as percentage) Connective | Popular | Non-popular
---|---|---
| Median/Mean/Std. | Median/Mean/Std.
ADD | 3.17/3.25/1.59 | 3.14/3.25/1.82
CE | 3.52/3.63/1.82 | 3.7/3.86/2.29
EMP | 0.0/0.05/0.15 | 0.0/0.06/0.31
ILL | 0.0/0.03/0.1 | 0.0/0.02/0.12
SEQ | 0.54/0.68 /0.73 | 0.51/0.78/1.21

From table IV, we find the percentages of common words (i.e., 5000 and
10000-word lists) are similar in both review groups. The results indicate that
the presence of unfamiliar or less-known words does not show any relation to
the likeability of the reviews. We notice that the presence of unique words is
a bit higher in the popular group, which suggests that people tend to prefer
reviews with fewer repetitive words. Similarly, the percentage of capitalized
words (Table V) is higher in popular reviews. This observation suggests that
strong feelings are somehow related to the likability of reviews.

We observe the percentages of connectives are similar to both popular and non-
popular groups (Table VI). Among all the five types of connectives considered,
none of them show any differences between the two groups. Similar patterns are
observed for subordinate conjunctions (Table VII), which suggest discourse-
related features hardly provide any distinguishing signals to separate popular
and non-popular groups.

TABLE VII: The presence of subordinating conjunctions in the reviews of popular and non-popular groups Review Group | Median | Mean | Std.
---|---|---|---
Popular | 4.22% | 4.15% | 2.09%
Non-Popular | 4.17% | 4.15% | 1.76%

Regarding the presence of the named entity (Table VIII), although not
significant, we observe some differences in their presence in popular and non-
popular reviews. We find that popular reviews usually have more occurrences of
all three types of entities. The readability tests, which indicate how easy it
is to read a particular piece of text, show minimal differences between the
two groups of reviews (Table IX), which implies that the readability measures
are not a good indicator of the likeability of the review.

TABLE VIII: The presence of named entities in the reviews of popular and non-popular groups Entity | Review Group | Median | Mean | Std.
---|---|---|---|---
PER | Popular | 0.8% | 1.02% | 1.0%
| Non-Popular | 0.6% | 0.9% | 1.0%
ORG | Popular | 0.6% | 0.77% | 1.14%
| Non-popular | 0.3% | 0.68% | 1.02%
GE | Popular | 0.1% | 0.32% | 0.55%
| Non-popular | 0.0% | 0.28% | 0.68%
TABLE IX: Statistics of FRE and FKG scores in popular and non-popular groups Readability | Group | Median | Mean | Std.
---|---|---|---|---
FRE | Popular | 59.29 | 58.70 | 13.18
| Non-popular | 61.32 | 60.32 | 13.52
FKG | Popular | 9.47 | 9.78 | 3.58
| Non-popular | 9.04 | 9.29 | 3.50

Table X provides the precision, recall, macro F1 scores, and accuracy of
various traditional ML classifiers and transformers-based BERT and RoBERTa
models. We find that classical ML algorithms provide macro F1 scores in the
range of 0.56-0.63. All classifiers yield poor performance, the best F1 score
obtained at 0.65 by transformer-based language models.

TABLE X: Performance of various classifiers for determining review popularity using word n-gram features Classifier | Precision | Recall | Macro F1 | Accuracy
---|---|---|---|---
LR | 0.63 | 0.63 | 0.63 | 0.63
SVM | 0.63 | 0.63 | 0.63 | 0.63
XGBoost | 0.60 | 0.60 | 0.60 | 0.60
K-NN | 0.56 | 0.55 | 0.56 | 0.54
RoBERTa | 0.65 | 0.64 | 0.65 | 0.66
BERT | 0.65 | 0.64 | 0.64 | 0.66

The analysis regarding the connection between various textual features and
review likeability suggests that most of the features are unrelated to review
likeability. In addition, the poor performance of word n-gram-based features
reveals the limitation of word-based content features.

###  VII-A Statistical Significance

We further investigate which features exhibit significant differences between
the two groups of book reviews. We employ the Mann-Whitney U test with a
probabilistic p-value of 0.05 to discover the significant attributes. Based on
the significance test, we notice the review length yields significant
differences in the two groups having a p-value less than 0.05. Regarding the
sentiment lexicon, we find that negative sentiment shows significant
differences in the two groups, while positive sentiment does not.
Additionally, the number of unique words differs significantly in the two
groups. Regarding the discourse and semantic features, we observe no
significant differences.

###  VII-B Practical Implications

The study investigates how user-generated text in book reviews becomes
popular. Understanding these factors provides valuable insights into user
behavior and preferences, benefiting authors, publishers, and online platforms
in content creation. By incorporating these insights into marketing
strategies, authors and publishers can effectively promote their books.
Additionally, the study contributes to existing research with a comprehensive
analysis of popularity factors, making it valuable for scholars interested in
user-generated content, text analytics, and online communities. Overall, the
study enhances our understanding of user behavior and content creation in the
digital age.

##  VIII Summary and Future Work

In this study, we explore various linguistic aspects of book reviews to
understand how they are related to the attained popularity. We utilize a
dataset of around 9000 book reviews to extract varied types of features and
investigate their influences on the popularity of the reviews. We find, except
for a few attributes like review length and the presence of capitalized and
unique words, most of the other linguistic features exhibit similar patterns
in popular and non-popular reviews. Besides, we observe that ML classifiers
demonstrate inadequate efficacy for automatic review popularity determination,
which reveals the challenges associated with this task. The results suggest
that determining the popularity of reviews in the creative domain is highly
challenging.

Despite predominantly negative results and findings, this study lays the
groundwork for future analyses by highlighting the challenges and limitations
of conventional attributes in distinguishing the likeness of user-generated
content in creative domains. The future work will incorporate additional
lexical, semantic, and grammatical features and an enlarged dataset for
identifying better signals for review popularity projection.

## References

  * [1] M. P. O’Mahony and B. Smyth, “Using readability tests to predict helpful product reviews,” in Paper presented at RIAO 2010 the 9th international conference on Adaptivity, Personalization and Fusion of Heterogeneous Information, Paris, France, April 28-30, 2010, 2010.
  * [2] H. Hong, D. Xu, G. A. Wang, and W. Fan, “Understanding the determinants of online review helpfulness: A meta-analytic investigation,” Decision Support Systems, vol. 102, pp. 1–11, 2017.
  * [3] B. Lutz, N. Pröllochs, and D. Neumann, “Are longer reviews always more helpful? disentangling the interplay between review length and line of argumentation,” Journal of Business Research, vol. 144, pp. 888–901, 2022.
  * [4] A. H. Huang, K. Chen, D. C. Yen, and T. P. Tran, “A study of factors that contribute to online review helpfulness,” Computers in Human Behavior, vol. 48, pp. 17–27, 2015.
  * [5] M. Malik and A. Hussain, “An analysis of review content and reviewer variables that contribute to review helpfulness,” Information Processing & Management, vol. 54, no. 1, pp. 88–104, 2018.
  * [6] S. Sazzed, “A hybrid approach of opinion mining and comparative linguistic analysis of restaurant reviews,” in Proceedings of the International Conference on Recent Advances in Natural Language Processing (RANLP 2021), pp. 1281–1288, 2021.
  * [7] S. Sazzed, “Stylometric and semantic analysis of demographically diverse non-native english review data,” in 2022 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM), pp. 470–476, IEEE, 2022.
  * [8] S. Sazzed and S. Jayarathna, “Ssentia: a self-supervised sentiment analyzer for classification from unlabeled data,” Machine Learning with Applications, vol. 4, p. 100026, 2021.
  * [9] S. Sazzed, “Understanding linguistic variations in neutral and strongly opinionated reviews,” in 2022 21st IEEE International Conference on Machine Learning and Applications (ICMLA), pp. 1512–1516, IEEE, 2022.
  * [10] J. Wu, “Review popularity and review helpfulness: A model for user review effectiveness,” Decision Support Systems, vol. 97, pp. 92–103, 2017.
  * [11] S. Krishnamoorthy, “Linguistic features for review helpfulness prediction,” Expert Systems with Applications, vol. 42, no. 7, pp. 3751–3759, 2015.
  * [12] S. P. Eslami, M. Ghasemaghaei, and K. Hassanein, “Which online reviews do consumers find most helpful? a multi-method investigation,” Decision Support Systems, vol. 113, pp. 32–42, 2018.
  * [13] A. Y. Chua and S. Banerjee, “Understanding review helpfulness as a function of reviewer reputation, review rating, and review depth,” Journal of the Association for Information Science and Technology, vol. 66, no. 2, pp. 354–362, 2015.
  * [14] H. S. Choi and S. Leon, “An empirical investigation of online review helpfulness: A big data perspective,” Decision Support Systems, vol. 139, p. 113403, 2020.
  * [15] A. Ghose and P. G. Ipeirotis, “Estimating the helpfulness and economic impact of product reviews: Mining text and reviewer characteristics,” IEEE transactions on knowledge and data engineering, vol. 23, no. 10, pp. 1498–1512, 2010.
  * [16] G. O. Diaz and V. Ng, “Modeling and prediction of online product review helpfulness: a survey,” in Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pp. 698–708, 2018.
  * [17] M. Malik and A. Hussain, “Helpfulness of product reviews as a function of discrete positive and negative emotions,” Computers in Human Behavior, vol. 73, pp. 290–302, 2017.
  * [18] M. Hu and B. Liu, “Mining and summarizing customer reviews,” in Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 168–177, 2004.
  * [19] C. Hutto and E. Gilbert, “Vader: A parsimonious rule-based model for sentiment analysis of social media text,” in Proceedings of the international AAAI conference on web and social media, vol. 8, pp. 216–225, 2014.
  * [20] S. Mohammad and P. Turney, “Emotions evoked by common words and phrases: Using mechanical turk to create an emotion lexicon,” in Proceedings of the NAACL HLT 2010 workshop on computational approaches to analysis and generation of emotion in text, pp. 26–34, 2010.
  * [21] I. Pak and P. L. Teh, “Value of expressions behind the letter capitalization in product reviews,” in Proceedings of the 2018 7th International Conference on Software and Computer Applications, pp. 147–152, 2018.
  * [22] S. Chan and A. Fyshe, “Social and emotional correlates of capitalization on twitter,” in Proceedings of the Second Workshop on Computational Modeling of People’s Opinions, Personality, and Emotions in Social Media, pp. 10–15, 2018.
  * [23] P. Durán, D. Malvern, B. Richards, and N. Chipere, “Developmental trends in lexical diversity,” Applied Linguistics, vol. 25, no. 2, pp. 220–242, 2004.
  * [24] A. C. Graesser, D. S. McNamara, M. M. Louwerse, and Z. Cai, “Coh-metrix: Analysis of text on cohesion and language,” Behavior research methods, instruments, & computers, vol. 36, no. 2, pp. 193–202, 2004.
  * [25] M. Honnibal and I. Montani, “spaCy 2: Natural language understanding with Bloom embeddings, convolutional neural networks and incremental parsing.” To appear, 2017.
  * [26] B. L. Zakaluk and S. J. Samuels, Readability: Its Past, Present, and Future. ERIC, 1988.
  * [27] Z. Liu and S. Park, “What makes a useful online review? implication for travel product websites,” Tourism management, vol. 47, pp. 140–151, 2015.
  * [28] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” arXiv preprint arXiv:1810.04805, 2018.
  * [29] Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis, L. Zettlemoyer, and V. Stoyanov, “Roberta: A robustly optimized bert pretraining approach,” arXiv preprint arXiv:1907.11692, 2019.
  * [30] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.


# Linear Regression, Covariate Selection and the Failure of Modelling.

Laurie Davies
University of Duisburg-Essen
<EMAIL_ADDRESS>
###### Abstract

It is argued that all model based approaches to the selection of covariates in
linear regression have failed. This applies to frequentist approaches based on
P-values and to Bayesian approaches although for different reasons. In the
first part of the paper 13 model based procedures are compared to the model-
free Gaussian covariate procedure in terms of the covariates selected and the
time required. The comparison is based on seven data sets and three
simulations. There is nothing special about these data sets which are often
used as examples in the literature. All the model based procedures failed.

In the second part of the paper it is argued that the cause of this failure is
the very use of a model. If the model involves all the available covariates
standard P-values can be used. The use of P-values in this situation is quite
straightforward. As soon as the model specifies only some unknown subset of
the covariates the problem being to identify this subset the situation changes
radically. There are many P-values, they are dependent and most of them are
invalid. The P-value based approach collapses. The Bayesian paradigm also
assumes a correct model but although there are no conceptual problems with a
large number of covariates there is a considerable overhead causing
computational and allocation problems even for moderately sized data sets.

The Gaussian covariate procedure is based on P-values which are defined as the
probability that a random Gaussian covariate is better than the covariate
being considered. These P-values are exact and valid whatever the situation.
The allocation requirements and the algorithmic complexity are both linear in
the size of the data making the procedure capable of handling large data sets.
It outperforms all the other procedures in every respect.

Keywords. linear regression, covariate selection, models, real data

##  1 Introduction

The claim that model based procedures have failed will be made in Section 2.
It is based on a comparison of the Gaussian covariate method Davies and
Dümbgen, (2022) with the following 13 selection procedures: lasso (Tibshirani,
(1996)), knockoff (Candès et al., (2018)), scaled sparse linear regression
(Sun and Zhang, (2012)), SIS (Sure Independence Screening) (Fan and Lv,
(2008)), desparsified lasso (Zhang and Zhang, (2014)), stability selection
(Meinshausen and Bühlmann, (2010)), ridge regression (Bühlmann, (2013)),
multiple splitting (Wasserman and Roeder, (2009)), EMVS (Expectation-
Maximization Approach to Bayesian Variable Selection) (Rockova and George,
(2014)) and Spike and Slab Regression (Scott, (2021)), Threshold Adaptive
Validation (Laszkiewicz et al., (2021)), graphical lasso (Friedman et al.,
(2008, 2019)) and huge (High-Dimensional Undirected Graph Estimation) (Jiang
et al., (2021). The latter three construct dependency graphs. A short
description of the Gaussian covariate method which is sufficient for this
paper is given in the Appendix.

The comparison is based on three simulations and the following seven data
sets: riboflavin Bühlmann et al., (2014), leukemia Golub et al., (1999),
lymphoma Alizadeh et al., (2000) and Dettling and Bühlmann, (2002),
osteoarthritis Cox and Battey, (2017), the Boston Housing data set Harrison
and Rubinfeld, (1978), sunspot data SILSO, (2020) and the American Business
Cycle data Gordon, (1986). Their dimensions $(n,q)$ are respectively
(71,4088), (72,3571), (62,4026),(129 48802), (506,176358), (3253, 3252) and
(240,22). The large number of covariates, 176358, for the Boston Housing data
is that these are all interactions of degree at most eight of the original 13
covariates.

The reasons for this failure are discussed in Section 3. Statistical models
are often successful so a case where they fail is of particular interest, the
issues involved are not limited to this particular failure.

##  2 Comparing procedures

Two Gaussian procedures f1st and f3st with $p0=0.01$ and $kmn=10$ (see Davies
and Dümbgen, (2022) and the Appendix), will be compared with the model based
procedures. All these are used with their default values, for example $B=100$
in multi.split, with the exceptions of stability which requires a value for
$EV$ set here to $EV=1$ and lm.spike which requires a value for $niter$ set
here to $niter=1000$.

###  2.1 Riboflavin and leukemia simulations

The simplest linear model for formulating the problem of covariate selection
is

$\boldsymbol{Y}=\sum_{\boldsymbol{x}_{i}\in{\mathcal{S}}_{0}}\beta_{i}\boldsymbol{x}_{i}+\sigma\boldsymbol{\varepsilon}$
(1)

where $\boldsymbol{\varepsilon}$ is standard Gaussian nose and
${\mathcal{S}}_{0}$ is a subset of the $q$ covariates, the so called active
set, with $\beta_{i}\neq 0$ for $\boldsymbol{x}_{i}\in{\mathcal{S}}_{0}$. The
problem is to determine ${\mathcal{S}}_{0}$ and to construct confidence
intervals and P-values for the corresponding values of the $\beta_{i}$ ( see
Dezeure et al., (2015)). Having postulated the model (1) the statistician will
now work within it, that is, treating the model as true.

The following simulations based on the model (1) were performed but using the
4088 covariates of the riboflavin data and the 3571 covariates of the leukemia
data . The covariates were first standardized to have mean zero and variance
one. Four covariates
$\\{\boldsymbol{x}_{i_{1}},\boldsymbol{x}_{i_{2}},\boldsymbol{x}_{i_{3}},\boldsymbol{x}_{i_{4}}\\}$
are chosen at random and the dependent variable given by

$\boldsymbol{y}=20\sum_{j=1}^{4}\boldsymbol{x}_{i_{j}}+\boldsymbol{\varepsilon}$

where $\boldsymbol{\varepsilon}$ is standard Gaussian noise. The reason for
taking just four covariates is that the simplest Gaussian covariate procedure
applied to the riboflavin data results in just the four covariates
$\boldsymbol{x}_{4003},\boldsymbol{x}_{2564},\boldsymbol{x}_{73},\boldsymbol{x}_{2034}$.
In the case of the leukemia data six covariates were chosen at random.

All ten of the model based procedures were considered using the functions
cv.glmnet for lasso (Friedman et al., (2021)), knockoff.filter for knockoff
(Barber et al., (2020)), scalreg for scaled sparse linear regression (Sun,
(2019)) , SIS for Sure Independence Screening (Feng et al., (2020)), stability
for stability selection, multi.split for multiple splitting , ridge.proj for
ridge regression, lasso.proj for desparsifed regression (the last four
available from Meier et al., (2021)), EMVS for EMVS variable selection
(Rockova, (2021)) and lm.spike for BoomSpikeSlab (Scott, (2021)). Two Gaussian
covariate functions were used, f1st and f3st with $kmn=10$ and with two
different values of the parameter $m=1,2$ of f3st.

Two, knockoff and desparse.lasso, required more than one hour for a single
simulation and were killed after this time. multi.split often halted with an
error, this occurred after 70 simulations in Table 1 and 38 simulations in
Table 2. BoomSpikeSlab in conjunction with the operating system Leap 15.1
occasionally shut the laptop down, this happened after 72 simulations with the
riboflavin data. ridge is very slow and was killed after 18 simulations. For
the remaining procedures 100 simulations were performed.

Riboflavin: 100 simulations; (*) 70, (**) 72, (***) 18 simulations
---
method | $fp$ | $fn$ | % correct | time
f1st | 0.77 | 0.72 | 75 | 1 (0.026)
f3st,m=1 | 0.18 | 0.17 | 93 | 5
f3st,m=2 | 0.07 | 0.05 | 98 | 24
lasso | 25.0 | 0.07 | 0 | 19
scalreg | 16.3 | 1.08 | 0 | 85
SIS | 13.5 | 2.45 | 3 | 150
stability | 0.24 | 2.16 | 8 | 96
multi.split(*) | 0.23 | 1.59 | 36 | 1570
BoomSpikeSlab(**) | 0.33 | 0.42 | 87 | 1540
EMVS | 0.00 | 4.00 | 0 | 27
knockoff | ? | ? | ? |  $>$150000
desparse.lasso | ? | ? | ? |  $>$150000
ridge(***) | 0.00 | 4.00 | 0.00 | 10000
Table 1: Columns 2-4 give the average number of false positives, the average
number of false negatives and the % of correct selections. Column 5 gives time
compared with Gaussian covariates which required on average 0.026 seconds per
simulation.  Leukemia: 100 simulations, (*) 38, (**) 18 simulations
---
method | $fp$ | $fn$ | % correct | time
f1st | 0.43 | 1.21 | 71 | 1 ( 0.025)
f3st,m=1 | 0.26 | 0.68 | 79 | 7
f3st,m=2 | 0.29 | 0.69 | 81 | 31
lasso | 19.8 | 0.11 | 0 | 18
scalreg | 10.8 | 1.55 | 0 | 144
SIS | 11.6 | 2.35 | 6 | 130
stability | 0.04 | 4.73 | 1 | 94
multi.split (*) | 0.10 | 4.7 | 3 | 1500
BoomSpikeSlab | 0.50 | 3.15 | 38 | 1340
EMVS | 0.00 | 6.00 | 0 | 19
knockoff | ? | ? | ? |  $>$150000
desparse.lasso | ? | ? | ? |  $>$150000
ridge(**) | 0.00 | 6.00 | 0.00 | 7500
Table 2: Columns as for Table 1.

It is clear from Tables 1 and 2 that the Gaussian procedures are
overwhelmingly superior. This is in spite of the fact that the data were
generated under the model (1). In terms of covariate selection BoomSpikeSlab
was the best of the model based procedures but worse than f3st with $kmn=10$
and $m=1$, it was also between 190 and 300 times slower. When applied to the
osteoarthritis and the Boston Housing datasets (see Tables 8 and 9 below)
BoomSpikeSlab fails completely.

###  2.2 The riboflavin, the leukemia, the lymphoma, the osteoarthritis, the
Boston Housing, the sunspot and the American Business Cycle data sets

Table 3 gives the detailed results for the first five data sets using the
default version f1st of the stepwise procedure with $kmn=10$. The first rows
give in order the data set followed by the selected covariates. The second
rows gives the time in seconds and the sum of squared residuals followed by
the Gaussian P-values of the selected covariates.

Gaussian covariate results for the first five data sets
---
riboflavin | 4003 | 2564 | 73 | 2034 | * | * | *
0.02, 8.448 | 1.97e-14 | 2.81e-13 | 4.1e-09 | 2.84e-05 | * | * | *
leukemia | 1182 | 1219 | 2888 | 183 | 3038 | 2558 | *
0.181, 0.931 | 2.7e-07 | 3.23e-07 | 4.90e-05 | 1.47e-04 | 4.47e-03 | 8.00e-03 | *
lymphoma | 2805 | 3727 | 714 | 2036 | 3201 | 1638 | 2057
0.024, 0.962 | 1.86e-06 | 2.30e-15 | 3.75e-07 | 4.29e-09 | 8.19e04 | 5.55e-03 | 9.96e-03
Osteoarthritis | 1149 | 31848 | 33321 | * | * | * | *
0.13, 6.767 | 3.80e-18 | 1.54e-05 | 6.25e-03 | * | * | * | *
Boston inter. | 7222 | 164793 | 106761 | 81454 | 29154 | 170455 | 67241
4.59, 6129.2 | 1.95e-143 | 2.15e-10 | 4.65e-07 | 8.13e-13 | 1.62e-16 | 1.18e-52 | 1.99e-04
Table 3: The Gaussian covariate results for the data sets riboflavin,
osteoarthritis and Boston Housing

Table 3 gives only one approximation for each data sets. More can be obtained
by using either f2st or f3st, the latter usually giving better results.
Applying f3st to the riboflavin data with $kmn=15,m=5$ results in 129
approximations. The first 10 are given in Table 4 in order of the sum of
squared residuals. The approximation in Table 3 is 32nd on the list.

$\boldsymbol{ss}$ | Riboflavin: Included covariates
---|---
3.72 | 4004 | 2564 | 73 | 315 | 2936 | 997 | 991 | 1661 | 3255
4.23 | 4004 | 2564 | 73 | 315 | 2936 | 997 | 1661 | 2048 | *
4.87 | 4004 | 2564 | 144 | 1131 | 3138 | 2186 | 1771 | * | *
5.43 | 1279 | 4004 | 2564 | 73 | 1131 | 2140 | * | * | *
5.47 | 4003 | 2564 | 69 | 1425 | 413 | 2484 | 1194 | * | *
5.53 | 1279 | 4004 | 2564 | 144 | 1131 | 2140 | * | * | *
5.57 | 1279 | 2564 | 73 | 1131 | 2140 | 4005 | * | * | *
5.71 | 1278 | 4004 | 2564 | 73 | 1131 | 2140 | * | * | *
5.71 | 4004 | 1304 | 2564 | 73 | 1131 | 2140 | * | * | *
5.80 | 1278 | 4004 | 2564 | 144 | 1131 | 2140 | * | * | *
Table 4: The first ten of the 129 approximations given by f3st with $kmn=15$
and $m=5$ in order of the sum of squared residuals $\boldsymbol{ss}$.

In the Tables 5, 6, 8, 9 and 11 the first column gives the method and the
second the numbers $k$ of covariates selected and the number $fp$ of selected
covariates which had Gaussian P-values exceeding 0.99. The third column
indicates whether the procedure returns P-values. The fourth gives the sum of
squared residuals. The fifth gives the factor $t/t0$ where $t$ is the time in
seconds of the procedure and $t0$ is the time for the Gaussian covariate
method f1st with $kmn=10$. The results show that all ten modelled based
procedures failed for all seven data sets. Note that f1st, f3st, scalreg and
EMVS are deterministic and give the same result when repeated. The others make
use of random choices and if repeated may give somewhat different results.

When applied to the riboflavin data Table 5 knockoff was killed after two
hours and had at this point selected zero covariates. For this reason it was
not applied to the remaining data sets. Similarly desparsified was very slow
and was applied only to the riboflavin, leukemia and the lagged sunspot and
American Business Cycle data sets. The osteoarthritis data set was too large
for EMVS and BoonSpikeSlab. This left only five model based procedures for the
Boston Housing data Table 9. Three failed with errors presumably because of
the size of the data. Lasso returned 63 covariates 61 of which had Gaussian
P-values of 1. In further applications of lasso it returned between four and
over 100 covariates. Stability returned just three covariates and took over 11
minutes to do so.

riboflavin (71,4088)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 4,0 | yes | 8.45 | 1 (0.024)
f3st,m=1 | 6,0 | yes | 6.21 | 4
f3st,m=2 | 6,0 | yes | 5.43 | 22
lasso | 32,30 | no | 2.05 | 25.7
knockoff | 0,0 | no | * |  $>$7e+05 (killed)
scalreg | 9,6 | no | 10.62 | 28.7
SIS | 4,0 | no | 11.49 | 89
desparsified lasso | 0,0 | yes | * | 130012
stability | 0,0 | no | * | 103
ridge.proj | 0,0 | yes | * | 12248
multi.split | 4,2 | yes | 17.45 | 1421
EMVS | 0,0 | no | * | 22
BoomSpikeSlab | (5,2) | no | 16.92 | 2290
Table 5: The results for the riboflavin data.  leukemia $(n,q)=(72,3571)$
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 6,0 | yes | 0.93 | 1 (0.0212)
f3st, m=1 | 6,0 | yes | 0.49 | 7
f3st, m=2 | 6,0 | yes | 0.49 | 30
lasso | 17,14 | no | 0.661 | 20
scalreg | 14,10 | no | 1.08 | 14
SIS | 3,0 | no | 1.89 | 85
desparsified lasso | 4,1 | yes | 1.56 | 104670
stability | 1,0 | no | 5.43 | 103
ridge.proj | 0,0 | yes | * | 11515
multi.split | 1,0 | yes | 4.42 | 1420
EMVS | 0,0 | no | * | 20
BoomSpikeSlab | (3,3) | no | 13.73 | 1600
Table 6: The results for the leukemia data.  Lymphoma data (62, 4026)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 7,0 | yes | 0.96 | 1 (0.024)
f3st,m=1 | 7,0 | yes | 0.96 | 8
lasso | 43,42 | no | 0.058 | 21
scalreg | 42,41 | no | 0.083 | 32
SIS | 5,0 | no | 1.75 | 430
desparsified lasso | 5,0 | yes | 1.62 | 111180
stability | 0,0 | no | * | 90
ridge.proj | 0,0 | yes | * | 11700
multi.split | 0,0 | yes | * | 1540
EMVS | 0,0 | no | * | 20
BoomSpikeSlab | (6,4) | no | 6.82 | 1650
Table 7: The results for the lymphoma data.  Osteoarthritis (129,48802)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 3,0 | yes | 6.77 | 1 (0.29)
f3st,m=1 | 3,0 | yes | 6.23 | 4
f3st,m=2 | 4,0 | yes | 5.02 | 16
lasso | 66,64 | no | 0.35 | 24.5
scalreg | 5,3 | no | 11.14 | 45
SIS | 0,0 | no | * | 80
stability | 1,0 | no | 13.3 | 116
multi.split | 3,1 | yes | 8.35 | 1470
EMVS | 0,0 | no | * | bad alloc.
BoomSpikeSlab | 0,0 | no | * | 17.7 Gb
Table 8: The results for the osteoarthritis data.  Boston data with
interactions (506,176358)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 7,0 | yes | 6130 | 1 (4.42)
f3st,m=1 | 9,0 | yes | 5589 | 8
f3st,m=3 | 10,0 | yes | 5557 | 68
lasso | 4,2 | no | 9950 | 30
SIS | 0,0 | no | * | max.iter
scalreg | 0,0 | no | * | self killed
stability | 4,1 | no | 14375 | 160
multi.split | 0,0 | yes | * | incorrect no. of P-values
Table 9: The results for the Boston Housing data.

The results for the sunspot data are given in Tables 10 and 11. The covariates
in Table 10 are $\boldsymbol{x}_{2j}=\sin(\pi jt)$ and
$\boldsymbol{x}_{2j-1}=\cos(\pi jt)$ for $j=1,\ldots,1626$ with
$t=(1:3253)/3253$.

Table 11 the covariates are the first 500 lags of the number of sunspots. The
lags selected by f1st were 1, 2 , 4, 6, 9, 27 and 117. We note that the
maximum likelihood version of ar of R selects the first 12 lags and has a
residual sum of squares of 1583885.

The American Business Cycle Gordon, (1986) data consists of quarterly data
from 1919-1941 of the variables GNP72, CPRATE, CORPYIELD, M1, M2, BASE,
CSTOCK, WRICE67, PRODUR72, NONRES72, IRES72, DBUSI72, CDUR72, CNDUR72, XPT72,
MPT72, GOVPUR72, NCSPDE72, NCSBS72, NCSCON72, CCSPDE72 and CCSBS72 available
from

    http://data.nber.org/data/abc/

In Table 12 the dependent variable is GNP72 (Gross National Product in 1972
dollar terms), the covariates are the lags 1:16 of all the 22 variables. These
are ordered as follows. The first 16 covariates are the 16 lags of the first
variable GNP. The second sixteen are the 16 lags of the second variable CPRATE
and so on. Thus covariate 180 is the fourth lag of the 11th variable IRES72,
$180=11*16+4$.

In contrast to all other data sets desparsified lasso and ridge.proj select a
large number of covariates, EMVS selects 223 which together with the intercept
give a perfect fit. Of the model based procedure SIS is the only one which
gives a reasonable selection. All the five selected covariates have Gaussian
P-values less then 0.05, but two have values exceeding 0.01. If we run f1st
with $p0=0.05$ it returns six covariates with a sum of squared residuals of
16870.

The analysis can be repeated for the remaining 21 variables giving a vector
autoregressive scheme. The results for f1st are given in Table 13 with
$p0=0.01$ replaced by $p0=0.01/22$.

Sunspot data (3253,3252)
---
method | $k,fp$ | ssq.resid | time
f1st | 54,0 | 2023014 | 1(2.41)
lasso | 787,649 | 632108 | 11
scalreg | 62 | 2095413 | 1285
SIS | * | * | *
stability | 34,3 | 2854149 | 15
multi.split | 57 | 2027120 | 400
BoomSpikeSlab | 52 | 4630928 | 24
EMVS | 0 | * | 20
Table 10: The results for the sunspot data.

Sunspot, lags (2753,500)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 7,0 | yes | 1507616 | 1 (0.436)
f3st,m=1 | 7,0 | yes | 1507616 | 7
lasso | 11,4 | no | 1522261 | 35
scalreg | 451,392 | no | 1255481 | 32
SIS | 314,286 | no | 1263136 | 825
desparsified lasso | 7,0 | yes | 1557744 | 14853
stability | 8,1 | no | 1545191 | 30
ridge.proj | 5,0 | yes | 1579760 | 82
multi.split | 7,0 | yes | 1546039 | 3403
EMVS | 500,0 | no | 1254916 | 6
BoomSpikeSlab | (7,0) | no | 1956351 | 20
Table 11: The results for the sunspot data with lags 1:500.  American Business
Cycle (224,352)
---
method | $k,fp$ | P-values | $\boldsymbol{ss}$ | time
f1st | 3,0 | yes | 18765 | 1 (0.039)
f3st,m=1 | 6,0 | yes | 18405 | 3
lasso | 4,2 | no | 24980 | 3
scalreg | 83,69 | no | 4960 | 19
SIS | 5,0 | no | 17854 | 16
desparsified lasso | 186,185 | yes | 1880 | 1000
stability | 2,0 | no | 25460 | 12
ridge.proj | 102,95 | yes | 8130 | 30
multi.split | 2,0 | yes | 25460 | 200
EMVS | 223, NaN | no | 0 | 2.8
BoomSpikeSlab | (4,0) | no | 48750 | 65
Table 12: The results for the American Business Cycle data with lags 1:16.

American Business Cycle
---
Vector autoregression (224,352)
variable | selected covariates
1 | 1 | 18 | 180 | *
2 | 17 | 161 | * | *
3 | 33 | 17 | * | *
4 | 49 | * | * | *
5 | 65 | * | * | *
6 | 81 | * | * | *
7 | 97 | 98 | * | *
8 | 113 | 117 | 69 | *
9 | 129 | 193 | * | *
10 | 145 | * | * | *
11 | 161 | 162 | * | *
12 | 177 | * | * | *
13 | 193 | 17 | 21 | *
14 | 209 | * | * | *
15 | 225 | 193 | 12 | *
16 | 241 | * | * | *
17 | 257 | 258 | * | *
18 | 273 | 129 | * | *
19 | 289 | 145 | 161 | *
20 | 305 | 193 | 17 | 26
21 | 321 | 323 | 322 | *
22 | 337 | 326 | * | *
Table 13: American Business Cycle, vector autoregressive scheme with lags
1:16.

###  2.3 Constructing models

The residuals of the riboflavin data look more or less normal so the data can
be modelled with the selected covariates plus Gaussian noise.

Data such as the leukemia data are usually analysed using the logit model.
This can be done by taking the covariates selected for the linear
approximations using the Gaussian step-wise procedure and running a logistic
regression using these covariates. If this is done for the 142 linear
approximations returned by f2st 32 of them give a perfect fit.

The sunspot data can be modelled as follow. Let $s(t)$ be the number of
sunspots at time $t$ and denote by $f_{1}$ the function returned by the step-
wise procedure. The residuals are defined as $r_{1}=s-\log(1+\exp(f_{1}))$
where the transformation is made to force the function to be non-negative.
They are reasonably described by an autoregressive process of order 1 with
coefficient 0.4. We write

$d(t)=s(t)-0.4s(t-1),t=2,\ldots,3253$

which removes the autocorrelation. We now apply the step-wise procedure to $d$
with the covariates $tr_{j}(t),t=2,\ldots,3253,j=1,\ldots,3252$ and denote the
resulting function by $f_{2}$. This eliminates the high frequency
trigonometric functions caused by the autocorrelations. We define a new
approximating function by

$f_{3}(t)=\sum_{j=0}^{6}0.4^{j}f_{2}(t-j)$

with residuals $r_{3}=s-\log(1+\exp(f_{3}))$. The size of the residuals
depends on $\log(1+\exp(f_{3}))$ and a simple linear regression shows that

$z=r_{3}/(1.28(9.4+0.125\log(1+\exp(f_{3})))$

can be well approximated by an autoregressive Gaussian process of order one
with parameter 0.5. This leads to the model

$S(t)=f_{3}(t)+1.28(9.4+0.125\log(1+\exp(f_{3})))Z(t)$

where $Z(t)\sim N(0,1)$ is a Gaussian autoregressive process of order with
parameter 0.5.

###  2.4 Graphs

There are several methods for calculating dependency graphs. The ones we
compare here are the Gaussian covariate graph function fgr1st, thav.glasso
(Laszkiewicz, (2021)), huge (Jiang et al., (2021)) and glasso (Friedman et
al., (2019)).

We consider firstly a random graph generated as in Meinshausen and Bühlmann,
(2006) but with the correction given in Davies and Dümbgen, (2022). The
undirected graph of the Gaussian procedure was used. The procedures
thav.glasso and glasso have regularization parameters C and rho respectively
with no indication as to how to choose them. They were chosen to equalize
approximately the numbers of false positives and false negatives. The values
used were C=0.4 and rho=0.236. Table 14 gives the results. The second column
gives the number of edges, the third the number of false positives $fp$, the
fourth the number of false negatives $fn$ and the fifth the time in seconds.

Random graph (1000,1000)
---
method | no. edges | $fp$ | $fn$ | time
fgr1st | 1821 | 1 | 3 | 25
thav.glasso | 1776 | 218 | 265 | 90
huge | 1830 | 26 | 19 | 30
glasso | 1840 | 293 | 276 | 16
Table 14: The results for one simulation of the random graph.

The results for the covariates of the riboflavin data were as follows. The
procedures thav.glasso and glasso were killed after one hour with no results.
There may be values of the regularization parameters which work but the user
is given no hint. The huge procedure took 35 seconds but returned zero edges.
The Gaussian covariate procedure with the default values took 16 seconds and
yielded a directed graph with 4491 edges and an undirected graph with 3882
edges.

###  2.5 Summary

The above demonstrates the overwhelming superiority of the Gaussian covariate
procedure. This is exemplified by the two simulations Tables 5 and 6 where
there is a true model which all but BoomSpikeSlab fail to find, Table 5 where
the best in terms of sum of squared residuals scalreg would be number 88 in
the complete list of Table 4, and the Boston Housing data Table 9 where only
five of the ten model based procedures got so far, three of these resulted in
errors and the remaining two remaining failed both with respect to the
selection of covariates and the times required. This may be compared with
Table 15 in the Appendix.

##  3 Reasons for the failure of model based procedures

###  3.1 Truth, P-values and confidence intervals

The decision as to whether a given covariate $\boldsymbol{x}_{j}$ lies in
${\mathcal{S}}_{0}$ is made by formulating and calculating a P-value for the
hypothesis $H_{j}:\beta_{j}=0$. Here is an example taken from Lockhart,
(2017). On page 5 the covariate $\boldsymbol{x}_{1278}$ is considered and
regressing the dependent riboflavin variable on $\boldsymbol{x}_{1278}$
results in the $P_{F}$-value 9.06e-10. This is the smallest value amongst all
covariates. It is not however included amongst the selected covariates of
Table 3. The reason is that when regressed together with the four other
covariates in Table 3 its $P_{G}$-values is 0.2474591 which exceeds the
threshold of $p0=0.01$. The covariate $\boldsymbol{x}_{73}$ exhibits the
opposite behaviour. Its $P_{F}$-value when regressed alone on the dependent
variable is 2.051118e-02 but its $P_{F}$-value when regressed with the
covariates $\boldsymbol{x}_{400},\boldsymbol{x}_{2564}$ and
$\boldsymbol{x}_{2034}$ is 9.968584e-13.

Which $P_{F}$-values are correct? Under the model the only correct
$P_{F}$-values are those for the subset ${\mathcal{S}}_{0}$, all others are
wrong. All P-value correction methods such as Bonferroni assume that the
P-values are correct. It is therefore not surprising that the corrected
P-values are of no help. But even if they are correct the corrected P-values
may still be of no help. It we take the first line of Table 4 the smallest
$P_{F}$-value is 1.83e-23 for the covariate $\boldsymbol{x}_{4004}$. There are
3.766288e+21 subsets of size seven leading to a corrected $P_{F}$-value of
0.068 (see page 19 of Lockhart, (2017)) which is not very convincing. The
corresponding $P_{G}$-value is 7.42e-20 which is correct as it stands.

All the above holds mutatis mutandis for confidence intervals. A confidence
interval for a parameter $\beta_{i}$ is an interval which contains the true
value, that is, the value of $\beta_{i}$ under which the dependent variable
$\boldsymbol{Y}$ was generated, with a specified probability. For this it is
necessary to know the true subset ${\mathcal{S}}_{0}$.

All P-values calculated for data generated under the model (1) have the same
source of randomness, namely ${\boldsymbol{\varepsilon}}$. Thus not only are
the resulting P-values in general false, they are also dependent further
undermining their usefulness data analysis.

The conclusion is that concepts such as F P-values and confidence intervals
derived from modelling with all covariates cannot be transferred to the
problem of determining the correct model even when such a model exists.

###  3.2 Data splitting

The above mentioned problems with P-values and confidence intervals are well
known. One method of trying to overcome them is data splitting. We follow here
Dezeure et al., (2015) but see also Wasserman and Roeder, (2009). Split the
sample into two $I_{1}$ and $I_{2}$ at random with $I_{1}$ of size $n/2$.
Based on $I_{1}$ select a subset $\widetilde{{\mathcal{S}}}_{1}$ of the
covariates and calculate the P-values of the covariates in
$\widetilde{{\mathcal{S}}}_{1}$ using $I_{2}$. For covariates not in
$\widetilde{{\mathcal{S}}}_{1}$ the P-value is set to 1. This is repeated
about 100 times and the P-value for each covariate is calculated from the 100
P-values, for example twice the median. For the procedure to work
${\mathcal{S}}_{0}\subset\widetilde{{\mathcal{S}}}_{1}$ must hold with high
probability. There are various suggestion as to how to choose the subset
$\widetilde{{\mathcal{S}}}_{1}$ (see Wasserman and Roeder, (2009)) but in the
end it does not work. Four applications of multi.split to the riboflavin data
resulted in the covariates $(\boldsymbol{x}_{4003})$,
$(\boldsymbol{x}_{1516},\boldsymbol{x}_{4003})$,
$(\boldsymbol{x}_{1278},\boldsymbol{x}_{4003})$ and $(\boldsymbol{x}_{4003})$
whose corrected P-values varied from 0.019 to 0.984. As a comparison the
$P_{G}$-value of $\boldsymbol{x}_{4003}$ alone is 7.88e-05.

###  3.3 Lasso

Lasso is one of most used selection procedures and is an integral part of
others, for example stability. Nevertheless lasso does not work and one reason
for this may well be that it shrinks the coefficients, one interpretation of
which is that it penalizes the best. We consider the riboflavin data. The
first nine covariates on the lasso path are in order (see Lockhart, (2017)
page 17)

$\boldsymbol{x}_{1278},\boldsymbol{x}_{4003},\boldsymbol{x}_{1516},\boldsymbol{x}_{2564},\boldsymbol{x}_{1588},\boldsymbol{x}_{624},\boldsymbol{x}_{1312},\boldsymbol{x}_{1502},\boldsymbol{x}_{1639}.$

The first nine selected by the Gaussian covariate method are, in order

$\boldsymbol{x}_{1278},\boldsymbol{x}_{4003},\boldsymbol{x}_{2564},\boldsymbol{x}_{73},\boldsymbol{x}_{2034},\boldsymbol{x}_{1131},\boldsymbol{x}_{1762},\boldsymbol{x}_{2186},\boldsymbol{x}_{3459}.$

If the dependent variable is regressed on both sets of covariates the sum of
squared residuals is 9.97 for the lasso covariates and 3.98 for the Gaussian
covariates. The Gaussian covariate function f1st selects the four covariates
$\boldsymbol{x}_{73},\boldsymbol{x}_{2034},\boldsymbol{x}_{2564}$ and
$\boldsymbol{x}_{4003}$, Table 3. The sum of squared residuals is 8.448 which
is less than that for all nine lasso covariates. The same occurs for the
Boston data. With $\lambda=1.2$ lasso returns 13 covariates with sum of
squared residuals 11819. The first thirteen covariates using the Gaussian
procedure have a sum of squared residuals of 5450. The nine covariates
selected by the procedure have a sum of squared residuals of 5589, Table 15.
Given this it is perhaps not surprising why lasso fails. The covariates on the
lasso path are so to speak, greatly inferior to those on the greedy forward
search path.

There have been attempts to provide valid P-values along the lasso path, see
Lockhart et al., (2014) and Tibshirani et al., (2016). There are however
problems. On page 19 of Lockhart, (2017) a significance test for the covariate
1278 of the riboflavin data is considered. The test statistic of Lockhart et
al., (2014) for the first covariate is

$T_{1}=\frac{\lambda_{1}(\lambda_{1}-\lambda_{2})}{\widehat{\sigma}^{2}}=24\text{
or }2.55$

with P-values respectively 0.078 or $3.7^{-11}$ depending on the choice of
$\widehat{\sigma}^{2}$. Lockhart comments ‘Estimation of $\sigma$ is crucial
and hard, I think’. On page 43 we read that one of the $\widehat{\sigma}$
values is the standard deviation of the dependent variable. The second is
‘based on a fit to some 30 predictors’. On page 19 we read ‘Lockhart et al.,
(2014) mostly consider a fictitious universe in which $\sigma$ is known’. The
same applies to Tibshirani et al., (2016). The problem of estimating $\sigma$
remains.

###  3.4 Bayesian procedures

The Bayesian paradigm carries over without change to the situation $q>n$. The
simplest model is

$\boldsymbol{Y}=\boldsymbol{X}{\boldsymbol{\beta}}+\sigma{\boldsymbol{\varepsilon}}$
(2)

which must be supplemented with priors for ${\boldsymbol{\beta}}$ and
$\sigma$. This has a certain elegance when compared with the collapse of the
P-value paradigm but it is not without problems.

In the riboflavin simulations of Section 2.1 there was one unique correct
model. BoomSpikeSlab was slow but in Table 1 it was the third best procedure
in terms of identifying the correct model . In the case of the riboflavin data
itself it was the second best of the model based procedures after SIS but much
worse than f1st and a factor 2290 slower. It failed completely for the
osteoarthritis (Table 8) and for the Boston Housing (Table 9) data sets.

All these failures are due to the problems with the prior over a high
dimensional space $\mathbb{R}^{q}$. Firstly there is the algorithmic problem
of calculating the posterior. BoomSpikeSlab uses MCMC making it exceptionally
slow. Secondly it is based on the covariance matrix of the covariates which in
the case of the osteoarthritis data has 1190842003 entries and a size of about
17Gb. This explains the errors for both Bayesian procedures in Table 8. They
would be completely overwhelmed by the Boston Housing data whose covariance
matrix is of size about 230Gb.

###  3.5 The mathematics of covariate selection

The mathematics of model based procedures is not simple, see for example
Tibshirani et al., (2016), the Appendix of Dezeure et al., (2015) and
Bühlmann, (2013). All the P-value based procedures make use of regularization
and in particular many make use of lasso. Such procedures requires a specific
value of the regularization parameter $\lambda$ and this value depends on the
data. All results on choosing $\lambda$ are of asymptotic nature and many
require a value for $\sigma$ in (1). None of this is any help in choosing
$\lambda$ for a specific data set with the consequence that the default method
for lasso is cross-validation as it was 26 years ago in Tibshirani, (1996),
and it fails.

Contrast this with the Gaussian method which only needs the threshold $p0$ for
the Gaussian P-values. It has an immediate interpretation, the default value
$p0=0.01$ was use for all data sets and simulations. By far the most
mathematically demanding sections of Davies and Dümbgen, (2022) are those on
the asymptotic behaviour for a given $p0$ ($\alpha$ in Davies and Dümbgen,
(2022)). They have nothing to do with the choice of $p0$, they only confirm
its interpretation.

###  3.6 Memory requirements and algorithmic complexity

Both EMVS and BoomSpikeSlab fail for the osteoarthritis and Boston data sets
because of allocation problems alone as shown in the previous Section. For the
same reason thav.glasso, huge and glasso also fail for these two data sets as
they also require the covariance matrix of the covariates. The memory
requirements of f1st are about $nq$ and those of fgr1st about $4nq$.

In Huber, (2011) we read ‘Pay attention to computational complexity; keep it
below $O(n^{3/2})$, or forget about the algorithm’. To which we can add, if
your algorithm requires the covariance matrix of the covariates, forget about
it.

The stepwise version of the Gaussian covariate procedure is very fast. If $k$
covariates are chosen it requires about $6nq(k+1)$ multiplications which is a
factor of $\sqrt{nq}$ smaller than the complexity given by Huber. This is
reflected in the times for the various procedures given above. The fastest of
these is EMVS which is a factor of about 20 slower than Gaussian covariates
but only seems to select zero covariates. The algorithmic complexity of the
other procedures seems not to be given.

###  3.7 Multiple approximations and models

The model based approach postulates a single model such as (1). The goal is to
identify the active subset ${\mathcal{S}}_{0}$. Consequently model based
procedures are designed to output just one set of selected covariates. In
practice repetitions will result in different selected covariates but this is
not intended, it is simply the use of some form of random data splitting. The
word true is part of the vocabulary, sometimes in inverted commas ‘true’ to
emphasize that the author does not believe it. Articles such as Benjamini and
Hochberg, (1995) are completely reliant on it. There is no concept of
approximation.

One advantage of the model approach may be thought to be that the problem of
covariate choice can be posed as an exact mathematical problem, essentially
that of finding the true model (Candès et al., (2018)). When thinking in terms
of approximations it makes little sense to talk of the true approximation. In
general there will be several adequate approximations and they may well be
equally valid. This can be posed as a precise mathematical problem using the
concept of Gaussian P-values as follows: find all subsets with the property
that no $P_{G}$-value exceeds $p0=0.01$. If $q$ is not too large the all
subset version of the Gaussian covariates procedure solves this problem.

For large $q$ the all subset version fails and there is no general strategy
for producing adequate approximations. The greedy forward selection of the
Gaussian stepwise procedure is one of the oldest selection procedures. The
success of the Gaussian version is due to the stopping rule and the
availability of exact and universally valid $P_{G}$-values. Its disadvantage
is that it yields only one approximation. The Gaussian covariate R package
offers two functions f2st and f3st to overcome this. They are described in the
Appendix.

This is not an argument against models, it is an argument against a vocabulary
and practice of treating models as true. A model can be explicitly treated as
an approximation but this requires a concept of approximation and the
definition of some form of adequate approximation. A concept such as
confidence regions is no longer defined as a region containing the true
parameter values with a given probability, but a region which contains those
parameter values which give an adequate approximation (see Davies, (2014) and
also Tukey, (1993).

Models can be obtained from the approximations by modelling the residuals. The
residuals of the riboflavin data look more or less like i.i.d. Gaussian with
of the model based procedures$\sigma$ set to the standard deviation. For the
osteoarthritis data a logit model can be used. The covariates for the leukemia
data can be used for a logit model. The Boston data are more difficult as the
residuals exhibit a structure and there are outliers. In Davies and Dümbgen,
(2022) the sunspot data are modelled using sines and cosines with the
residuals being an AR(1) Gaussian process with parameter 0.5.

###  3.8 Gaussian covariates

As has been repeatedly stated the Gaussian covariate P-values are exact and
valid no matter what the data and what the subset. The memory demands and the
algorithmic complexity are linear in the size of the data. These may be
necessary conditions for good performance but they are not sufficient. To
understand better why it works so well one can investigate conditions under
which it can be shown to work. This is done by means of asymptotics in Section
9 of Davies and Dümbgen, (2022) which gives sufficient conditions on $(n,q)$,
the size $k^{*}$ of the active set ${\mathcal{S}}_{k^{*}}$ and the sizes of
the $\beta_{j}$ for $\boldsymbol{x}_{j}\in{\mathcal{S}}_{k^{*}}$ for the
stepwise procedure to work.

The first step in the stepwise procedures consists of calculating $\max_{1\leq
i\leq q}|\boldsymbol{x}_{i}^{t}\boldsymbol{y}|$ and comparing this with
$\max_{1\leq i\leq q}|\boldsymbol{Z}_{i}^{t}\boldsymbol{y}|$ where the
$\boldsymbol{Z}_{i}$ are i.i.d. Gaussian covariates. Gaussian covariates can
be thought of as the most irrelevant of the universally irrelevant covariates.
If $\boldsymbol{y}=\boldsymbol{Z}$ is just Gaussian noise and the
$\boldsymbol{x}_{i}$ have norm 1 it follows from the Gaussian correlation
inequality (Royen, (2014)) that

$\boldsymbol{P}\bigl{(}\max_{\boldsymbol{x}_{i}}(\boldsymbol{x}_{i}^{t}\boldsymbol{Z})^{2}\geq
x\bigr{)}\ \leq\ \boldsymbol{P}\bigl{(}\max_{i}Z_{i}^{2}\geq x\bigr{)}$

with equality when the $\boldsymbol{x}_{i}$ are uncorrelated. This is used in
the proof of Theorem 3 of Davies and Dümbgen, (2022).

A similar approximate result for fixed $\boldsymbol{y}$ can be obtained by
considering uniform orthogonal rotations of the $\boldsymbol{x}_{i}$ (see
Dümbgen and Davies, (2020)). These retain the geometric structure of the
covariates. The message is that highly correlated covariates
$\boldsymbol{x}_{i}$ reduce the effective value of $q$. The extreme example is
when the $\boldsymbol{x}_{i}$ are all equal in which case the effective value
of $q$ is $q=1$. It may be possible to improve the performance of the Gaussian
covariate procedure by taking the covariance structure of the
$\boldsymbol{x}_{i}$ into account but it is not clear how to do this.

Correlation structure of covariates of data sets

Figure 1: From top left clockwise: riboflavin, osteoarthritis, sunspot, Boston
housing.

Figure 1 shows the correlation structure of the covariates of the four data
sets based on a random choice of 20000 pairs. The Gaussian covariate procedure
performs uniformly well over all data sets. The procedures scalreg, stability,
multi.split and BoomSpikeSlab performed best for the sunspot data of Table 10
in terms of covariates selection. Of these stability performed the worst as it
clearly underfitted. Again the Gaussian covariate procedure performed best in
that it selected fewer covariates but nevertheless had the smallest sum of
squared residuals. It was also by far the fastest procedure. The other three
all failed for the sunspot data: lasso performed worse for this data set than
for the other three, EMVS returned zero as for the other data sets with the
exception of Table 11 where it returned all 500, the initial screening in SIS
failed completely presumably because of the near orthogonality of the
covariates.

###  3.9 Publications and rejections

The model based procedures all of which have failed have an impressive record
of publication in leading statistical journals. Of those mentioned in this
article three were published in the AoS, five in JRSSB, two in JASA, one in
Statistical Science and two in Biometrika.

The Gaussian procedure, the only one which works, also has an impressive
record with rejections from AoS (twice), JRSSB as a discussion paper, JASA,
EJS and AoAS. Clearly statistical journals prefer mathematically demanding and
complicated methods that don’t work to simple ones which do.

###  3.10 Simplicity

The Gaussian covariates procedures are based on (3) which is a special case of
Theorem 1 of Davies and Dümbgen, (2022). The Gaussian $P_{G}$ values follow
directly from this. All inherit the properties of (3), namely that they are
exact and valid no matter what the data, the subset or the covariate. The
calculations require only least squares without needing either a QR or
Singular Value Decomposition. The only exceptions are the subroutines for the
Beta probabilities which are taken from Press et al., (2003). The algorithms
are of linear complexity as are the memory requirements. The result is
overwhelmingly superior to all model based procedures. In short, (3) replaces
about 60 years of failed model based research on covariate selection for
linear regression. There is no reason why the method could not have been
developed 60 years ago, apart from the reluctance of statisticians to consider
anything outside the statistical canon (see Section 3.9).

##  4 Appendix

This gives a short and simple introduction to the Gaussian covariate procedure
of Davies and Dümbgen, (2022). In spite of its simplicity it represents a
complete break with the current statistical paradigm.

Given a dependent variable $\boldsymbol{y}$ of size $n$, covariates
$\boldsymbol{x}_{i},i=1,\ldots,q$ and subset $\mathcal{S}_{k}$ of the
$\boldsymbol{x}_{i}$ of size $k<\min(n,q)$ regress $\boldsymbol{y}$ on the
covariates in $\mathcal{S}_{k}$ to give a sum $\mathrm{rss}_{k}$ of squared
residuals. Choose an $\boldsymbol{x}_{i}\in\mathcal{S}_{k}$ and replace it
with a Gaussian covariate $Z$ to give a random sum of squared residual
$\mathrm{RSS}$. Then

$RSS/rss_{k,-i}\sim\text{Beta}((n-k)/2,1/2))$ (3)

where $\mathrm{rss}_{k,-i}$ is the sum of squared residuals when
$\boldsymbol{y}$ is regressed on
$\mathcal{S}_{k}\setminus\\{\boldsymbol{x}_{i}\\}$. From this and the identity

$\text{F}_{\nu_{1},\nu_{2}}(x)=\text{Beta}_{\nu_{1}/2,\nu_{2}/2}((\nu_{1}/\nu_{2})x/((\nu_{1}/\nu_{2})x+1))$
(4)

it follows after some manipulation that the standard P-value of
$\boldsymbol{Z}$ is uniformly distributed over $(0,1)$. Equation (3) is a
special case of Theorem 1 of Davies and Dümbgen, (2022). It holds whatever the
data and whatever the subset.

It follows from (3) and (4) that the probability that the Gaussian covariate
$\boldsymbol{Z}$ is better than $\boldsymbol{x}_{i}$, that is if
$\mathrm{RSS}<\mathrm{rss}_{k}$, is given by

$\boldsymbol{P}(\mathrm{RSS}<\mathrm{rss}_{k})=\text{Beta}_{(n-k)/2,1/2}(\mathrm{rss}_{k}/\mathrm{rss}_{k,-i})=P_{F}(\boldsymbol{x}_{i})$
(5)

where $P_{F}(\boldsymbol{x}_{i})$ is the standard F P-value based on the
subset $\mathcal{S}_{k}$.

Any subset $\mathcal{S}_{k},k<n$ of the covariates specifies a selection of
covariates and defines a linear approximation to $\boldsymbol{y}$ obtained by
regressing $\boldsymbol{y}$ on the covariates in $\mathcal{S}_{k}$. When
selecting covariates the problem is to avoid including ones which are
irrelevant, defined here as being no better than Gaussian noise. This will be
done by defining a Gaussian P-value as follows. Given a subset
$\mathcal{S}_{k}$ size $k$ and a covariate
$\boldsymbol{x}_{i}\in\mathcal{S}_{k}$ we replace $\boldsymbol{x}_{i}$ and the
remaining $q-k$ covariates by $q-k+1$ by i.i.d. Gaussian covariates
$\boldsymbol{Z}_{j},j=1,\ldots,q-k+1$ . For each $\boldsymbol{Z}_{j}$ we
regress $\boldsymbol{y}$ on
$\mathcal{S}_{k}\setminus\\{\boldsymbol{x}_{i}\\}\cup\\{\boldsymbol{Z}_{j}\\}$
to give a sum of squared residuals $\mathrm{RSS}_{k-1,+j}$. The Gaussian
P-values of $\boldsymbol{x}_{i}$ should measure the degree to which
$\boldsymbol{x}_{i}$ is better than the best of the Gaussian covariates
$\boldsymbol{Z}_{j}$. This is accomplished by defining the the P-value of each
$\boldsymbol{x}_{i}$ by

$\displaystyle P_{G}(\boldsymbol{x}_{i})$ $\displaystyle=$
$\displaystyle\boldsymbol{P}(\min_{j}\mathrm{RSS}_{j}<\mathrm{rss}_{k})=\text{Beta}_{1,q-k+1}(\text{Beta}_{(n-k)/2,1/2}(\mathrm{rss}_{k}/\mathrm{rss}_{k,-i}))$
$\displaystyle=$ $\displaystyle 1-(1-P_{F}(\boldsymbol{x}_{i}))^{q-k+1}.$

The second equality follows form the independence of the
$\mathrm{RSS}_{k-1,+j}$ and (5).

The all subset version of Gaussian covariate selection retains only those
subsets for which each covariate $\boldsymbol{x}_{i}$ in the subset has a
$P_{G}$ value (4) $P_{G}(\boldsymbol{x}_{i})<p0$. The procedure in gausscov
further eliminates all retained subsets which are subsets of some other
retained subset. Those left are ordered according to the sum of squared
residuals. These subsets are maximal in the sense that it is not possible to
include another covariate and still retain the upper bound $p0$. Clearly the
all subset procedure is only possible if $q$ is not too large, say $q\leq 25$.

The stepwise selection procedure based on Gaussian covariates makes use of a
slightly different $P_{G}$-value. If a subset $\mathcal{S}_{k}$ of size $k$
with sum of squared residuals $\mathrm{rss}_{k}$ has already been selected the
next candidate for selection is $\boldsymbol{x}_{b}$ the best of the remaining
$q-k$ covariates. Its $P_{G}$-value is given by

$\displaystyle P_{G}(\boldsymbol{x}_{b})$ $\displaystyle=$
$\displaystyle\boldsymbol{P}(\min_{\ell}\mathrm{RSS}_{k,+\ell}<\mathrm{rss}_{k,+b})=\text{Beta}_{1,q-k}(\text{Beta}_{(n-k-1)/2,1/2}(\mathrm{rss}_{b}/\mathrm{rss}_{k}))$
(7) $\displaystyle=$ $\displaystyle 1-(1-P_{F}(\boldsymbol{x}_{b}))^{q-k}$

where $\mathrm{rss}_{k,+b}$ is the sum of squared residuals when
$\boldsymbol{y}$ is regressed om $\mathcal{S}_{k}\cup\\{\boldsymbol{x}_{b}\\}$
and $\mathrm{RSS}_{k,+\ell}$ is the sum of squared residuals when
$\boldsymbol{y}$ is regressed om
$\mathcal{S}_{k}\cup\\{\boldsymbol{Z}_{\ell}\\}$

The stepwise selection procedure uses the $P_{G}$-values (7). If a subset
$\mathcal{S}_{k}$ has already been selected the best of the remaining
covariates is considered and if its $P_{G}$-value is less than $p0$ it is
selected and the procedure continues. Otherwise the procedure terminates. If
the size of the selected subset does not exceed an upper bound specified by
the user all subsets of the selected subset are considered as above and the
best subset in the sense of least squares is returned. The default maximum
size for all subsets is 20. If this is exceeded the selected subset is
returned with the $P_{G}$-values.

If the intercept is specifically included (the default case) it is denoted
$\boldsymbol{x}_{0}$ in the output,it is always the last covariate and its
$P_{G}$-value is the standard $P_{F}$-value.

It can happen with the stepwise procedure that it fails at the first step as
the first $P_{G}$-value exceeds $p0$. This is the case for the dental data of
Seheult and Tukey, (2001). This can be avoided to some extent by specifying a
minimum size $kmn$ for the initial selected subset. We give some examples
below.

Two further functions are f2st and f3st. The first f2st makes an initial
selection with f1st, eliminates the selected covariates and applies f1st to
those remaining. This is continued until f1st selects no more covariates. If
it is applied to the riboflavin data with kmn=10 it results in 44
approximations involving 132 covariates. The time required is about one
second.

The second f3st again the initial selection is made with f1st but now each of
the selected covariates is eliminated in turn whilst maintaining the others.
f1st is applied again and this it iterated $m$ times where $m$ is specifies by
the user. The initial selection for the riboflavin data consists of the
covariates
$\boldsymbol{x}_{73},\boldsymbol{x}_{2034},\boldsymbol{x}_{2564},\boldsymbol{x}_{4003}$
with a sum of squared residuals 8.45, Table 3. Eliminating
$\boldsymbol{x}_{73}$ results in
$\boldsymbol{x}_{4003},\boldsymbol{x}_{2564},\boldsymbol{x}_{143},\boldsymbol{x}_{2034}$
with sum of squared residuals 9.41, eliminating $\boldsymbol{x}_{4003}$
results in
$\boldsymbol{x}_{73},\boldsymbol{x}_{1131},\boldsymbol{x}_{1278},\boldsymbol{x}_{2140},\boldsymbol{x}_{2564},\boldsymbol{x}_{4006}$
with a sum of the squared residuals 6.21.

We illustrate the stepwise procedure f1st using the Boston Housing data where
the covariates are all interactions of order eight and less. The dimensions
are $(n,q)=(506,176358)$. The results are given in Table 15 for various values
of $kmn$.

Boston Housing data (506,176358)
---
$kmn$ | covariates | $\mathrm{rss}$ | time
0 | 7222 | 170037 | 164793 | 106761 | 82567 | 81454 |  |  |  | 6566 | 2.94
10 | 7222 | 164793 | 106761 | 81454 | 29154 | 170455 | 67241 |  |  | 6130 | 4.59
15 | 7222 | 164793 | 106761 | 81454 | 29154 | 170455 | 67241 | 92 | 23339 | 5589 | 6.50
17 | 7222 | 164793 | 106761 | 81454 | 29154 | 170455 | 92 | 23339 | 170246 | 4711 | 10.1
| 5862 | 5927 | 169947 |  |  |  |  |  |  |  |
Table 15: The stepwise procedure with $p0=0.01$ and for various values of the
size $kmn$ of the initial selection. The time is given in seconds.

##  5 Acknowledgements

Thanks to John Hodgson for pertinent comments on earlier versions. He is not
responsible for any nonsense still standing.

## References

  * Alizadeh et al., (2000) Alizadeh, A., Eisen, M., Davis, R., Ma, C., Lossos, I., Rosenwald, A., Boldrick, J., Sabet, H., Tran, T., Yu, X., Powell, J., Yang, L., Marti, G., Moore, T., Hudson, J. J., Lu, L., Lewis, D., Tibshirani, R., Sherlock, G., Chan, W., Greiner, T., Weisenburger, D., Armitage, J., Warnke, R., Levy, R., Wilson, W., Grever, M., Byrd, J., Botstein, D., Brown, P., and Staudt, L. (2000).  Distinct types of diffuse large b-cell lymphoma identified by gene expression profiling.  Nature, 403:503–511.
  * Barber et al., (2020) Barber, R., Candes, E., Janson, L., Patterson, E., and Sesia, M. (2020).  The knockoff filter for controlled variable selection.  https://CRAN.R-project.org/package=knockoff.
  * Benjamini and Hochberg, (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: A practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society: Series B, 57:289–300.
  * Bühlmann, (2013) Bühlmann, P. (2013).  Statistical inference in high dimensional models.  Bernoulli, 19:1212–1242.
  * Bühlmann et al., (2014) Bühlmann, P., Kalisch, M., and Meier, L. (2014).  High-dimensional statistics with a view toward applications in biology.  Annual Review of Statistics and Its Applications, 1(1):255–278.
  * Candès et al., (2018) Candès, E., Fan, Y., Janson, L., and Lv, J. (2018).  Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection.  JRSSB, 80(3):551–577.
  * Cox and Battey, (2017) Cox, D. R. and Battey, H. S. (2017).  Large numbers of explanatory variables, a semi-descriptive analysis.  Proc. Natl. Acad. Sci. USA, 114(32):8592 8595.
  * Davies, (2014) Davies, L. (2014).  Data Analysis and Approximate Models.  Monographs on Statistics and Applied Probability 133. CRC Press.
  * Davies and Dümbgen, (2022) Davies, L. and Dümbgen, L. (2022).  Covariate selection based on an model-free approach to linear regression with exact probabilities.  arxiv.org/2202.01553.
  * Dettling and Bühlmann, (2002) Dettling, M. and Bühlmann, P. (2002).  Supervised clustering of genes.  Genome Biology, 3(2):1–15.
  * Dezeure et al., (2015) Dezeure, R., Bühlmann, P., Meier, L., and Meinshausen, N. (2015).  High-dimensional inference: confidence intervals, p-values and R-software hdi.  Statistical Science, 30(4):533–558.
  * Dümbgen and Davies, (2020) Dümbgen, L. and Davies, L. (2020).  Connecting model-based and model-free approaches to linear least squares regression.  arxiv.org/abs/1807.09633.
  * Fan and Lv, (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultrahigh dimensional feature space.  Journal of the Royal Statistical Society: Series B, 70(5):849–911.
  * Feng et al., (2020) Feng, Y., Fan, J., Saldana, D., and Samworth, R. (2020).  Sis: Sure independence screening.  https://CRAN.R-project.org/package=SIS.
  * Friedman et al., (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation withthe graphical lasso.  Biostatistics, 9(3):432–441.
  * Friedman et al., (2019) Friedman, J., Hastie, T., and Tibshirani, R. (2019).  Graphical lasso: Estimation of gaussian graphical models.  https://CRAN.R-project.org/package=glasso.
  * Friedman et al., (2021) Friedman, J., Hastie, T., Tibshirani, R., Narasimhan, B., Tay, K., Simon, N., and Qian, J. (2021).  Lasso and elastic-net regularized generalized linear models.  https://CRAN.R-project.org/package=glmnet.
  * Golub et al., (1999) Golub, T., Slonim, D., P., T., Huard, C., Gaasenbeek, M., Mesirov, J., Coller, H., Loh, M., Downing, J., Caligiuri, M., Bloomfield, C., and Lander, E. (1999).  Molecular classification of cancer: class discovery and class prediction by gene expression monitoring.  Science, 286(15):531–537.
  * Gordon, (1986) Gordon, R., editor (1986).  The American Business Cycle: Continuity and Change, volume 25 of National Bureau of Economic Research Studies in Business Cycles.  Univerisity of Chicago Press.
  * Harrison and Rubinfeld, (1978) Harrison, D. and Rubinfeld, D. (1978).  Hedonic prices and the demand for clean air.  J. Environ. Economics and Management, 5:81–102.
  * Huber, (2011) Huber, P. J. (2011).  Data Analysis.  Wiley, New Jersey.
  * Jiang et al., (2021) Jiang, H., Fei, X., Liu, H., Roeder, K., Lafferty, J., Wasserman, L., Li, X., and Zhao, T. (2021).  huge: High-dimensional undirected graph estimation.  https://CRAN.R-project.org/package=huge.
  * Laszkiewicz, (2021) Laszkiewicz, M. (2021).  Thresholded adaptive validation: Tuning the graphical lasso for graph recovery.  https://github.com/MikeLasz/thav.glasso.
  * Laszkiewicz et al., (2021) Laszkiewicz, M., Fischer, A., and Lederer, J. (2021).  Thresholded adaptive validation:tuning the graphical lasso for graph recovery.  https://arxiv.org/pdf/2005.00466.pdf.
  * Lockhart, (2017) Lockhart, R. (2017).  Inference in high-dimensional linear models course notes.  httpDimensionals://www.sfu.ca/ lockhart/richard/Cambridge/Notes.pdf.
  * Lockhart et al., (2014) Lockhart, R., Taylor, J., Tibshirani, R. J., and Tibshirani, R. (2014).  A significance test for the lasso.  Ann. Statist., 42(2):413–468.
  * Meier et al., (2021) Meier, L., Dezeure, R., Meinshausen, N., Maechler, M., and Buehlmann, P. (2021).  hdi: High-dimensional inference.  https://CRAN.R-project.org/package=hdi.
  * Meinshausen and Bühlmann, (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the lasso.  Annals of Statistics, 34(3):1436–1462.
  * Meinshausen and Bühlmann, (2010) Meinshausen, N. and Bühlmann, P. (2010).  Stability selection.  Journal of the Royal Statistical Society: Series B, 72:1436–146.
  * Press et al., (2003) Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2003).  Numerical Recipes in Fortran 77: The Art of Scientific Computing, volume 1.  Cambridge University Press, second edition.
  * Rockova, (2021) Rockova, V.and Moran, G. (2021).  Emvs: The expectation-maximization approach to bayesian variable selection.  https://CRAN.R-project.org/package=EMVS.
  * Rockova and George, (2014) Rockova, V. and George, E. (2014).  The em approch to Bayesian variable selection.  Journal of the American Statistical Association, 109(506).
  * Royen, (2014) Royen, T. (2014).  A simple proof of the Gaussian correlation conjecture extended to some multivariate gamma distributions.  Far East Journal of Theoretical Statistics, 48(2):139–145.
  * Scott, (2021) Scott, S. (2021).  Boomspikeslab: MCMC for spike and slab regression.  https://cran.r-project.org/web/packages/BoomSpikeSlab/index.html.
  * Seheult and Tukey, (2001) Seheult, A. H. and Tukey, J. W. (2001).  Towards robust analysis of variance.  In Saleh, A. K. M. E., editor, Data Analysis from Statistical Foundations: A Festschrift in Honor of the 75th Birthday of D.A.S. Fraser, pages 217–244. Nova Science, New York.
  * SILSO, (2020) SILSO (2020).  The international sunspot number.  International Sunspot Number Monthly Bulletin and online catalogue.  Royal Observatory of Belgium, avenue Circulaire 3, 1180 Brussels, Belgium.
  * Sun, (2019) Sun, T. (2019).  Scaled sparse linear regression.  https://CRAN.R-project.org/package=scalreg.
  * Sun and Zhang, (2012) Sun, T. and Zhang, C.-H. (2012).  Scaled sparse linear regression.  Biometrika, 99:879–898.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  J. Royal. Statist. Soc B., 58(1):267–288.
  * Tibshirani et al., (2016) Tibshirani, R., Taylor, J., Lockhardt, R., and Tibshiranis, R. (2016).  Exact post-selection inference for sequential regression procedures.  Journal of the American Statistical Association, 111.
  * Tukey, (1993) Tukey, J. W. (1993).  Issues relevant to an honest account of data-based inference, partially in the light of Laurie Davies’s paper.  Princeton University, Princeton.
  * Wasserman and Roeder, (2009) Wasserman, L. and Roeder, K. (2009).  High dimensional variable selection.  Annals of Statistics, 37(5A):2178–2201.
  * Zhang and Zhang, (2014) Zhang, C.-H. and Zhang, S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  J. R. Stat. Soc. Ser. B. Stat. Methodol., 76:217–242.