{"id":855,"date":"2024-12-28T07:03:39","date_gmt":"2024-12-28T07:03:39","guid":{"rendered":"https:\/\/mailitics.com\/index.php\/2024\/12\/28\/introduction-to-the-finite-normal-mixtures-in-regression-with-6a884810a692\/"},"modified":"2024-12-28T07:03:39","modified_gmt":"2024-12-28T07:03:39","slug":"introduction-to-the-finite-normal-mixtures-in-regression-with-6a884810a692","status":"publish","type":"post","link":"https:\/\/mailitics.com\/index.php\/2024\/12\/28\/introduction-to-the-finite-normal-mixtures-in-regression-with-6a884810a692\/","title":{"rendered":"Introduction to the Finite Normal Mixtures in Regression with"},"content":{"rendered":"

Introduction to the Finite Normal Mixtures in Regression with
\n \t

\n
<\/BR>
\n
\n \t

\n
<\/BR><\/p>\n

\n

Introduction to the Finite Normal Mixtures in Regression with\u00a0R<\/h3>\n

How to make linear regression flexible enough for non-linear data<\/h4>\n

The linear regression is usually considered not flexible enough to tackle the nonlinear data. From theoretical viewpoint it is not capable to dealing with them. However, we can make it work for us with any dataset by using finite normal mixtures in a regression model. This way it becomes a very powerful machine learning tool which can be applied to virtually any dataset, even highly non-normal with non-linear dependencies across the variables.<\/p>\n

What makes this approach particularly interesting comes with interpretability. Despite an extremely high level of flexibility all the detected relations can be directly interpreted. The model is as general as neural network, still it does not become a black-box. You can read the relations and understand the impact of individual variables.<\/p>\n

In this post, we demonstrate how to simulate a finite mixture model for regression using Markov Chain Monte Carlo (MCMC) sampling. We will generate data with multiple components (groups) and fit a mixture model to recover these components using Bayesian inference. This process involves regression models and mixture models, combining them with MCMC techniques for parameter estimation.<\/p>\n

\"\"
Data simulated as a mixtures of three linear regressions<\/figcaption><\/figure>\n

Loading Required Libraries<\/h3>\n

We begin by loading the necessary libraries to work with regression models, MCMC, and multivariate distributions<\/p>\n

# Loading the required libraries for various functions
library(\"pscl\")         # For pscl specific functions, like regression models
library(\"MCMCpack\")     # For MCMC sampling functions, including posterior distributions
library(mvtnorm)        # For multivariate normal distribution functio<\/pre>\n
    \n
  • \npscl<\/strong>: Used for various statistical functions like regression models.<\/li>\n
  • \nMCMCpack<\/strong>: Contains functions for Bayesian inference, particularly MCMC sampling.<\/li>\n
  • \nmvtnorm<\/strong>: Provides tools for working with multivariate normal distributions.<\/li>\n<\/ul>\n
    Data Generation<\/h3>\n
    We simulate a dataset where each observation belongs to one of several groups (components of the mixture model), and the response variable is generated using a regression model with random coefficients.<\/p>\n
    We consider a general setup for a regression model using G Normal mixture components.<\/p>\n
    \"\"<\/figure>\n
    ## Generate the observations
    # Set the length of the time series (number of observations per group)
    N <- 1000
    # Set the number of simulations (iterations of the MCMC process)
    nSim <- 200
    # Set the number of components in the mixture model (G is the number of groups)
    G <- 3<\/pre>\n
      \n
    • \nN<\/strong>: The number of observations per\u00a0group.<\/li>\n
    • \nnSim<\/strong>: The number of MCMC iterations.<\/li>\n
    • \nG<\/strong>: The number of components (groups) in our mixture\u00a0model.<\/li>\n<\/ul>\n
      Simulating Data<\/h4>\n
      Each group is modeled using a univariate regression model, where the explanatory variables (X) and the response variable (y) are simulated from normal distributions. The betas represent the regression coefficients for each group, and sigmas represent the variance for each\u00a0group.<\/p>\n
      # Set the values for the regression coefficients (betas) for each group
      betas <- 1:sum(dimG) * 2.5  # Generating sequential betas with a multiplier of 2.5
      # Define the variance (sigma) for each component (group) in the mixture
      sigmas <- rep(1, G) \/ 1  # Set variance to 1 for each component, with a fixed divisor of 1<\/pre>\n
        \n
      • \nbetas<\/strong>: These are the regression coefficients. Each group\u2019s coefficient is sequentially assigned.<\/li>\n
      • \nsigmas<\/strong>: Represents the variance for each group in the mixture\u00a0model.<\/li>\n<\/ul>\n
        In this model we allow each mixture component to possess its own variance paraameter and set of regression parameters.<\/p>\n
        Group Assignment and\u00a0Mixing<\/h4>\n
        We then simulate the group assignment of each observation using a random assignment and mix the data for all components.<\/p>\n
        We augment the model with a set of component label vectors\u00a0for<\/p>\n
        \"\"<\/figure>\n
        where<\/p>\n
        \"\"<\/figure>\n
        and thus z_gi=1<\/em> implies that the i-<\/em>th individual is drawn from the g-<\/em>th component of the\u00a0mixture.<\/p>\n
        This random assignment forms the z_original vector, representing the true group each observation belongs\u00a0to.<\/p>\n
        # Initialize the original group assignments (z_original)
        z_original <- matrix(NA, N * G, 1)
        # Repeat each group label N times (assign labels to each observation per group)
        z_original <- rep(1:G, rep(N, G))
        # Resample the data rows by random order
        sampled_order <- sample(nrow(data))
        # Apply the resampled order to the data
        data <- data[sampled_order,]<\/pre>\n
        Bayesian Inference: Priors and Initialization<\/h3>\n
        We set prior distributions for the regression coefficients and variances. These priors will guide our Bayesian estimation.<\/p>\n
        \"\"<\/figure>\n
        ## Define Priors for Bayesian estimation# Define the prior mean (muBeta) for the regression coefficients
        muBeta <- matrix(0, G, 1)# Define the prior variance (VBeta) for the regression coefficients
        VBeta <- 100 * diag(G)  # Large variance (100) as a prior for the beta coefficients# Prior for the sigma parameters (variance of each component)
        ag <- 3  # Shape parameter
        bg <- 1\/2  # Rate parameter for the prior on sigma
        shSigma <- ag
        raSigma <- bg^(-1)<\/pre>\n
          \n
        • \nmuBeta<\/strong>: The prior mean for the regression coefficients. We set it to 0 for all components.<\/li>\n
        • \nVBeta<\/strong>: The prior variance, which is large (100) to allow flexibility in the coefficients.<\/li>\n
        • \nshSigma<\/strong> and raSigma<\/strong>: Shape and rate parameters for the prior on the variance (sigma) of each\u00a0group.<\/li>\n<\/ul>\n
          For the component indicators and component probabilities we consider following prior assignment<\/p>\n
          \"\"<\/figure>\n
          The multinomial prior M is the multivariate generalizations of the binomial, and the Dirichlet prior D is a multivariate generalization of the beta distribution.<\/p>\n
          MCMC Initialization<\/h3>\n
          In this section, we initialize the MCMC process by setting up matrices to store the samples of the regression coefficients, variances, and mixing proportions.<\/p>\n
          ## Initialize MCMC sampling# Initialize matrix to store the samples for beta
          mBeta <- matrix(NA, nSim, G)# Assign the first value of beta using a random normal distribution
          for (g in 1:G) {
            mBeta[1, g] <- rnorm(1, muBeta[g, 1], VBeta[g, g])
          }# Initialize the sigma^2 values (variance for each component)
          mSigma2 <- matrix(NA, nSim, G)
          mSigma2[1, ] <- rigamma(1, shSigma, raSigma)# Initialize the mixing proportions (pi), using a Dirichlet distribution
          mPi <- matrix(NA, nSim, G)
          alphaPrior <- rep(N\/G, G)  # Prior for the mixing proportions, uniform across groups
          mPi[1, ] <- rdirichlet(1, alphaPrior)<\/pre>\n
            \n
          • \nmBeta<\/strong>: Matrix to store samples of the regression coefficients.<\/li>\n
          • \nmSigma2<\/strong>: Matrix to store the variances (sigma squared) for each component.<\/li>\n
          • \nmPi<\/strong>: Matrix to store the mixing proportions, initialized using a Dirichlet distribution.<\/li>\n<\/ul>\n
            MCMC Sampling: Posterior Updates<\/h3>\n
            If we condition on the values of the component indicator variables z, the conditional likelihood can be expressed as<\/p>\n
            \"\"<\/figure>\n
            In the MCMC sampling loop, we update the group assignments (z), regression coefficients (beta), and variances (sigma) based on the posterior distributions. The likelihood of each group assignment is calculated, and the group with the highest posterior probability is selected.<\/p>\n
            The following complete posterior conditionals can be obtained:<\/p>\n
            \"\"<\/figure>\n
            where<\/p>\n
            \"\"<\/figure>\n
            denotes all the parameters in our posterior other than\u00a0x<\/em>.<\/p>\n
            \"\"<\/figure>\n
            and where n_g<\/em> denotes the number of observations in the g<\/em>-th component of the\u00a0mixture.<\/p>\n
            \"\"<\/figure>\n
            and<\/p>\n
            \"\"<\/figure>\n
            Algorithm below draws from the series of posterior distributions above in a sequential order.<\/p>\n
            ## Start the MCMC iterations for posterior sampling# Loop over the number of simulations
            for (i in 2:nSim) {
              print(i)  # Print the current iteration number
              
              # For each observation, update the group assignment (z)
              for (t in 1:(N*G)) {
                fig <- NULL
                for (g in 1:G) {
                  # Calculate the likelihood of each group and the corresponding posterior probability
                  fig[g] <- dnorm(y[t, 1], X[t, ] %*% mBeta[i-1, g], sqrt(mSigma2[i-1, g])) * mPi[i-1, g]
                }
                # Avoid zero likelihood and adjust it
                if (all(fig) == 0) {
                  fig <- fig + 1\/G
                }
                
                # Sample a new group assignment based on the posterior probabilities
                z[i, t] <- which(rmultinom(1, 1, fig\/sum(fig)) == 1)
              }
              
              # Update the regression coefficients for each group
              for (g in 1:G) {
                # Compute the posterior mean and variance for beta (using the data for group g)
                DBeta <- solve(t(X[z[i, ] == g, ]) %*% X[z[i, ] == g, ] \/ mSigma2[i-1, g] + solve(VBeta[g, g]))
                dBeta <- t(X[z[i, ] == g, ]) %*% y[z[i, ] == g, 1] \/ mSigma2[i-1, g] + solve(VBeta[g, g]) %*% muBeta[g, 1]
                
                # Sample a new value for beta from the multivariate normal distribution
                mBeta[i, g] <- rmvnorm(1, DBeta %*% dBeta, DBeta)
                
                # Update the number of observations in group g
                ng[i, g] <- sum(z[i, ] == g)
                
                # Update the variance (sigma^2) for each group
                mSigma2[i, g] <- rigamma(1, ng[i, g]\/2 + shSigma, raSigma + 1\/2 * sum((y[z[i, ] == g, 1] - (X[z[i, ] == g, ] * mBeta[i, g]))^2))
              }
              
              # Reorder the group labels to maintain consistency
              reorderWay <- order(mBeta[i, ])
              mBeta[i, ] <- mBeta[i, reorderWay]
              ng[i, ] <- ng[i, reorderWay]
              mSigma2[i, ] <- mSigma2[i, reorderWay]
              
              # Update the mixing proportions (pi) based on the number of observations in each group
              mPi[i, ] <- rdirichlet(1, alphaPrior + ng[i, ])
            }<\/pre>\n
            This block of code performs the key steps in\u00a0MCMC:<\/p>\n
              \n
            • \nGroup Assignment Update<\/strong>: For each observation, we calculate the likelihood of the data belonging to each group and update the group assignment accordingly.<\/li>\n
            • \nRegression Coefficient Update<\/strong>: The regression coefficients for each group are updated using the posterior mean and variance, which are calculated based on the observed\u00a0data.<\/li>\n
            • \nVariance Update<\/strong>: The variance of the response variable for each group is updated using the inverse gamma distribution.<\/li>\n<\/ul>\n
              Visualizing the\u00a0Results<\/h3>\n
              Finally, we visualize the results of the MCMC sampling. We plot the posterior distributions for each regression coefficient, compare them to the true values, and plot the most likely group assignments.<\/p>\n
              # Plot the posterior distributions for each beta coefficient
              par(mfrow=c(G,1))
              for (g in 1:G) {
                plot(density(mBeta[5:nSim, g]), main = 'True parameter (vertical) and the distribution of the samples')  # Plot the density for the beta estimates
                abline(v = betas[g])  # Add a vertical line at the true value of beta for comparison
              }<\/pre>\n
              This plot shows how the MCMC samples (posterior distribution) for the regression coefficients converge to the true values\u00a0(betas).<\/p>\n
              Conclusion<\/h3>\n
              Through this process, we demonstrated how finite normal mixtures can be used in a regression context, combined with MCMC for parameter estimation. By simulating data with known groupings and recovering the parameters through Bayesian inference, we can assess how well our model captures the underlying structure of the\u00a0data.<\/p>\n