The package hdpGLM makes it easy to estimate semi-parametric regression models, and summarize and visualize the results. The package is useful for many purposes:

- Find clusters using semi-parametric Bayesian methods (Dirichlet Process)
- Find a set of generalized linear models (GLM) that apply to different latent subpopulations
- Estimate latent heterogeneity in the marginal effect of regression coefficients
- Investigate if Simpson’s Paradox occurs due to latent or omitted features
- Cluster the data points based on marginal effect heterogeneity
- All the above, but with hierarchical data from many contexts (e.g., schools, hospitals, cities, states, countries, etc.)
- Find the effect of context-level features (upper-level covariates) on the latent heterogeneity of within-context features (lower-level covariates)

The `hdpGLM`

works similarly to linear regression models. It estimates coefficients of linear regressions, including generalized linear models, such as logit coefficients. But it simultaneously searches for latent clusters in the data and estimates the linear coefficients for those clusters. The result of the estimation is \(K\) vectors of linear coefficients, one vector for each cluster. If there are no clusters in the data, it returns the estimated coefficients similar to classical regression models.

The clustering procedure is based on a hierarchical semi-parametric Bayesian model proposed in Ferrari (2020). The model, called Hierarchical Dirichlet process Generalized Linear Models (hdpGLM), can be used to deal with latent heterogeneity in different situations, including those that emerge due to unobserved variables. It deals with the latent heterogeneity in two ways: (1) it finds latent clusters which can be better described by different regression models and (2) estimate the coefficient of those models. The hdpGLM can also be used with hierarchical data to estimate latent heterogeneity in multiple contexts and check if the clusters are context-dependent (see an example in section Estimating Context-dependent latent heterogeneity).

The linear model is estimated by sampling from the posterior distribution using a Gibbs sampler. Non-linear models (e.g., logit with binary outcomes) use Hamiltonian Monte Carlo within Gibbs. The algorithms are presented in Ferrari (2020).

Why should we estimate clusters of linear regressions instead of fitting a single regression model?

First, it improves the predictive performance of the regression model and keeps the interpretability of the regression coefficients. `hdpGLM`

is much more flexible than traditional regression and produces monotonically lower root mean square error (see Ferrari (2020) for details).

Second, latent heterogeneity emerges when there are omitted variables in the estimation of regression models. The `hdpGLM`

can be used to estimate marginal effects even when interactions were omitted. It recovers the linear coefficients of each latent group.

The function `hdpGLM`

estimates a semi-parametric Bayesian regression model. The syntax is similar to other R functions such as `lm()`

, `glm()`

, and `lmer()`

.

Here is a toy example. Suppose we are studying how income inequality affects support policies that help alleviate poverty in a given country A. Yet, suppose further that (1) the effect of inequality varies between groups of people; for some people, inequality increases support for welfare policies, but for others, it decreases welfare policy support; (2) we don’t know which individual belongs to which group. The data set `welfare`

contains simulated data for this example.

```
## loading and looking at the data
welfare = read.csv2('welfare.csv')
head(welfare)
#> support inequality income ideology
#> 1 -18.5649610 0.3392724 0.1425111 1.9225985
#> 2 -9.3905812 -0.9906646 -0.5117102 0.2483346
#> 3 0.9276234 -2.2318510 -0.3856288 -1.3619216
#> 4 -12.3594498 -3.0079501 -0.9440585 -0.2088675
#> 5 -2.4834411 0.1000455 0.8322192 0.1321378
#> 6 -11.4187853 -0.9543883 -0.8810503 0.2916444
```

Now, suppose that inequality increases support for welfare only among women, but it decreases support among men. We didn’t collect data on gender (male versus female). We could estimate the `hdpGLM`

and recover the coefficients even if gender wasn’t observed. The package provides a function called `hdpGLM`

, which estimates a semi-parametric Bayesian generalized linear model using a Dirichlet mixture. Let’s estimate the model. The example uses few iterations in the MCMC, but in real applications, one should use a much larger number.

```
library(hdpGLM)
#>
#> ## =======================================================================
#> ## Hierarchial Dirichlet Process Generalized Linear Model (hdpGLM)
#> ## =======================================================================
#>
#> Attaching package: 'hdpGLM'
#> The following object is masked _by_ '.GlobalEnv':
#>
#> welfare
```

```
## estimating the model
mcmc = list(burn.in=10, ## MCMC burn-in period
n.iter =500) ## number of MCMC iterations to keep
mod = hdpGLM(support ~ inequality + income + ideology, data=welfare,
mcmc=mcmc)
#> Warning in model.matrix.default(mt1, reg.matrix, stats::contrasts): non-list
#> contrasts argument ignored
```

