<figure><img src="image_1.png" /></figure>

HAL
open science

# Predicting sex-based thermal sensation differences: a multilevel ordinal regression approach

Matteo Favero, Marika Vellei, Andrew Sonta

<figure><img src="image_2.png" /></figure>

## To cite this version:

Matteo Favero, Marika Vellei, Andrew Sonta. Predicting sex-based thermal sensation differences: a multilevel ordinal regression approach. 2025 Building Simulation Conference, Aug 2025, Brisbane, Australia. 10.26868/25222708.2025.1385. hal-05411129

HAL Id: hal-05411129

https://hal.science/hal-05411129v1

Submitted on 11 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Predicting sex-based thermal sensation differences: a multilevel ordinal regression approach

Matteo Favero¹*, Marika Vellei², Andrew Sonta¹

¹ETHOS Lab, École Polytechnique Fédérale de Lausanne (EPFL), Fribourg, Switzerland

²Institute of Mechanics and Engineering (I2M), University of Bordeaux, Talence, France

*Corresponding Author: matteo.favero@epfl.ch

## Abstract

Building energy consumption plays a major role in global energy use, with a significant portion allocated to maintaining thermal comfort. Yet, many occupants remain dissatisfied with their thermal environment, partly due to the limitations of the ‘one-size-fits-all’ approach, such as Fanger’s PMV model, which does not account for individual differences. This study addresses sex-based differences by re-analysing existing data involving 1,600 participants. Using a multilevel ordinal regression within a Bayesian framework, the results reveal that females exhibit a greater response to temperature than males, suggesting that accounting for sex differences could improve thermal comfort. This research contributes to developing more inclusive thermal comfort standards for occupant well-being and energy efficiency.

## Key Innovations

* This paper uses a multilevel ordinal regression to predict sex-based differences in the thermal sensation of 160 groups of people (five males and five females in each group) exposed to 160 combinations of air temperature and relative humidity.

* The use of Bayesian inference allows for (i) the quantification of uncertainty surrounding parameter estimates and its propagation within the model, (ii) the integration of prior knowledge, (iii) the ability to make direct probability statements, and (iv) the avoidance of the dichotomisation of evidence (‘significant’ vs ‘non-significant’) inherent in null hypothesis significance testing.

## Practical Implications

This research shows that despite high variability, females exhibit, on average, a greater response to cold and hot environments than males, and this is more pronounced at lower and higher temperatures. As a result, it is essential not to overlook sex-based differences when designing indoor thermal environments.

## Introduction

Energy consumption in buildings constitutes a significant portion of global energy use, where a substantial share of this energy is used to provide thermally comfortable indoor environments. Despite these considerable energy investments, many occupants remain dissatisfied with their thermal environment (Graham et al., 2021). This dissatisfaction could be attributed to the prevalent ‘one-size-fits-all’ approach (e.g., Fanger’s PMV model), which overlooks individual differences. Among them, sex-based

differences are disputed, where well-known variations in morphology, body composition and metabolism are neglected in practice. Understanding and integrating these differences for the provision of comfortable spaces is essential for developing inclusive standards and enhancing energy efficiency while ensuring occupant well-being. Although some statistical analyses support these differences, the thermal comfort literature does not provide clear and consistent results (Wang et al., 2018). Using a frequentist statistical framework, subjective comfort ratings (e.g., thermal sensation vote, TSV) measured on an ordinal scale are often analysed with linear regression, which can lead to inaccurate estimates and misleading conclusions (Favero et al., 2023; Liddell & Kruschke, 2018). Additionally, conclusions often rely on classifying results as either ‘significant’ or ‘nonsignificant’ (i.e., *p*-values ≤ 0.05 or *p*-values > 0.05) to determine the presence or absence of an ‘effect.’ This approach is detrimental as it disregards the variability in human response and the uncertainty associated with statistical inference.

In this study, we re-analyse existing experimental data (extracted from a technical report written by Rohles (1970)) that were originally analysed using linear regression and were part of the basis on which Fanger (1970) built his thermal comfort equations. Our objective is to predict the thermal sensation votes of males and females exposed to a wide range of steady-state thermal environments and compare it with the PMV in the same conditions. To perform the re-analysis, we used Bayesian statistical inference. Bayesian inference is an example of inductive inference (or inverse probability), where the probability of a hypothesis is estimated conditional on the observed data and can mathematically be written as $Pr(H|D)$ (where *H* stands for hypothesis and *D* for data). As a result, it allows a direct probability statement (i.e., the probability that a hypothesis is true) conditionally on the model, the data and the priors. Additionally, Bayesian data analysis allows for (i) the quantification of uncertainty surrounding parameter estimates and its propagation within the model, (ii) the integration of prior knowledge, and (iii) the avoidance of the dichotomisation of results (‘significant’ vs ‘non-significant’).

## Methods

