Abstract

This study employs machine learning techniques (AI), specifically multiple linear regression with nonlinear covariates, to model rental prices per square meter in Berlin, Germany. The research investigates major quantitative and qualitative variables influencing rent dynamics by leveraging a robust dataset comprising over 2.6 million apartments with 59 variables from 2007 to 2020, sourced from FDZ Ruhr and ImmobilienScout24. Drawing from over 99,000 rental records (2015 and 2019 datasets) and 31 variables (9 quantitative and 22 qualitative), the analysis evaluates the influence of factors such as furnishing quality, modernization year, apartment size, and energy efficiency on rent pricing. Polynomial and log transformations were applied to improve model robustness. The use of log-transformed rent as the response variable, combined with nonlinear covariates, yielded the best model performance, with the highest adjusted Rsquared values of 0.3645 and 0.492 in the 2015 and 2019 models, respectively, among the tested models. The results suggest that both quantitative and qualitative variables significantly influence rent sqm, with influential predictors varying in linearity and significance across the two years. In 2015, rent was influenced by nonlinear trends in living space and construction years, while in 2019, heat cost and modernization showed linear increases in rent. Apartments with upscale furnishings, high energy efficiency, elevators, and guest toilets consistently commanded higher rents, whereas a pet allowance was associated with lower rents. Residual analysis, variance inflation factors (VIF), and the Akaike Information Criterion (AIC) confirmed the model’s statistical validity. Results indicate significant shifts in rent patterns and offer predictive insights relevant to policymakers, urban planners, and educational leaders addressing housing affordability and equity. This research builds prior literature in AI-supported urban analytics and contributes a replicable data-driven framework for strategic housing decisions in high-demand university cities.

Keywords

AI Rent Forecasting, Berlin Housing Market, Rental Price Modeling, Multiple Linear Regression, Social Equity, Policy Innovation

Share and Cite:

1. Background and Study Framework

1.1. Introduction

In this paper, artificial intelligence (AI) refers to the application of predictive modeling techniques, particularly regression analysis, within decision systems that use data and digital tools to support urban planning. Advanced statistical methods are used to model the rent per square meter (rent_sqm) in Berlin’s housing market. Berlin’s 2024 housing market exhibits a significant imbalance, with demand outpacing supply. In 2023, the population increased by 0.7% to 3.878 million, driven by 187,971 new arrivals, making it the third-highest year of immigration since 1991 [1]. Rental demand is high in Germany, particularly in Berlin and Munich, where prices have risen significantly, presenting challenges like those faced in other high-income countries, such as the UK, France, the US, and Canada [2]. Germany has a significantly higher share of renters compared to its international counterparts. In 2018, Germany’s homeownership rate stood at 51.5%, while the rates in the UK, Italy, and Romania were 65.1%, 72.4%, and 96.4%, respectively [2]. Although AI often involves complex algorithms, it also encompasses predictive modeling for pattern recognition and analysis. This study uses a transformed multiple linear regression with nonlinear variables to forecast rent prices in Berlin, contributing to urban analysis efforts that support thoughtful policy decisions and promote fairness in housing [3].

1.2. Objective

This paper employs multiple linear regression to model rent_sqm in Berlin, examining major market trends, assessing the need for transformation of the response variable, and identifying significant predictors. The results aim to support the development of housing policy and promote equity, while also assisting educational leaders in strategic and equitable housing planning.

1.3. Literature Review

Skewness and variability in housing data are often addressed using log-linear regression, which highlights furnishing quality, energy efficiency, and modernization as main influencing factors [4] [5]. Germany’s public housing contrasts with the U.S.’s system based on market behavior [6]. Energy-efficient design has garnered increasing attention to sustainability [7]. This study employs AI-driven modeling to analyze Berlin’s rental market and connect statistical insights to educational policy, with a focus on policy innovation and social equity in housing [8] [9].

1.4. Research Questions

• RQ1: Does a relationship exist between the response variable (rent_sqm) and the selected predictor variables?

• RQ2: Is a transformation of response variable (rent_sqm) necessary to meet the assumptions of linear regression?

• RQ3: Which predictor variables significantly affect the rental price per square meter in Berlin’s housing market?

2. Methodology

This section explores the methodology and required mathematical and statistical background for this paper. The study uses multiple regression analysis to model rent trends in Berlin. Although regression is a classical statistical method, its application in this context serves as a robust AI tool (machine learning algorithms) for predictive modeling in housing policy. We look at linear models, their formulations, assumptions, estimations, validation, predictions, and hypothesis testing. More details can be found in [10]. The data analysis was executed using the R programming language.

2.1. Research Design

This research uses a quantitative approach and applies AI-informed predictive modeling, focusing on multiple linear regression with nonlinear covariates, to analyze Berlin’s rental prices.

2.2. Data Collection

A secondary source of data collection was used for this study. The data was provided by the FDZ Ruhr at RWI (and ImmobilienScout24) institution. The Immobilien-Scout24 GmbH, founded in 1998, deals with real estate properties in Germany. The data set contains 2,651,885 observations and 59 attributes from 2007 to 2020. We selected Berlin City as it had the highest number of rental transactions in Germany. Thereafter, we chose the years 2015 and 2019, which had 49,724 and 49,536 records, respectively, based on the significant impact observed in the plotted scatter of years with rent prices as shown in the data description in Table 3.

2.3. Data Analysis: Multiple Linear Regression with Nonlinear Covariates

We cleaned the data and removed outliers using the interquartile range (IQR) method. The missing values recorded and the NAs were part of the labels for most categorical variables, as shown in Table 2. Exploratory data analysis (EDA) was used to visualize the behavior of rent per square meter and its covariates. The distribution of rent prices was first explored through histograms, overlaid with fitted Normal and Log-Normal density functions to assess its shape as shown in Figure 1 and Figure 2. Scatter and box plots were used to observe the influential variables of rent price and how each variable enters the model. Eight models were fitted using multiple linear regression, and residual analysis, variance inflation factors (VIF), the Akaike information criterion (AIC), and the adjusted R-squared were used for model evaluation. Among the eight models, multiple linear regression with nonlinear covariates was used for the prediction of rent per square meter based on model evaluation.

