\section*{Acknowledgements}
\vspace{-7mm}
We acknowledge the financial support by the Ministry of Science and Technology (Grants No. MOST-107-2112-M006-025-MY3, No. MOST-108-2638-M-006-002-MY2,
No. MOST 110-2112-M-006-020-, and No. MOST 110-2112- M-003-008-MY3) and National Center for Theoretical Sciences in Taiwan.
C. O. acknowledges support from a VIDI grant (Project 680-47-543) financed by the Netherlands Organization for Scientific Research (NWO).
C.-H. C. acknowledges support from the Yushan Young Scholar Program under the Ministry of Education (MOE) in Taiwan.
We acknowledge technical assistance by Botsz Huang and You-Ting Huang.
The authors gratefully acknowledge many helpful discussions with Szu-Chao Chen, Hsiu-Chuan Hsu, Ion Cosma Fulga, Chao-Cheng Kaun, Chao-Ping Hsu, and Chong-Der Hu.

\section*{References}


\section*{Acknowledgement}

This research is funded by University of Information Technology-Vietnam National University HoChiMinh City under grant number D1-2021-15.

\bibliographystyle{plain}

\section{Introduction}
In the era of technology development, the trend of online shopping on e-commerce sites is rapidly increasing. Consumers can easily express their opinions or feelings by leaving comments to resolve the complaint issues. Manually processing responses takes a significant amount of time and effort, hence identifying complaints should be automated. According to Olshtain and Weinbach
1987 \cite{olshtain198710}, complaining is a basic speech act used to express a negative mismatch between the expectation and reality towards a state of affairs, products, organizations, or events.

Datasets used for complaint detection are now not so much, especially in Vietnam; there is no dataset about that so far. Therefore, in this paper, we present a human-annotated dataset for this task to automatically identify complaints and contribute to a more robust Vietnamese dataset. After building the dataset, we implement machine learning-based methods as baseline systems to evaluate our dataset.

Our contributions in this paper include:

\begin{itemize}
    \item We introduce a Vietnamese dataset named UIT-ViOCD, which contains 5,485 human-annotated reviews for identifying complaints on the open-domain such as fashion, cosmetics, applications, and phones.

    \item We conduct multiple experiments on different machine learning-based models, including traditional machine learning, deep learning, transfer learning models with low to high complexity and analyze the quantitative dataset. The highest performance is reached by PhoBERT with 92.16\% F1-score for complaint speech detection. We also implement different experiments to figure out the role of Vietnamese structure ingredients.
\end{itemize}

We organize the rest of the paper as follows. Section 2 presents the works related to our task. The process of building and analyze the dataset is described in Section 3. In Section 4, we propose our methodology and discuss experimental results. Finally, we draw a conclusion and future work in Section 5.

\section{Related Work}
In 1987, Olshain et al. \cite{olshtain198710} stated that complaints are the result when a speaker is adversely affected or expects a favorable event to occur. Complaints can also be described in a positive sense as a report from a consumer documenting a problem with a product or service \cite{felix2016effectiveness}. Complaints usually co-occur with other speech acts such as warnings, threats, suggestions, or advice. Moreover, Cohen et al. \cite{cohen1993production} have found that the complaint often occurs as "speech act set".

In the world, there are various studies about customer complaints in different languages. In 2002, Pang et al. \cite{pang2002thumbs} used supervised machine learning methods to categorize film reviews with datasets available from the rec.arts.movies.reviews archive. In 2004, Hu and Liu \cite{hu2004mining} did the research with an aim to mine and summarize all the customer reviews of products, namely, they mined the features of the products collected from the reviews on Amazon.com and C|net.com, on which the customers have expressed their opinions positively or negatively. Moreover, Gaman and N. Aletras used Distant Supervision, deep learning methods to identify Twitter complaints in English with nine domains \cite{preotiuc2019automatically}.

