Random scaling factors in Bayesian distributional regression models with an application to real estate data

Abstract

Distributional structured additive regression provides a flexible framework for modelling each parameter of a potentially complex response distribution in dependence of covariates. Structured additive predictors allow for an additive decomposition of covariate effects with non-linear effects and time trends, unit- or cluster-specific heterogeneity, spatial heterogeneity and complex interactions between covariates of different type. Within this framework, we present a simultaneous estimation approach for multiplicative random effects that allow for cluster-specific heterogeneity with respect to the scaling of a covariate′s effect. More specifically, a possibly non-linear function f(z) of a covariate z may be scaled by a multiplicative and possibly spatially correlated cluster-specific random effect (1+α_c). Inference is fully Bayesian and is based on highly efficient Markov Chain Monte Carlo (MCMC) algorithms.

We investigate the statistical properties of our approach within extensive simulation experiments for different response distributions. Furthermore, we apply the methodology to German real estate data where we identify significant district-specific scaling factors. According to the deviance information criterion, the models incorporating these factors perform significantly better than standard models without (spatially correlated) random scaling factors.

Keywords

Bayesian P-splines Iteratively weighted least squares proposals Markov random fields MCMC multiplicative random effects structured additive predictors

1 Introduction

Classical regression models, such as generalized linear models (GLMs, see McCullagh and Nelder, 1989), generalized additive models (GAMs, see Hastie and Tibshirani, 1990; Wood, 2006) or structured additive regression models (STAR models, see Brezger and Lang, 2006, or Fahrmeir et al., 2013), assume that the distribution of a response variable y belongs to an exponential family and relate its mean to a number of covariates. A potential dependence of other moments of the response distribution, however, is neglected.

Generalized additive models for location, scale and shape (GAMLSS; see Rigby and Stasinopoulos, 2005 and Stasinopoulos et al., 2017) and its Bayesian version of distributional regression (see Klein et al., 2015) provide a more flexible framework. On the one hand, it is no longer restricted to the exponential family and on the other hand, it allows for modelling each parameter of the response distribution—and not only the mean—in dependence of a set of covariates. An alternative to distributional regression is quantile regression. A comparison between the two approaches with particular focus on real estate data is given in Razen et al. (2017). Here, distributional regression turned out to be favourable compared to quantile regression.

In many applications, the data consists of a number of different clusters. In general, it is not guaranteed that a covariate′s effect on a parameter of the response distribution—be it its mean or another parameter—is homogeneous across these clusters. In real estate data, for example, a frequently observed phenomenon is that the price effects of covariates vary from one spatial unit to another. However, completely different functional forms in each cluster are not common.

In order to deal with this challenge, Wechselberger et al. (2008) suggest the use of cluster-specific random scaling factors. In doing so, one still assumes homogeneity for the functional form of the response function but allows for heterogeneity with respect to its scaling. Lang et al. (2015) and Weber et al. (2015) have successfully applied this approach to store sales models, considerably improving the predictive validity of the models. In a real estate context, Brunauer et al. (2010) introduced spatial scaling factors to account for district-specific heterogeneity in rent prices.

Nevertheless, the present method has several drawbacks. First, the analyses so far are restricted to modelling the mean in Gaussian response distributions, ignoring both alternative response distributions and covariate effects on other moments. Second, due to identifiability reasons, one currently has to assume monotonicity for the response functions of the covariates, which in many applications is not justified a priori. Third, the present method only works for independent scaling factors, ignoring feasible correlations that particularly can be observed in applications with spatial clustering.

The aim of this article is to provide a general framework that allows the use of random scaling factors for arbitrary response functions in applications that go beyond the modelling of the response distribution′s mean. For this purpose, we extend the idea of random scaling factors to distributional regression models and develop identifiability constraints that no longer restrict the response functions to be monotone. Furthermore, we provide an alternative estimation technique for spatially correlated scaling factors based on Markov random fields, see, for example, Besag et al. (1991), Fahrmeir and Lang (2001) or Rue and Held (2005).

We investigate the properties of our approach in comprehensive simulation scenarios for different response distributions. We then apply the method to a German real estate dataset with almost 100 000 observations. We consider and compare different distributional regression models where house-specific attributes flexibly are estimated using P-splines that at the same time are scaled according to correlated, district-specific random scaling factors. The results are compared to standard models without scaling factors.

The remainder of the article is structured as follows: In Section 2, we present the methodology that subsequently is tested in different simulation scenarios in Section 3. Section 4 attends to the real estate data and the model specification before we present the results in Section 5. We conclude in Section 6 with an outlook on future research perspectives.

2 Methodology

2.1 Distributional regression models

Suppose we are given data on n observations in the form $(y_{i}, z_{i}, x_{i})$ , $i = 1, \dots, n$ , with response $y$ and a number of covariates $z$ and $x$ . Bayesian distributional regression as introduced in Klein et al. (2015) now assumes an $L$ -parametric distribution of the response $y$ , given the covariates, and links its parameters $θ_{1}, \dots, θ_{L}$ to structured additive predictors $η_{l}$ via known response functions $h_{l}$ , $θ_{l} = h_{l} (η_{l})$ , $l = 1, \dots, L$ . The predictors are defined in terms of (potentially different) subsets of the covariates,

\begin{matrix} η_{l} = f_{1 l} (z_{1 l}) + \dots + f_{q_{l} l} (z_{q_{l} l}) + X_{l} γ_{l}, \end{matrix}

(2.1)

where the functions $f_{jl}$ are possibly non-linear functions of the covariates $z_{jl}$ and the term $X_{l} γ_{l}$ comprises the linear effects of the model. For the sake of simplicity, we will suppress the index discriminating between the $L$ parameters in the following whenever possible.

Using known basis functions $B_{k}$ , a particular function $f$ can be approximated by

f (z) = \sum_{k = 1}^{K} β_{k} B_{k} (z),