```
## printing the outcome
summary(mod)
#> Warning: `as_data_frame()` is deprecated as of tibble 2.0.0.
#> Please use `as_tibble()` instead.
#> The signature and semantics have changed, see `?as_tibble`.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> # A tibble: 10 x 8
#> k Parameter term Mean Median SD HPD.lower HPD.upper
#> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 beta1 (Intercept) -3.87 -3.87 0.0378 -3.94 -3.80
#> 2 1 beta2 inequality 2.00 2.00 0.0358 1.93 2.07
#> 3 1 beta3 income 3.85 3.85 0.0374 3.78 3.92
#> 4 1 beta4 ideology -8.31 -8.31 0.0387 -8.39 -8.24
#> 5 1 sigma sigma 0.999 0.999 0.0444 0.914 1.08
#> 6 10 beta1 (Intercept) -3.82 -3.82 0.0373 -3.90 -3.75
#> 7 10 beta2 inequality -1.53 -1.53 0.0353 -1.59 -1.45
#> 8 10 beta3 income 3.88 3.88 0.0349 3.81 3.95
#> 9 10 beta4 ideology -8.26 -8.26 0.0375 -8.33 -8.18
#> 10 10 sigma sigma 0.968 0.969 0.0606 0.854 1.08
```

The summary function prints the result in a tidy format. The column `k`

in the summary shows the label of the estimated clusters. The column `Mean`

is the average of the posterior distribution for each linear coefficient in each cluster.

The function `hdpGLM_classify`

can be used to classify the data points into clusters based on the estimation.

```
welfare_clustered = hdpGLM_classify(welfare, mod)
head(welfare_clustered)
#> Cluster support inequality income ideology
#> 1 1 -18.5649610 0.3392724 0.1425111 1.9225985
#> 2 1 -9.3905812 -0.9906646 -0.5117102 0.2483346
#> 3 1 0.9276234 -2.2318510 -0.3856288 -1.3619216
#> 4 1 -12.3594498 -3.0079501 -0.9440585 -0.2088675
#> 5 10 -2.4834411 0.1000455 0.8322192 0.1321378
#> 6 1 -11.4187853 -0.9543883 -0.8810503 0.2916444
tail(welfare_clustered)
#> Cluster support inequality income ideology
#> 1995 10 -1.5230053 1.055855140 -0.7295937 -0.7067871
#> 1996 10 0.4814892 0.582588091 2.0051082 0.3090389
#> 1997 10 -14.1929956 0.391164197 -0.9607449 0.7765482
#> 1998 10 -8.2396789 0.074437376 1.2020300 1.0874928
#> 1999 10 -23.1583753 0.434223018 -0.6176438 2.0387294
#> 2000 1 -7.2075582 0.008355317 -0.4538951 0.2268072
```

There are a series of built-in functions, with various options, to plot the results. In the example below, you see two of those options. The `separate`

parameter plot the posterior samples for each cluster separately, and the option `ncols`

controls how many columns to use for the panels in the figure (to see more, run `help(plot.hdpGLM)`

and `help(plot.dpGLM)`

).

To continue the previous toy example, suppose that we are analyzing data from many countries, and we suspect that the latent heterogeneity is different in each country. The effect of inequality on support for welfare may be gender-specific only in some countries (contexts). Or maybe the way it is gender-specific varies from country to country. Suppose we didn’t have data on gender, but we collect information on countries’ gender gap in welfare provision. Let’s look at this new data set.

```
## loading and looking at the data
welfare = read.csv2('welfare2.csv')
head(welfare)
#> support inequality income ideology country gap
#> 1 -18.5649610 0.3392724 0.1425111 1.9225985 0 0.1
#> 2 -9.3905812 -0.9906646 -0.5117102 0.2483346 0 0.1
#> 3 0.9276234 -2.2318510 -0.3856288 -1.3619216 0 0.1
#> 4 -12.3594498 -3.0079501 -0.9440585 -0.2088675 0 0.1
#> 5 -2.4834411 0.1000455 0.8322192 0.1321378 0 0.1
#> 6 -11.4187853 -0.9543883 -0.8810503 0.2916444 0 0.1
tail(welfare)
#> support inequality income ideology country gap
#> 3195 0.3190583 -0.7504798 -0.7839583 0.92300705 4 -0.8280808
#> 3196 -1.3837239 0.6620435 -1.5566268 0.05634618 4 -0.8280808
#> 3197 -1.3820016 -0.4298706 -1.0945688 0.71559078 4 -0.8280808
#> 3198 0.6878775 0.5450604 2.6175887 -1.94844469 4 -0.8280808
#> 3199 -7.9282930 1.7846004 1.6755823 1.29160208 4 -0.8280808
#> 3200 -1.7472485 0.5030992 -0.5395479 0.20109879 4 -0.8280808
```