In Vietnamese, the current datasets mostly focus on the emotional recognition problem. Minh et al. researched aspect-based emotional analysis and published a dataset including 7,828 restaurant reviews for solving two sub-tasks: aspect detection and polarity detection \cite{nguyen2019corpus}. In 2019, Vong et al. built a standard Vietnamese Social Media Emotion Corpus (UIT-VSMEC) consisting of 6,927 human-annotated comments with six emotion labels \cite{ho2019emotion}. In 2021, Nguyen et al. published a Vietnamese dataset about constructive and toxic speech detection named UIT-ViCTSD \cite{nguyen2021constructive}.  Nonetheless, there is no dataset for identifying complaints so far. Therefore, we decide to build a dataset about customer complaints about automatically classifying comments into complaints or non-complaints.
\section{Dataset}
We collect customer product reviews from e-commerce sites and Google Play Store\footnote{https://play.google.com/store} with 5,485 reviews belonging to four areas: fashion, cosmetics, applications, and phones.
\subsection{Task Definition And Guidelines}
Based on the complaint definition proposed by Olshtain and Weinbach \cite{olshtain198710}, we build annotation guidelines for our task.

\textbf{Complaint (Label 1):} Complaining reviews express customers of dissatisfaction between reality and expectations, often accompanied by the build-up and desire to resolve. Complaints usually co-occur with other speech acts such as warnings, threats, suggestions, or advice \cite{preotiuc2019automatically}. Besides, complaining reviews sometimes contain compliments in addition to complaints. Complaints are constructive, not showing hateful emotions or heavy criticism regarding the service or product. Punctuation and special characters, icons can also be physical signs that show complaints.

\textbf{Non-complaint (Label 0):} Non-complaining reviews are compliments, praise comments, satisfaction about products and services.

\begin{table}[hbt!]
\centering
\caption{Examples of reviews annotated for complaint task}
\label{tab:examples_label}
\begin{tabular}{clcc}
\hline
\multicolumn{1}{l}{\textbf{No}} & \multicolumn{1}{c}{\textbf{Comment}}                                                                                & \textbf{Label} \\ \hline
1                                & \begin{tabular}[c]{@{}l@{}}Phải chi có tiếng việt nữa thì tuyệt vời.\\ \textit{(It's great to have Vietnamese.)}\end{tabular}                                           & 1\\
2                                & \begin{tabular}[c]{@{}l@{}}Chất lượng rất tốt, thời gian xử lý đơn hàng và giao hàng quá lâu.\\ \textit{(Very good quality, but take too much time for processing and delivery.)}\end{tabular}                                           & 1\\
3                                & \begin{tabular}[c]{@{}l@{}}Không thấy quà tặng theo là như nào lừa đảo à ???\\ \textit{(Don't see gifts, scam ???)}\end{tabular}                                        & 1 \\
4                                & \begin{tabular}[c]{@{}l@{}}Game như L** đừng có tải nha\\ \textit{(Game likes P*ssy don't download it)}\end{tabular}                                                           & 0\\
5                                & \begin{tabular}[c]{@{}l@{}}Mình đã sd đến chai t5 của srm senka r rất  ưnggg\\ (\textit{I have used the 5th bottle of senka cleanser so satisfied)}\end{tabular} &                   0\\ \hline
\end{tabular}
\end{table}

Table \ref{tab:examples_label} presents examples of reviews highlighting differences between complaints and non-complaints in our dataset. In the first sentence, the assessment sentence, which sets out unmet customer wishes but does not show negative sentiment, is constructive. The second sentence has both praise and dissatisfaction and is classified as a complaint. Meanwhile, the third sentence shows mixed negative sentiment (angry emotion) and expressed a complaint. In the fourth sentence, we see an insult that implies negative sentiment but is not constructive for the business organizations; hence, this sentence can not be labeled as a complaint. The fifth sentence shows a positive and satisfied attitude to the product.

\subsection{Annotation Process}
It is necessary to know the task for building a high-quality and reliable dataset, thus we create a guideline for annotators before labeling data. The annotator team consists of 3 people who have enough knowledge about natural language processing and labeling data. For evaluating the inter-agreement among annotators, we use the Fleiss' Kappa \cite{bhowmick2008agreement}, a statistical technique according to the formula (1).
\textit{}\begin{center}$$\mathrm{A}_{\mathrm{m}}=\frac{\mathrm{P}_{o}-\mathrm{P}_{e}}{1-\mathrm{P}_{e}} \quad\quad(1)$$  \end{center}

Where, the observed agreement ${P_o}$ is the proportion of sentences in which both of the annotators agreed on the class pairs. The chance agreement ${P_e}$ is the proportion of items for which agreement is expected by chance when the sentences are seen randomly. First, we label a part of the dataset. Then, we measure the ${A}_{m}$ and update guidelines. We repeat this process until the inter-agreement of annotators reaches over 80\%. After getting that score, annotators have enough knowledge about the task, and we begin labeling the rest of the data consecutively.

\subsection{Dataset Analysis}
There is no significant difference between 2,854 complaining reviews (52.03\%) and 2,631 non-complaining reviews (47.97\%).

To have a deeper understanding of the UIT-ViOCD dataset, we analyze the distribution of the complaints dataset according to the length of the sentence, and the results are shown in Table \ref{tab:len_of_sentence}. The number of sentences whose length is greater than 30 and makes up most of the dataset, especially in Label 1, occupies up to 46.45\%. It shows that customers aim to express their desires and opinions in detail and exact.

\begin{table}[hbt!]
\centering
\caption{Distribution of the complaint-annotated dataset according to the length of the sentence (\%).}
\label{tab:len_of_sentence}
\begin{tabular}{cccc}
\hline
\textbf{Length}           & \textbf{Label 0} & \textbf{Label 1} & \textbf{Overall} \\ \hline
\textbf{1-5}              & 0.40              & 0.04             & 0.44           \\
\textbf{6-10}             & 1.42             & 0.27             & 1.69           \\
\textbf{11-15}            & 4.41             & 0.86             & 5.27           \\
\textbf{16-20}            & 3.94             & 1.22             & 5.16           \\
\textbf{21-25}            & 2.94             & 1.22             & 4.16           \\
\textbf{26-30}            & 3.43             & 1.97             & 5.40            \\
\textbf{\textgreater{}30} & \textbf{31.43}            & \textbf{46.45 }           & \textbf{77.88}          \\
\textbf{Total}            & 47.97             & 52.03             & 100.00            \\
\hline
\end{tabular}
\end{table}

Most customers tend to complain about quality, delivery service, and price. The application field is slightly different from other areas, and customers often show a positive attitude to improve the application. The rest of the sectors tend to be straight to offer dissatisfaction.

We perform a vocabulary-based analysis experiment to understand linguistic features between complaints and non-complaints. This experiment is Part-of-speech, Word Frequency Dictionary, and Embedding Extraction from the identifying complaint task. We assign part-of-speech to each word in every review of our dataset, such as nouns, verbs, adjectives, using ViPosTagger (Vietnamese pos tagging F1\_score = 92.5\%) of Pyvi library\footnote{https://pypi.org/project/pyvi/}. We create a frequency dictionary to map the word and the class it appeared to in (complaint and non-complaint) to the number of times that word appeared in its corresponding class. Top part-of-speech features that belong to complaints and non-complaint are analyzed about word feature similarity using Cosine similarity \cite{li2013distance} (2). Where, $\vec{a}$ and $\vec{b}$ are two vectors of two words extracted from the embedding word model of the this task. Table  \ref{big_table} shows the results of this experiment.

$$\operatorname{Cos} (\vec{a},  \vec{b})=\frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}=\frac{\sum_{1}^{n} a_{i} b_{i}}{\sqrt{\sum_{1}^{n} a_{i}^{2}} \sqrt{\sum_{1}^{n} b_{i}^{2}}}     \quad\quad(2)$$

With the results in Table \ref{big_table}, linguistic features from part-of-speech are distinctive of complaints. In detail, several part-of-speech patterns such as coordinating conjunction and adverb are important Vietnamese complaints classification features. Linguistics research identified that a significant proportion of complaints tended to juxtapose overall negative evaluation with some positive appraisal \cite{vasquez2011complaints}. Definitely, "nhưng" (but) is the word that appears the most in complaining reviews and goes with the good and the bad of the product (e.g. "giao hàng chậm nhưng sản phẩm ok, nguyên seal.",  means: "slow delivery but the product is ok, fully sealed."). Moreover, "không" (not, no) often accompanies words that describe the service or product issue. Furthermore, "rất" (very) comes with good characteristics.

Besides, we realize that keywords with complaint meanings (top words) usually go with negative words. Likewise, words that go with non-complaint keywords often have not negative sense.
\begin{table}[hbt!]
\caption{Vietnamese abbreviations in the UIT-ViOCD dataset.}
\centering
\begin{tabular}{|c|c|c|}
\hline
\textbf{Abbreviations} & \textbf{Normalization} & \textbf{English meaning} \\ \hline
"dx", "dc", "dk" & "được" &"ok" \\
\hline
"dt" & "điện thoại" & "phone" \\
\hline
"god" & "tốt" & "good" \\
\hline
"ship" & "giao hàng" & "ship" \\ \hline
\end{tabular}
\end{table}
\begin{table}[]

\caption{The experiment results: Part-of-speech (POS), Word Frequency Dictionary, and Embedding Extraction from Identifying Complaint Task.}
\centering
\label{big_table}
\begin{tabular}{|c|c|l|c|}
\hline
\multicolumn{4}{|c|}{\textbf{Complaint}}                                                                                \\ \hline
\textbf{POS}                        & \textbf{Top words}                                                   & \multicolumn{1}{c|}{\textbf{Words have high similarity scores to top-word}}                                                                                & \textbf{Similarity score}                                        \\ \hline
\multirow{3}{*}{\textbf{Verb}}      & \begin{tabular}[c]{@{}c@{}}cập nhật\\ (update)\end{tabular}         & \begin{tabular}[c]{@{}l@{}}hôi cực (smell), trễ (late), lựa chọn (choose), xấu (bad),\\càng (more)\end{tabular}                                  & \begin{tabular}[c]{@{}c@{}}[0.294,  0.384]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}khắc phục\\ (fix)\end{tabular}           & \begin{tabular}[c]{@{}l@{}}tẩy trang (cleansing), không thể (can’t), trôi (drift), thanh\\toán (payment), lâu (slow)\end{tabular}                                         & \begin{tabular}[c]{@{}c@{}}[0.070,  0.104]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}thất vọng\\    (disappoint)\end{tabular} & \begin{tabular}[c]{@{}l@{}}mỏng (thin), nhưng (but), trừ (minus), đề nghị (suggest),\\tẩy trang (cleansing)\end{tabular}              & \begin{tabular}[c]{@{}c@{}}[0.089,  0.119]\end{tabular}  \\ \hline
\multirow{2}{*}{\textbf{Adjective}} & \begin{tabular}[c]{@{}c@{}}tệ \\    (bad)\end{tabular}               & \begin{tabular}[c]{@{}l@{}}sơ sài (incomplete), mỏng (thin), cồn (alcohol), rách (torn),\\chậm (slow), nhưng (but)\end{tabular}                               & \begin{tabular}[c]{@{}c@{}}[0.029,  0.082]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}sai \\    (wrong)\end{tabular}            & \begin{tabular}[c]{@{}l@{}}lộn (mismatch), lừa đảo (cheat), tệ (bad), mỏng (thin), sơ\\sài (imcomplete), hư hỏng (broken), kém, xấu (bad),\\không (not)\end{tabular}    & \begin{tabular}[c]{@{}c@{}}[0.077,  0.125]\end{tabular}  \\ \hline
\multirow{2}{*}{\textbf{Noun}}      & \begin{tabular}[c]{@{}c@{}}lỗi \\    (error)\end{tabular}            & \begin{tabular}[c]{@{}l@{}}chậm (slow), gửi (send), móp (dented), bắt (must), thanh\\toán (payment), nhu cầu (need), chán (bored), khắc phục\\(fix)\end{tabular}         & \begin{tabular}[c]{@{}c@{}}[0.105,  0.124]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}trò chơi\\ (game)\end{tabular}           & \begin{tabular}[c]{@{}l@{}}đề nghị (suggest), đéo (fuck), thất vọng (disappoint), không\\(not)\end{tabular}                                                               & \begin{tabular}[c]{@{}c@{}}[0.315,  0.377]\end{tabular}  \\ \hline
\multirow{2}{*}{\textbf{Others}}  & \begin{tabular}[c]{@{}c@{}}không\\ (not, no)\end{tabular}            & \begin{tabular}[c]{@{}l@{}}chậm (slow), gửi (send), móp (dented), bắt (must), thanh\\toán (payment), nhu cầu (need), chán (bored), khắc phục\\(fix)\end{tabular}         & \begin{tabular}[c]{@{}c@{}}[0.107,  0.174]\end{tabular} \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}nhưng\\ (but)\end{tabular}                & \begin{tabular}[c]{@{}l@{}}sơ sài (incomplete), rách (torn), gửi (send), mỏng (thin),\\thất vọng (disappoint), cồn (alcohol), chậm (slow), tệ\\(bad)\end{tabular}         & \begin{tabular}[c]{@{}c@{}}[0.082,  0.115]\end{tabular} \\ \hline
\multicolumn{4}{|c|}{\textbf{Non-Complaint}}                                                                                \\ \hline
\textbf{POS}                        & \textbf{Top words}                                                   & \multicolumn{1}{c|}{\textbf{Words have high similarity scores to top-word}}                                                                                & \textbf{Similarity score}                                               \\ \hline
\multirow{2}{*}{\textbf{Verb}}      & \begin{tabular}[c]{@{}c@{}}thích \\    (like)\end{tabular}           & \begin{tabular}[c]{@{}l@{}}êm (smooth), kỹ lưỡng (elaborate), cảm ơn (thank), tiện\\dụng (comfortable)\end{tabular}                                & \begin{tabular}[c]{@{}c@{}}[0.294,  0.384]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}cảm ơn\\ (thank)\end{tabular}            & \begin{tabular}[c]{@{}l@{}}xinh , dễ thương (lovely), nhờ (thanks to), yêu (love), rẻ\\(cheap), nhanh chóng (fast)\end{tabular}                & \begin{tabular}[c]{@{}c@{}}[0.071,  0.119]\end{tabular}  \\ \hline
\multirow{3}{*}{\textbf{Adjective}} & \begin{tabular}[c]{@{}c@{}}đẹp\\ (beautiful)\end{tabular}            & \begin{tabular}[c]{@{}l@{}}rẻ (cheap), biết (know), nhờ (thanks to), mua (buy), tiện\\(convenient), tốt, xịn (good)\end{tabular}                                            & \begin{tabular}[c]{@{}c@{}}[0.102,  0.137]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}nhanh \\    (fast)\end{tabular}           & \begin{tabular}[c]{@{}l@{}}hài lòng (sastified), dễ thương (lovely), thích (like), thương\\(love), xuất sắc (excellent)\end{tabular}                                     & \begin{tabular}[c]{@{}c@{}}[0.079,  0.117]\end{tabular}  \\ \cline{2-4}
                                    & \begin{tabular}[c]{@{}c@{}}đúng \\    (right)\end{tabular}           & \begin{tabular}[c]{@{}l@{}}vợ (wife), ưng ý (like), kỹ càng (elaborate), dễ thương\\(lovely), cảm ơn (thank)\end{tabular}                           & \begin{tabular}[c]{@{}c@{}}[0.189, 0.245]\end{tabular}  \\ \hline
\textbf{Noun}                       & \begin{tabular}[c]{@{}c@{}}sao \\    (star)\end{tabular}             & vui vẻ (happy), kĩ càng (elaborate)                                                                                & \begin{tabular}[c]{@{}c@{}}[0.156,  0.166]\end{tabular}  \\ \hline
\textbf{Others}                   & \begin{tabular}[c]{@{}c@{}}rất \\    (very)\end{tabular}             & \begin{tabular}[c]{@{}l@{}}mình (I), êm (smooth), đi (working), điện thoại\\(smartphone), tuyệt vời (excellent), rẻ (cheap), xài (use),\\đầy đủ (enough)\end{tabular} & \begin{tabular}[c]{@{}c@{}}[0.270, 0.312]\end{tabular}  \\ \hline
\end{tabular}
\end{table}

\subsection{Comparison with Other Vietnamese Datasets}
Table \ref{tab:datasetcomp} presents many Vietnamese datasets for different tasks, which is enough to build information recognition features on the social listening system. Compared to other Vietnamese datasets, our dataset (UIT-ViOCD) is the first Vietnamese dataset for complaint classification.