where $β = {(β_{1}, \dots, β_{K})}^{'}$ is a vector of unknown regression coefficients to be estimated. A standard choice for continuous covariates are B-spline basis functions, see further.

Defining the $n \times K$ design matrix $Z$ with elements $Z [i, k] = B_{k} (z_{i})$ , the vector $f = {(f (z_{1}), \dots, f (z_{n}))}^{'}$ of function evaluations can be written in matrix notation as $f = Z β$ . Accordingly, the predictors in (2.1) can be written as

\begin{matrix} η = Z_{1} β_{1} + \dots + Z_{q} β_{q} + X γ . \end{matrix}

In a Bayesian framework, overfitting of a particular function $f$ is usually avoided by employing a suitable smoothness prior for the regression coefficients $β$ , see, for example, Fahrmeir et al. (2013). A standard choice is a (possibly improper) Gaussian prior of the form

\begin{matrix} p (β | τ^{2}) \propto {(\frac{1}{τ^{2}})}^{rk (K) / 2} \exp (- \frac{1}{2 τ^{2}} β^{'} K β) \cdot I (A β = 0), \end{matrix}

(2.2)

where $I (\cdot)$ is the indicator function. The key components of the prior are the penalty matrix $K$ , the variance parameter $τ^{2}$ and the constraint $A β = 0$ . Usually the penalty matrix is rank deficient, that is, $rk (K) < K$ , resulting in a partially improper prior. The specific structure of $K$ depends on the covariate type and on prior assumptions about the smoothness of $f$ . The term $I (A β = 0)$ imposes required identifiability constraints on the parameter vector. As a standard, we use $a_{1 k} = \sum_{i = 1}^{n} B_{k} (z_{i})$ , $k = 1, \dots, K$ resulting in the constraint $\sum_{i = 1}^{n} f (z_{i}) = 0$ .

We apply, for example, a Bayesian version of P-splines while modelling a smooth function $f$ that depends on a continuous covariate $z$ , see Eilers and Marx (1996), Lang and Brezger (2004) and Eilers et al. (2015). Here, the columns of the design matrix $Z$ are given by B-spline basis functions evaluated at the observations $z_{i}$ and we use first or second order random walks as smoothness priors for the regression coefficients, that is, $β_{k} = β_{k - 1} + u_{k}, or β_{k} = 2 β_{k - 1} - β_{k - 2} + u_{k},$ with Gaussian errors $u_{k} \sim N (0, τ^{2})$ and diffuse priors $p (β_{1}) \propto const$ , or $p (β_{1})$ and $p (β_{2}) \propto const$ , for initial values. This prior is of the form (2.2) with penalty matrix given by $K = D^{'} D$ , where $D$ is a first or second order difference matrix.

The amount of smoothness is governed by the variance parameter $τ^{2}$ . A conjugate inverse-gamma prior is employed for $τ^{2}$ , that is, $τ^{2} \sim IG (a, b)$ with small values for the hyperparameters $a$ and $b$ resulting in an uninformative prior on the log scale. As a default, we choose $a = b = 0.001$ .

2.2 Multiplicative random effects

As outlined in the introduction, in many applications, the data is clustered. Real estate data, for example, typically is clustered in spatial units (e.g., districts, states, etc.). Usually, there is no economic reason to assume homogeneous covariate effects across these units. In contrast, different consumer price sensitivities originating from varying levels of income, diverse value of land or different ways of construction suggest spatial heterogeneity in price response. Indeed, it is reasonable to assume the effects to have the same functional form but to vary with respect to the scaling of the function. Thus, in order to account for this kind of heterogeneity, we allow for cluster-specific random scaling factors for some or all of the non-linear functions $f_{j}$ in (2.1). This leads to predictors of the form

\begin{matrix} η_{i} = (1 + α_{1 c_{i}}) f_{1} (z_{1 i}) + \dots + (1 + α_{{qc}_{i}}) f_{q} (z_{qi}) + x_{i}^{'} γ, \end{matrix}

(2.3)

where $i = 1, \dots n$ , $c_{i} \in {1, \dots, C}$ is the cluster index of the respective observation and $(1 + α_{{jc}_{i}}), j = 1, \dots, q$ are the corresponding scaling factors. Obviously, a positive random effect $α_{{jc}_{i}} > 0$ leads to a scaling up of the function $f_{j}$ , indicating an increased price sensitivity, while a negative random effect $α_{{jc}_{i}} < 0$ refers to weaker price sensitivity.

2.2.1 Independent scaling factors

In the most basic case of independent scaling factors, the random effects $α_{jc}$ , $j = 1, \dots, q$ , are assumed to be independent and normally distributed with mean $0$ and variance ${\tilde{τ}}_{j}^{2}$ , that is,

\begin{matrix} α_{jc} | {\tilde{τ}}_{j}^{2} \sim N (0, {\tilde{τ}}_{j}^{2}), c = 1, \dots, C . \end{matrix}

For the variance parameters ${\tilde{τ}}_{j}^{2}$ usual inverse-gamma priors ${\tilde{τ}}_{j}^{2} \sim IG ({\tilde{a}}_{j}, {\tilde{b}}_{j})$ are assigned.

A priori, the parameters are not identifiable since there is an arbitrary multiplicative constant for the functions $f_{j}$ . Thus, previous works (e.g., Lang et al., 2015, or Weber et al., 2015) assumed the response functions to be monotone and restricted their spread by assuming

\sum_{k = 1}^{K} β_{j k}^{2} = d_{j} .

Typically, the constant $d_{j}$ is chosen such that the squared sum of the coefficients is identical to that of the corresponding model without scaling factors.

Instead of making assumptions about the coefficients $β_{jk}$ , we propose to assume

\sum_{c = 1}^{C} α_{j c} = 0.

This assumption is preferable for at least two reasons: First, we no longer need monotonicity constraints for the response functions $f_{j}$ . Second, the unscaled functions now can be interpreted as the average effect over all clusters.

Then, the prior for the random effects is of the form (2.2) too.

2.2.2 Spatially correlated scaling factors

If we assume the scaling factors to be spatially correlated, we instead propose the use of a Markov random field prior for the random effects $α_{jc}$ in (2.3):

α_{j c} | α_{j r}, {\tilde{τ}}^{2}, r \in N (c) ~ N (\frac{1}{| N (c) |} \sum_{r \in N (c)} α_{j r}, \frac{{\tilde{τ}}^{2}}{| N (c) |}),

where $| N (c) |$ is the number of the neighbours $N (c)$ of region $c$ . The joint distribution of all random effects $α_{j} = (α_{j 1}, \dots, α_{jC})^{'}$ then is given by

p (α_{j} | {\tilde{τ}}_{j}^{2}) \propto {(\frac{1}{{\tilde{τ}}_{j}^{2}})}^{(C - 1) / 2} \exp (- \frac{1}{2 {\tilde{τ}}_{j}^{2}} α_{j}^{'} {\tilde{K}}_{j} α_{j}),

with the precision matrix ${\tilde{K}}_{j}$ being defined by

{\tilde{K}}_{j} [c, r] = {\begin{array}{l} - 1 & if c \neq r and r \in N (c), \\ 0 & if c \neq r and r \notin N (c), \\ | N (c) | & if c = r . \end{array}

(2.4)

For the variance parameter ${\tilde{τ}}_{j}^{2}$ , we assign the usual inverse-gamma prior ${\tilde{τ}}_{j}^{2} \sim IG ({\tilde{a}}_{j}, {\tilde{b}}_{j})$ and again centre the random effects around zero,

\sum_{c = 1}^{C} α_{j c} = 0.

Then, the prior for the correlated random effects again is a smoothing prior of the form (2.2).

2.3 Inference

For the sake of illustration, we first consider a Gaussian model with a single predictor for the mean parameter of the form (2.3).

The description of posterior inference is facilitated by rewriting the model equation in matrix notation. We obtain

y = D_{1} Z_{1} β_{1} + \dots + D_{q} Z_{q} β_{q} + X γ + ε,

(2.5)

where $Z_{j}$ , $j = 1, \dots, q$ is the usual design matrix for the $j$ th non-parametric term (e.g., P-spline), $D_{j} = diag (1 + α_{{jc}_{1}}, \dots, 1 + α_{{jc}_{n}})$ is an $n \times n$ diagonal matrix with the random scaling factors in the main diagonal and $β_{j}$ is the $j$ th vector of regression coefficients.

Each of the $q$ terms $D_{j} Z_{j} β_{j}$ is formally in the form of a varying coefficient term (Hastie and Tibshirani, 1993) with effect modifier matrix $Z_{j}$ and (pseudo) values of the interacting variable stored in the diagonal matrix $D_{j}$ .

An alternative formulation in terms of the scaling parameter vectors $α_{j} = (α_{j 1}, \dots, α_{jC})^{'}$ is given by

y = \dots + f_{j} + \tilde{D} j {\tilde{Z}}_{j} α_{j} + \dots,

(2.6)

where $f_{j}$ is a vector of function evaluations at the observed covariate values, ${\tilde{D}}_{j} = d i a g (f (z_{j 1}), \dots, f (z_{j n}))$ is a n n diagonal matrix with diagonal elements now given by the function evaluations of the non-linear effects, and ${\tilde{Z}}_{j}$ is a n C matrix indicating if observation i belongs to cluster c (in this case, ${\tilde{Z}}_{j} (i, c) = 1$ , otherwise it equals 0). Again, the expressions ${\tilde{D}}_{j} {\tilde{Z}}_{j} α_{j}$ in (2.6) are in the form of a varying coefficient term now with effect modifier matrix ${\tilde{Z}}_{j}$ and values of the interacting variable given in $\tilde{D} j$ .

The varying coefficient representation (2.6) suggests the following (heuristic) two-stage estimation procedure:

Assume the covariate effects to be homogeneous over all clusters, that is, set the random effects $α_{j}$ to zero, and estimate the model $y = f_{1} + \dots + f_{q} + X γ + ε$ as usual.

Treat the estimated functions ${\hat{f}}_{j}$ from stage 1 as a fixed offset and estimate the random effects from the varying coefficient representation (2.6).

Conceptually, this two-stage procedure works fine. However, simulations show that ignoring the random effects in stage $1$ can raise difficulties in properly estimating the functions $f_{j}$ , particularly if the variances of the random effects are large, see Section 5.2 for details. Inappropriate estimations of the functions $f_{j}$ may then lead to peculiar results for the random effects in stage $2$ . Thus, instead of this two-stage procedure, we propose a simultaneous estimation approach, where Gibbs updates for the random scaling terms are employed by alternately obeying the two different varying coefficient representations (2.5) and (2.6).

The full conditionals of the regression parameters $β_{j}$ are easily derived from (2.5) and are multivariate Gaussian $β_{j} | \cdot \sim N (μ_{j}, Σ_{j})$ with

Σ_{j}^{- 1} = \frac{1}{σ^{2}} (Z_{j}^{'} D_{j}^{2} Z_{j} + \frac{σ^{2}}{τ_{j}^{2}} K_{j}), Σ_{j}^{- 1} μ_{j} = \frac{1}{σ^{2}} Z_{j}^{'} D_{j} (y - η_{j}),

(2.7)

where $η_{j}$ contains the current predictor except the $j$ th term.

The full conditionals of the scaling factors are derived from (2.6) with $α_{j} | . ~ N ({\tilde{μ}}_{j}, {\tilde{Σ}}_{j})$ and

{\tilde{Σ}}_{j}^{- 1} = \frac{1}{σ^{2}} ({\tilde{Z}}_{j}^{'} {\tilde{D}}_{j}^{2} {\tilde{Z}}_{j} + \frac{σ^{2}}{{\tilde{τ}}_{j}^{2}} {\tilde{K}}_{j}), {\tilde{Σ}}_{j}^{- 1} {\tilde{μ}}_{j} = \frac{1}{σ^{2}} {\tilde{Z}}_{j}^{'} {\tilde{D}}_{j} (y - f_{j} - η_{j}) .

(2.8)

Here, the penalty matrix ${\tilde{K}}_{j}$ is either the identity matrix in the case of independent scaling factors or the precision matrix defined in (2.4) in the case of spatially correlated scaling factors.

In contrast to ‘usual’ varying coefficients terms, the diagonal matrices $D_{j}$ and ${\tilde{D}}_{j}$ of the (pseudo) interacting variables are not constant during the Gibbs sampler. Hence, cross products and other quantities can not be computed and stored in advance and numerical efficient updating is considerably complicated. However, the methodology for highly efficient Gibbs updates described in Lang et al. (2014), can be applied.

Finally, the full conditionals of the variance parameters are inverse-gamma and are given by

\begin{array}{l} τ_{j}^{2} ~ I G (a', b'), a' = a + 0.5 rk (K_{j}), b' = b + 0.5 β_{j}^{'} K_{j} β_{j}, \\ {\tilde{τ}}_{j}^{2} ~ I G (\tilde{a}', \tilde{b}'), \tilde{a}' = \tilde{a} + 0.5 rk ({\tilde{K}}_{j}), \tilde{b}' = \tilde{b} + 0.5 α_{j}^{'} {\tilde{K}}_{j} α_{j} . \end{array}

(2.9)

We finally arrive at the following:

Algorithm: Update of a particular term $f = D Z β$

Update an unconstrained version of $β$ for given $α$ by drawing from full conditional (2.7). Correct the unconstrained vector (at negligible computational cost) using steps 5–9 of Algorithm 2.6 in Rue and Held (2005).

Update $α$ for given $β$ by drawing from full conditional (2.8). Correct the unconstrained vector again using Algorithm 2.6 of Rue and Held (2005).

Update the variance parameters $τ_{j}^{2}$ and ${\tilde{τ}}_{j}^{2}$ by drawing from the inverse-gamma full conditionals (2.9).

For most updates of regression parameters $β_{j}$ and $α_{j}$ in distributional regression, Gibbs steps are not available because the full conditionals are no longer Gaussian. Then we rely on the Metropolis–Hastings updates with iteratively weighted least squares (IWLS) proposals as described in detail in Klein et al. (2015) and Klein et al. (2014). The key for updating $β_{j}$ and $α_{j}$ is again the varying coefficient type notation of the random scaling terms given in (2.5) and (2.6).

The basic idea of IWLS proposals is to determine an approximation of the full conditional that leads to a Gaussian proposal density with expectation and covariance matrix corresponding to the mode and the curvature at the mode at the current state of the Markov chain.

The resulting Metropolis–Hastings sampler has the advantage that the IWLS proposal automatically adapts to the form of the full conditional and thereby avoids manual tuning. Let $l (η)$ be the log-likelihood depending on the predictor $η$ . Then the full conditional for a parameter block $β_{j}$ is

\log (p (β_{j} | \cdot)) \propto l (η) - \frac{1}{2 τ_{j}^{2}} β_{j}^{'} K_{j} β_{j}

Taylor expansion leads to a Newton′s method type approximation and we can deduce a Gaussian working model, see Klein et al. (2015) for more details. We finally obtain the Gaussian proposal densities $N (μ_{j}, P_{j}^{- 1})$ for $β_{j}$ with expectation and precision matrix

μ_{j} = P_{j}^{- 1} Z_{j}^{'} D_{j} W (z - η_{j}) P_{j} = Z_{j}^{'} D_{j} {W D}_{j} Z_{j} + \frac{1}{τ_{j}^{2}} K_{j}

(2.10)

where $z = η + W^{- 1} \frac{\partial l}{\partial η}$ and $W = diag (w_{1}, \dots, w_{n})$ , $w_{i} = E (- \frac{\partial^{2} l}{\partial^{2} η_{i}})$ are the usual working observations and weights in GLMs. The parameters for the Gaussian proposal densities $N (̃ {\tilde{μ}}_{j}, ̃ {\tilde{P}}_{j}^{- 1})$ for $α_{j}$ are given by

{\tilde{μ}}_{j} = {\tilde{P}}_{j}^{- 1} {\tilde{Z}}_{j}^{'} {\tilde{D}}_{j} W (z - f_{j} - η_{j}) {\tilde{P}}_{j} = Z_{j}^{'} {\tilde{D}}_{j} W {\tilde{D}}_{j} Z_{j} + \frac{1}{{\tilde{τ}}_{j}^{2}} {\tilde{K}}_{j}

(2.11)

Finally, we replace the Gibbs steps of the previous algorithm for updating a particular function $f = D Z β$ by MH steps with proposals given by (2.10) and (2.11).

3 Simulation

3.1 Simulation design

An extensive simulation study for four distributions (Binomial, Gaussian, gamma and Student′s-t) has been carried out. To keep the article in reasonable length, we present exemplary results for models with gamma responses.

The parameters $μ$ and $σ$ of the gamma distribution are linked to separate predictors

η_{1} = (1 + α_{c 1}) f_{1} (z), η_{2} = (1 + α_{c 2}) f_{2} (z),

via the exponential link function. Here, $f_{1} (z) = f_{2} (z) = \sin (z)$ with $z$ being randomly drawn from a uniform distribution on the interval $[- π, π]$ . These non-linear functions are modified by random scaling factors $(1 + α_{cl})$ , $l = 1, 2$ , for the districts of Western Germany. The random effects $α_{cl}$ have mean $0$ and variance ${\tilde{τ}}_{l}^{2}$ and are either assumed to be independent and normally distributed (which we will call uncorrelated setting) or they are assumed to be spatially correlated (correlated setting). In the latter case, they are defined by the normalized centroid coordinates $x_{c}$ and $y_{c}$ of the individual districts $c$ as

α_{c 1} = 0.4 \cdot (x_{c} + y_{c}), α_{c 2} = 2.5 \cdot x_{c} + 0.5 \cdot y_{c} .

We additionally scale the random effects $α_{cl}$ in three different ways, so that their variance ${\tilde{τ}}_{l}^{2}$ is either $0 . 1^{2}$ , $0 . 5^{2}$ or $1 . 0^{2}$ , in order to evaluate the influence of ${\tilde{τ}}_{l}^{2}$ . Furthermore, we vary the number of observations per district. In total, we analyse five different models whose specifications are summarized in Table 1.

Table 1:
Model specifications

Model Observation per district Variance of $α_{cl}$

Model 1 10 $0 . 5^{2}$

Model 2 25 $0 . 5^{2}$

Model 3 50 $0 . 5^{2}$

Model 4 25 $0 . 1^{2}$

Model 5 25 $1 . 0^{2}$

Model	Observation per district	Variance of $α_{cl}$
Model 1	10	$0 . 5^{2}$
Model 2	25	$0 . 5^{2}$
Model 3	50	$0 . 5^{2}$
Model 4	25	$0 . 1^{2}$
Model 5	25	$1 . 0^{2}$

3.2 Results

For each model, we generate $250$ replications and carry out the respective estimation procedure for uncorrelated or correlated random scaling factors described in Section 2 based on a final Markov chain Monte Carlo (MCMC) run with $120 000$ iterations and a burn in period of $20 000$ iterations. We store every $100$ th iteration in order to obtain a sample of $1 000$ (uncorrelated) draws from the posterior. We then calculate the arithmetic mean from the $250$ replications.

Figures 1 and 2 show the average estimates of the effects $f_{l} (z)$ as well as of the district-specific effects $(1 + α_{cl}) f_{l} (z)$ for the smallest and largest random effects $α_{cl}$ (solid) of the gamma models 1, 3, 4 and 5. The true effects also are plotted (dashed) in order to facilitate comparison. The first column refers to the uncorrelated setting, the second column to the correlated setting. Figure 1 shows the results for the $μ$ -parameter, Figure 2 those for the $σ$ -parameter.

For both parameters, the scaled effects almost perfectly can be estimated in the correlated settings even with only $10$ observations per district. In the uncorrelated setting, we are still able to estimate the scaled effects, but here we need much more observations per district to get adequate results, especially for the $σ$ -parameter. If the variance of the random effects is too small, the effects can hardly be separated from the overall noise in the data in the uncorrelated setting. In the correlated setting, however, we still get reasonable results even for a small variance. If the variance gets larger, the estimation results get better in both settings. Thus, the estimates are more biased if the random effect is weak, which is a well-known feature of random effects estimators, see, for example, Gelman and Hill (2007).

Figure 1:

The average estimates of the functions $f_{1} (z)$ as well as of the smallest and largest district-specific effects $(1 + α_{c 1}) f_{1} (z)$ (solid) and the respective true effects (dashed) for the gamma models $1$ , $3$ , $4$ and $5$ . The first column shows the effects in the uncorrelated setting and the second column those in the correlated setting

Figure 2:

The average estimates of the functions $f_{2} (z)$ as well as of the smallest and largest district-specific effects $(1 + α_{c 2}) f_{2} (z)$ (solid) and the respective true effects (dashed) for the gamma models $1$ , $3$ , $4$ and $5$ . The first column shows the effects in the uncorrelated setting and the second column those in the correlated setting

For illustration, Figure 3 shows the geographical maps of the estimated random effects $\hat{α_{c 1}}$ for the $μ$ -parameter of the gamma model 2. The first row refers to the uncorrelated setting, the second row to the correlated setting. For comparison, the true effects $α_{c 1}$ are plotted in the middle column. We also depict the biases, defined as $\hat{α_{c 1}} - α_{c 1}$ , in the right column. In both settings, the biases are rather small, reflecting the good performance of our methods in identifying both independent and correlated random scaling factors.

Figure 3:

Maps of the random effects $α_{c 1}$ of the gamma model 2 for the uncorrelated (first row) and the correlated setting (second row). The first column shows the estimated effects, the second column shows the true effects and the third column shows the biases

In order to evaluate the performance of our simultaneous estimation approach, we refer to mean squared errors (MSEs). In each replication of our models, we calculate the MSE for the $l$ different parameters as follows:

{MSE}_{l} = \frac{1}{n} {(\hat{η_{l}} - η_{l})}^{'} (\hat{η_{l}} - η_{l}) .

We then reestimate all replications using the two-stage estimation procedure sketched in Section 2.3 and again calculate the corresponding MSEs.

For illustration, Figure 4 shows boxplots of the MSEs of the 250 replications in the Gamma models 1–3 in the correlated setting. As we can see, the simultaneous estimation approach (1) yields lower MSEs than the two-stage procedure (2) for all these models. Especially for the $σ$ -parameter, the estimation results are considerably better with our simultaneous approach.

Figure 4:

MSEs for the $μ$ - and $σ$ -parameter of the gamma models in the correlated setting with ${\tilde{τ}}_{l}^{2} = 0 . 5^{2}$ from the simultaneous estimation approach (‘1’) and the two-stage procedure (‘2’)

With respect to the variance of the random scaling factors, we find that the performance of the simultaneous estimation approach compared to the two-stage procedure substantially gets better the higher the variance is (results not shown).

The same findings hold for the gamma models in the uncorrelated seeting. The main reason is that ignoring the scaling factors in the first step of the two-stage procedure may lead to peculiar results for the average effect, in particular if the scaling factors have a large variance. Then, the scaled effects estimated in the second step obviously are biased as well.

4 Data description and model specification

We apply our methodology to a dataset of almost $100 000$ single family homes all over Germany. The data was provided by F+B Research and Consulting for Habitation, Real Estate and Environment Ltd, a business consultancy in Hamburg, Germany.

Since the raw data initially was supply data, prices obviously were upward biased. Thus, F+B adjusted these prices by a transaction discount estimated from a regression model in order to provide realistic purchase prices. Dividing these adjusted prices by the floor area of the houses leads to the prices per square meter ( $p_{qm}$ ), which we use as the dependent variable in our analysis. The set of explanatory variables include continuous and categorical covariates that characterize the building and its location:

Continuous covariates: The floor area of the building ( $area$ ) is expected to have a decreasing effect on the price per square meter due to the law of diminishing marginal utility, while the plot area where the house is built on ( $plot_a rea$ ) should have a positive effect. Due to depreciation over time, the year of construction ( $year$ ) in general should also have an increasing effect on house prices per square meter. Besides from these structural covariates, an expert rating ( $rating$ ) is included that characterizes the area in which the building is situated. Since the rating ranges from $1$ (excellent) to $9$ (bad), we expect a negative effect on the prices.

Categorical covariates: The equipment of the house ( $equipment$ ) is classified by four categories. Obviously, we expect an increasing effect for the better equipment. In order to control for state-specific price differences, we include dummy variables for the states in which the buildings are located.

The standard model for explaining house prices is a log-Gaussian model with location parameter $μ$ and scale parameter $σ$ . We link these parameters to predictors $η_{1}$ and $η_{2}$ via $μ = η_{1}$ and $σ = \exp (η_{2})$ . As an alternative, we additionally consider a gamma model with mean parameter $μ$ and shape parameter $σ$ , which we link to the predictors via $μ = \exp (η_{1})$ and $σ = \exp (η_{2})$ .

In both models, the predictors $η_{}$ , $l = 1, 2$ , are constructed as follows:

\begin{matrix} \begin{matrix} η_{} & = & D_{1 l} f_{1 l} (area) + D_{2 l} f_{2 l} (plot_a rea) \\ + D_{3 l} f_{3 l} (year) + f_{4 l} (rating) + {X γ}_{l} . \end{matrix} \end{matrix}

(4.1)

$f_{1 l}, \dots, f_{4 l}$ are possibly non-linear functions of the continuous covariates and will be modelled with P-splines. $D_{jl} = diag (1 + α_{{jc}_{1} l}, \dots, 1 + α_{{jc}_{n} l})$ , $c_{i} \in {1, \dots, C}$ , contains random scaling factors for the $C$ districts that allow for regional heterogeneity in the respective price response functions. Here, we particularly expect spatial differences between Eastern and Western Germany. A random intercept as well as the dummy variables for the equipment of the houses and the states where the buildings are located in are subsumed in the design matrix $X$ with parameters $γ_{l}$ .

Since we assume the random scaling factors to be spatially correlated, we apply the respective estimation procedure based on Markov random field priors described in Section 2.

5 Results

The estimation results for the models described in the previous section are based on a final MCMC run with $270 000$ iterations and a burn in period of $20 000$ iterations. We stored every $250$ th iteration, leading to a sample of $1 000$ practically independent draws from the posterior.

When analysing the results, one has to be aware that the parameters of the distributions do not exactly correspond to each other. For example, the mean of the gamma model immediately is given by the respective $μ$ -parameter, while it is given by $\exp (μ + σ^{2} / 2)$ in the log-Gaussian model. Therefore, a simple comparison of the parameters (or of the involved scaling factors) is not reasonable. Instead, we always derive the mean effects of the two models and compare them to each other. However, we additionally present the results of the individual scaling factors at least for the best model (the gamma model, see Subsection 5.2).

5.1 Effect estimates

For the sake of illustration, we restrict the presentation of the results to the covariates being multiplied by random scaling factors, that is, the floor area, the plot area and the year of construction.

5.1.1 Mean effects

As expected, the mean effect of the floor area on house prices per square meter, averaged over all districts, is monotonically decreasing in both models (see first row of Figure 5, Panel [a]). Here, in order to get an impression of the magnitude of the effects and to make the results comparable, the other continuous covariates are held constant at mean level of attributes and the categorical variables are held at their mode level (which we will call the mean level). As we can see, the results of the two models almost coincide.

With respect to the scaling, we find that the effect of the floor area considerably differs between the districts. Panels (b) and (c) in the first row of Figure 5 show the scaled effects of the different districts for the two models together with the respective average effect. In the gamma model, for example, the effect of the floor area accounts for a variation between $300$ euro per square meter in the district with the smallest scaling factor and $950$ euro per square meter in the district with the largest scaling factor.

Regarding the geographic distribution, the scaling factors of the mean parameter tend to be higher in Western Germany, while they are considerably lower in Eastern Germany (Panel [a] of Figure 6). In order to verify the significance of these differences, we refer to the posterior probabilities based on a nominal level of $80$ %. If the $80$ % credible interval of a random effect $α$ is strictly positive or negative, we assign the corresponding district a value of $1$ or $- 1$ , respectively. If the credible interval contains zero, we record a value of $0$ for this district. As Panel (d) of Figure 6 shows, the scaled effects of the floor area significantly differ from the average effect in about one third of the districts.

For the plot area, the average mean effect over all districts is monotonically increasing up to an area of $1 900$ square meters and then slightly reverses for both models (Plot [a] of second row in Figure 5). However, this minor decrease for very large plots is not significant according to the simultaneous $95$ % confidence bands (not depicted in the figure). There is again considerable variation in the effects for the different districts (second row of Figure 5, Panels [b] and [c]). In the gamma model, for example, the respective scaling factors of the mean parameter range from $0.06$ to $1.86$ .

According to these scaling factors, the magnitude of the mean effect of the plot area is much more pronounced in the Southwestern states than in the remaining parts of Germany (Panel [b] of Figure 6). This coincides with the regions that have a higher population density and where land therefore is a scarcer resource. The deviation from the average effect is significant in almost half of the districts, see Panel (e) of Figure 6.

In general, the results for the year of construction show an increasing effect on the mean of the house prices per square meter (third row of Figure 5, Panel[a]). However, there are three periods where the effect is negative: During and shortly after the First World War, construction naturally was cheaper. The same holds for the time of the Second world war. The third period showing a negative effect is the time after the latest financial crisis when the real estate market in Germany was in trouble. Indeed, a recovery started in recent years, but this is not yet covered by our data ending in $2012$ . Again, the scaling factors, which for the gamma model vary between $0.16$ and $1.70$ , reveal considerable spatial heterogeneity in the magnitude of the effects (third row of Figure 5, Panels [b] and [c]).

With the exception of the surrounding of Berlin, the mean effect of the year of construction is more pronounced in Eastern Germany than in most parts of Western Germany, see Panel (c) of Figure 6. The main reason is that the construction industry in Eastern Germany was much more affected by the reunification of $1990$ , causing at first a huge boom in the Eastern states followed by a severe downturn afterwards. Thus, the year of construction has more influence on house prices in Eastern than in Western Germany. Overall, the scaled effects significantly differ from the average effect in almost three-fourth of the districts (Panel [e] of Figure 6).

Figure 5:

Mean effect of the floor area, plot area and year of construction evaluated at the mean level. [a]: Average effect over all districts for the two different models. [b], [c]: Scaled effects of the individual districts for the two models together with the respective average effect

Figure 6:

Random scaling factors and $80$ % posterior probabilities of the floor area (Panels [a] and [d]), plot area (Panels [b] and [e]) and year of construction (Panels [c] and [f]) for the mean parameter in the gamma model

5.1.2 Shape parameters

For the shape parameter $σ$ , the marginal effect of the floor area, averaged over all districts, is increasing up to an area of about $120$ square meters and then levels off (see Figure 7, Panel [a]). However, as can be seen from the simultaneous $95$ % confidence bands, the effect is not significant for a wide range of the floor area. Accordingly, the majority of the respective scaling factors are not significant. Only in the Southern part of Bavaria, the scaling factors are consistently lower than in the rest of Germany (see Figure 8, Panels [a] and [c]).

For the shape parameter of the gamma model, the average marginal effect of the plot area is slightly decreasing, see Panel (b) of Figure 7. Again, there are considerable differences in the scaling of the effect over the districts with the corresponding scaling factors ranging from $- 1.96$ to $4.01$ . The largest scaling factors can be found in the Southern and Western states of Germany, the lowest scaling factors are located in Central Germany (Panels [b] and [d] of Figure 8). In Eastern Germany, there are hardly any districts where the scaled effects significantly differ from the average effect.

Except for the last few years, the effect of the year of construction for the shape parameter in the gamma model looks similar to the one for the mean parameter (Panel [c] of Figure 7). The corresponding scaling factors range from $- 0.24$ to $3.16$ .

The effect is more pronounced in Central and Eastern Germany than in the Western parts, see Panels (c) and (e) of Figure 8, with significant differences from the average effect in about half of the districts.

Figure 7:

Average marginal effect of the floor area, plot area and year of construction for the shape parameter in the gamma model together with simultaneous $95$ % confidence bands

Figure 8:

Random scaling factors and corresponding posterior probabilities of the floor area (Panels [a] and [d]), plot area (Panels [b] and [e]) and year of construction (Panels [c] and [f]) for the shape parameter in the gamma model

5.2 Model comparison

In order to evaluate the performance of the models, we refer to the well known deviance information criterion (DIC) of Spiegelhalter et al. (2002), see also Klein et al. (2015), who have shown that the DIC also can be used to discriminate between different types of response distributions in distributional regression. As we can see from Table 2, the DIC of the models with scaling factors are by far lower than the DIC of models without scaling factors, and the gamma distribution clearly outperforms the log-Gaussian model. The findings based on the DIC are further confirmed by proper scoring rules as proposed originally by Gneiting and Raftery (2007) and transferred to distributional regression in Klein et al. (2015). We used the logarithmic score, quadratic score, spherical score and the continuous ranked probability score which again clearly favour the gamma model (results not shown).

Table 2:
DIC of the standard models and the extended models with random scaling factors

DIC Log-Gaussian Gamma

Standard model $50 882$ $50 556$

Model with scaling factors $45 303$ $45 081$

DIC	Log-Gaussian	Gamma
Standard model	$50 882$	$50 556$
Model with scaling factors	$45 303$	$45 081$

6 Conclusion

This article presents a simultaneous estimation approach for Bayesian distributional regression models with either independent or spatially correlated random scaling factors and provides a comprehensive simulation study showing the accuracy of this approach. A comparison to an ordinary two-stage estimation procedure reveals considerable improvement particularly for models where the variance of the random scaling factors is large.

We apply our methodology to real estate data from Germany and identify district-specific random scaling factors that are significant in up to three-fourth of the districts. The results confirm expected spatial heterogeneity in the covariates’ effects and provide new insights into the valuation of house prices. Furthermore, allowing for district-specific random scaling factors significantly improves the performance of the models with respect to their DIC.

Several directions for further research are conceivable. First, we are already working on providing other spatial smoothing techniques (e.g. Kriging or tensor product splines) for estimating correlated scaling factors. Second, we aim to implement structured additive predictors for the scaling factors themselves, which would allow us to explain them, for example, by cluster-specific covariates. Third, in the application, we have obtained negative scaling factors for some districts. Notwithstanding that they have occurred only occasionally, negative scaling factors seem rather unlikely in most applications from an economic perspective. Thus, we intend to provide alternative scaling factors that are positive by default, for example by defining them by $\exp (α_{c})$ or by implementing truncated factors. Finally, the choice of the link function is more important in distributional regression then in pure mean regression models. Methodology for simultaneous choice of the link function(s) could be particularly useful.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

This work was supported by funds of the Oesterreichische Nationalbank (Oesterreichische Nationalbank, Anniversary Fund, project number: 15309).

References

Besag

York

Mollie

(1991) Bayesian image restoration with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics , 43, 1–59.

Brezger

Lang

(2006) Generalized structured additive regression based on Bayesian P-splines. Computational Statistics and Data Analysis , 50, 967–91.

Brunauer

Lang

Wechselberger

Bienert

(2010) Additive Hedonic regression models with spatial scaling factors: An application for rents in Vienna. Journal of Real Estate Finance and Economics , 40, 390–410.

Eilers

PHC

Marx

(1996) Flexible Smoothing with B-splines and Penalties. Statistical Science , 11, 89–121.

Eilers

Marx

Durban

(2015) Twenty years of P-splines. Statistics and Operations Research Transactions , 39, 149–86.

Fahrmeir

Lang

(2001) Bayesian inference for generalized additive mixed models based on Markov random field priors. Applied Statistics , 50, 201–20.

Fahrmeir

Kneib

Lang

Marx

(2013) Regression: Models, Methods and Applications. Berlin: Springer.

Gelman

Hill

(2007) Data Analysis Using Regression and Multi-level/Hierarc- hical Models. Cambridge: Cambridge University Press.

Gneiting

Raftery

A E

(2007) Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association , 102, 359–78.

10.

Hastie

T J

Tibshirani

(1990) Generalized Additive Models . London: Chapman & Hall.

11.

Hastie

T J

Tibshirani

(1993) Varying-coefficient models. Journal of the Royal Statistical Society B , 55, 757–96.

12.

Klein

Kneib

Lang

(2014) Bayesian generalized additive models for location, scale and shape for zero-inated and overdispersed count data. Journal of the American Statistical Association , 10, 405–19.

13.

Klein

Kneib

Lang

Sohn

(2015) Bayesian structured additive distributional regression with an application to regional income inequality in Germany. The Annals of Applied Statistics , 9, 1024–52.

14.

Lang

Brezger

(2004) Bayesian P-splines. Journal of Computational and Graphical Statistics , 13, 183–212.

15.

Lang S, Umlauf

Wechselberger

Harttgen

Kneib

(2014) Multilevel structured additive regression. Statistics and Computing , 24, 223–38.

16.

Lang

Steiner

Weber

Wechselberger

(2015) Accommodating heterogeneity and nonlinearity in price effects for predicting brand sales and profits. European Journal of Operational Research , 246, 232–41.

17.

McCullagh

Nelder

(1989) Generalized Linear Models. Boca Raton, FL: Chapman & Hall.

18.

Razen

Brunauer

Klein

Kneib

Lang S Umlauf

(2017) Statistical risk analysis for real estate collateral valuation using Bayesian distributional and quantile regression. Technical report .

19.

Rigby

Stasinopoulos

(2005) Generalized additive models for location, scale and shape. Applied Statistics , 54, 507–54.

20.

Rue

Held

(2005) Gaussian Markov Random Fields. Boca Raton, FL: Chapman & Hall/CRC.

21.

Spiegelhalter

, Best

Carlin

van der Linde

(2002) Bayesian measures of model complexity and t (with discussion). Journal of the Royal Statistical Society Series B , 64, 583–639.

22.

Stasinopoulos

Rigby

Heller

Voudouris

De Bastiani

(2017) Flexible Regression and Smoothing: Using GAMLSS in R. Boca Raton, FL: Champman and Hall.

23.

Weber

Steiner

Lang

(2015) A comparison of semiparametric and heterogeneous store sales models for optimal category pricing OR Spectrum, 39, 403–45.

24.

Wechselberger

Lang

Steiner

(2008) Additive models with random scaling factors: Applications to modeling price response functions. Austrian Journal of Statistics , 37, 255–70.

25.

Wood

(2006) Generalized Additive Models: An Introduction with R. Boca Raton, FL: Chapman & Hall.

Random scaling factors in Bayesian distributional regression models with an application to real estate data

Abstract

Keywords

1 Introduction

2 Methodology

2.1 Distributional regression models

2.2.2 Spatially correlated scaling factors

3.1 Simulation design

Table 1: Model specifications Model Observation per district Variance of α cl Model 1 10 0 . 5 2 Model 2 25 0 . 5 2 Model 3 50 0 . 5 2 Model 4 25 0 . 1 2 Model 5 25 1 . 0 2

Figure 1:

Maps of the random effects α c 1 of the gamma model 2 for the uncorrelated (first row) and the correlated setting (second row). The first column shows the estimated effects, the second column shows the true effects and the third column shows the biases

MSEs for the μ - and σ -parameter of the gamma models in the correlated setting with τ ̃ l 2 = 0 . 5 2 from the simultaneous estimation approach (‘1’) and the two-stage procedure (‘2’)

5.1 Effect estimates

5.1.1 Mean effects

Figure 5:

Mean effect of the floor area, plot area and year of construction evaluated at the mean level. [a]: Average effect over all districts for the two different models. [b], [c]: Scaled effects of the individual districts for the two models together with the respective average effect

Random scaling factors and 80 % posterior probabilities of the floor area (Panels [a] and [d]), plot area (Panels [b] and [e]) and year of construction (Panels [c] and [f]) for the mean parameter in the gamma model

Figure 7:

Average marginal effect of the floor area, plot area and year of construction for the shape parameter in the gamma model together with simultaneous 95 % confidence bands

Random scaling factors and corresponding posterior probabilities of the floor area (Panels [a] and [d]), plot area (Panels [b] and [e]) and year of construction (Panels [c] and [f]) for the shape parameter in the gamma model

Table 2: DIC of the standard models and the extended models with random scaling factors DIC Log-Gaussian Gamma Standard model 50 882 50 556 Model with scaling factors 45 303 45 081

Declaration of Conflicting Interests

Funding

References

Table 1:
Model specifications

Model Observation per district Variance of $α_{cl}$

Model 1 10 $0 . 5^{2}$

Model 2 25 $0 . 5^{2}$

Model 3 50 $0 . 5^{2}$

Model 4 25 $0 . 1^{2}$

Model 5 25 $1 . 0^{2}$

Maps of the random effects $α_{c 1}$ of the gamma model 2 for the uncorrelated (first row) and the correlated setting (second row). The first column shows the estimated effects, the second column shows the true effects and the third column shows the biases

MSEs for the $μ$ - and $σ$ -parameter of the gamma models in the correlated setting with ${\tilde{τ}}_{l}^{2} = 0 . 5^{2}$ from the simultaneous estimation approach (‘1’) and the two-stage procedure (‘2’)

Random scaling factors and $80$ % posterior probabilities of the floor area (Panels [a] and [d]), plot area (Panels [b] and [e]) and year of construction (Panels [c] and [f]) for the mean parameter in the gamma model

Average marginal effect of the floor area, plot area and year of construction for the shape parameter in the gamma model together with simultaneous $95$ % confidence bands

Table 2:
DIC of the standard models and the extended models with random scaling factors

DIC Log-Gaussian Gamma

Standard model $50 882$ $50 556$

Model with scaling factors $45 303$ $45 081$