We analysed data derived from a remarkably large laboratory experiment performed by Rohles (1971) in the climate chamber of the Institute for Environmental Research, Department of Mechanical Engineering, Kansas State University (KSU). A total of 1,600 college-
age students (800 males and 800 females), divided into groups of ten (five males and five females), were exposed to 160 different test conditions for three hours. Each participant experienced only one condition. These conditions result from the combinations of twenty dry-bulb temperatures (ranging from 15.6 to 36.7°C in increments of 1.1°C) and eight relative humidities (ranging from 15% to 85% in increments of 10%). Airspeed was maintained at 0.1 m/s, the mean radiant temperature matched the air temperature, and half of the tests were conducted in the afternoon (1 pm to 4 pm) and half in the evening (7 pm to 10 pm). Participants remained sedentary and, dressed in standard clothing with a thermal resistance of 0.6 clo (i.e., cotton twill shirts, trousers and sweat socks provided on the day of the experiment), were asked to describe their thermal sensations using a seven-point ordinal scale ranging from 1 to 7 (1 ‘cold,’ 2 ‘cool,’ 3 ‘slightly cool,’ 4 ‘comfortable,’ 5 ‘slightly warm,’ 6 ‘warm’ and 7 ‘hot’) every half-hour after the first hour.

## Data preparation

Table 1 describes the lists of variables we considered during the analysis. We first converted indoor air temperature from degrees Fahrenheit to Celsius and then centred it by subtracting 25.35 °C (we chose this value for air temperature (ta) because it corresponds to a PMV ≈ 0.00, assuming mean radiant temperature (tr) = 25.35 °C, airspeed (va) = 0 m/s, relative humidity (rh) = 50%, metabolic rate (M) = 1 met, and clothing level (Icl) = 0.6 clo). We first scaled relative humidity by dividing by 5 (so that 1-unit change corresponds to a 5% RH change) and then centred it by subtracting 10.

<table><caption>Table 1: List of variables</caption><thead><tr><th>Variable</th><th>Code</th><th>Unit</th></tr></thead><tbody><tr><td>Time of day</td><td>time.day</td><td>Afternoon (reference) / Evening</td></tr><tr><td>Air temperature</td><td>temp_c</td><td>°C (centred)</td></tr><tr><td>Relative humidity</td><td>rh_c</td><td>% (scaled and centred)</td></tr><tr><td>Sex (male reference)</td><td>sex</td><td>Male (reference) / Female</td></tr><tr><td>Experimental condition</td><td>COND</td><td>—</td></tr></tbody></table>

We extracted the data from the tabulation of the thermal sensation assessments provided in a technical report written by Rohles (1970). Unfortunately, the tabulated data made it impossible for us to recover the subjects' identifiers and determine how each participant voted over time.

## Statistical modelling

In this study, we selected a multilevel ordinal regression model to account for the data's clustering and the ordinal nature of the response variable (i.e., thermal sensation). In the literature, ordinal regression models refer to a broad class of models, each with a different rationale and, consequently, a different application. In this study, we used the ordered multinomial model. Like other generalised linear models, the ordered multinomial model requires the specification of a link function. We used the

probit link in this paper since it was used by Fanger (1970) and Ballantyne et al. (1977) when re-analysing the KSU data; however, other link functions are possible (e.g., logit link). We refer to this model as a multilevel ordered multinomial model with a probit link. This model can be generally written as:

$$ 
\begin{array}{l}
Y_i \sim \text{Multinomial}(n, \pi_i) \\
\text{Probit}(\gamma_{ik}) = \frac{\tau_k - \eta_i}{\sigma} \\
\eta_i = \alpha_0 + \sum_{l=1}^{l} \beta_l x_{l,i} + \alpha_{COND[i]}
\end{array}
\qquad (1) $$

In Eq. (1), the conditional distribution of the ordinal response variable $Y_i$ (with $K$ ordered levels) for the $i^{th}$ observation is assumed to follow a multinomial distribution with probability vector $\pi_i = \{\pi_{i1}, \dots, \pi_{ik}\}$, where $\pi_{ik} = \text{Pr}(Y_i = k)$. The term $n$ indicates the number of trials and is fixed to 1 since each of the $i^{th}$ observations was a single entry in the model. Since we assumed a probit link, the cumulative probability for $\pi_{ik}$ is the cumulative distribution function of the standard normal distribution $\Phi$, hence $\gamma_{ik} = \text{Pr}(Y_i \le k) = \Phi((\tau_k - \eta_i)/\sigma)$. The terms $\{\tau_k\}$ are the strictly ordered threshold parameters (with $k \in \{1, \dots, K-1\}$). The term $\eta_i$ is the linear predictor composed of a constant intercept $\alpha_0$ (fixed to 0 for identification), the coefficients $\beta_l$ for the $l^{th}$ predictor of observation $i$ (i.e., $x_{l,i}$) and a varying intercept $\alpha_{COND[i]}$ (this term is generally called ‘random effect;’ we avoid this term following the recommendations of Gelman and Hill (2006, pp. 245–256)). The term $\sigma$ is the standard deviation (fixed to 1 for identification). The model in Eq. (1) can be simplified as:

$$ 
\begin{array}{l}
Y_i \sim \text{Multinomial}(1, \pi_i) \\
\text{Probit}(\gamma_{ik}) = \tau_k - \eta_i \\
\eta_i = \sum_{l=1}^{l} \beta_l x_{l,i} + \alpha_{COND[i]}
\end{array}
\qquad (2) $$