2.3.1. Concept of a Multiple Linear Regression Model

Researchers frequently explore whether a relationship exists between variables and how they are connected if one is found. Regression models the relationship between a response variable and one or more independent variables (covariates); see [11]. For example, we might examine the relationship between apartment size (independent variable) and monthly rent (dependent or response variable). Regression analysis aims to estimate the parameters of the linear function that best describes the joint distribution of the dependent variable and the covariates [12]. Relationships between variables can be linear, nonlinear (e.g., quadratic or cubic), or absent. Exploratory Data Analysis (EDA) helps identify suitable model structures before fitting regression models. When multiple independent variables are used to predict a continuous response, the approach is called multiple linear regression. In this study, we aim to investigate the relationship between the rent per square meter in Berlin charged for an apartment, characterized by both continuous and discrete covariates.

2.3.2. Concept of a Multiple Linear Regression Model

Researchers frequently explore whether a relationship exists between variables and how they are connected if one is found. Regression models the relationship between a response variable and one or more independent variables (covariates); see [11]. For example, we might examine the relationship between apartment size (independent variable) and monthly rent (dependent or response variable). Regression analysis aims to estimate the parameters of the linear function that best describes the joint distribution of the dependent variable and the covariates [12]. Relationships between variables can be linear, nonlinear (e.g., quadratic or cubic), or absent. Exploratory Data Analysis (EDA) helps identify suitable model structures before fitting regression models. When multiple independent variables are used to predict a continuous response, the approach is called multiple linear regression. In this study, we aim to investigate the relationship between the rent per square meter in Berlin charged for an apartment, characterized by both continuous and discrete covariates.

2.3.3. Model Formulation

In a regression analysis with a continuous response variable $Y_i$ and $p$ covariates or predictors $X_{i1},X_{i2},\cdots ,X_{ik}$ which may be continuous or qualitative (ordinal or nominal) with $n$ observations, let $\left(y_i,x_i^{\top }\right):={\left(y_i,x_{i1},\cdots ,x_{ik}\right)}^{\top }$ , $i=1,\cdots ,n$ , $k=p-1$ , be a pair of the $i\text{th}$ observation $\left(y_i,x_i^{\top }\right)$ of the random vector $\left(Y_i,x_i^{\top }\right)$ , where $x_i={\left(x_{i1},x_{i2},\cdots ,x_{ik}\right)}^{\top }$ , then our objective is to analyze the effects of the covariates on the mean value of the response variable (${\mu }_i\equiv E\left[Y_i\right]$ ). The linear model models the response as a linear function of the predictors together plus an error term, i.e.

$Y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_k{x}_{ik}+{ϵ}_i={\beta }_0+{\displaystyle \sum }_{j=1}^k{\beta }_j{x}_j+{ϵ}_i$ (2.1)

with mean $E\left[Y_i\right]={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_k{x}_{ik}$ .

Multiple linear regression model: The multiple linear regression model is defined as

$Y_i={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_k{x}_{ik}+{\epsilon }_i,i=1,\cdots ,n,$ (2.2)

where ${\epsilon }_i$ is the random error variable, ${\beta }_0$ is the intercept, and the $k$ parameters ${\beta }_1,\cdots ,{\beta }_k$ are the unknown regression parameters to be estimated from $n$ observations $\left(y_i,x_{i1},\cdots ,x_{ik}\right)$ , for $i=1,\cdots ,n$ .

2.3.4. Matrix Representation

It is very easy to represent our linear model in a matrix form [13]. The four different model components below, are defined in order to represent the multiple linear regression model of (2.2) in the matrix-vector notation for our model formulation and calculation.

a) Let $Y=\left(\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right)\in {ℝ}^n$ , be the vector of the response variables.

b) Let $X=\left(\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1k}\\ 1& x_{21}& x_{22}& \cdots & x_{2k}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& x_{n2}& \cdots & x_{nk}\end{array}\right)\in {ℝ}^{n\times p}$ , be the design matrix that contains

$p=k+1$ predictors with their $n$ observations in its rows. The columns correspond to the $p$ unknown regression parameters. Note that the first column which corresponds to the intercept ${\beta }_0$ equals 1 for all $n$ entries.

Denote by $x_i\in {ℝ}^p$ the $i\text{th}$ row of the design matrix.

c) Let $\beta =\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right)\in {ℝ}^p$ , be the unknown vector of the regression coefficients.

d) Let $\epsilon =\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_n\end{array}\right)\in {ℝ}^n$ , be the vector of random error variables.

(Linear regression model in matrix-vector notation) With these notations model (2.2) can be expressed as

$\underset{n\times 1}{\underbrace{\left(\begin{array}{c}{Y}_1\\ Y_2\\ ⋮\\ Y_n\end{array}\right)}}=\underset{n\times p}{\underbrace{\left(\begin{array}{ccccc}1& x_{11}& x_{12}& \cdots & x_{1k}\\ 1& x_{21}& x_{22}& \cdots & x_{2k}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& x_{n2}& \cdots & x_{nk}\end{array}\right)}}\underset{p\times 1}{\underbrace{\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_k\end{array}\right)}}+\underset{n\times 1}{\underbrace{\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ ⋮\\ {ϵ}_n\end{array}\right)}},$

Thus, the multiple linear regression (2.2) can now be written as

$Y=X\beta +{ϵ}_i,\text{with}{ϵ}_i~N_n\left(0,{\sigma }^2{I}_n\right)$ (2.3)

where $I_n$ is the $n\times n$ identity matrix and $N_m\left(\mu ,\Sigma \right)$ denotes the m-dimensional multivariate normal distribution with the mean vector $\mu$ and covariance matrix $\Sigma$ .

2.3.5. Assumptions of the Linear Model

a) Linearity in the covariates: In (2.1), we introduced that the relationship between the covariate vector $x_i$ and the random response $Y_i$ has the form

$Y_i={\beta }_0+{\displaystyle \sum }_{j=1}^p{\beta }_j{x}_{ij}+{\epsilon }_i$

