The coronavirus disease 2019 (COVID-19), caused by a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has become an unprecedented public health crisis. Around the world, many governments issued a call to researchers in machine learning (ML) and artificial intelligence (AI) to address high-priority questions related to COVID-19. This call was not unusual because ML methods are finding many uses in medical diagnosis applications. The ML field is rich with examples where based on predictive models one can, for example, estimate disease severity and consequently, the state of a patient's health. These models employ data-driven algorithms that can extract features and discover complicated patterns that could have not been recognized or interpreted by humans [1] .
Much of the medical data are of multivariate categorical type, and typically they represent examinations and test results of patients. Besides, they include random errors and systematic biases, and sometimes they are missing. The development of accurate models and efficient inference algorithms for overcoming the challenges that clinical data bring in and getting the best out of them is of great importance.
Pregnant women are a particularly important patient population to study due to their vulnerability to disease and the often underrepresentation of the population in clinical research [2] . There has been a relative sparcity of data in regards to SARS-CoV-2 and its effect in pregnancy. Utilizing ML techniques to study this population during the pandemic can help build pregnancy-specific evidence to guide clinical recommendations.
The key to have successful predictive methods largely depends on feature selection and data representation. A common approach is to have a clinical doctor specify the variables and label the clinical data to be used as training sets. Then the ML method will find the mapping from data and will be tested on new data sets. Although appropriate in many situations, a supervised definition of the features contributes to losing an opportunity to learn latent patterns and features [3] . In combating the subjectivity of defining the features, an unsupervised learning approach can be used to extract useful information from data. One might argue that a limitation of this method is that patients are often represented in a low-dimensional space that summarizes the information available in the data but with losses in information [3] .
In the ML literature, Gaussian processes (GPs) provide data efficient and powerful Bayesian framework for learning latent functions or patterns [4] . It can be shown that a GP is equivalent to a single layer of fully connected neural network with an i.i.d. prior over its parameters and with infinite network width [5] . GPs have been successfully applied in supervised and unsupervised task of learning from medical data [6] . In particular, GP latent variable models (GP-LVMs) are used for Bayesian unsupervised learning. They are endowed with the ability to automatically learn latent structures and dependencies in the data [7] .
When the data consist of categorical variables, i.e., all the possible examinations, co-morbidities and patient's symptoms in medical data, we face with long vectors of categorical variables which lend to huge numbers of possible realizations. This in turn creates very sparse spaces when we deal with a limited number of data.
The problem of sparsity usually occurs when the dataset diversity is poor and the number of patients is relatively small. The recent models [8] - [9] , has achieved significant improvements in performance. For example, the authors in [9] , has non-linearly embedded the observations in the continuous space using Gaussian processes.
In this paper, we explore modeling vectors of categorical variables within the Bayesian nonparametric framework, and in particular, we work with Gaussian process latent variable models (GP-LVMs). We focus on categorical latent GPs and use them in the context of GP-LVM models. We study data related to test-positive COVID-19 pregnant women who tested positive between March 13, 2020, and August 08, 2020 at Stony Brook University Hospital (SBUH).
Our contribution in the paper is in formulating the problem of density estimation of categorical data in test-positive COVID-19 pregnant women. The purpose of this effort was to discover hidden patterns among the women that could be used for prediction of outcomes based on a number of categorical variables. For example, we could use the found latent space to find similarities among patients and use them to predict the severity of the disease for a particular patient given the "location" of the patient in the latent space. Further, such studies could allow for easy visualization of the cohort of patients under consideration.
The remainder of this paper is organized as follows. In the next section, we provide a brief overview of GP-LVMs. In Section 3, we present the details of the model used for our multivariate categorical data analysis and its inference process. In Section 4, we show our results from both simulations and real COVID-19 data. Finally, with Section 5, we conclude the paper with final remarks.
Gaussian process latent variable models (GP-LVMs) are Bayesian nonparametric frameworks that offer unsupervised learning [7] . GP-LVMs can be seen as multi-output Gaussian process regression when the inputs are unobserved. To be more specific, if Y ∈ R N ×D is the observed data (where N is the number of observations and D is the dimensionality of each data point), we can say that these data are latently associated with the inputs X ∈ R N ×Q through D different GPs. If we assume that these GPs are independent, then the likelihood function can be written as
where y d represents the d-th feature of Y and
The covariance matrix KNN is of size N × N and is defined by the covariance function (also called kernel), k (x, x ). In order to automatically learn the dimensionality of the latent space, we will use Automatic Relevance Determination (ARD) with the kernel
Since in GP-LVMs, X is a latent variable, we should assign it a prior density and then try to find the MAP estimate of it. For example, if we use the standard Gaussian distribution, we will have
where the xn's are the rows of X. Now that we have the distribution on X, we can write the joint probability of X and Y as
Next we need to find the MAP estimate of X by optimizing the kernel hyperparameters θ = σ 2 f , α1, . . . , αQ and the inverse variance parameter σ [10] .
In this section we present the method based on GP-LVM for distribution estimation of multivariate categorical data, first proposed in [9] . We first present the underlying generative model, and then we show how to infer the hidden model variables given the observed data.
We consider a generative model for a dataset Y with N observations and D categorical variables. We show the d-th categorical variable in the n-th observation by y nd . Each categorical variable y nd can take values, for example, from 1 to K. Each of them are samples from their corresponding probability mass function with weights f nd = (f nd1 , . . . , f ndK ). We assume each random variable f ndk is a nonlinear function of the input variable xn ∈ R Q , where Q is the dimension of xn. Therefore, f ndk = F dk (xn) .
Since x and F are latent we need to assign prior distributions for them. We assume a Gaussian distribution prior with standard deviation σ 2
x for the latent variables x, and a Gaussian process (GP) prior for each of the functions F. We also consider a set of M variational inducing points. For each vector of latent function values f dk , we introduce a separate set of M inducing variables u dk , evaluated at a set of inducing input locations given by Z. It is assumed that all u dk s are evaluated at the same inducing locations. The inducing variables are function points drawn from the GP prior and lie in same latent space as the F variables (Fig. 1 ). The generative model can be summarized as
. The pictorial description of the generative model is displayed in Fig.2 . 
The marginal log-likelihood is intractable because of the covariance function of the GP and the Softmax likelihood. We will consider a variational approximation to the posterior distribution of X, F and U factorized as, q(X, F, U) = q(X)q(U)p(F|X, U).
By applying Jensen's inequality, we can write a lower bound of the log-evidence (ELBO) as
where,
The lower bound is still intractable because of the softmax likelihood, log p (y nd | f nd ). Therefore, we will compute the lower bound L and its derivatives with the Monte Carlo method. We draw samples of xn, U d and f nd from q (xn) , q (U d ) , and p (f nd | xn, U d ) respectively and estimate L with the sample average. We consider mean field variational approximation for the latent points q(X) and a joint Gaussian distribution for q(U) as,
where the covariance matrix Σ d is shared for the same categorical variable d. The KL divergence in L can be computed analytically with the given variational distributions. We need to optimize the hyperparameters of each GP (parameters of K d ), parameters of the variational random variables u dk , µ dk , Σ d , mean mnq and variance σ 2 nq of the latent inputs. The graphical model of the variational distributions for the inference part is shown in Fig. 3 . For further details about the inference and learning hyperparameters refer to [9] .
In this section, we present our results with synthetic data. We first generated a dataset Y for N = 100 patients and 10 categorical variables as shown in Fig. 2 . These patients belonged to two clusters. We assumed each categorical variable can take values 0 and 1. In this way, we only needed to estimate the probability of one category, since Pr{k = 0} = 1 − Pr{k = 1}. The probability of the d-th category to be 1 for patient xn was proportional to the output of the function f d (xn), Fig. 4 (a) . We generated samples xn from a mixture of two Gaussian distributions in order to see if the model could detect the underlying cluster of patients. These two clusters are illustrated in Fig. 4 (a) . Note that the dimension of the latent space was one. However, we initialized the dimension at a higher value and let the model learn the lower dimension of the latent space. Figure 4 (b) shows the inferred latent space for the first two dimensions. Figure 4 (c) plotted the train error. We can infer the dimensionality of the latent space by comparing the maximized ELBO for different initialization of Q. 
We used data collected at SBUH of test-positive COVID-19 pregnant women. The dataset was composed of categorical data of 89 patients. All the variables are listed in Table 1 . The COVID-19 patients may carry severe or mild symptoms, and some of them end up in an ICU. There have been studies that utilized ML algorithms to build a detection model for the severeness of COVID-19. The authors in [11] have investigated the binary classification problem between severe and mild cases. The problem with these methods is that they do not provide uncertainty about the prediction. Moreover, they need subjective definition of features and outcomes. Most of the time, disease severity levels are more than just described by two classes. Thus, it is important to study a range of classes of disease severeness, from mild, such as common cough, various in-between levels of disease, to harsh that entail ending up in ICUs. This is why it is important to employ methods like categorical latent GPs to decide on their own on the number of categories of patients and based on the patterns they discover in the patients' data.
In this analysis, we mapped the patient population into a twodimensional space using only the symptoms at the time of diagnosis. Then we labeled the points by severity of the outcomes. The latent space representing outcomes is shown in Fig. 5 . By looking at the mapped data, we can easily see that there is a population of patients at x = 0.1 most of which did not need hospitalization. Instead, all the patients that have been admitted to the hospital are concentrated at the upper left corner of the space. Among them, three patients were admitted to ICU after developing severe outcomes. Here we have also observed one patient who is an outlier because that patient has been in ICU, and in the latent space she was lying in the far right of the space. Note that the pregnant women under investigation could have been in the hospital for reasons other than COVID. Many of them tested positive when they were completely asymptomatic.
Finally we mapped patients using all the categorical variables from Table 1 into a two-dimensional latent space. The latent points and the clustering of the patients are shown in Fig. 6 . The intensity of green color in the figure reflects the density of the population in the latent space. The highest density is around the rightmost group of patients. One might argue that there are three clusters of patients, of which the one in the middle is the least concentrated. These clusters allow us to study the relationships among the patients and in particular, understand why patients grouped in a cluster are similar to other patients in the same cluster.
In this paper, we analyzed multivariate categorical data with models based on categorical latent Gaussian processes. With these models, we can discover much lower dimensional latent spaces that can facilitate classification, prediction and visualization. More specifically, we used a data-efficient Bayesian framework for clustering of highdimensional categorical data. Our tests with synthetic data showed that the method is capable of finding latent structures of the data. Further, we applied the method to data obtained from test-positive COVID-19 pregnant women. There, too, the method discovered latent structures. These structures can be useful for many purposes including gaining important insights of high interest to physicians. 