All statistical analyses were performed using R (R Core Team, 2023) with the RStudio integrated development environment (Posit team, 2024). Regression analyses were performed with the *brms* package (Bürkner, 2017), and models were compared with the *loo* package (Vehtari et al., 2017). Data wrangling and visualisation were carried out via the *tidyverse* package (Wickham et al., 2019) and *tidybayes* package (Kay, 2023).

## Prior specification

The prior distributions for the different model parameters are:

$$ 
\begin{array}{l}
\tau_k \sim \text{Normal}(0, 4) \\
\beta_D \sim \text{Normal}(0, 0.4) \\
\beta_T \sim \text{Normal}(0.5, 0.3) \\
\beta_H \sim \text{Normal}(0.25, 0.15) \\
\beta_F \sim \text{Normal}(0, 0.4) \\
\beta_{TH} \sim \text{Normal}(0, 0.06) \\
\beta_{TF} \sim \text{Normal}(0, 0.15) \\
\beta_{HF} \sim \text{Normal}(0, 0.08) \\
\beta_{THF} \sim \text{Normal}(0, 0.08)
\end{array}
\qquad (3) $$
$$ 
\begin{aligned}
\alpha_j &\sim \text{Normal}(0, \delta), \text{ for } j = 1 \dots 160 \\
\delta &\sim \text{Exponential}(1)
\end{aligned}
 $$

where $\tau_k$ are the threshold parameters, $\beta_D$ is the coefficient for time of day, $\beta_T$ is the coefficient for air temperature, $\beta_H$ is the coefficient for relative humidity, $\beta_F$ is the coefficient for females, $\beta_{TH}$, $\beta_{TF}$, $\beta_{HF}$ are the coefficients for the two-way interactions and $\beta_{THF}$ is the coefficient for the three-way interaction. The term $\alpha_j$ represents the varying intercept in the model, one for each of the 160 conditions.

### Non-centred parameterisation

In Eq. (3), the term $\alpha_j$ is normally distributed with zero mean and standard deviation $\delta$ to be estimated from the data, assuming an exponential prior. The normal distribution for $\alpha_j$ is a hierarchical prior (i.e., a prior inside a prior) using a centred parametrisation. However, for multilevel models, *brms* does not support this parametrisation and internally re-write a statistically equivalent model using a non-centred parameterisation:

$$ 
\begin{aligned}
\alpha_j &= z_{\alpha,j} \times \delta, \text{ for } j = 1 \dots 160 \\
z_{\alpha,j} &\sim \text{Normal}(0, 1) \\
\delta &\sim \text{Exponential}(1)
\end{aligned}
 \quad (4) $$

### Computing predictions using a multilevel model

In a multilevel ordinal model, two types of prediction can be computed depending on how the varying intercepts are handled. The cluster-specific predictions were obtained by fixing the varying intercept $\alpha_{COND[i]}$ equal to 0 and then calculating the probabilities for specific values of the predictors $x$. The sample-averaged predictions were calculated by drawing each posterior sample from the posterior draws of a randomly chosen existing intercept $\alpha_{COND[i]}$ and then calculating the probabilities for the selected $x$-values. Therefore, the cluster-specific predictions can be interpreted as the predictions for the typical (i.e., average) condition, whereas the sample-averaged predictions can be interpreted as the prediction across conditions on average (i.e., averaged over the 160 conditions).

### Model comparison and selection

Ten different models were fitted and compared to select the one with the best predictive performance. Table 2 summarises the variables included in the ten models. These ten models were fitted separately for the five

different exposure times (see Table 3). For example, five instances of *mod.0* were fitted to data from exposure at times 1, 2, 3, 4, and 5 separately. As previously mentioned, fitting a model including all observations was inappropriate because the subjects' identifiers from the tabulated data were unavailable. Fitting such a model will ignore that the data were clustered in groups (i.e., 160 conditions).

<table><caption>Table 2: Models summary.</caption><thead><tr><th>Model →<br>Variable code ↓</th><th>mod.0</th><th>mod.1</th><th>mod.2</th><th>mod.3</th><th>mod.4</th><th>mod.5</th><th>mod.6</th><th>mod.7</th><th>mod.8</th><th>mod.9</th></tr></thead><tbody><tr><td>time.day</td><td>•</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>temp_c</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>rh_c</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>sex</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>temp_c × rh_c</td><td></td><td></td><td>•</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>temp_c × sex</td><td></td><td></td><td></td><td>•</td><td></td><td>•</td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>rh_c × sex</td><td></td><td></td><td></td><td></td><td>•</td><td></td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>temp_c × rh_c × sex</td><td></td><td></td><td></td><td></td><td></td><td>•</td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>COND</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr></tbody></table>

