Determining conditional independence relationships through undirected graphical models is a key component in the statistical analysis of complex obervational data in a wide variety of disciplines. In many situations one seeks to estimate the underlying graphical model of a dataset that includes variables of different domains. As an example, take a typical dataset in the social, behavioral and medical sciences, where one is interested in interactions, for example between gender or country (categorical), frequencies of behaviors or experiences (count) and the dose of a drug (continuous). Other examples are Internet-scale marketing data or high-throughput sequencing data. There are methods available to estimate mixed graphical models from mixed continuous data, however, these usually have two drawbacks: first, there is a possible information loss due to necessary transformations and second, they cannot incorporate (nominal) categorical variables (for an overview see here). A new method implemented in the R-package mgm addresses these limitations. In the following, we use the mgm-package to estimate the conditional independence network in a dataset of questionnaire responses of individuals diagnosed with Autism Spectrum Disoder. This dataset includes variables of different domains, such as age (continuous), type of housing (categorical) and number of treatments (count). The dataset consists of responses of 3521 individuals to a questionnaire including 28 variables of domains continuous, count and categorical. We now use our knowledge about the variables to specify the domain (or type) of each variable and the number of categories for categorical variables (for non-categorical variables we choose 1). “c”, “g”, “p” stands for categorical, Gaussian and Poisson (count), respectively. The estimation algorithm requires us to make an assumption about the highest order interaction in the true graph. Here we assume that there are at most pairwise interactions in the true graph and set d = 2. The algorithm includes an L1-penalty to obtain a sparse estimate. We can select the regularization parameter lambda using cross validation (CV) or the Extended Bayesian Information Criterion (EBIC). Here, we choose the EBIC, which is known to be a bit more conservative than CV but is computationally faster. The fit function returns all estimated parameters and a weighted and unweighted (binarized) adjacency matrix. Here we use the qgraph package to visualize the graph: A reproducible example can be found in the examples of the package or more elaboratly explained in the corresponding paper. Here is a paper explaining the theory behind the implemented algorithm. Computationally efficient methods for Gaussian data are implemented in the huge package and the glasso package. For binary data, there is the IsingFit package. Great free resources about graphical models are Chapter 17 in the freely available book The Elements of Statistical Learning and the Coursera course Probabilistic Graphical Models.