with random error variable ${\epsilon }_i$ satisfying $E\left[{\epsilon }_i\right]=0$ , so that

$E\left[Y_i\right]={\beta }_0+{\displaystyle \sum }_{j=1}^p{\beta }_j{x}_{ij}$

$i=1,\cdots ,n.$

In matrix notation: $E\left[Y\right]=X\beta$ .

b) Homoscedasticity: The error variables ${\epsilon }_i$ have constant variance

$\text{Var}\left[Y_i\right]=\text{Var}\left[{\epsilon }_i\right]={\sigma }^2,i=1,\cdots ,n.$

c) Independence of the random errors: We assume that the error variables ${\epsilon }_i$ are independent and identically distributed (i.i.d.). Then it follows

$\text{Cov}\left(Y_j,Y_{j^{\prime }}\right)=\text{Cov}\left({\epsilon }_j,{\epsilon }_{j^{\prime }}\right)=0,\forall j\ne j^{\prime }$

d) Normality: The random error variables ${\epsilon }_i$ follow a normal distribution.

$Y~N_n\left(X\beta ,{\sigma }^2{I}_n\right)$ (2.4)

where $I_n$ is the $n\times n$ identity matrix and $N_m\left(\mu ,\Sigma \right)$ denotes the m-dimensional multivariate normal distribution with the mean vector $\mu$ and covariance matrix $\Sigma$ , respectively.

In general, we assume a Gaussian error ${\epsilon }_i~N_n\left(0,{\sigma }^2{I}_n\right)$ . This allows us to construct confidence intervals and conduct statistical tests.

Here, $I_n$ is the n-dimensional identity matrix and $N_m\left(\mu ,\Sigma \right)$ denotes the m-dimensional multivariate normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$ , respectively.

2.3.6. Polynomial Regression

Polynomial regression is often appropriate when there exists a relationship between the response and the covariates.

(Polynomial regression). Given a continuous covariate $V_i$ with observations $v_i$ that has a polynomial effect of degree $d$ on the response, then the model $Y_i={\beta }_0+{\beta }_1{V}_i+{\beta }_2{V}_i^2+\cdots +{\beta }_d{V}_i^d+\cdots +{\epsilon }_i$ can be used. Note, it is a linear regression model of the form (2.2) with $x_{ij}=v_i^j$ , $j=1,\cdots ,d$ [14] and [15].

In order to increase numerical stability, we orthonomalize the corresponding design matrix $X=\left(\begin{array}{ccc}1& v_1& v_1^d\\ ⋮& ⋮& ⋮\\ 1& v_n& v_n^d\end{array}\right)$ to $X^{\text{*}}$ , where all columns have unit norms and are orthogonal. In $R$ , this is achieved by $\text{poly}\left(v,d\right)$ , see [16].

2.3.7. Transformations of the Response Variable

Sometimes, the transformation of the response variable is appropriate when non-normality and/or unequal error variances are present in the data. We will consider three different transformations of the response variable in this paper.

(Logarithmic and inverse transformation) Given a response variable $Y$ that has an exponential relationship with the covariates. Let $Y_i^{ln}:=\text{ln}\left(Y_i\right)$ , let $Y_i^{lnln}:=\text{ln}\left(\text{ln}\left(Y_i\right)\right)$ , let $Y_i^{inv}:=\frac{1}{Y_i}$ , then the formulated model

$Y_i=\text{exp}\left({\beta }_0+{\beta }_1{x}_{i1},\cdots ,{\beta }_k{x}_{ik}+{\epsilon }_i\right)$ can be expressed in the form of the linear regression model (2.2) as