\begin{table}[H]
\centering
\caption{Vietnamese datasets for evaluating different text classification tasks.}
\label{tab:datasetcomp}
\begin{tabular}{|l|r|l|r|l|}
\hline
\multicolumn{1}{|c|}{\textbf{Name}} & \multicolumn{1}{c|}{\textbf{Year}} & \multicolumn{1}{c|}{\textbf{Task}}             & \multicolumn{1}{c|}{\textbf{Size}} & \multicolumn{1}{c|}{\textbf{Domain}} \\ \hline
SA-VLSP2016 \cite{nguyen2018vlsp}                         & 2016                               & Sentiment Analysis                             & 12,196                             & Open-domain                          \\ \hline
SA-VLSP2018 \cite{nguyen2018vlsp}                         & 2018                               & Aspect-based Sentiment Analysis                & 10.351                             & Hotel and Restaurant                 \\ \hline
UIT-VSFC \cite{van2018uit}                            & 2018                               & Sentiment Analysis                             & 16,175                             & Education                            \\ \hline
UIT-VSMEC \cite{ho2019emotion}                           & 2019                               & Emotion Recognition                            & 6,927                              & Open-domain                          \\ \hline
ReINTEL2020 \cite{le2020reintel}                         & 2020                               & Responsible Information Identification         & 10,007                             & Open-domain                          \\ \hline
UIT-ViCTSD \cite{nguyen2021constructive}                          & 2021                               & Constructive and Toxic Detection & 10,000                             & Open-domain                          \\ \hline
ViHSD \cite{luu2021large}                           & 2021                               & Hate Speech Detection                          & 33,400                             & Open-domain                          \\ \hline
UIT-ViSFD \cite{phan2021sa2sl}                           & 2021                               & Aspect-based Sentiment Analysis                & 11,122                             & Smartphone                           \\ \hline
UIT-ViOCD (Ours)                    & 2021                               & Complaint Detection                            & 5,485                              & Open-domain                          \\ \hline
\end{tabular}
\end{table}
\section{Experiments}
We implement various methods on the UIT-ViOCD dataset and then evaluate its quality by four evaluation metrics such as macro-averaged f1 score, macro-averaged precision, macro-averaged recall, and accuracy \cite{scikit-learn}. We use macro for every metric because our dataset is not imbalanced.
\subsection{Data Pre-processing and Preparation}
We conduct the pre-processing phase to normalize the raw data: lowercasing all texts, removing the HTML tags, extracting white spaces, and applying Unicode standardization. There are many freestyle words in the reviews, like teen code, extending the last sound or acronyms, error typing, or Vietnamese sentences mixed with a few English words. For that reason, we normalize them by accurate words (e.g., "rồii": "rồi" (already)). The old type and the new type are transferred to a unified form (e.g., "òa": "oà"). Moreover, we convert emotional icons to equivalent meaning in Vietnamese. After that, we replace abbreviations with complete words and correct the misspelling. The table \ref{tab:abbreviations} shows some examples. Unlike English, white space between two Vietnamese words is the sign for separating syllables, not words. We tokenize every review using Pyvi library for Vietnamese.

After pre-processing data, we split the dataset into training, validation, and test sets with an 80:10:10 ratio by train\_test\_split function of scikit-learn \cite{scikit-learn}. The results shown in this paper are reported on the test set.

\begin{table}[H]
\centering
\caption{Vietnamese abbreviations in the UIT-ViOCD dataset.}
\label{tab:abbreviations}
\begin{tabular}{ccc}
\hline
\textbf{Abbreviations} & \textbf{Normalization} & \textbf{English meaning} \\ \hline
"dx", "dc", "dk"       & "được"                 & "ok"                     \\
"dt"                   & "điện thoại"           & "phone"                  \\
"god"                  & "tốt"                  & "good"                   \\
"ship"                 & "giao hàng"            & "ship"                   \\ \hline
\end{tabular}
\end{table}
\subsection{Experimental Methods}
\subsubsection{Experimental Methods}We conduct various experiments on different models, including traditional machine learning, deep learning, and transfer learning.

After we convert every review to one hot vector, we encode them mapped into embeddings. In particular, each word is represented by a 64-dimensional word embedding. The embedded representation feeds into three models: Logistic Regression, MLP, and BiLSTM.
\paragraph{Traditional Machine Learning Systems}
      \begin{itemize}
          \item  \textbf{Logistic Regression:} Logistic Regression \cite{hilbe2009logistic} is a binary classification method, as a one-layer neural network with a sigmoid activation function for the prediction task.
       \end{itemize}
\paragraph{Deep Learning Systems}
        \begin{itemize}
             \item  \textbf{MLP (Multi-layer Perceptron):}
 MLP \cite{hornik1989multilayer} has two dense hidden layers (D=512, D=64) with ReLU activation function and the last dense hidden layer (D=1) with sigmoid activation function.
            \item  \textbf{BiLSTM (Bidirectional Long Short-term Memory):} LSTM was developed and achieved superior performance (Hochreiter and Schmidhuber, 1997) \cite{hochreiter1997lstm}. LSTM is a variation of a recurrent neural network, which has an input gate, an output gate, a forget gate, and a cell. In 2005, Graves et al. \cite{graves2005bidirectional,graves2005framewise} studied BiLSTM in various classification tasks and received good results compared to other neural network architectures. Besides, BiLSTM which consists of two LSTMs has two inputs from forward direction and backward direction. Therefore, BiLSTM can learn contextual information extracted from two directions. Moreover, BiLSTM can learn more features and perform slightly better than LSTM in emotion classification task \cite{zhou2018nlp}. In this paper, we use BiLSTM having 128 units for each.

        \end{itemize}
        In addition, we use two pre-training word embedding models (PhoW2V and fastText):
        \begin{itemize}

          \item \textbf{PhoW2V:} PhoW2V is published by Nguyen et al. \cite{nguyen2020phobert}, that were pre-trained on a dataset of 20GB Vietnamese texts. In this paper, we use PhoW2V syllables 300 dims and PhoW2V words 300 dims.
          \item \textbf{fastText:} fastText is released by Grave et al. \cite{grave2018learning}. This model is known for achieving good performance in word performance and sentence classification because they use character-level representations.
           \end{itemize}
\paragraph{Transfer Learning Systems}
        \begin{itemize}
          \item  \textbf{PhoBERT:} Pre-trained PhoBERT models are the first public large-scale monolingual language models pre-trained for Vietnamese (Nguyen et al. 2020) \cite{nguyen2020phobert}. This model approach is based on the same RoBERTa \cite{liu2019roberta}, which eliminates the Next Sentence Prediction (NSP) task of the pre-trained model BERT. There are two PhoBERT versions: base and large. In this paper, we use PhoBERT base and employ RDRSegmenter of VnCoreNLP \cite{vu2018vncorenlp} to pre-process the pre-training data and Transformers\footnote{https://huggingface.co/transformers/} (Hugging Face) for loading pre-trained BERT.
        \end{itemize}
We train the networks using the Adam optimizer by minimizing the binary cross-entropy.
\subsubsection{Results and Discussion}
As the results in Table \ref{result_model}, most models achieve the art performance state with upper 85\% of the F1-score. The best performing model is PhoBERT with 92.16\% F1-score. These results prove that the dataset has featured enough to help the model classify well.
\begin{table}[H]
\centering
\caption{Results of implemented models on the UIT-ViOCD dataset (\%).}
\label{result_model}
\begin{tabular}{llrrr}
\cline{2-5}
 & \multicolumn{1}{c}{\textbf{Model}} & \textbf{F1-score} & \textbf{Recall} & \textbf{Precision} \\ \cline{2-5}
 & Logistic Regression                & 89.62             & 88.17           & 91.11              \\ \cline{2-5}
 & MLP                                & 90.97             & 90.32           & 91.64              \\
 & BiLSTM                             & 89.55             & 92.11           & 87.12              \\
 & BiLSTM + PhoW2V\_words             & 88.46             & 90.68           & 86.35              \\
 & BiLSTM + PhoW2V\_syllables         & 86.36             & 89.61           & 83.33              \\
 & BiLSTM + fasttext                  & 86.33             & 92.83           & 80.69              \\ \cline{2-5}
 & \textbf{PhoBERT}                   & \textbf{92.16}    & \textbf{90.63}  & \textbf{94.20}     \\ \cline{2-5}
\end{tabular}
\end{table}
\subsection{Ablation Tests}
We use the ablation technique eliminating structure ingredients for a Vietnamese sentence to evaluate their role in the complaint classification task.
\subsubsection{Punctuation}With the desire to exploit the importance of punctuation marks in Vietnamese complaints, we proceed to include in the dataset model with punctuation and without punctuation to assess. The results are shown in Table \ref{punc_result}, we notice that removing the punctuation marks has affected the performance of the model, significantly reduced results. This experiment proves that punctuation marks contribute to a better classification model.
\begin{table}[!http]
\centering
\caption{The experiment results about punctuation.}
\label{punc_result}
\begin{tabular}{lrr}
\hline
\multicolumn{1}{c}{\multirow{2}{*}{\textbf{Model}}} & \multicolumn{2}{c}{\textbf{F1-score}}                                                     \\ \cline{2-3}
\multicolumn{1}{c}{}                                & \multicolumn{1}{c}{\textbf{Punctation}} & \multicolumn{1}{c}{\textbf{Without punctation}} \\ \hline
Logistic Regression                                 & 89.62                                   & 89.17                                           \\
MLP                                                 & 90.97                                   & 89.55                                           \\
BiLSTM                                              & 89.55                                   & 88.57                                           \\
BiLSTM + PhoW2V\_syllables                          & 86.36                                   & 86.36                                           \\
BiLSTM + PhoW2V\_words                              & 88.46                                   & 88.06                                           \\
BiLSTM + fastText                                   & 86.33                                   & 85.13                                           \\
\textbf{PhoBERT}                                    & \textbf{92.16}                          & \textbf{91.47}                                  \\ \hline
\end{tabular}
\end{table}

Punctuation is the grammar medium used in writing. It shows the sentence's intonation, expressing nuances about the sentence's meaning and the writer's feelings and attitudes. In specific, question marks differ significantly between complaints and non-complaints comments and mostly appear in complaining reviews. Customers often express their inquiries, the emphasis on the content by using question marks. They desire businesses to explain why the actual product has not been as expected. For example, the review “hàng chính hãng mà bốc hơi quà km là khăn giấy tẩy trang mà không có một lời giải thích hay báo trước là sao?” (“The official item for which the promotional gift is lost is the makeup remover without an explanation, why?”). Besides, users also often use two dots, bullet marks in their complaints. They list the dissatisfaction, express comments with the desire for businesses to improve products and services. Moreover, punctuation is combined to form emotional icons. For example, the most common icons in complaining reviews are ": (" or ": ((" expressing a sad mood, their frustration. Therefore, the models with dataset retaining punctuation have better performance.
\subsubsection{Part-of-speech (POS)}
In this experiment, we use Logistic Regression algorithm with input POS features from the UIT-ViOCD dataset. Then, we conduct repeated experiments eliminating one of the part-of-speech (noun, verb, adjective, and others) from reviews to test the classification performance versus the result of the experiment with input full of POS features. As result shown in Table \ref{pos}, the performance of model for classifying decreased slightly by 1-4\% with F1-score. This result proves that each part-of-speech is a significant feature for predicting complaining reviews.

\begin{table}[H]
\centering
\caption{The result of removing part-of-speech (POS) on the dataset (\%).}
\label{pos}
\begin{tabular}{lrrrr}
\cline{1-3}
 & \multicolumn{1}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}F1-score \end{tabular}}} & \multicolumn{1}{c}{\textbf{\begin{tabular}[c]{@{}c@{}}Accuracy\end{tabular}}} & \multicolumn{1}{c}{\textbf{}} & \multicolumn{1}{c}{\textbf{}} \\ \cline{1-3}