Table 3 shows the comparison among models and ranks them according to their expected log predictive density (*elpd*), which was estimated by the approximate leave-one-out (LOO) cross-validation using Pareto smoothed importance sampling (PSIS). More precisely, the models are compared in terms of the difference in the expected log predictive density (i.e., *elpd<sub>diff</sub>*). This difference is computed by making a pairwise comparison between each model and the model with the largest *elpd* (i.e., the model ranking first); the associated *se<sub>diff</sub>* is the standard error of this difference. As a 'rule of thumb' if |*elpd<sub>diff</sub>*| < 4, the difference is small, and the models compared have very similar predictive performance. If |*elpd<sub>diff</sub>*| > 4, then one can compare this difference to the standard error of the difference (i.e., *se<sub>diff</sub>*). (Sivula et al., 2020). As a metric for comparison, we used |*elpd<sub>diff</sub>*| > 2 *se<sub>diff</sub>*. We selected the value of 2 times the *se<sub>diff</sub>* for the following reasons: assuming that the distribution of *elpd<sub>diff</sub>* is approximately normal, then 95 % of the *elpd<sub>diff</sub>* values will be within 2 *se<sub>diff</sub>* (i.e., two standard errors of the difference), and zero will not be included in this 95 % interval (where zero implies equal *elpd*). Given these assumptions, any values of |*elpd<sub>diff</sub>*| > 2 *se<sub>diff</sub>*

<table><caption>Table 3: Models' comparison (difference in the expected log predictive density, *elpd<sub>diff</sub>* and its standard error, *se<sub>diff</sub>*).</caption><thead><tr><th>Model order given elpd</th><th>Obs 1 (1.0 h) [elpd<sub>diff</sub>, se<sub>diff</sub>]</th><th>Obs 2 (1.5 h) [elpd<sub>diff</sub>, se<sub>diff</sub>]</th><th>Obs 3 (2.0 h) [elpd<sub>diff</sub>, se<sub>diff</sub>]</th><th>Obs 4 (2.5 h) [elpd<sub>diff</sub>, se<sub>diff</sub>]</th><th>Obs 5 (3.0 h) [elpd<sub>diff</sub>, se<sub>diff</sub>]</th></tr></thead><tbody><tr><td>1</td><td>mod.5 [0.0, 0.0]</td><td>mod.8 [0.0, 0.0]</td><td>mod.5 [0.0, 0.0]</td><td>mod.5 [0.0, 0.0]</td><td>mod.5 [0.0, 0.0]</td></tr><tr><td>2</td><td>mod.3 [-0.5, 1.4]</td><td>mod.5 [0.0, 1.6]</td><td>mod.8 [-1.0, 0.2]</td><td>mod.8 [-0.7, 0.2]</td><td>mod.8 [-0.6, 0.3]</td></tr><tr><td>3</td><td>mod.8 [-0.7, 1.3]</td><td>mod.9 [-1.1, 0.5]</td><td>mod.9 [-1.3, 1.2]</td><td>mod.9 [-1.6, 0.4]</td><td>mod.9 [-2.1, 0.5]</td></tr><tr><td>4</td><td>mod.7 [-0.8, 1.9]</td><td>mod.3 [-1.5, 3.1]</td><td>mod.3 [-2.3, 2.7]</td><td>mod.3 [-3.5, 3.2]</td><td>mod.3 [-3.5, 3.2]</td></tr><tr><td>5</td><td>mod.9 [-1.6, 1.4]</td><td>mod.7 [-1.6, 2.6]</td><td>mod.7 [-3.2, 2.7]</td><td>mod.7 [-4.4, 3.2]</td><td>mod.7 [-4.6, 3.2]</td></tr><tr><td>6</td><td>mod.1 [-47.8, 10.1]</td><td>mod.2 [-50.8, 10.4]</td><td>mod.2 [-50.6, 10.2]</td><td>mod.2 [-67.5, 12.0]</td><td>mod.2 [-53.4, 10.4]</td></tr><tr><td>7</td><td>mod.2 [-48.2, 10.1]</td><td>mod.6 [-51.7, 10.4]</td><td>mod.6 [-51.4, 10.1]</td><td>mod.6 [-68.7, 12.0]</td><td>mod.6 [-54.6, 10.4]</td></tr><tr><td>8</td><td>mod.0 [-48.5, 10.0]</td><td>mod.1 [-52.4, 10.7]</td><td>mod.1 [-53.4, 10.5]</td><td>mod.1 [-71.7, 12.4]</td><td>mod.1 [-57.2, 10.7]</td></tr><tr><td>9</td><td>mod.6 [-49.1, 10.1]</td><td>mod.0 [-52.7, 10.7]</td><td>mod.0 [-53.6, 10.6]</td><td>mod.0 [-72.3, 12.4]</td><td>mod.0 [-57.9, 10.7]</td></tr><tr><td>10</td><td>mod.4 [-49.1, 10.1]</td><td>mod.4 [-53.7, 10.6]</td><td>mod.4 [-54.0, 10.5]</td><td>mod.4 [-72.7, 12.3]</td><td>mod.4 [-58.4, 10.7]</td></tr></tbody></table>
will support the two models of being different (different in terms of predictive performance).

In Table 3, the models in the first five positions are the same across all the exposure times (i.e., *mod.3*, *mod.5*, *mod.7*, *mod.8* and *mod.9*). All these models have in common an interaction between air temperature and sex (see Table 2). This provides a clear indication that this interaction should be considered from a predictive point of view. However, since $|elpd_{diff}| < 4$ overall, we consider these five models to have equal predictive performance. All models in the last five positions (from 6 to 10) have at least $|elpd_{diff}| > 4 se_{diff}$; this suggests that all these models differ from the reference models in predictive ability.