$Y_i^{ln}={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_k{x}_{ik}+{\epsilon }_i,i=1,\cdots ,n$ (2.5)

$Y_i^{lnln}={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_k{x}_{ik}+{\epsilon }_i,i=1,\cdots ,n$ (2.6)

$Y_i^{inv}={\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_k{x}_{ik}+{\epsilon }_i,i=1,\cdots ,n$ (2.7)

2.3.8. Estimation of Model Parameters

In this section, we will examine the methods for estimating the unknown parameters in the linear regression model defined in Equation (2.2). Our goal is to determine estimates

$\widehat{\beta }={\left({\widehat{\beta }}_0,\cdots ,{\widehat{\beta }}_k\right)}^{\top }\in {ℝ}^p$ (2.8)

and the error variance $\sigma$ based on $n$ observations. Here $\beta$ is the unknown regression parameter vector.

Note that parameter estimators, which are random quantities are different from their realizations called estimates, which are determined by the values of the observations. We will consider two approaches, least squares (LS) estimation, and maximum likelihood (ML) estimation. These two estimation methods yield the same estimator if the assumptions of independence, homoscedasticity, and normality of errors are satisfied.

2.3.9. Least Squares Estimation Method

Let the fitted values of the Model (2.2) be given as

$\begin{array}{c}{\widehat{Y}}_i={\widehat{\beta }}_0+{\widehat{\beta }}_1{x}_{i1}+\cdots +{\widehat{\beta }}_k{x}_{ik},i=1,\cdots ,n\\ ={x^{\prime }}_i\widehat{\beta }\end{array}$ (2.9)

Also, let the residual denoted by $\widehat{\epsilon }={\left({\widehat{\epsilon }}_1,\cdots ,{\widehat{\epsilon }}_n\right)}^{\prime }\in {ℝ}^n$ , which is the difference between the observed response values $y_i$ and the corresponding fitted values of (2.11), be given as

$\widehat{\epsilon }=y-\widehat{y}=Y-X\widehat{\beta },$ (2.10)

where $\widehat{y}={\left({\widehat{y}}_1,\cdots ,{\widehat{y}}_n\right)}^{\top }\in {ℝ}^n$ in the vector notation. Then, least squares minimizes the residual sum of squares (the sum of the squared deviations) of Equation (2.12).

(Sum of squared deviations) Given the data $\left(y_i,x_i\right),i=1,2,\cdots ,n$ , the sum of the squared deviations which is used in obtaining the estimates $\widehat{\beta }$ of Equation (2.10) for the unknown regression parameters $\beta$ is given as

$Q_{LS}\left(\beta \right)={\displaystyle \sum }_{i=1}^n{\left(y_i-x_i^{\text{T}}\beta \right)}^2={\displaystyle \sum }_{i=1}^n{\widehat{\epsilon }}_i^2={\widehat{\epsilon }}^{\text{T}}\widehat{\epsilon }$ (2.11)

In order to minimize $Q_{LS}\left(\beta \right)$ (2.13), we take the partial derivative of $Q_{LS}\left(\beta \right)$ with respect to $\beta$ and set the result to zero. Then, it follows

$\frac{\partial \left(Q_{LS}\left(\beta \right)\right)}{\partial \beta }=0\iff -2X^{\text{T}}y+2X^{\text{T}}X\beta =0\iff X^{\top }X\beta =X^{\top }y$ (2.12)

We are now interested in solving the least squares normal equations given in (2.14). If the matrix $X$ has a full rank $p$ , then $X^{\text{T}}X$ will be positive definite and will have a unique solution. Thus, the minimum of $Q_{LS}\left(\beta \right)$ is attained at

${\widehat{\beta }}_{LS}={\left(X^{\top }X\right)}^{-1}{X}^{\top }y$ (2.13)

which is the least squares estimate from the normal equations.

2.3.10. Maximum Likelihood Estimation Method

The method of maximum likelihood estimation is based on specifying the distribution we are sampling from and writing the joint density of our sample, unlike in the least squares method where we do not specify the distribution of the response variable $Y_i$ . Considering the assumptions of our linear model, we assumed in Equation (2.4) that the random variables $Y_i$ are normally distributed ($Y~N_n\left(X\beta ,{\sigma }^2{I}_n\right)$ ). Thus, it follows that the likelihood of the vector $\left(\beta ,\sigma \right)$ given the data values $y$ is

$L\left(\beta ,\sigma |y\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{\frac{n}{2}}}\text{exp}\left(-\frac{1}{2{\sigma }^2}{\left(y-X\beta \right)}^{\text{T}}\left(y-X\beta \right)\right)$ (2.14)

Therefore, the corresponding log likelihood is given by

$l\left(\beta ,\sigma |y\right)=-\frac{n}{2}\text{log}\left(2\pi \right)-\frac{n}{2}\text{log}\left({\sigma }^2\right)-\frac{1}{2{\sigma }^2}{\left(y-X\beta \right)}^{\text{T}}\left(y-X\beta \right)$ (2.15)

To maximize this log-likelihood (2.17) with respect to $\beta$ , we differentiate Equation (2.17) with respect to $\beta$ and set it equal to zero [17]. Thus, we have

$\frac{\partial \left(l\left(\beta ,\sigma |y\right)\right)}{\partial \beta }=0\iff -\frac{1}{2{\sigma }^2}\left(-2X^{\text{T}}y+2X^{\text{T}}X\beta \right)=0\iff X^{\top }X\beta =X^{\top }y$ (2.16)

This shows that ${\widehat{\beta }}_{ML}={\widehat{\beta }}_{LS}$ .

Also, differentiating Equation (2.17) with respect to ${\sigma }^2$ and maximizing over ${\sigma }^2$ , we have

${\widehat{\sigma }}^2:=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2=\frac{1}{n}{\Vert \widehat{\epsilon }\Vert }^2$ (2.17)

2.3.11. Distribution of the Estimators

($\widehat{Y}$ and $H$ ) We define the vector of the fitted random values $\widehat{Y}$ as

$\widehat{Y}=X\widehat{\beta }=X{\left(X^{\top }X\right)}^{-1}{X}^{\top }Y$ (2.18)

Also, we define the hat matrix $H$ which gives the projection of the vector $Y$ onto the space that is spanned by the columns of the design matrix $X$ as

$H:=X{\left(X^{\top }X\right)}^{-1}{X}^{\top }\in {ℝ}^{n\times n}$ (2.19)

It can be easily shown in Lemma (2.7.3) that $H$ is both symmetric ($H^{\text{T}}=H$ ) and idempotent ($H^2=H$ ) using the fact that ${\left(X^{\prime }X\right)}^{-1}$ is symmetric (${\left[{\left(X^{\prime }X\right)}^{-1}\right]}^{\prime }={\left(X^{\prime }X\right)}^{-1}$ ).

Lemma 2.1 (Symmetry and Idempotence of H) Let $H=X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }$ . Then $H$ is symmetric and idempotent, i.e., $H^{\prime }=H$ and $H^2=H$ .

Proof.

$H^{\prime }={\left(X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }\right)}^{\prime }=X{\left[{\left(X^{\prime }X\right)}^{-1}\right]}^{\prime }{X}^{\prime }=X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }=H$ , using ${\left(AB\right)}^{\prime }=B^{\prime }{A}^{\prime }$ .

$HH=X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }=H$ .

Since the estimators $\widehat{\beta }$ of the regression coefficients $\beta$ , the fitted values $\widehat{Y}$ and the raw residuals $\widehat{\epsilon }$ are all linear functions of the vector of random variables $Y_i$ , we can apply the transformation rules for expectation and variance-covariance matrix respectively, to show that

$\begin{array}{ll}E\left[\widehat{\beta }\right]=\beta ,\hfill & \text{Var}\left[\widehat{\beta }\right]={\sigma }^2{\left(X^{\top }X\right)}^{-1},\hfill \\ E\left[\widehat{Y}\right]=X\beta ,\hfill & \text{Var}\left[\widehat{Y}\right]={\sigma }^{2}H,\hfill \\ E\left[\widehat{\epsilon }\right]=0,\hfill & \text{Var}\left[\widehat{\epsilon }\right]={\sigma }^2\left(I_n-H\right)\hfill \end{array}$ (2.20)

Considering the normality assumption since $\widehat{\beta }$ , $\widehat{Y}$ , and $\widehat{\epsilon }$ are linear functions of $Y$ , we have