\textbf{Normal} & \textbf{70.09} & 71.83 &  &  \\ \cline{1-3}
\textbf{Removing Verb} & \textbf{67.38} & 69.28 &  &  \\
\textbf{Removing Adjective} & \textbf{66.41} & 68.64 &  &  \\
\textbf{Removing Noun} & \textbf{69.43} & 71.10 &  &  \\
\textbf{Removing Others} & \textbf{69.25} & 71.01 &  &  \\ \cline{1-3}
 & \multicolumn{1}{l}{\textbf{}} &  &  &  \\
 & \multicolumn{1}{l}{\textbf{}} &  &  &  \\
 & \multicolumn{1}{l}{\textbf{}} &  &  &  \\
 & \multicolumn{1}{l}{\textbf{}} &  &  &
\end{tabular}
\end{table}
\section{Conclusion and Future work}
In this paper, we created the UIT-ViOCD dataset about customer complaints in Vietnamese on open-domain. After building the dataset with 5,485 human-annotated comments on four domains, we implemented various methods for the first evaluations on the dataset.  As a result, we achieved the highest performance by the fine-tuned PhoBERT classifier, a Vietnamese transfer learning model, for identifying complaints with the F1-score of 92.16\%.

In future, we are going to extend the dataset with new comments on different domains besides the current ones. From that results, we plan to create a tool for identifying complaints automatically. Moreover, the UIT-ViOCD dataset can be extended into many other tasks, such as aspect-based complaint analysis (inspired by the task of aspect-based sentiment analysis \cite{phan2021sa2sl}) or complaint span detection (inspired by the task of toxic span detection \cite{pavlopoulos2021semeval}). Social listening systems can integrate a complaint classification function.


\section{Introduction}
\label{sec:intro}

Finite mixture models have seen increasing use over the last decades \citep[for a recent survey, see][]{mcnicholas16b}.
Because of their flexibility, they are a suitable statistical tool for modeling a wide range of phenomena characterized by unobserved heterogeneity, and constitute a powerful device for clustering and classification.
Specifically, in their ``direct application'', each mixture component represents a group (or cluster) within data, and the scope is to identify these groups and estimate the parameters of the conditional-group distributions \citep{mclachlan2000finite, mcnicholas2016mixture}.
If no exogenous variables explain the means and the variances of each component, they are also called unconditional mixture models.
However, when there is a linear relationship between some variables, important insight can be gained by accounting for functional dependencies between them.
For this reason, finite mixtures of regression models with fixed covariates (FMR) have been proposed in the literature \citep[see][for examples]{desarbo1988, fruhwirth2006finite}.
Finite mixtures of regression models with concomitant variables (FMRC; \citealp{dayton1988concomitant}) are an extension of FMR where the mixing weights depend on some concomitant variables (which are often the same covariates) and are usually modeled by a multinomial logistic model (see \citealp{ingrassia2016decision, ingrassia2019cluster} and \citealp{Mazz:Batt:Ingr:Punz:Mode:2019} for details).
Unfortunately, these methodologies do not explicitly use the distribution of the covariates for clustering, i.e.~the assignment of data points to clusters does not directly utilize any information from the distribution of the covariates.
Differently from these approaches, finite mixtures of regressions with random covariates (\citealp{gershenfeld1997, gershenfeld1999}), also known as cluster-weighted models (CWMs), allow for such functional dependency.
This occurs because, for each mixture component, CWMs decompose the joint distribution of responses and covariates into the product between the marginal distribution of the covariates and the conditional distribution of the responses given the covariates.

Several CWMs have been introduced in the univariate and multivariate literature.
Most of them consider a univariate response variable, along with a set of covariates, modeled by a univariate and a multivariate distribution, respectively \citep[see][for examples]{ingrassia2012local, ingrassia2014model, punzo2014flexible}.
Instead, fewer CWMs exist in the case of a multivariate response \citep{punzo2017robust, dang2017multivariate}.
However, over the years, there has been an increasing interest in applications involving matrix-variate (three-way) data.
This data structure can occur in several and different application domains, such as longitudinal data on multiple response variables, spatial multivariate data, multivariate repeated measures or spatio-temporal data.
Nevertheless, there exists a limited number of contributions involving matrix-variate regression models.
First introduced by \citet{viroli2012matrix}, these models are a natural generalization of the multivariate-multiple regression (see also \citealp{anderlucci2014matrix}).
In the mixture framework, finite mixtures of matrix-variate regressions (MN-FMR) have been recently considered by \citet{melnykov2019studying}.
There are no matrix-variate CWMs in the literature and this paper aims to fill this gap by introducing and discussing a matrix-variate CWM in which the cluster-specific marginal distribution of the covariates, and the cluster-specific conditional distribution of the responses given the covariates, are assumed to be matrix normal.

The remainder of the paper is organized as follows.
In Section~\ref{sec:back}, some preliminary aspects are described.
The matrix normal cluster-weighted model (MN-CWM) and the expectation-conditional maximization (ECM) algorithm for parameter estimation are discussed in Section~\ref{sec:model}.
Computational and operational aspects are laid out in Section~\ref{sec:comp}.
In the simulation study outlined in Section~\ref{sec:artif}, the parameter recovery and the classification performance of the MN-CWM are investigated as well as the capability of the Bayesian information criterion (BIC; \citealp{schwarz1978estimating}) to detect the underlying group structure.
The MN-CWM is also therein compared to the MN-FMR and to the multivariate-multiple normal CWM (MMN-CWM).
The application of the MN-CWM to two real datasets concerning educational indicators and the Italian non-life insurance market is therefore analyzed in Section~\ref{sec:data}, whereas some conclusions and ideas for future developments are drawn in Section~\ref{sec:conclude}.

\section{Background}
\label{sec:back}

\subsection{Matrix normal distribution}
\label{subsec:MXVN}

A $p \times r$ continuous random matrix $\bY \in \real^{p \times r}$ has a matrix normal (MN) distribution, denoted by $\pazocal{N}_{p \times r} \left(\bM,\bPhi,\bPsi\right)$, if its density is
\begin{equation}
\resizebox{.94\textwidth}{!}{
\bphi_{p \times r}\left(\bY;\bM,\bPhi,\bPsi\right)=\left(2\pi\right)^{-\frac{pr}{2}}\left|\bPhi\right|^{-\frac{r}{2}}\left|\bPsi\right|^{-\frac{p}{2}}\exp\left\{-\frac{1}{2}\tr\left[\bPhi^{-1}\left(\bY-\bM\right)\bPsi^{-1}\left(\bY-\bM\right)^\top\right]\right\},
$}
\label{eq:MatNorm}
\end{equation}
where $\bM$ is the $p\times r$ mean matrix, and $\bPhi$ and $\bPsi$ are the $p \times p$ and $r \times r$ covariance matrices associated with the $p$ rows and $r$ columns, respectively.
An equivalent definition specifies the $\left(p\times r\right)$-matrix normal distribution as a special case of the $pr$-multivariate normal distribution.
Specifically,
\begin{equation}
\bY \sim \pazocal{N}_{p \times r} \left(\bM,\bPhi,\bPsi\right)\Leftrightarrow \vect\left(\bY\right) \sim \pazocal{N}_{pr}\left(\vect\left(\bM\right),\bPsi \otimes \bPhi \right),
\label{eq:vectorize}
\end{equation}
where $\pazocal{N}_{pr}\left(\cdot\right)$ denotes the $pr$-variate normal distribution, $\vect\left(\cdot\right)$ is the vectorization operator and $\otimes$ denotes the Kronecker product.
However, a MN distribution has the desirable feature of simultaneously modeling and identify the between and within-variables variabilities as well as reducing the number of free parameters from $pr\left(pr+1\right)/2$ to $\left[r\left(r+1\right)/2\right]+\left[p\left(p+1\right)/2\right]-1$.

\subsection{The matrix-variate regression model}
\label{subsec:MatReg}

Let $\bY \in \real^{p \times r}$ be a continuous random matrix of dimension $p \times r$, containing $p$ responses measured in $r$ occasions.
Suppose we observe a set of $q$ covariates for each occasion, inserted in a matrix $\bX$ of dimension $q \times r$.
A generic matrix-variate regression model for $\bY$ has the form
\begin{equation}
\bY = \bbeta\bw^{\top} + \bB \bX + \bU ,
\label{eq:regModel1}
\end{equation}
where $\bbeta$ is the $p \times 1 $ vector consisting in the parameters related with the intercept, $\bw$ is a $r \times 1$ vector of ones, $\bB$ is the $p \times q$ matrix containing the parameters related to the $q$ covariates and $\bU$
is the $p \times r$ error term matrix.
Model~\eqref{eq:regModel1} can be expressed in compact notation as
\begin{equation}
\bY = \bB^{*} \bX^{*} + \bU ,
\label{eq:regModel2}
\end{equation}
where $\bB^{*}$ is the $p \times \left(1+q\right)$ matrix involving all the parameters to be estimated and $\bX^{*}$ is the $\left(1+q\right) \times r$ matrix containing the information about the intercept and $q$ covariates \citep{viroli2012matrix}.
If we assume $\bU \sim \pazocal{N}_{p \times r}\left(\b0,\bPhi,\bPsi\right)$, then $\bY|\bX^{*} \sim \pazocal{N}_{p \times r}\left(\bB^{*} \bX^{*},\bPhi,\bPsi\right)$.
Therefore, a matrix-variate regression can be viewed as an encompassing framework containing as special cases the multivariate-multiple regression, when $r = 1$, and the univariate-multiple regression when $r = 1$ and $p = 1$.

\section{Methodology}
\label{sec:model}

\subsection{The matrix normal CWM}
\label{subsec:MatCWM}