Among the models in the first five positions, we have selected *mod.5* for the following two reasons:

1. This model was selected four out of five times as the reference model.

2. Considering the interaction between air temperature and relative humidity is a sensible choice because as temperature increases, relative humidity will have a greater effect on the evaporative heat loss necessary to maintain body temperature.

The results are presented and discussed only for the data from exposure at time 5 (i.e., at 3.0 h), where it can be assumed that the subjects were at or near thermal equilibrium with their surroundings.

Family: cumulative
Links: mu = probit; disc = identity
Formula: vote = 1 + c.air.temp + c_rel.hum + female + c.air.temp:c_rel.hum + c.air.temp:female + (1 | condition_id)
Data: df.experiment (Number of observations: 1600)
Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
total post-warmup draws = 12000

Multilevel Hyperparameters:
~condition_id (Number of levels: 160)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.42 0.05 0.34 0.51 1.00 4507 7293

Regression Coefficients:

<table><thead><tr><th></th><th colspan="2">Estimate</th><th colspan="2">Est Error</th><th colspan="2">l-95% CI</th><th colspan="2">u-95% CI</th><th colspan="2">Rhat</th><th colspan="2">Bulk ESS</th><th colspan="2">Tail ESS</th></tr></thead><tbody><tr><td>Intercept[1]</td><td>-2.85</td><td>0.09</td><td>-3.04</td><td>-2.67</td><td>1.00</td><td>6651</td><td>8435</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept[2]</td><td>-1.83</td><td>0.08</td><td>-2.00</td><td>-1.67</td><td>1.00</td><td>8180</td><td>9243</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept[3]</td><td>-0.78</td><td>0.07</td><td>-0.92</td><td>-0.64</td><td>1.00</td><td>10973</td><td>10004</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept[4]</td><td>0.93</td><td>0.08</td><td>0.79</td><td>1.08</td><td>1.00</td><td>11244</td><td>9493</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept[5]</td><td>2.00</td><td>0.09</td><td>1.84</td><td>2.17</td><td>1.00</td><td>10460</td><td>9797</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept[6]</td><td>3.14</td><td>0.10</td><td>2.94</td><td>3.34</td><td>1.00</td><td>10058</td><td>9354</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>c.air.temp</td><td>0.32</td><td>0.01</td><td>0.30</td><td>0.34</td><td>1.00</td><td>7644</td><td>8011</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>c_rel.hum</td><td>0.07</td><td>0.01</td><td>0.05</td><td>0.09</td><td>1.00</td><td>6982</td><td>8041</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>female</td><td>-0.20</td><td>0.06</td><td>-0.31</td><td>-0.09</td><td>1.00</td><td>16115</td><td>9011</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>c.air.temp:c_rel.hum</td><td>0.01</td><td>0.00</td><td>0.00</td><td>0.01</td><td>1.00</td><td>9332</td><td>9217</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>c-air.temp:female</td><td>0.10</td><td>0.01</td><td>0.08</td><td>0.12</td><td>1.00</td><td>17123</td><td>8610</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Further Distributional Parameters:

Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
disc 1.00 0.00 1.00 1.00 NA NA NA
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS and Tail_ESS are effective sample size measures, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1).

<caption>Figure 1: *brms* output for *mod.5*.</caption>

The *brms* output for *mod.5* is presented in Figure 1. For all the model's parameters, *Rhat* is less than 1.01, and ESS (*Bluk_ESS* and *Tail_ESS*) exceeds 400, indicating good MCMC convergence and stable uncertainty estimates, respectively (Aki Vehtari et al., 2021). These checks were also performed for the other models but were not reported here for brevity.

## Results

In this section, the results are presented and discussed.

## Posterior predictive check

Figure 2 shows the posterior predictive check, where many simulated data are generated under the fitted model and compared with the observed ones. This check is used to look for ‘systematic discrepancies that would be meaningful to address’ (Kruschke, 2015). This method has the disadvantage of assessing a model using the same data employed to determine its parameters. Unsurprisingly, the model can predict the data used to fit the parameters, but this basic test is helpful because it will fail when the model’s assumptions are severely violated. Figure 2 shows that males’ and females’ simulated thermal sensation votes agree with the observed ones, implying that the model fits the data well.

<figure><img src="image_4.png" /><figcaption>Figure 2: Posterior predictive checks comparing the distribution of TSV for males and females (yellow and green lines) to predictive replications using all the drawn from the posterior distribution (grey bands) of replicated data from the fitted model.</figcaption></figure>

## Mean probability and comparison with PMV

An ordinal model estimates a vector of probabilities, one for each category of the TSV (i.e., seven in this case). Given these probabilities, we calculate its expected value (or mean) as follows:

$$ \text{MeanPr} = \sum_{1}^{K} \pi_{k} k \quad (5) $$

where $\pi_k$ is the probability of a specific category $k$, $k \in \{1, \dots, K\}$. Although discussing the mean in an ordinal distribution often lacks relevance, the expected TSV is helpful since it can be compared directly with Fanger’s PMV.