$\begin{array}{l}\widehat{\beta }~N_p\left(\beta ,{\sigma }^2{\left(X^{\top }X\right)}^{-1}\right)\\ \widehat{Y}~N_n\left(X\beta ,{\sigma }^{2}H\right)\\ \widehat{\epsilon }~N_n\left(0,{\sigma }^2\left(I_n-H\right)\right)\end{array}$ (2.21)

It can be shown that the variance estimator ${\widehat{\sigma }}^2$ given in $\left(2.19\right)$ is given by

$\text{E}\left({\widehat{\sigma }}^2\right)=\frac{n-p}{n}{\sigma }^2.$

and an unbiased estimator $s^2$ of ${\sigma }^2$ is given by

$s^2:=\frac{1}{n-p}{\displaystyle \sum }_{i=1}^n{\left(Y_i-{\widehat{Y}}_i\right)}^2=\frac{n}{n-p}{\widehat{\sigma }}^2=\frac{1}{n-p}{\Vert \widehat{\epsilon }\Vert }^2.$ (2.22)

2.4. Goodness of Fit and Model Selection

It is of great importance to know the goodness of the fitted model after estimating the parameters of the linear regression model of (2.2). Thus, we need suitable measures of the goodness of fit. Therefore, we will introduce one of the appropriate measures of the goodness of fit called the coefficient of determination (R²), which determines the proportion of variation of the response variable that is explained by the covariates.

2.4.1. Sum of Squares

(Sum of squares) We define the sum of squares SST (total sum of squares), SSR (regression sum of squares) and SSE (error sum of squares) to quantify the amount of variability explained by the regression model as follows