Let $\left(\bX,\bY\right)$ be a pair of random matrices, defined as in Section~\ref{subsec:MatReg}, with joint distribution $p\left(\bX,\bY\right)$.
Then, a general matrix CWM has the following joint distribution
\begin{equation}
p\left(\bX,\bY\right) = \sum\limits_{g=1}^{G} p_{g}\left(\bY|\bX^{*}\right) p_{g}\left(\bX\right) \pi_{g},
\label{eq:CWM}
\end{equation}
where $p_{g}\left(\bY|\bX^{*}\right)$ is the cluster-specific conditional distribution of the responses, $p_{g}\left(\bX\right)$ is the cluster-specific marginal distribution of the covariates, and $\pi_{g}$ is the mixing weight (with $\pi_{g} > 0 $ and $\sum\limits_{g=1}^{G} \pi_{g} =1$).
Furthermore, we assume that in each group the conditional expectation $\E\left(\bY|\bX^{*}\right)$ is a linear function of $\bX^{*}$ depending on some parameters.

In this paper, we focus on model~\eqref{eq:CWM} by assuming that both $p_{g}\left(\bY|\bX^{*}\right)$ and $p_{g}\left(\bX\right)$ are matrix normal densities, and $\E\left(\bY|\bX^{*}\right) = \bB^{*} \bX^{*}$, as described in Section~\ref{subsec:MatReg}.
Thus, model~\eqref{eq:CWM} can be rewritten as
\begin{equation}
\resizebox{.94\textwidth}{!}{
p\left(\bX,\bY;\bTheta\right) = \sum\limits_{g=1}^{G} \bphi_{p \times r}\left(\bY|\bX^{*};\bB_{g}^{*}\;\bX^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX;\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g},
$}
\label{eq:NOR_CWM}
\end{equation}
where $\bTheta$ denotes the set of all parameters.
For ease of understanding, a subscript with the variable name is added to the respective covariance matrices.
Notice that there is an identifiability issue concerning the covariance matrices.
Indeed, for a MN distribution, $\bPsi \otimes \bPhi = \bPsi^{*} \otimes \bPhi^{*}$ if $\bPhi^{*}= a\bPhi$ and $\bPsi^{*}= a^{-1}\bPsi$.
As a result, matrices $\bPhi$ and $\bPsi$ are identifiable up to a multiplicative constant $a$ \citep{melnykov2019studying}.
According to \citet{gallaugher17,gallaugher2018finite}, a way to obtain a unique solution is to fix the first diagonal element of the row covariance matrix to 1.
Therefore, in model~\eqref{eq:NOR_CWM} we adopt this approach by setting the first diagonal element of $\bPhi_{\bY_{g}}$ and $\bPhi_{\bX_{g}}$ to $1$.

If the MN-CWM was not available, a possible approach would be to vectorize the matrices and consider the MMN-CWM, of which it is a special case; see Equation~\eqref{eq:vectorize}.
However, such a procedure leads to two main concerns.
The first one is the overparameterization of the vectorized model.
Secondly, this increased number of parameters may have direct effects on model selection, as will be better analyzed in Section~\ref{subsec:simu4}.

\subsection{Parameter estimation}
\label{sec:EM}

Parameter estimation is carried out via the expectation-conditional maximization (ECM) algorithm \citep{meng1997algorithm}.
The only difference with respect to the expectation-maximization (EM) algorithm \citep{dempster1977maximum} is that the M-step is replaced by a sequence of simpler and computationally convenient CM-steps.
The EM algorithm cannot be directly implemented since there is no closed form solution for the covariance matrices of the MN distribution, i.e.~one of the two depends on the value of the other at the previous iteration \citep{dutilleul1999mle}.

Let $\Sa = \left\{\left(\bX_{i},\bY_{i}\right)\right\}_{i=1}^{N}$ be a sample of $N$ independent observations from model~\eqref{eq:NOR_CWM}.
Then, the incomplete-data likelihood function is
\begin{equation}
\resizebox{.94\textwidth}{!}{
  \begin{aligned}
	L\left(\bTheta|\Sa\right) = & \prod_{i=1}^{N} p\left(\bX_{i},\bY_{i};\bTheta\right)\\
                            = & \prod_{i=1}^{N}\left[\sum\limits_{g=1}^{G} \bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*};\bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g}\right].
\end{aligned}
}
\label{eq:lik}
\end{equation}
Within the formulation of mixture models, $\Sa$ is viewed as being incomplete because, for each observation, we do not know its component membership.
Let $\bz_{i} = \left( z_{i1},\ldots,z_{iG}\right)^{\top}$ be the component membership vector such that $z_{ig}= 1$ if $\left(\bX_{i},\bY_{i}\right)$ comes from group $g$ and $z_{ig}= 0$ otherwise.
Now, the complete-data are $\Sa_{c} = \left\{\left(\bX_{i},\bY_{i},\bz_{i}\right)\right\}_{i=1}^{N}$, and the complete-data likelihood is
\begin{equation}
\resizebox{.94\textwidth}{!}{
\begin{aligned}
L_{c}\left(\bTheta|\Sa_{c}\right) = \prod_{i=1}^{N} \prod_{g=1}^{G} \left[\bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*}; \bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right) \bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right) \pi_{g}\right]^{z_{ig}}.
\end{aligned}
}
\label{eq:likcomp}
\end{equation}
Therefore, the corresponding complete-data log-likelihood can be written as
\begin{align}
l_{c}\left(\bTheta|\Sa_{c}\right) = & \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left(\pi_{g}\right) + \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left[\bphi_{p \times r}\left(\bY_{i}|\bX_{i}^{*};\bB_{g}^{*}\;\bX_{i}^{*},\bPhi_{\bY_{g}},\bPsi_{\bY_{g}}\right)\right]\nonumber\\
& + \sum\limits_{i=1}^{N} \sum\limits_{g=1}^{G} z_{ig} \ln\left[\bphi_{q \times r}\left(\bX_{i};\bM_{g},\bPhi_{\bX_{g}},\bPsi_{\bX_{g}}\right)\right].\label{eq:llk}
\end{align
In the following, by adopting the notation used in \citet{tomarchio2020two}, the parameters marked with one dot correspond to the updates at the previous iteration and those marked with
two dots represent the updates at the current iteration.

\paragraph{E-Step}
\label{subsec:E-step}

The E-step requires calculation of the conditional expectation of \eqref{eq:llk}, given the observed data and the current estimate of the parameters $\dot{\bTheta}$.
To do this, we need to calculate
\begin{align}
\ddot{z}_{ig} & = \E_{\dot{\bTheta}} \left[Z_{ig}|\bX_{i},\bY_{i}\right]\nonumber\\
                          &= \frac{\dot{\pi}_{g} \bphi_{p \times r}(\bY_{i}|\bX_{i}^{*};\dot{\bB}_{g}^{*}\bX_{i}^{*},\dot{\bPhi}_{\bY_{g}},\dot{\bPsi}_{\bY_{g}}) \bphi_{q \times r}(\bX_{i};\dot{\bM}_{g},\dot{\bPhi}_{\bX_{g}},\dot{\bPsi}_{\bX_{g}})}  {\sum\limits_{j=1}^{G} \dot{\pi}_{j} \bphi_{p \times r}(\bY_{i}|\bX_{i}^{*};\dot{\bB}_{j}^{*}\bX_{i}^{*},\dot{\bPhi}_{\bY_{j}},\dot{\bPsi}_{\bY_{j}}) \bphi_{q \times r}(\bX_{i};\dot{\bM}_{j},\dot{\bPhi}_{\bX_{j}},\dot{\bPsi}_{\bX_{j}})}, \label{eq:E-step}
\end{align}
which corresponds to the posterior probability that the unlabeled observation $\left(\bX_{i},\bY_{i}\right)$ belongs to the $g$th component of the mixture.

\paragraph{CM-Step 1}
\label{subsec:cm1}

At the first CM-step, we maximize the expectation of the complete-data log-likelihood with respect to $\bTheta_{1}=\left\{\pi_g,\bM_{g},\bPhi_{\bX_{g}},\bB_{g}^{*},\bPhi_{\bY_{g}}\right\}_{g=1}^{G}$, fixing $\bTheta_{2}=\left\{\bPsi_{\bX_{g}},\bPsi_{\bY_{g}}\right\}_{g=1}^{G}$ at $\dot\bTheta_{2}$.
In particular, we obtain
\begin{equation}
\ddot{\pi}_g   = \frac{\sum_{i=1}^N \ddot{z}_{ig}}{N},
\label{eq:cwm_cm1.1}
\end{equation}
\begin{equation}
\ddot{\bM}_g    = \frac{\sum_{i=1}^N \ddot{z}_{ig}\bX_i}{\sum_{i=1}^N \ddot{z}_{ig}},
\label{eq:cwm_cm1.2}
\end{equation}
\vspace{-8mm}
\begin{align}
\ddot{\bB}_{g}^{*} & = \left[\sum_{i=1}^N \ddot{z}_{ig} \bY_{i} (\dot{\bPsi}_{\bY_{g}})^{-1} \bX_{i}^{*\top}\right] \left[\sum_{i=1}^N \ddot{z}_{ig} \bX_{i}^{*} (\dot{\bPsi}_{\bY_{g}})^{-1} \bX_{i}^{*\top}\right]^{-1}, \label{eq:cwm_cm1.3} \\
\ddot{\bPhi}_{\bX_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bX_i-\ddot{\bM}_g\right)\left(\dot\bPsi_{\bX_{g}}\right)^{-1}\left(\bX_i-\ddot{\bM}_g\right)^{\top}}{r \sum_{i=1}^N \ddot{z}_{ig}}, \label{eq:cwm_cm1.4} \\
\ddot{\bPhi}_{\bY_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)\left(\dot\bPsi_{\bY_{g}}\right)^{-1}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)^{\top}}{r \sum_{i=1}^N \ddot{z}_{ig}}. \label{eq:cwm_cm1.5}
\end{align}

\paragraph{CM-Step 2}

At the second CM-step, we maximize the expectation of the complete-data log-likelihood with respect to $\bTheta_{2}$, keeping fixed $\bTheta_{1}$ at $\ddot\bTheta_{1}$.
Therefore, we have
\begin{align}
\ddot{\bPsi}_{\bX_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bX_i-\ddot{\bM}_g\right)^{\top}\left(\ddot\bPhi_{\bX_{g}}\right)^{-1}\left(\bX_i-\ddot{\bM}_g\right)}{q \sum_{i=1}^N \ddot{z}_{ig}},\label{eq:cwm_cm2.1} \\
\ddot{\bPsi}_{\bY_{g}} & = \frac{\sum_{i=1}^N \ddot{z}_{ig}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)^{\top}\left(\ddot\bPhi_{\bY_{g}}\right)^{-1}\left(\bY_i-\ddot{\bB}^{*}_{g}\bX_{i}^{*}\right)}{p \sum_{i=1}^N \ddot{z}_{ig}}.\label{eq:cwm_cm2.2}
\end{align}