Figure 3 describes the posterior distribution (summarised by its mean and 95 % credible intervals) for the expected
<figure><img src="image_5.png" /><figcaption>Figure 3: Posterior distribution (summarised by its mean and 95 % credible intervals) for the expected TSV. The black lines represent the PMV (assuming $ta = tr, rh = \{30; 50; 70\}, va = 0.1, Icl = 0.6, \text{and } M = \{1.0; 1.1; 1.2\}$).</figcaption></figure>

TSV given a range of temperatures (16–36°C) and three relative humidity levels for males and females. The black line represents the PMV calculated assuming $ta = tr, rh = \{30; 50; 70\}, va = 0.1, Icl = 0.6, \text{and } M = \{1.0; 1.1; 1.2\}$. As mentioned previously, in a multilevel ordinal model, two types of prediction can be computed depending on how the varying intercepts are handled. Considering the cluster-specific approach (Figure 3 top row), we can see that the expected TSV for males and females is different and that this difference is larger at lower temperatures. A comparison with the calculated PMV (black lines) shows that at lower temperatures, the expected TSV for males is better predicted using a value of 1.1 met while for females using 1.0 met. In contrast, at higher temperatures, the calculated PMV overpredicts the TSV for both sexes independently of the met value, but this difference attenuates at higher relative humidity and is lower for females across humidity. The overall picture changes when considering the sample-averaged approach (Figure 3 bottom row). As previously mentioned, the posterior distribution of the sample-averaged approach can be interpreted as the expected TSV across conditions (i.e., groups) on average. The expected TSV for males and females maintains the trend previously described with the cluster-specific approach, but now the posterior distributions are broader since they represent the variation across conditions. Considering these posterior distributions, it is possible to calculate the probability that the expected TSV is above or below the PMV or within a range of practical equivalence. In this study, we defined this region considering a ± 0.25 scale unit around the calculated PMV. Humphreys and Nicol (2002) proposed this value to evaluate the difference between observed thermal sensation and PMV.

Figure 4 shows these probabilities for males and females across three temperatures (i.e., 20, 26 and 32 °C) with relative humidity equal to 50 % and two metabolic rate values (i.e., 1.0 and 1.1 met). At 20 °C and 1.0 met, males

<figure><img src="image_6.png" /><figcaption>Figure 4: Posterior distributions for the expected TSV for males and females calculated with the sample-averaged approach. The dashed black lines represent the PMV (assuming $ta = tr, rh = 50, va = 0.1, Icl = 0.6, \text{and } M = \{1.0; 1.1\}$) while the red region around it is the range of practical equivalence (± 0.25 PMV).</figcaption></figure>
are more likely than females to have their expected TSVs exceed the PMV, while females are more likely to fall within the ± 0.25 PMV range. At 1.1 met, males still have a higher probability of exceeding the PMV, but their TSVs are more likely to fall within the ± 0.25 PMV range than females. Consequently, considering the same met value for males and females at lower temperatures can result in higher inaccuracies in the calculated PMV, where using a higher metabolic rate increases the underprediction of the expected TSV across conditions for females. At 26 °C, independently of the metabolic rate used, the probability of exceeding the PMV decreases for both sexes, with males more likely to fall within the ± 0.25 PMV range. At 32 °C, males are less likely than females to exceed the PMV, and males' probability of falling within the ± 0.25 PMV range is also lower. In the PMV model, the sweat evaporative heat loss only depends on metabolic rate and is assumed to be zero for 1.0 met; consequently, at higher temperatures, the probability that the expected TSV across conditions is within the ± 0.25 PMV for different metabolic rates changes slightly for both males and females. However, the PMV better predicts the expected TSV for females.

<figure><img src="image_7.png" /><figcaption>Figure 5: Contrast between the posterior distribution for the expected TSV for males and females.</figcaption></figure>

In Figure 4 the posterior distribution for the expected TSV for males and females overlaps substantially. However, asserting no difference between males and females because their posterior distributions overlap would be incorrect. It is never legitimate to compare overlaps in distribution; one must compute the contrast (i.e., difference) between them. Figure 5 shows the contrast between the posterior distributions for the expected TSV for males and females (the respective posterior distributions are shown in Figure 4). For example, although the posterior distributions for the expected TSV for males and females at 26 °C overlap almost completely (see Figure 4), their contrast is not centred at zero. Only 1% of the expected TSV across conditions is expected to be higher for females than males, with a mean difference of about +0.1 scale unit. Similar behaviour can also be observed for the other two temperatures. Although the posterior distributions for the expected TSV for males and females overlap at 20 and 32 °C, their contrast is far from

zero. At 20 °C, 0% of the expected TSV for females is expected to be higher than that of males (mean difference of about +0.60 scale unit). At 32 °C, 100% of the expected TSV for females is expected to be higher than that of males (mean difference of about -0.35 scale unit). Therefore, asserting no difference between males and females because of the overlap in distribution would have resulted in an incorrect conclusion. In this case, it is not appropriate to compare overlap in distribution because of how the posterior distribution for the expected TSV for the cluster-specific approach is calculated. In this approach, the posterior distribution is calculated by fixing the varying intercept $\alpha_{COND[i]}$ equal to 0. In Figure 3 (top row), the expected TSV for males and females are different, and this difference changes based on temperature and relative humidity. For example, considering 50% relative humidity, the expected TSV for males is higher than that of females at 20 °C, similar at 26 °C, and lower at 32 °C. Calculating the posterior distribution for the expected TSV for a different fixed value of $\alpha_{COND[i]}$ will result in shifting the lines in Figure 3 (top row) upwards or downwards by a constant amount (where its magnitude and direction depend on the specific value of $\alpha_{COND[i]}$). However, the relative difference between males and females does not change. The sample-averaged approach, average over the posterior distribution for the expected TSV calculated using different $\alpha_{COND[i]}$, but this process does not change the relative difference between males and females for each $\alpha_{COND[i]}$.