$\begin{array}{l}\text{SST}:={\displaystyle \sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2\iff \left(\text{total sum of squares}\right)\\ \text{SSR}:={\displaystyle \sum }_{i=1}^n{\left({\widehat{y}}_i-\overline{y}\right)}^2\iff \left(\text{regression sum ofsquares}\right)\\ \text{SSE}:={\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2\iff \left(\text{error sum of squares}\right)\end{array}$ (2.23)

where $\overline{y}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{y}_i}$ . Thus, we can have the decomposition as

${\displaystyle \sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2={\displaystyle \sum }_{i=1}^n{\left({\widehat{y}}_i-\overline{y}\right)}^2+{\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2$ (2.24)

and using the fact that ${\displaystyle {\sum }_{i=1}^n\left({\widehat{y}}_i-\overline{y}\right)\left(y_i-{\widehat{y}}_i\right)}=0$ , it follows from (2.26) that

$\text{SST}=\text{SSR}+\text{SSE}$ (2.25)

2.4.2. Selection of Model (R² and Adjusted R²)

The multiple coefficient of determination R² is a measure of goodness of fit. It measures how well the covariates in the model explain the variance in the response variable, see [18].

(Multiple coefficient of determination ) We define the multiple coefficient of determination R² as

${\text{R}}^2:=\frac{\text{SSR}}{\text{SST}}=1-\frac{\text{SSE}}{\text{SST}}$ (2.26)

We also define the adjusted multiple coefficient of determination ${\text{R}}_{\text{adj}}^2$ as

${\text{R}}_{\text{adj}}^2:=1-\frac{\text{SSE}/\left(n-p\right)}{\text{SST}/\left(n-1\right)}$ (2.27)

The values of the multiple coefficient of determination range from zero to one ($0\le {\text{R}}^2\le 1$ ). Our model accounts for a larger amount of variation of the response when the R² is closer to 1. However, the weakness of R² is that, it always increases when we add more covariates to our model, and therefore cannot be used to compare the goodness of fit for models with different numbers of covariates, see [19]. Thus, the need to establish an appropriate measure ${\text{R}}_{\text{adj}}^2$ which compares models with different numbers of covariates. We will therefore make use of the adjusted multiple coefficient of determination (${\text{R}}_{\text{adj}}^2$ ) as a measure of our model selection in this paper.

2.4.3. Correlation Analysis

To measure the strength and direction of the linear relationship between two continuous variables, we use the correlation analysis. The most commonly used metric is the Pearson correlation coefficient, denoted by $\rho$ for the population and $r$ for the sample. It ranges from −1 to 1, where values close to 1 or −1 indicate strong positive or negative linear relationships, respectively, and values near 0 suggest no linear relationship.

The sample Pearson correlation coefficient between two variables $X$ and $Y$ is given by:

$r=\frac{{{\displaystyle \sum }}_{i=1}^n\left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(X_i-\overline{X}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(Y_i-\overline{Y}\right)}^2}},$

where $\overline{X}$ and $\overline{Y}$ are the sample means of $X$ and $Y$ , respectively. This metric provides a preliminary indication of potential multicollinearity when applied to predictor variables.

2.5. Hypothesis Testing

A statistical hypothesis is an assumption about the form of a population, which based on sample information from the population, seeks to support or reject this assumption. If there is evidence that the null hypothesis (hypothesis of no difference) denoted by $H_0$ is not true, then it is rejected and its alternative denoted by $H_1$ is accepted. Thus, a test of hypothesis is a rule or a procedure used for deciding whether to accept or reject $H_0$ or to determine whether the observed sample differs significantly from expected results under $H_0$ [20]. This concept can be extended in statistical inference for the model parameters of linear regression [21]. For instance, we may want to know if the response variable is significantly influenced by a particular set of covariate variables, which can be expressed in terms of linear combinations of the unknown regression parameters $\beta ={\left({\beta }_0,\cdots ,{\beta }_k\right)}^{\top }$ . We will use the Chi-square, F and the univariate t-distribution since the t-test and the F-test rely on quantities of these distributions.

(Chi-square distribution) A continuous random variable X is said to have a Chi-square distribution with parameter, $\nu$ , if its probability density function is given by

$f_X\left(x|\nu \right)=\frac{2^{-\nu /2}}{\Gamma \left(\nu /2\right)}{x}^{\nu /2-1}{\text{e}}^{-x/2},\nu >0,x>0$

Here, $\nu$ is the degree of freedom, $\text{E}\left(X\right)=\nu$ , $\text{Var}\left(X\right)=2\nu$ . Thus, we say that X follows a Chi-square distribution with $\nu$ degree of freedom ($X~{\chi }_{\nu }^2$ ).

(F-distribution) A continuous random variable X is said to have an F-distribution with degrees of freedom (df) ${\nu }_1$ and ${\nu }_2$ , if its pdf is given by

$f\left(x\right)=\frac{\Gamma \left(\frac{{\nu }_1+{\nu }_2}{2}\right){\left(\frac{{\nu }_1}{{\nu }_2}\right)}^{\frac{{\nu }_1}{2}}{x}^{\frac{{\nu }_1}{2}-1}}{\Gamma \left(\frac{{\nu }_1}{2}\right)\Gamma \left(\frac{{\nu }_2}{2}\right){\left(1+\frac{{\nu }_{1}x}{{\nu }_2}\right)}^{\frac{{\nu }_1+{\nu }_2}{2}}},x\ge 0.$ (2.28)

If $X_1~{\chi }_{v_1}^2$ and $X_2~{\chi }_{{\nu }_2}^2$ and are independent, it follows in (2.30) that X is F-distributed with ${\nu }_1$ and ${\nu }_1$ df.

$X=\frac{X_1/{\nu }_1}{X_2/{\nu }_2}~F_{{\nu }_1,{\nu }_2}$ (2.29)

(Univariate t-distribution) A continuous random variable X is said to have a Univariate t-distribution with degree of freedom df $\nu$ , if its pdf is given by

$f_{\nu }\left(x;\mu ,{\sigma }^2\right):=\frac{\Gamma \left(\frac{\nu +1}{2}\right)}{\Gamma \left(\frac{\nu }{2}\right)\sqrt{\left(\pi \nu \right)}\sigma }{\left\{1+{\left(\frac{x-\mu }{\sigma }\right)}^2\frac{1}{\nu }\right\}}^{-\frac{\nu +1}{2}},\nu \ge 1$ (2.30)

$\text{E}\left(X\right)=\mu \text{and}\text{Var}\left(X\right)=\frac{\nu }{\nu -2}{\sigma }^2.$

If $X_1~N\left(0,1\right)$ and $X_2~{\chi }_n^2$ and are independent, it can be shown in (2.32) that T has a t-distribution with $\nu$ df.

$T=\frac{X_1}{\sqrt{\frac{X_2}{\nu }}}~t_{\nu }.$ (2.31)

2.6. F-Test

(General testing problem) We define the general testing problem as

$H_0:C\beta =d\text{versus}{H}_1:C\beta \ne d$ (2.32)

where matrix $C\in {ℝ}^{q\times p}$ with $\text{rank}\left(C\right)=q$ and $d\in {ℝ}^q$ is called the general linear hypothesis. Using the distribution of the estimated regression coefficients $\widehat{\beta }$ given in (2.23), if $H_0$ is true, it follows that

$\widehat{\vartheta }=C\widehat{\beta }-d\stackrel{H_0}{~}{N}_q\left(0,{\sigma }^{2}C{\left(X^{\top }X\right)}^{-1}{C}^{\top }\right).$ (2.33)

We used the fact that

$\begin{array}{l}\text{E}\left(\widehat{\vartheta }\right)=C\text{E}\left(\widehat{\beta }\right)-d=C\text{E}\left(\beta \right)-d=0\\ \text{Var}\left(\widehat{\vartheta }\right)=C\text{Var}\widehat{\beta }{C}^{\top }={\sigma }^{2}C{\left(X^{\top }X\right)}^{-1}{C}^{\top }.\end{array}$ (2.34)

Also, using the spectral decomposition for a specific covariance and considering the definition of ${\chi }^2$ -distribution, it can be shown that

$\frac{1}{{\sigma }^2}{\widehat{\vartheta }}^{\top }{\left(C{\left(X^{\top }X\right)}^{-1}{C}^{\top }\right)}^{-1}\widehat{\vartheta }\stackrel{H_0}{~}{\chi }_q^2$ (2.35)

One can also show that $\frac{1}{{\sigma }^2}{\widehat{\vartheta }}^{\top }{\left(C{\left(X^{\top }X\right)}^{-1}{C}^{\top }\right)}^{-1}\widehat{\vartheta }~{\chi }_q^2$ and $\text{SSE}/{\sigma }^2~{\chi }_{n-p}^2$ , and are independent ${\chi }^2$ distributed. We therefore define the statistic under $H_0$ as

$F=\frac{{\widehat{\vartheta }}^{\top }{\left(C{\left(X^{\top }X\right)}^{-1}{C}^{\top }\right)}^{-1}\widehat{\vartheta }/q}{\text{SSE}/\left(n-p\right)}\stackrel{H_0}{~}{F}_{q,n-p}.$ (2.36)

If the null hypothesis $H_0:C\beta =d$ holds (ie, $C\beta -d=0$ ), then we will reject $H_0$ for large values of $F$ since small value of $C\widehat{\beta }-d$ is expected.

Let the least square estimate among those $\beta$ vectors which satisfy $C\beta =d$ be denoted as ${\beta }_{H_0}$ , i.e. it minimizes (2.13) under the condition $C\beta =d$ . We define the corresponding sum of squares error for the LS fit under $H_0$ as

${\text{SSE}}_{H_0}:={\Vert y-X{\widehat{\beta }}_{H_0}\Vert }^2$

It can also be shown that ${\widehat{\vartheta }}^{\top }{\left(C{\left(X^{\top }X\right)}^{-1}{C}^{\top }\right)}^{-1}\widehat{\vartheta }={\text{SSE}}_{H_0}-\text{SSE}$ is true which allows us to give the general F-test.

(General F test in linear regression) We define the test statistic F under the null hypothesis $H_0:C\beta =d$ versus $H_1:C\beta \ne d$ in the linear regression model of (2.2) as

$F=\frac{\left({\text{SSE}}_{H_0}-\text{SSE}\right)/q}{\text{SSE}/\left(n-p\right)}\stackrel{H_0}{~}{F}_{q,n-p}$ (2.37)

We reject $H_0$ against $H_1$ at level $\alpha$ if

$F>F_{\left(1-\alpha \right),q,n-p}$ (2.38)

Here, $F_{\left(1-\alpha \right),q,n-p}$ is the $\left(1-\alpha \right)$ quantile of an $\text{F}$ distribution with $q$ and $n-p$ df. The quantity $n-p$ is also called the residual degree of freedom. Thus, we can now summarize the F-test procedure for our model (2.2) as follows:

Hypothesis

$\begin{array}{l}{H}_0:\left({\beta }_1,\cdots ,{\beta }_k\right)=0\\ H_1:\left({\beta }_1,\cdots ,{\beta }_k\right)\ne 0\end{array}$

$H_0:{\beta }_1={\beta }_2=\cdots ={\beta }_k=0$

(No significant relationship exists between rent sqm and the predictors.)

$H_1:\text{At}\text{least}\text{one}{\beta }_j\ne 0,\text{for}\text{some}j\in \left\{1,2,\cdots ,k\right\}$

(A significant relationship exists between rent sqm and at least one predictor.)

Test statistic = F, defined in (2.37)

Rejection Rule: Reject $H_0$ at level $\alpha$ , if $F>F_{\left(1-\alpha \right),q,n-p}$

2.7. T-Test

(t-test) We define the t-test procedure for our model (2.2) as follows, since in a t-test, the test statistic is computed for each ${\beta }_j$ , see [22].

Hypotheses:

$H_0:{\beta }_j=0\text{versus}{H}_1:{\beta }_j\ne 0$

Test statistic:

$T_j=\frac{{\widehat{\beta }}_j}{\widehat{se}\left({\widehat{\beta }}_j\right)}~t_{n-p},\text{under}{H}_0$ (2.39)

Here, $\widehat{se}\left({\widehat{\beta }}_j\right):=s\sqrt{{\left({\left(X^{\top }X\right)}^{-1}\right)}_{jj}}$ is the estimated standard error of ${\widehat{\beta }}_j$ and $s=\sqrt{s^2}$ defined in Equation (2.22)

Rejection Rule: Reject $H_0$ at level $\alpha$ , if $\left|T_j\right|>t_{n-p,1-\alpha /2}$

Note that t-test is a special case of the F-test, in particular we have $F_j=T_j^2$ since if

$\begin{array}{l}{F}_j:=\frac{{\widehat{\beta }}_j^2}{{\left({\left(X^{\top }X\right)}^{-1}\right)}_{jj}\text{SSE}/\left(n-p\right)}\stackrel{H_0}{~}{F}_{1,n-p}\\ T_j:=\frac{{\widehat{\beta }}_j}{\widehat{se}\left({\widehat{\beta }}_j\right)}\stackrel{H_0}{~}{t}_{n-p}\end{array}$ (2.40)

2.8. Analysis of Variance (ANOVA)

ANOVA is mostly used to summarize the hypothesis tests results in linear models in a tabular form. Given two models $M_{\text{reduced}}$ and $M_{\text{full}}$ which are nested: $M_{\text{reduced}}\subset M_{\text{full}}$ , that is, all covariates of the reduced model are contained in the full model, we define the ANOVA-test ratio for the comparison of $M_{\text{reduced}}$ and $M_{\text{full}}$ as follows

$F=\frac{\left({\text{SSE}}_{\text{reduced}}-{\text{SSE}}_{\text{full}}\right)/\left(n-p_{\text{full}}\right)}{{\text{SSE}}_{\text{reduced}}/\left(p_{\text{full}}-p_{\text{reduced}}\right)}~F_{n-p_{\text{full}},p_{\text{full}}-p_{\text{reduced}}}$ (2.41)

Hypotheses

$H_0:{\beta }_i=0\text{versus}{H}_1:{\beta }_i\ne 0$

Test statistic: $F$ , defined in Equation (2.41)

Rejection Rule: Reject $H_0$ at level $\alpha$ , if $F>F_{\left(1-\alpha \right),n-p_{\text{full}},p_{\text{full}}-p_{\text{reduced}}}$

2.9. Analysis of Residuals

After estimating the model parameters, the credibility of the assumptions of linearity, normality of errors, and homoscedasticity for the given data can be assessed using residuals. It is therefore of importance to study the residual in order to examine in what extent our model assumptions may be violated. Therefore, taking a look at the patterns in the residual plots could help us understand if our model assumptions are violated or not. This is called analysis of residual. Residual plots can equally help us to decide whether to transform any of the covariates which we may want to include in the model or not. We will introduce three types of residuals for the $i\text{th}$ observation namely: raw residuals, internally studentized residuals and externally studentized residuals.

2.10. Raw Residuals Check

The raw residual vector $\widehat{\epsilon }$ was defined in (2.12), and its distribution in (2.22). Thus, we have that

$\widehat{\epsilon }=y-\widehat{y}=Y-X\widehat{\beta }=Y-\underset{H}{\underbrace{X{\left(X^{\prime }X\right)}^{-1}{X}^{\prime }}}Y=Y-HY=\left(I_n-H\right)Y$

and

$\text{Var}\left({\widehat{\epsilon }}_i\right)={\sigma }^2\left(1-h_{ii}\right),$ (2.42)

where $h_{ii}=x_i^{\top }{\left(X^{\top }X\right)}^{-1}{x}_i$ is the $i$ th diagonal element of the hat-matrix $H$ defined in (2.21). Since $\text{Var}\left({\widehat{\epsilon }}_i\right)$ still changes under the homoscedasticity condition, we need to standardize the raw residuals. This standardization gives rise to both internally studentized residuals and externally studentized residuals.

Internally studentized residuals

The internally studentized residuals are defined by

$s_i:=\frac{{\widehat{\epsilon }}_i}{\sqrt{1-h_{ii}}s},$ (2.43)

Here, $s^2$ is the estimate of ${\sigma }^2$ defined in (2.19). One can now analyze the variances and conclude whether the assumption of homoscedasticity is violated or not using the standardized residuals by plotting the standardized residuals versus the predicted values ${\widehat{y}}_i$ . However, the deficiency of the internally studentized residuals is that it is not robust against outliers since all data is used to estimate $\sigma$ and Equation (2.43) is t-distributed. This leads to the definition of externally studentized residuals which is more robust against outliers.

2.11. Externally Studentized Residuals or Jackknifed Residuals

To solve the problem of the non robustness of the internally studentized residuals, we define a new model just like the model of (2.3), but it is based on “drop-one-observation” of the data, which contains all observations except the $i_{th}$ observation as

$Y_{-i}=X_{-i}{\beta }_{-i}+{\epsilon }_{-i},$ (2.44)

Here, $X_{-i}$ denotes the design matrix without the row $i\text{th}$ and $Y_{-i}$ is the response vector $Y$ with the $i\text{th}$ observation $y_i$ removed. We define the fitted values corresponding to the model given in (2.44) as

${\widehat{y}}_{i,-i}:=x_i^{\top }{\widehat{\beta }}_{-i}$ (2.45)

where ${\widehat{\beta }}_{-i}$ is the associated least squares estimates of ${\beta }_{-i}$ . We define the corresponding residual also called the $i\text{th}$ predictive residual as

${\widehat{\epsilon }}_{i,-1}:=y_i-{\widehat{y}}_{i,-i}$ (2.46)

Using (2.46), we obtain the estimate $s_{-i}^2$ of the error variance ${\sigma }^2$ which does not include the $i\text{th}$ observation as

$s_{-i}^2:=\frac{{{\displaystyle \sum }}_{j=1,j\ne i}^n{\left(y_j-x_j^{\top }{\widehat{\beta }}_{-i}\right)}^2}{n-p-1}$ (2.47)

Finally, we now define the externally studentized residuals which is also called jackknifed residual. The externally studentized residuals $t_i$ which is based on the “drop-one-observation” of the data, are defined as

$t_i:=\frac{r_{i,-i}}{\sqrt{1-h_{ii}}{s}_{-i}}$ (2.48)

2.12. Statistical Checks for the Plausibility of the Linear Model Assumptions

2.12.1. Linearity

The check we are going to use is the residuals versus the fitted values plot. If this plot has no trend, then we assume the linearity assumption as plausible [23].

2.12.2. Homoscedasticity

We are interested in checking if $\text{Var}\left[Y_i\right]=\text{Var}\left[{\epsilon }_i\right]={\sigma }^2$ holds. To check this, we use again the standardized residual versus the residual plots. If the standardized residuals are not spread equally along the range of the fitted values, then we interpret the homoscedasticity assumption as not plausible, see [24].

2.12.3. Independence

To check if $\text{Cov}\left({\epsilon }_j,{{\epsilon }^{\prime }}_j\right)=\rho =0$ holds, we plot the residuals versus the covariates to see if the residuals are randomly and symmetrically distributed around zero. If this is true, we assume that the independence assumption is plausible [23].

2.12.4. Normality

To check for ${\epsilon }_i~N_n\left(0,{\sigma }^2{I}_n\right)$ , we use the Quantile versus Quantile plot (QQPlot). If we do not have a straight line on the QQ plots of our variable versus the theoretical normal quantile, then we assume that the normality assumption is not plausible [25].

2.12.5. Multicollinearity

To check for multicollinearity among explanatory variables $X_1,X_2,\cdots ,X_p$ , we assess whether there is a strong linear relationship between them, which can inflate the standard errors of the estimated coefficients ${\widehat{\beta }}_j$ . This is commonly evaluated using the Variance Inflation Factor (VIF), defined as

${\text{VIF}}_j=\frac{1}{1-R_j^2},$

where $R_j^2$ is the coefficient of determination from regressing $X_j$ on the remaining predictors. A ${\text{VIF}}_j>5$ suggests a potentially problematic level of multicollinearity. Variables exceeding this threshold must be examined and removed if necessary to enhance model stability and interpretability. [26] [27].

3. Data Description and Management

The data set contains 2,651,885 observations and 59 attributes from 2007 to 2020. The data has both quantitative and qualitative covariates with rent per square meter (rent_sqm) as the response variable. We focus on the most relevant 31 variables such as “the additional cost”, “heat cost”, “construction year”, etc. The quantitative covariates are summarized as follows: Min = Minimum, 25% = 1st quartile, 50% = Median, $\overline{X}$ = Mean, 75% = 3rd quartile, Max = Maximum and Not available = NA. On the other hand, the qualitative covariates are summarized with their respective categories. Note that costs are expressed in EUR and rounded to two decimal digits and the following data summaries in Table 1 and Table 2 represent the whole data set.

Table 1. Description of quantitative variables.

The descriptive statistics reveal notable variation across key quantitative variables relevant to rent modeling. Rent per square meter varies significantly, indicating differences in property value influenced by factors such as location, condition, and amenities. Additional cost (adcost) and heating costs (heatcost) also show a wide spread, with some properties incurring substantially higher energy costs, possibly due to inefficiency. The construction (conyear) and modernization years (lmod) extend across several decades, which capture both historic and newly renovated buildings. Living (lspace) and floor space (fspace) differ widely, which shows a mix of compact and spacious apartments in the dataset. Energy consumption (energycon) values highlight a range of efficiency levels, which likely impact rental prices. The duration of the advertisement ranges from zero to several weeks, suggesting varying market demand and listing strategies. These differences confirm the complexity of accurately predicting rental prices. As a result, the study proposes the use of a log-normal model incorporating nonlinear covariates to better capture the heterogeneity and improve prediction accuracy in rent data.

Table 2. Description of qualitative variables.

3.1. Data Sets

We split the date set described in Table 1 and Table 2 into two sub data sets: Berlin 2015 and Berlin 2019. The number of rental properties contained in each data set is given in Table 3. The summaries of the response variable and the quantitative covariates are given in Table 4, while in Table 5, we give the summary of each qualitative variable followed by their percentage.

Table 3. Number of rental properties in the two data sets.

Table 4 Univariate data summaries of quantitative covariates: first row = Berlin 2015, second row = Berlin 2019.

Table 5. Univariate data summaries of qualitative covariates: first row = Berlin 2015, second row = Berlin 2019.