\section{Computational and operational aspects}
\label{sec:comp}

The code for the ECM algorithm previously described, and for the analysis in the following sections, is written within the R computing environment \citep{R_software}.

\subsection{ECM initialization}
\label{subsec:init}

The choice of the starting values is an important aspect for EM-based algorithms (see, e.g., \citealp{maitra2010simulating, michael2016effective}).
A common way to start an EM-based algorithm consists in providing initial values for the $z_{ig}$ in~\eqref{eq:E-step}, during the first E-step of the algorithm \citep{mclachlan2000finite}.
In our case, we need to provide an initial value also for $\bPsi_{\bX_{g}}$ and $\bPsi_{\bY_{g}}$ in~\eqref{eq:cwm_cm1.4} and~\eqref{eq:cwm_cm1.5}, respectively, during the first CM-step 1 of the algorithm.
Therefore, the following initialization strategy is implemented:
\begin{enumerate}
\item Generate $G$ random positive definite matrices for both $\bPsi_{\bX_{g}}$ and $\bPsi_{\bY_{g}}$.
This is done via the \texttt{genPositiveDefMat()} function of the \textbf{clusterGeneration} package, by using the ``eigen'' method. For further details see \citet{clusterG}.
\item generate $N$ random vectors $z_{i}=\left(z_{i1},\ldots,z_{iG}\right)^{\top}$, $i=1,\ldots,N$.
This is done by using the following three approaches:
\begin{enumerate}[label*=\arabic*.]
    \item in a ``soft'' way, by generating $G$ positive random values from a uniform distribution on [0,1] for each observation, that are subsequently normalized to have a unitary sum.
		Being purely random, this procedure is repeated 15 times, and the solution maximizing the observed-data log-likelihood among these runs is considered;
		\item in a ``hard'' way, by using the classification produced by the $k$-means algorithm on the vectorized and merged data.
		Specifically, after computing $\left\{\vect\left(\bX_{i}\right)\right\}_{i=1}^{N}$ and $\left\{\vect\left(\bY_{i}\right)\right\}_{i=1}^{N}$, the data are merged so that for each
		observation we have a vector of dimension $\left(pr+qr\right) \times 1$;
		\item in a ``hard'' way, by using the classification produced by mixtures of matrix-normal distributions, computed on the merged data.
		In detail, for each observation we have a $\left(p+q\right) \times r$ matrix.
  \end{enumerate}
\end{enumerate}
The approach producing the highest observed-data log-likelihood is finally selected.

\subsection{Spurious clusters}
\label{subsec:spur}

A well-known issue in mixture models is related to the possibility for EM-based algorithms to converge to spurious local maximizers.
This is not a failing of these algorithms, and such solutions reflect a random pattern in the data rather than an underlying group structure.
They can be usually detected by the presence of data groups with very few observations or small eigenvalues compared to the other ones.
Consequently, these solutions often have a high likelihood, but are of little practical use or real-world interpretation \citep{mclachlan2000finite}.
Various approaches for mitigate this problem have been proposed in the literature.
For example, \citet{leisch2004flexmix} and \citet{dang2015families} impose a minimum size for the clusters of $\pi_{g}=0.05$.
\citet{melnykov2019studying} fought with this issue by excluding all the solutions that involved clusters consisting of less than 5 points.
We faced this problem by removing the solutions that included clusters with estimated ${\pi}_{g}\approx 0.05$ or less, and with eigenvalues close to zero, in one or more of its estimated covariance matrices.

\subsection{Model selection and performance evaluation}
\label{subsec:classperf}

It is often the case that the number of groups $G$ is not known in advance, and model selection criteria are commonly used for estimating it.
Among them, the BIC is one of the most popular, and will be used in the following.
Furthermore, when the true classification of the data is known, the adjusted Rand index (ARI; \citealp{hubert1985comparing}) can be adopted to evaluate the clustering performance of a model.
Specifically, the ARI calculates the agreement between the true classification and the one predicted by the model.
An ARI of 1 indicates perfect agreement between the two partitions, whereas the expected value of the ARI under a random classification is 0.
It will be used in the manuscript along with the misclassification percentage rate $\eta$, which is the percentage of units classified in the wrong classes.

\section{Simulation studies}
\label{sec:artif}

\subsection{Simulation 1: A focus on the matrix-normal CWM}
\label{subsec:simu1}

In this study, several aspects related to our model are analyzed.
First of all, since the ECM algorithm is used to fit the model, it is desirable to evaluate its parameter recovery, i.e.~whether it can recover the generating parameters accurately.
For this reason, data are generated from a four-component MN-CWM with $p=q=r=3$.
Two scenarios are then evaluated, according to the different level of overlap of the mixture components.
In the first scenario (labeled as ``Scenario A$_{1}$''), the mixture components are well-separated both in $\bX$, by assuming relatively distant mean matrices, and in $\bY|\bX^{*}$, by using different intercepts and slopes.
On the contrary, in the second scenario (labeled as ``Scenario B$_{1}$''), there is a certain amount of overlap, since the intercepts are all equal among the mixture components, while the slopes and the mean matrices assume approximately the same values among the mixture components.
The parameters used for Scenario A$_{1}$ are displayed in~


\section{Introduction}
\label{sec:intro}
The discrepancy between the ratio $\mu_p G_E / G_M$ of the the proton's electromagnetic form factors extracted from polarization
asymmetry measurements, and the ratio extracted from unpolarized cross section measurements, leaves the field of form factor physics
in an uncomfortable state (see \cite{Afanasev:2017gsk} for a recent review).
On the one hand, there is a consistent and viable hypothesis that the discrepancy is caused by non-negligible hard two-photon exchange (TPE) \cite{Guichon:2003qm,Blunden:2003sp}, the one radiative correction omitted from the standard radiative correction prescriptions \cite{Mo:1968cg,Maximon:2000hm}.
On the other hand, three recent measurements of hard TPE (at VEPP-3, at CLAS, and with OLYMPUS) found that the effect of TPE is small in the region of $Q^2 < 2$~GeV$^2/c^2$ \cite{Rachek:2014fam,Adikaram:2014ykv,Rimal:2016toz,Henderson:2016dea}.
The TPE hypothesis is still viable; it is possible that hard TPE contributes more substantially at higher momentum transfers, and can fully resolve the form factor discrepancy.
But the lack of a definitive conclusion from this recent set of measurements is an indication that alternative approaches are needed to illuminate the situation, and it may be prudent to concentrate experimental effort on constraining and validating model-dependent theoretical calculations of TPE.
There are multiple theoretical approaches, with different assumptions and different regimes of validity~\cite{Chen:2004tw,Afanasev:2005mp,Tomalak:2014sva,Blunden:2017nby,Kuraev:2007dn}.
If new experimental data could validate and solidify confidence in one or more theoretical approaches, then hard TPE could be treated in the future like any of the other standard radiative corrections, i.e., a correction that is calculated, applied, and trusted.

VEPP-3, CLAS, and OLYMPUS all looked for hard TPE through measurements of the $e^+p$ to $e^-p$ elastic scattering cross section ratio.
After applying radiative corrections, any deviation in this ratio from unity indicates a contribution from hard TPE.
However, this is not the only experimental signature one could use. Hard TPE can also appear in a number of polarization asymmetries.
Having constraints from many orthogonal directions, i.e., from both cross section ratios and various polarization asymmetries, would be valuable for testing and validating theories of hard TPE.
As with unpolarized cross sections, seeing an opposite effect for electrons and positrons is a clear signature of TPE.

Here, we propose one such polarization measurement that could both be feasibly accomplished with a positron beam at Jefferson Lab, and contribute new information about two photon exchange that could be used to constrain theoretical models.
We propose to measure the polarization transfer (PT) from a polarized positron beam scattering elastically from a proton target, for which no data currently exist.
The proposed experiment uses the Hall A Super Big-Bite Spectrometer (SBS) to measure the polarization of recoiling protons and a lead-glass calorimeter for detecting scattered positrons in coincidence.
In the following sections, we review polarization transfer, sketch the proposed measurement, and discuss possible systematic uncertainties.

\section{Polarization Transfer}
\label{sec:pt}
In the Born approximation (i.e.~one-photon exchange), the polarization transferred from a polarized lepton to the recoiling proton is
\begin{eqnarray}
    P_t &= -hP_e\sqrt{\frac{2\epsilon (1-\epsilon)}{\tau}} \frac{G_E G_M}{G_M^2 + \frac{\epsilon}{\tau}G_E^2}, \nonumber \\
    P_l &= hP_e\sqrt{1 - \epsilon^2} \frac{G_M^2}{G_M^2 + \frac{\epsilon}{\tau}G_E^2}, \label{PTeqns}
\end{eqnarray}
where $P_t$ is the polarization transverse to the momentum transfer 3-vector (in the reaction plane), $P_l$ is the longitudinal polarization, $P_e$ is the initial lepton polarization, $h$ is the lepton helicity, $\tau \equiv \frac{Q^2}{4M^2}$ is the dimensionless 4-momentum transfer squared, $\epsilon \equiv \left[1+2(1+\tau)\tan^2 \left(\frac{\theta_e}{2}\right)\right]^{-1}$, with $\theta_e$ the electron scattering angle in the nucleon rest (laboratory) frame, is the virtual photon polarization parameter, and $G_E$ and $G_M$ are the proton's electromagnetic form factors.
The strength of the polarization transfer technique is to measure $P_t/P_l$, thereby cancelling some systematics associated with polarimetry, and isolating the ratio of the proton's form factors:
\begin{equation}
\frac{P_t}{P_l}= -\sqrt{\frac{2\epsilon}{\tau(1+\epsilon)}}\frac{G_E}{G_M}. \label{ratio}
\end{equation}
This technique has several advantages over the traditional Rosenbluth separation technique for determining form factors.
This polarization ratio can be measured at a single kinematic setting, avoiding the systematics associated with comparing data taken from different spectrometer settings.
This technique allows the relative sign of the form factors to be determined, rather than simply their magnitudes.
And furthermore, whereas the sensitivity of Rosenbluth separation to $G_E^2$ diminishes at large momentum transfer, polarization transfer retains sensitivity to $G_E$ even when $Q^2$ becomes large.
When used in combination at high $Q^2$, Rosenbluth separation can determine $G_M^2$, while polarization transfer can determine $G_E/G_M$, allowing the form factors to be separately determined.