## Cumulative probabilities and transition temperatures

Figure 6 shows the cumulative probabilities for both males and females in the sample-averaged approach. The cumulative probabilities for males span a wider temperature range than those of females, indicating that females have more stringent conditions than males,

<figure><img src="image_8.png" /><figcaption>Figure 6: Posterior distribution (summarised by its mean and 95% credible intervals) for cumulative probabilities for males and females using the sample-averaged approach.</figcaption></figure>
particularly for lower temperatures. The horizontal dashed black line intersects the cumulative probabilities in their middle points (black dots), and the projections of these points onto the x-axis (vertical dashed black lines) indicate the temperatures at which voting lower and equal or higher to a category *k* is equally likely (i.e., Pr(*Y*<sub>*i*</sub> ≤ *k*) = Pr(*Y*<sub>*i*</sub> > *k*) = 50%). For example, the green line in Figure 6 represents the probability of voting higher than category 3 ‘slightly cool’ (i.e., Pr(*Y*<sub>*i*</sub> > 3)). The temperatures corresponding to the intersection of this curve with the 50% probability are 22.9°C and 24.0°C for males and females, respectively. Those temperatures represent the point estimates of the transition temperatures, that is, the temperature at which it is equally probable to vote equal or lower than category 3 ‘slightly cool’ (i.e., Pr(*Y*<sub>*i*</sub> ≤ 3)) or vote higher (i.e., Pr(*Y*<sub>*i*</sub> > 3)).

<figure><img src="image_9.png" /><figcaption>Figure 7: Posterior distribution of the transition temperatures for males and females using the sample-averaged approach.</figcaption></figure>

Figure 7 shows the posterior distribution of the transition temperatures for males and females. Although these posterior distributions are highly variable, females exhibit a greater response to cold and hot environments than males.

## Discussion

Ballantyne et al. (1977) re-analysed the KSU data using probit analysis, a more appropriate approach than Rohles’s (1971) linear regression. Ballantyne et al. found that males exhibit a weaker response to cold and hot environments than females and consider a broader range of temperatures acceptable. However, since probit analysis can be applied only to binary response variables, they collapsed the thermal sensation into two groups of adjacent categories (i.e., category 1 vs the others, category 1 and 2 vs the others, etc.) to analyse the data. Collapsing categories suffer from loss of information, particularly when dichotomising. Additionally, this approach does not consider the clustering in the data, which could lead to an underestimate of the standard errors and, consequently, an overstatement of statistical

significance. Our re-analysis used a multilevel ordered multinomial model to overcome these limitations.

A comparison of the expected TSV with the PMV shows that the difference between these two quantities varies depending on the metabolic rate used as an input for calculating the PMV. At lower temperatures, the expected TSV for males is better predicted using a value of 1.1 met, while for females it is better predicted using 1.0 met. In contrast, the calculated PMV overpredicts the expected TSV for both sexes at higher temperatures independently of the met value, but this difference attenuates at higher relative humidity and is lower for females across humidity. Instead, if a metabolic rate of 1.2 met is used as an input for calculating the PMV, the difference between the expected TSV and the PMV would increase on the cooler sensations (especially for females) while not affecting the difference in the warmer ones. This result is consistent with previous studies that found a larger discrepancy between the observed TSV and the PMV on the cooler sensations than the warmer ones (Cheung et al., 2019; Humphreys & Fergus Nicol, 2002). Nevertheless, it cannot be generally concluded that calculating the PMV with different metabolic activities for males and females will result in more accurate predictions in cold and warm environments, since sex-based differences are not limited to different metabolic rates. Additionally, the PMV model does not consider the human thermoregulatory system; in such environments, automatic thermoregulatory responses will take place and should be accounted for.

Although this study leverages experimental data with a large sample size and a wide range of conditions, it was impossible to define the best model from a predictive point of view since five candidates’ models showed an equally good predictive performance. Additionally, the sample-averaged approach shows high variability in the expected TSV. This highlights how the mean thermal sensation can vary substantially (up to 1 scale unit) even in controlled conditions. This should be taken as an indication of the inherent challenge of detecting differences and interactions. Given the presence of even more uncontrolled factors in field studies, this task would inevitably be more challenging. Defining a range of practical equivalence would be necessary to properly evaluate differences (e.g., sex-based differences), especially if statistical significance is used since the latter is not related to practical significance. If the Bayesian approach is used, it will be possible to use this range to calculate the probability that the difference of interest is within in, given some conditions (e.g., temperature, relative humidity, etc.).

## Conclusion