The variable `country`

indicates the country (context) of the observation, and the variable `gap`

the gender gap in welfare provision in the respective country. The estimation is similar to the previous example, but now there is a second `formula`

for the context-level variables. Again, the example below uses few iterations in the MCMC, but in practical applications, one needs to increase that).

```
## estimating the model
mcmc = list(burn.in=1, ## MCMC burn-in period
n.iter =50) ## number of MCMC iterations to keep
mod = hdpGLM(support ~ inequality + income + ideology,
support ~ gap,
data=welfare, mcmc=mcmc)
#> Warning in model.matrix.default(mt1, reg.matrix, stats::contrasts): non-list
#> contrasts argument ignored
#> Warning in model.matrix.default(mt1, reg.matrix, stats::contrasts): non-list
#> contrasts argument ignored
```

```
summary(mod)
#>
#> Generating summary for beta...
#> Warning: `data_frame()` is deprecated as of tibble 1.1.0.
#> Please use `tibble()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#>
#> Generating summary for tau...
#> $beta
#> # A tibble: 110 x 9
#> k j Parameter term Mean Median SD HPD.lower HPD.upper
#> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 1 beta[0] (Intercept) -3.87 -3.87 0.0275 -3.93 -3.81
#> 2 1 1 beta[1] inequality 1.71 1.96 0.590 0.168 2.05
#> 3 1 1 beta[2] income 3.86 3.86 0.0415 3.78 3.96
#> 4 1 1 beta[3] ideology -8.34 -8.34 0.0456 -8.44 -8.25
#> 5 1 1 sigma sigma 1.00 0.999 0.0391 0.952 1.11
#> 6 1 2 beta[0] (Intercept) -0.128 0.0136 0.379 -0.669 0.395
#> 7 1 2 beta[1] inequality -0.776 -0.725 0.443 -1.05 0.567
#> 8 1 2 beta[2] income 0.432 0.378 0.441 -0.697 0.993
#> 9 1 2 beta[3] ideology -3.16 -3.08 1.09 -3.77 0.236
#> 10 1 2 sigma sigma 0.0923 0.0910 0.0369 0.0497 0.108
#> # … with 100 more rows
#>
#> $tau
#> # A tibble: 8 x 8
#> w beta Parameter Description Mean SD HPD.lower HPD.upper
#> <int> <int> <fct> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 0 0 tau[0][0] Intercept of beta[0] -0.566 1.58 -3.97 1.98
#> 2 0 1 tau[0][1] Intercept of beta[1] -0.207 1.37 -2.97 2.24
#> 3 0 2 tau[0][2] Intercept of beta[2] -0.269 1.52 -3.59 2.12
#> 4 0 3 tau[0][3] Intercept of beta[3] -1.04 1.27 -3.32 1.30
#> 5 1 0 tau[1][0] Effect of gap on beta… 0.137 1.20 -2.16 2.03
#> 6 1 1 tau[1][1] Effect of gap on beta… 0.0860 1.30 -2.50 2.32
#> 7 1 2 tau[1][2] Effect of gap on beta… -0.935 1.20 -3.23 0.973
#> 8 1 3 tau[1][3] Effect of gap on beta… 0.857 1.01 -1.06 2.64
```

The `summary`

contains more information now. As before, the column `k`

indicates the estimated clusters. The column `j`

indicates the country (context) of the estimated value for the respective cluster’s coefficient. The second summary (`$tau`

) shows the marginal effect of the context-level feature (`gap`

). Details of the interpretation can be found in Ferrari (2020).

There are a series of built-in functions to visualize the output. The function `plot_tau()`

displays the estimation of the effect of the context-level variables.

The function `plot_pexp_beta()`

displays the association between the context-level features and the latent heterogeneity in the effect of the linear coefficients in each context. The paramter ‘smooth.line’ plots a line representing the linear association between the context-level feature (`gap`

) and the posterior averages of the marginal effects in each cluster. The parameter `ncol.beta`

controls the number of columns in the figure for the panels. For more options, see `help(plot_pexp_beta)`

Ferrari, Diogo. 2020. “Modeling Context-Dependent Latent Effect Heterogeneity.” *Political Analysis* 28 (1): 20–46.