Polarization transfer using electron scattering has been used extensively to map out the proton's form factor ratio over a wide range of $Q^2$, with experiments conducted at MIT Bates~\cite{Milbrath:1997de}, Mainz~\cite{Pospischil:2001pp}, and Jefferson Lab~\cite{Gayou:2001qt,MacLachlan:2006vw,Ron:2011rd,Paolone:2010qc,Zhan:2011ji}, including three experiments, GEp-I~\cite{Jones:1999rz,Punjabi:2005wq}, GEp-II~\cite{Gayou:2001qd,Puckett:2011xg}, and GEp-III~\cite{Puckett:2010ac,Puckett:2017flj} that pushed to high momentum transfer.
Another experiment, GEp-2$\gamma$, looked for hints of TPE in the $\epsilon$-dependence in polarization transfer \cite{Meziane:2010xc,Puckett:2011xg}.
Two other experiments made equivalent measurements by polarizing the proton target instead of measuring recoil polarization \cite{Jones:2006kf,Crawford:2006rz}.

While polarization transfer is less sensitive to the effects of hard TPE, it is not immune. Following the formalism of Ref.~\cite{Carlson:2007sp}, one finds that
\begin{eqnarray}
\frac{P_t}{P_l} &=& -\sqrt{\frac{2\epsilon}{\tau (1+\epsilon)}}\frac{G_E}{G_M} \times \Bigg[ 1 \pm
                    \textrm{Re}\left(\frac{\delta \widetilde{G}_M}{G_M}\right)  \nonumber \\
                & & \pm \frac{1}{G_E}\textrm{Re}\left(\delta\widetilde{G}_E + \frac{\nu}{M^2}\widetilde{F}_3\right) \nonumber \\
                & & \mp \frac{2}{G_M}\textrm{Re}\left(\delta\widetilde{G}_M + \frac{\epsilon\nu}{(1+\epsilon)M^2} \widetilde{F}_3\right) + \mathcal{O}(\alpha^2) \Bigg], \label{eq3}
\end{eqnarray}
with $\nu \equiv (p_e + p_{e'})_\mu (p_p + p_{p'})^\mu$, and where $\delta \widetilde{G}_E$,  $\delta \widetilde{G}_M$, and $\delta \widetilde{F}_3$ are additional form factors that become non-zero when moving beyond the one-photon exchange approximation and, crucially, depend on both $Q^2$ and $\epsilon$, whereas the one-photon-exchange form factors depend only on $Q^2$. The correction terms $\delta \widetilde{G}_E$, $\delta\widetilde{G}_M$, and $\widetilde{F}_3$ are $\mathcal{O}(\alpha)$ relative to the one-photon-exchange form factors $G_E, G_M$. The $\pm$/$\mp$ symbols in Eq.~\eqref{eq3} indicate the sign with which the two-photon-exchange amplitudes enter the observable $P_t/P_l$ depending on the lepton charge, with the upper (lower) symbol indicating the appropriate sign for $e^-(e^+)$ beams.
This particular dependence on new form factors is slightly different than what one finds when taking a positron to electron cross section ratio:
\begin{eqnarray}
\frac{\sigma_{e^+p}}{\sigma_{e^-p}} &=& 1 + 4G_M\textrm{Re}\left(\delta \widetilde{G}_M + \frac{\epsilon \nu}{M^2}\widetilde{F}_3\right) \nonumber \\
  & & - \frac{4\epsilon}{\tau}G_E\textrm{Re}\left(\delta \widetilde{G}_E + \frac{\nu}{M^2} \widetilde{F}_3\right) + \mathcal{O}(\alpha^2).
\end{eqnarray}
A measurement of the difference in polarization transfer between electron and positron scattering therefore adds information about TPE in addition to what can be learned from cross section ratios alone.

The GEp-2$\gamma$ experiment looked for the effects of TPE in polarization transfer by making measurements at three kinematic points with varying values of $\epsilon$, but with $Q^2$ fixed at 2.5~GeV$^2/c^2$~\cite{Meziane:2010xc}.
Since in the absence of hard TPE the ratio $G_E/G_M$ has no $\epsilon$-dependence, any variation with $\epsilon$ is a sign of hard TPE.
The GEp-2$\gamma$ measurement was statistically consistent with no $\epsilon$-dependence, though its measurement of the relative variation with $\epsilon$ of the longitudinal component $P_\ell$ showed deviations from the one-photon exchange expectation at the level of 1.4\%, with a statistical significance of roughly 2$\sigma$.

A measurement with positron scattering for several $Q^2$ values where the discrepancy between cross section and polarization data is large, and where the $e^-p$ polarization transfer observables have already been measured precisely, will be useful for constraining TPE effects, because deviations from the Born approximation should have the opposite sign from those in electron scattering. This will help determine if deviations are truly caused by TPE, or if they arise from systematic effects.
As the largest systematic uncertainties in polarization transfer are associated with proton polarimetry, a measurement with positrons would have largely the same systematics as an experiment with electrons.

\section{Proposed Measurement}

We propose a first set of measurements of $e^+p$ polarization transfer observables in the $Q^2$ region where the Rosenbluth-polarization discrepancy is large, using the newly constructed Super BigBite Spectrometer (SBS), which was designed to measure $G_E^p/G_M^p$ to $Q^2 \approx 12\text{ GeV}^2$ using the polarization transfer technique. Details of the planned high-$Q^2$ measurements of $\vec{e}^- p$ polarization transfer, hereafter referred to as the SBS GEP experiment, can be found in the original and jeopardy proposals to Jefferson Lab's Program Advisory Committee \cite{GEP5,GEP5_PAC47}. Despite the lower expected figure-of-merit $P^2 I$ of polarized positron beams compared to polarized electron beams, these measurements can be accomplished in a reasonable time frame owing to the large solid-angle acceptance of the new SBS apparatus.
\begin{figure*}[ht]
  \begin{center}
    \includegraphics[width=0.75\textwidth]{GEPscreenshot_June24_2020_labeled.png}

  \end{center}
  \caption{\label{fig:SBS} Screenshot from the GEANT4-based Monte Carlo simulation of the SBS-GEP apparatus, illustrating one elastic $e^+ p$ event generated within the 40-cm liquid hydrogen target, with the electron detected in the lead-glass calorimeter (located on beam left) and the outgoing polarized proton detected in the SBS on beam right.}
\end{figure*}
Figure~\ref{fig:SBS} shows the layout of the proposed experiment in \textit{g4sbs}, the SBS GEANT4-based Monte Carlo simulation package. Polarized positrons are elastically scattered from free protons at rest in a 40-cm liquid hydrogen target. Scattered positrons are detected in a lead-glass calorimeter (ECAL) and a ``coordinate detector'' (CDET), consisting of two planes of scintillator strips with high segmentation in the vertical direction. The combination of CDET and ECAL provides a highly efficient and selective trigger for elastically scattered positrons and precise measurement of the positron's scattering angles, for a clean selection of the elastic $e^+p$ channel in the presence of higher-rate inelastic background processes, predominantly $\pi^0$ photoproduction.

Elastically scattered protons are detected in the SBS, which consists of a large dipole magnet with a transverse field integral along the direction of particle motion of up to 2.5 T$\cdot$m, a proton polarimeter with Gas Electron Multiplier (GEM)-based tracking and CH$_2$ as analyzer material, and a large hadron calorimeter. The role of the dipole magnet is for momentum analysis and to precess the longitudinal polarization of the recoiling proton into a transverse component that can be measured by the secondary analyzing scattering in the CH$_2$. The tracking in SBS relies on the relatively recently invented technology of Gas Electron Multipliers (GEMs) ~\cite{Gnanvo:2014hpa}, which can operate with stable gain at very high charged particle fluxes. The SBS front tracker, made of six GEM layers of area $40 \times 150 \text{ cm}^2$, is used for reconstruction of the proton's momentum, scattering angles, and interaction vertex, and also to define the proton's incident trajectory on the polarimeter, for subsequent measurement of the angular distribution of the secondary scattering. The spin-orbit coupling in the $\vec{p} + \text{CH}_2 \rightarrow \vec{p} + X$ scattering gives rise to an azimuthal asymmetry in the distribution of scattered protons that is proportional to their initial transverse polarization. Each of the two CH$_2$ analyzer blocks has a thickness of approximately one nuclear interaction length, and is followed by a tracker assembled from five GEM layers of area $60 \times 200\text{ cm}^2$, to measure the angular distribution of the polarization-analyzing scattering. Finally, a large iron-scintillator sampling hadronic calorimeter (HCAL), absorbs the energy of the protons and provides for efficient triggering on the events of interest, which are those in which the proton undergoes forward-angle elastic scattering in either (or both) of the two analyzers~\cite{Basilev:2019sno}.

\begin{table*}
  \caption{\label{tab:prop} Summary of proposed measurements. $E_{e}$ is the incident lepton energy, $\left<Q^2\right>$ is the acceptance averaged $Q^2$, $\theta_{e}$ is the central lepton scattering angle, $\left<\epsilon\right>$ is the acceptance averaged $\epsilon$ value, $\theta_p$ is the central proton scattering angle, and $p_p$ is the central proton momentum. The expected event rate is based on the assumption of a 200 nA (30 $\mu$A) positron (electron) beam, and $\Delta R$ is the projected absolute statistical uncertainty for the indicated number of beam days in the ratio $R \equiv -\mu_p \frac{P_t}{P_\ell} \sqrt{\frac{\tau(1+\epsilon)}{2\epsilon}}$, which equals $\mu_p G_E^p/G_M^p$ in the one-photon approximation, assuming 60\% (85\%) positron (electron) polarization. On the third line, we depict an ancillary $e^- p$ measurement at kinematics identical to the higher $Q^2\ e^+p$ measurement, that could achieve 1\% statistical precision in 24 hours (not including any time required to change CEBAF from $e^+$ to $e^-$ running).}
  \begin{center}
    \begin{tabular}{cccccccccc}
    \hline\hline\noalign{\smallskip}
      Lepton & $E_{e}$ & $\left<Q^2\right>$ & $\theta_{e}$ & $\left<\epsilon\right>$ & $\theta_p$ & $p_p$ & Event rate & Days & $\Delta R$  \\
      & GeV & GeV$^2$ & deg. & & deg. & GeV & Hz & & (absolute) \\
      \noalign{\smallskip}\hline\noalign{\smallskip}
      $e^+$ & 2.2 & 2.5 & 69.8 & 0.39 & 23.2 & 2.04 & 11 & 30 & 0.015 \\
      $e^{+}$ & 4.4  & 2.6 & 27.0 & 0.84 & 36.2 & 2.15 & 16 & 30 & 0.021 \\
      $e^{+}$ & 4.4  & 3.4 & 32.5 & 0.76 & 31.1 & 2.56 & 7 & 60 & 0.023 \\ \noalign{\smallskip}\hline
      $e^{-}$ & 4.4 & 3.4 & 32.5 & 0.76 & 31.1 & 2.56 & 1,050 & 1 & 0.01 \\ \noalign{\smallskip}\hline\hline
    \end{tabular}
  \end{center}