Understanding and considering sex-based differences in perception of the indoor (thermal) environment is essential to developing inclusive standards and enhancing energy efficiency while ensuring occupant well-being. This is particularly relevant considering that, to achieve net zero by 2050, the IEA (2021) report identified an average space heating and cooling temperature of 19–20 °C and 24–25 °C, respectively, as required behavioural
changes for the building sector by 2030. In this study, we re-analysed existing experimental data, comprising 1,600 college-age students (800 males and 800 females) exposed to a wide range of steady-state thermal environments for three hours. The results show that, from a predictive point of view, despite high variability, considering the interaction between sex and temperature better predicts the expected thermal sensation vote. Although it is inappropriate to generalise outside the data sample (e.g., older age) and other experimental conditions (e.g., different clothing levels, non-steady-state), for the conditions described in this study, it appears likely that taking measures to minimise the probability of ‘cold’ and ‘hot’ thermal sensations for females will result in a comfortable environment for all.

## References

Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, & Paul-Christian Bürkner. (2021). Rank-Normalization, Folding, and Localization: An Improved $\hat{R}$ for Assessing Convergence of MCMC (with Discussion). *Bayesian Analysis*, 16(2), 667–718.
https://doi.org/10.1214/20-BA1221

Ballantyne, E. R., Hill, R. K., & Spencer, J. W. (1977). Probit analysis of thermal sensation assessments. *International Journal of Biometeorology*, 21(1), 29–43.
https://doi.org/10.1007/BF01552964

Bürkner, P.-C. (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. *Journal of Statistical Software*, 80(1), 1–28.
https://doi.org/10.18637/jss.v080.i01

Cheung, T., Schiavon, S., Parkinson, T., Li, P., & Brager, G. (2019). Analysis of the accuracy on PMV – PPD model using the ASHRAE Global Thermal Comfort Database II. *Building and Environment*, 153, 205–217.
https://doi.org/10.1016/j.buildenv.2019.01.055

Fanger, P. O. (1970). *Thermal comfort: Analysis and applications in environmental engineering*. Danish Technical Press.

Favero, M., Luparelli, A., & Carlucci, S. (2023). Analysis of subjective thermal comfort data: A statistical point of view. *Energy and Buildings*, 281, 112755.
https://doi.org/10.1016/j.enbuild.2022.112755

Gelman, A., & Hill, J. (2006). *Data Analysis Using Regression and Multilevel/Hierarchical Models*. Cambridge University Press; Cambridge Core.
https://www.cambridge.org/core/product/32A29531C7FD730C3A68951A17C9D983

Graham, L. T., Parkinson, T., & Schiavon, S. (2021). Lessons learned from 20 years of CBE’s occupant surveys. *Buildings and Cities*.
https://doi.org/10.5334/bc.76

Humphreys, M. A., & Fergus Nicol, J. (2002). The validity of ISO-PMV for predicting comfort votes in

every-day thermal environments. *Special Issue on Thermal Comfort Standards*, 34(6), 667–684.
https://doi.org/10.1016/S0378-7788(02)00018-X

IEA. (2021). *Net Zero by 2050*. https://www.iea.org/reports/net-zero-by-2050

Kay, M. (2023). *tidybayes: Tidy Data and Geoms for Bayesian Models*. https://doi.org/10.5281/zenodo.1308151

Kruschke, J. K. (Ed.). (2015). *Doing Bayesian Data Analysis* (2nd ed.). Academic Press.

Liddell, T. M., & Kruschke, J. K. (2018). Analyzing ordinal data with metric models: What could possibly go wrong? *Journal of Experimental Social Psychology*, 79, 328–348.
https://doi.org/10.1016/j.jesp.2018.08.009

Posit team. (2024). *RStudio: Integrated Development Environment for R* [Computer software]. Posit Software, PBC. http://www.posit.co/

R Core Team. (2023). *R: A Language and Environment for Statistical Computing* [Computer software]. R Foundation for Statistical Computing. https://www.R-project.org/

Rohles, F. H. (1970). *Thermal Sensations of Sedentary Man in Moderate Temperatures* (Technical Report AD756492; p. 131).

Rohles, F. H. (1971). Thermal Sensations of Sedentary Man in Moderate Temperatures. *Human Factors*, 13(6), 553–560.
https://doi.org/10.1177/001872087101300606

Sivula, T., Magnusson, M., Alonzo Matamoros, A., & Vehtari, A. (2020). Uncertainty in Bayesian Leave-One-Out Cross-Validation Based Model Comparison. *arXiv e-Prints*, arXiv:2008.10296.
https://doi.org/10.48550/arXiv.2008.10296

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. *Statistics and Computing*, 27(5), 1413–1432.
https://doi.org/10.1007/s11222-016-9696-4

Wang, Z., de Dear, R., Luo, M., Lin, B., He, Y., Ghahramani, A., & Zhu, Y. (2018). Individual difference in thermal comfort: A literature review. *Building and Environment*, 138, 181–193.
https://doi.org/10.1016/j.buildenv.2018.04.040

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., ... Yutani, H. (2019). Welcome to the tidyverse. *Journal of Open Source Software*, 4(43), 1686.
https://doi.org/10.21105/joss.01686