\end{table*}

Optimizing a measurement of $R_{e^+p} \equiv -\mu_p \frac{P_t}{P_\ell}\sqrt{\frac{\tau(1+\epsilon)}{2\epsilon}}$ (which equals $\mu_p G_E^p/G_M^p$ in the one-photon approximation) requires a choice of $Q^2$ and beam energy that maximizes the product of the asymmetry magnitude squared and the event rate. Merely maximizing the electron differential cross section $d\sigma/d\Omega_e$ by choosing the highest available beam energy does not always lead to the highest figure of merit (FOM) at a given $Q^2$, due to the $\epsilon$ dependence of $P_t$ and $P_\ell$, which both vanish in the limit $\epsilon \rightarrow 1$ (see equations~\eqref{PTeqns}), and also the diminishing reaction Jacobian at forward angles of the electron, where the solid angle $\Delta \Omega_e$ corresponding to the fixed proton solid angle $\Delta \Omega_p$ becomes small. The uncertainty of the ratio $R$ is typically dominated by the uncertainty of the transverse component $P_t$ of the transferred polarization, which reaches a maximum at $\epsilon \approx 0.5$, which usually occurs around $\theta_e \approx 45^\circ$. Generally speaking, for a fixed proton solid angle acceptance $\Delta \Omega_p$, the FOM at a constant $Q^2$ has a broad central maximum in the region $0.3 \le \epsilon \le 0.9$, in which it does not vary strongly.

In planning an optimized program of exploratory $R_{e^+p}$ measurements, $Q^2$ should be chosen large enough that significant TPE corrections to this observable might reasonably be expected, but small enough that useful precision can be achieved in a ``reasonable'' amount of beam time. For an initial, exploratory set of measurements, a reasonable precision goal is $\approx 2\%$ absolute statistical uncertainty in $R_{e^+p}$ at each kinematic point (or less). Such precision would be sufficient to set a useful upper limit on the size of TPE corrections in this observable and sensitive enough to detect any TPE signal as large as, e.g., the correction needed to account for the discrepancy between cross section and polarization data. Indeed, the TPE correction to the reduced cross section (at $\epsilon = 0$ relative to $\epsilon = 1$) for elastic $ep$ scattering required to explain the discrepancy with polarization data is estimated at $\approx 4\%$ for $Q^2 \gtrsim 2$ GeV$^2$, depending on the assumed $Q^2$ and/or $\epsilon$ dependence of the TPE effect~\cite{Christy:2021snt}.

It is also desirable to choose a $Q^2$ for which $R_{e^-p}$ is already precisely known. $Q^2 = 2.5$ GeV$^2$ is an obvious choice, being close to the most precise existing measurements~\cite{Puckett:2017flj} in the $Q^2$ region where the discrepancy is significant. At this $Q^2$, it is also possible to measure the $\epsilon$ dependence of $R_{e^+p}$ in a reasonable amount of beam time, as shown in Fig.~\ref{fig:SBS2gamma}. In 60 days' beam time, $R_{e^+p}$ can be measured over a significant range of $\epsilon$ with an absolute statistical uncertainty of 2\% or less. Combining the two measurements in a weighted average would give $\Delta R(e^+p) \approx 0.012$, which is competitive in precision with the best existing electron measurements at this $Q^2$. Another measurement at a meaningfully larger $Q^2 \approx 3.5$ GeV$^2$ would also be attractive, as it would be very close to two existing measurements from the GEp-I~\cite{Punjabi:2005wq} and GEp-II~\cite{Gayou:2001qd} experiments which, however, are significantly less precise. Going significantly higher in $Q^2$ than 3.5 GeV$^2$ would most likely require prohibitive beam time to reach a precision goal of 1-2\%, given the low maximum current for polarized positrons.

Table~\ref{tab:prop} shows the basic parameters of a plausible program of $R_{e^+p}$ measurements using the SBS GEP apparatus. To estimate the precision of these measurements, elastic $e^+p$ scattering events were generated in a range of angles sufficient to populate the combined acceptance of SBS and ECAL, and tracked through the GEANT4 simulation of the experiment, including the transport of the outgoing proton's spin through the SBS magnetic field. The event rates and the figure of merit for polarimetry were evaluated using the methods described in Ref.~\cite{GEP5_PAC47}, assuming a 200 nA, 60\% polarized $e^+$ beam~\cite{Abb16} on a 40-cm liquid hydrogen target. Since the precision of the existing data at the higher $Q^2$ is only about 4\% (absolute), it would be desirable to include an additional measurement of $e^-p$ scattering in identical kinematics. This could be accomplished in a tiny fraction of the total beam time, as shown in Tab.~\ref{tab:prop}, plus any time that would be required to change CEBAF from positron mode to electron mode and back again. Alternatively, such a measurement could be added to the planned SBS GEP experiment with very little additional beam time, in anticipation of a future comparison to positron measurements.
\begin{figure}
    \centering
    \includegraphics[width=0.98\columnwidth]{GEp2gamma_positron.pdf}
    \caption{Projected statistical precision of a measurement of the $\epsilon$ dependence of $R(e^+p)$ at $Q^2 = 2.5$ GeV$^2$ in a hypothetical 60-day run using the SBS. The black filled circles show the existing $e^- p$ data from the Hall C GEp-2$\gamma$ experiment~\cite{Puckett:2017flj,Meziane:2010xc}. The empty square shows the GEp-I result~\cite{Punjabi:2005wq} for $R(e^-p)$ at this $Q^2$, offset by +0.03 in $\epsilon$ for clarity. The blue triangles show the projected statistical precision of the two 30-day SBS measurements depicted in Tab.~\ref{tab:prop}, plotted at the weighted average $R$ of the Hall C data. The red square shows the projected precision of the proposed SBS $e^+p$ measurements when combined in a weighted average, plotted arbitrarily at $\epsilon = 0.95$. The pink inverted triangle shows the weighted average of the GEp-2$\gamma$ data for $e^-p$, plotted arbitrarily at $\epsilon = 0.975$.}
    \label{fig:SBS2gamma}
\end{figure}

The systematic uncertainties of the polarization transfer method are typically extremely small. Because both $P_t$ and $P_\ell$ are measured simultaneously in a single kinematic configuration, a number of sources of systematic uncertainty, such as beam polarization and analyzing power, cancel in the ratio $P_t/P_\ell$. The luminosity also doesn't need to be known precisely. Moreover, the $ep$ reaction is self-calibrating with respect to the analyzing power, and the rapid beam helicity reversal cancels the effects of false or instrumental asymmetries in the polarimeter. In previous experiments of this type, a dominant source of systematic error was the calculation of the proton spin precession in focusing magnetic spectrometers with several quadrupole magnets in addition to the main, momentum-analyzing dipole. In the SBS case, the spin precession calculation is much simpler, as the SBS is a single, simple dipole magnet which is non-focusing. It is therefore anticipated that any measurement of $e^+ p$ polarization transfer observables will be statistics-limited in terms of accuracy. In addition, the relatively low luminosity of the proposed $e^+p$ measurements means that the event reconstruction in the SBS detectors will be extremely clean, and far \emph{less} challenging than in the approved high-$Q^2$ measurements from the SBS GEP experiment~\cite{GEP5_PAC47}.

\begin{figure}[ht]
  \begin{center}

    \includegraphics[width=0.98\columnwidth]{PositronFig_new.pdf}
  \end{center}
  \caption{\label{fig:SBSeplus} Projected results of the proposed future measurements of $R_{e^+p}$ using the polarization transfer method with the SBS GEP apparatus, compared to selected existing data and the projected results of the SBS GEP experiment~\cite{GEP5_PAC47}. Selected Rosenbluth separation data are from Refs.~\cite{Andivahis:1994rq,Qattan:2004ht,Christy:2004rc}. Selected polarization transfer data are from Refs.~\cite{Punjabi:2005wq} (GEp-I), ~\cite{Puckett:2011xg} (GEp-II), and~\cite{Puckett:2017flj} (GEp-III and GEp-2$\gamma$). Projections of future experiments, including the proposed measurements from Tab.~\ref{tab:prop}, are shown at values of $R$ from the global fit described in Ref.~\cite{Puckett:2017flj}, and the projected precision of the SBS result at 2.5 GeV$^2$ is the weighted average of the two measurements at different $\epsilon$ values (see Fig.~\ref{fig:SBS2gamma} and Tab.~\ref{tab:prop}).}
\end{figure}
Figure~\ref{fig:SBSeplus} shows what could be accomplished in the 120-day experiment at ``standard'' CEBAF 1st-pass and 2nd-pass beam energies of 2.2 and 4.4 GeV, under the perhaps somewhat optimistic assumption that a positron beam of 200 nA current and 60\% polarization could be realized at CEBAF. This would be the first measurement of polarization transfer in $e^+p$ scattering, reaching very respectable precision in the $Q^2$ regime where the discrepancy between cross sections and polarization observables is large, and where $R_p$ is falling most rapidly as a function of $Q^2$. Such data would provide important model-independent constraints on hard TPEX amplitudes, toward the goal of finding a conclusive explanation of the discrepancy and a model-independent, data-constrained theoretical prescription for applying hard TPEX corrections to elastic $e^\pm p$ scattering observables.

\begin{acknowledgements}
This work was supported in part by the US Department of Energy Office of Science, Office of Nuclear Physics, Award ID DE-SC0021200 and in part by the National Science Foundation grant number 2012114
\end{acknowledgements}

\iffalse
\begin{figure}
  \includegraphics{example.eps}
\caption{Please write your figure caption here}
\label{fig:1}
\end{figure}
\begin{figure*}
  \includegraphics[width=0.75\textwidth]{example.eps}
\caption{Please write your figure caption here}
\label{fig:2}
\end{figure*}
\fi

\bibliographystyle{spphys}