The most important impacts on occupant energy behaviors are environmental factors such as indoor air temperature (Parsons, 2002; Cui et al., 2013; ANSI/ASHRAE, 2017). Several studies have reported that indoor air temperature has a substantial impact on occupant energy behaviors like thermostat adjustments for heating/cooling (Bae and Chun, 2009; Lin et al., 2016), turning on/off AC (Bae and Chun, 2009; Ren et al., 2014), window opening/closing (Haldi and Robinson, 2008, 2009; Yun and Steemers, 2008; Yun et al., 2008, 2009), drinking cold/hot beverages (Haldi and Robinson, 2008, 2011), desk fan usage (Haldi and Robinson, 2008), adjusting body posture (Schweiker et al., 2017), and adjusting clothing (Haldi and Robinson, 2008). Occupant thermal sensations under different indoor air temperatures also affect their energy behavior (Bonte et al., 2014). Individual factors of occupants, which are further classified into psychological and social factors, significantly impact occupant energy behavior and thermal comfort, and they have been thoroughly researched. Individual factors affecting occupant energy behavior include users' knowledge of energy matters, gender, employment, family size, and ethnicity (Wei et al., 2013; Martinaitis et al., 2015). People's education levels and psychology, such as their attitudes toward energy consumption, also impacted their energy related behaviors (Janda, 2011). However, the above studies followed conventional modeling approaches, and the results are interpreted as correlations rather than causation.

Furthermore, to analyze the factors affecting human behaviors, various theories such as the Theory of Planned Behavior (TPB) (Ajzen, 1991) have been applied widely to address individuals' energy behavior (Vlek and Steg, 2007; Tetlow et al., 2015). The Theory of Planned Behavior is a well-known theory that predicts the intention of a person to engage in a behavior. The three major elements of this theory consist of the attitudes toward performing the behavior (i.e., the point to which a behavior is regarded favorably or unfavorably), the social norms or pressure to perform the behavior (i.e., social pressure coming from people's immediate social network), and the perceived behavioral control (i.e., beliefs about the existence of factors that can aid or hinder the success of the behavior). A combination of these elements determines the behavioral intention, as well as the final action. For example, in building energy behavior research, TPB factors were used to explain the differences in the intentions of employees' to perform pro-environmental activities like recycling, switching off electrical appliances while leaving their desks, and using teleconference instead of traveling (Greaves et al., 2013). Another study indicated that viewpoints toward conserving energy together with perceived behavioral control strongly determine the propensity to conserve energy within low-income houses throughout the USA (Chen et al., 2017). Therefore, in this study, the factors of TPB will also be applied to analyze occupant energy behavioral intentions as a proximal antecedent of actual energy behavior in IVEs.

Researchers analyzed occupant thermal states since they also significantly affect their energy behavior (Hong et al., 2015). The thermal states of the occupants usually include sensation, acceptability, and comfort. Thermal sensation is defined as the apparent sensation that people have in relation to their surrounding environment and is usually measured using a 7-point Likert scale (−3-cold sensation, +3-hot sensation) of ASHRAE Standard 55 (de Dear and Brager, 1998; ASHRAE, 2013). Whereas thermal comfort refers to “the condition of mind which expresses satisfaction with the thermal environment” (ASHRAE, 2013), and thermal acceptability defines whether a particular environmental condition is acceptable or not. Comfort and acceptability are also measured using surveys such as 6-point Likert scales (−3-totally unacceptable, +3-totally acceptable; −3-very uncomfortable, +3-very comfortable) (Zhang, 2003; Goto et al., 2006; Zhang and Zhao, 2008; Zhang et al., 2016). Similar to the factors affecting occupant energy behavior, the occupant thermal states are also affected by environmental, individual, and physiological factors. Factors such as indoor air temperature and overall quality of environment were correlated with thermal sensation and comfort (Clausen et al., 1993; Alm et al., 1999). Some studies have indicated that there are variations between different groups of people, i.e., gender, age, and body mass index (BMI) on thermal preferences (Lan et al., 2008; Schweiker and Shukuya, 2009, 2010; Chow et al., 2010; Indraganti and Rao, 2010; Indraganti et al., 2015; Rupp et al., 2018). Furthermore, physiological factors such as skin temperatures and heart rates are the most important indicators of thermoregulation of the human body (Arens and Hui, 2006; Yao et al., 2008; Choi et al., 2012; Choi and Yeom, 2017, 2019; Nkurikiyeyezu et al., 2018). These physiological factors will help maintain body temperature equilibrium in different thermal environments and play a vital role in determining thermal states, which in turn affects occupants' energy behaviors (Bulcao et al., 2000).

The virtual experience factors are a collection of characteristics that indicate the quality of experience in immersive virtual environments (IVEs). These characteristics are usually grouped into three main factors, i.e., presence, cybersickness, and immersive tendencies. Presence is described as being the “psychological state of 'being there' mediated by an environment that activates our senses, captures our attention, and promotes our active participation” (Ijsselsteijn et al., 2000; Witmer et al., 2005). Igroup Presence Questionnaire (IPQ) (Schubert et al., 1999) is a frequently used questionnaire to quantify the users' presence. This 13-item questionnaire has four factors, i.e., general presence, spatial presence, involvement, and experienced realism. The questions are gauged on a 5-point Likert scale that ranges between 1-strongly disagree; 5-strongly agree. On the other hand, cybersickness is a collection of unpleasant sensations that impair users' interaction with a virtual world (LaViola, 2000). The Simulator Sickness Questionnaire (SSQ) (Kennedy et al., 1993) has been extensively utilized for determining motion sickness. This questionnaire contains 16 items that include four factors, i.e., nausea, oculomotor, disorientation, and total cybersickness. The total cybersickness is the summation of the responses of the remaining three factors. All the questions in the SSQ are gauged on a 4-point scale that ranges between 0-none; 3-severe. Immersive tendencies are defined as a users' proclivity to become absorbed in a simulated world (Slater and Wilbur, 1997). They are thought to be one of the crucial factors affecting a users' sensation of presence and cybersickness, with a high level of immersion predicting a high level of presence (Lessiter et al., 2001; Banerjee et al., 2002) and a higher level of cybersickness (Bangay and Preston, 1998). Immersive tendencies are usually measured using the ITQ questionnaire (Witmer and Singer, 1998). It is a questionnaire of 18 items that consists of four factors, i.e., focus (indicates the level of focus in the virtual world), involvement (indicates the level of involvement in the virtual world), emotion (indicates the level of emotion during and after exposure to the virtual world), and entertainment (indicates the level of entertainment while in the virtual world). The responses to these factors are on a 7-point Likert scale (1 = never; 4 = occasionally; 7 = often). The total number of points represents a users' immersive tendency.

Past evidence suggests that the aforementioned virtual experience factors will impact the users' physiological factors. For example, various studies have reported that presence influences the degree of stress in a virtual world (Bouchard et al., 2008), and cybersickness influences heart rates (Nalivaiko et al., 2015), respiratory rates (Dennison et al., 2016), skin conductance (Gavgani et al., 2017), and skin temperatures (Kim et al., 2005). Individual factors like age, gender, and ethnicity also influence presence and cybersickness (Klosterhalfen et al., 2005; Schuemie et al., 2005; Felnhofer et al., 2012; Paillard et al., 2013). Furthermore, the literature presented in Key Factors Affecting Occupant Energy Behavior and Key Factors Affecting Thermal States sections suggests that the physiological and individual factors influence occupant energy behaviors and thermal states. Consequently, a pathway may exist from virtual experiences to occupant energy behaviors and thermal states through individual and physiological factors. Therefore, when using IVEs for occupant energy behavior and thermal state experiments, it is necessary to investigate whether the virtual experience factors would influence those experiments outcomes.

The first step of a causal inference process is to select and describe the factors of a study. Two types of factors are usually included in a causality study. They are the treatment factor (or cause) and the outcome factor (or effect) (Tsapeli and Musolesi, 2015). A treatment factor is an influential or independent factor that impacts the values of the outcome factor; whereas, the outcome factor or effect is a dependent factor that can be modified by varying the treatment factor. A series of questions can be asked to identify an initial set of potential treatment factors, such as (1) what circumstances lead to the occurrence of the problem or the outcome? (2) What other problems exist in the context of the core problem's occurrence? On the basis of such questions, a literature search can be performed to uncover an initial set of treatment factors for analyzing the outcome (Petersen and van der Laan, 2014).

Once the treatment factors are identified, typically, they are coded in binary format, indicating whether or not a research item is subjected to treatment. However, if various levels of treatments are considered, a treatment factor may also be categorical. Continuous treatment or outcome factors are often discretized into categorical factors because the algorithms performing approximate causal inference function well on categorical data (Spirtes et al., 2000; Jensen et al., 2009; Nanda et al., 2016; Tosun et al., 2017). Discretization also makes the findings easier to interpret.

A causal model can define a formal representation of potential causal relationships between the key factors. A causal model is an arithmetic concept that explains the causal relationships between factors in quantitative terms (Guo et al., 2018). Different types of causal models exist, such as Granger causality (Granger, 1969) and potential outcome framework (Rubin, 1974; Maldonado and Greenland, 2002), among which structural causal models (SCMs) or causal graphs have demonstrated their ability to identify the causal relationships between various factors of a study (Pearl, 1988; Judea and Dana, 2018). SCMs provide a visualization of the interactions between the factors by representing the factors as nodes and connecting those nodes with arcs to show the interactions (Glymour et al., 2019). They can naturally deal with missing data. They offer a strategy for avoiding data overfitting, and they enable data to be combined with prior literature knowledge to create more accurate models (Heckerman, 1995; Scutari et al., 2017). Furthermore, SCMs can readily be coupled with decision-making tools to assist management (Uusitalo, 2007). As such, SCMs were chosen in this study to identify the factors causing energy behavioral intentions and thermal states.

An SCM G = (F, A) is represented by a directed acyclic graph (DAG) (i.e., graphs with no directed paths that begin and end at the same node) (Pearl, 2009b). An SCM (G) depicts the causal relationships between factors, with F being a set of nodes (or vertices) and A ⊆ F × F is the set of arcs (or edges). Each node in an SCM represents a treatment factor, outcome factor, and other observed or unobserved factors. The arcs represent the causal links between them. Figure 1 shows an example of a structural causal model representing the causal dependencies between OUTCOME and a collection of factors (FAC or Treatment) that may influence OUTCOME. The nodes have a hierarchical structure which is determined by the way the arcs are directed. For example, in Figure 1, X₁ and X₃ are connected by an arc, then node X₁ is the ancestor or immediate parent of X₃ and X₃ is the descendant or child of X₁. A directed arc X₃ → OUTCOME means that if the distribution of X₃ is altered, it will alter the distribution of the OUTCOME. However, if there is an absence of arc between X₃ and OUTCOME, it indicates that any change in the distribution of X₃ will not change the distribution of OUTCOME.

One of the advantages of SCMs is identifying back-door paths and confounding factors that are useful in blocking non-causal paths during causal effect estimation. The following sub-sections provide a brief description of confounding and learning of SCMs.

The causal effect exercised by a treatment factor X over an outcome factor Y refers to the magnitude of change in the probability distribution associated with the factor Y with respect to change in the probability distribution associated with X. In other words, it is a measure of how different the outcome is when the subject received the treatment compared to that when the subject did not receive the treatment (Rubin, 1974; Maldonado and Greenland, 2002). For example, in Figure 1, assume X₂ as indoor room temperature with two treatment levels, i.e., hot and cold temperature and OUTCOME as occupant's skin temperature. We would like to know what is the causal effect of hot indoor temperature versus cold indoor temperature on an occupant's skin temperature. Since there are two treatment levels, there will be two outcomes, i.e., the occupant's skin temperature under hot indoor temperature and their skin temperature under cold indoor temperature. The causal effect of the indoor room temperature on the occupant's skin temperature is then determined by taking the difference between those two outcomes. The above definition is known as individual causal effect (ICE), which applies to a single individual. However, the average causal effect (ACE) is used based on ICE to find the causal effects on a target population. ACE can be defined as the approximation of ICE over the entire target population.

There are two types of causal effects in an SCM; direct and indirect causal effects. Direct causal effects are those that occur when they go directly from one factor to another. Indirect causal effects occur when one or more factors mediate the relationship between two factors. For example, in Figure 1, X₂ affects OUTCOME both in a direct and indirect way via its direct impact on X₃. The magnitude of the direct causal effect can be determined by computing the coefficient between X₂ and OUTCOME. The indirect effects can be determined by computing the product of the coefficients between X₂ and X₃ and between X3 and OUTCOME. The combined effect is obtained by adding the direct and indirect impacts of X₂ on OUTCOME. Regression analysis is one way of determining the coefficients as an average causal effect estimate by training the model with the outcome as dependent factor and treatments as independent factors (Bellemare et al., 2017; Luo et al., 2017). The coefficients can be either positive or negative. The positive coefficients imply that the event is more likely to occur at a treatment factor's group level than at its reference level; whereas, negative coefficients imply that the event is more likely to occur at the reference level of the treatment factor than at its group level.

Furthermore, if there are back-door paths and confounding factors between treatment and outcome, those back-door paths must be blocked by conditioning on the confounding factors. Common methods of conditioning confounding factors and estimating causal effects for SCMs include Parametric G-computation (Zhang et al., 2018), Targeted Maximum Likelihood Estimation (TMLE) (van der Laan and Rubin, 2006), and Inverse Probability Weighting (IPW) (Rosenbaum and Rubin, 1983).

The performance and the accuracy of the causal models can be validated using several metrics such as receiver operating characteristics, error rates and scoring rules, and cross-validation (Marcot et al., 2006). In this study, cross-validation (Boyce et al., 2002) is chosen as the method for validating the SCMs (Barthelmes et al., 2017). Cross-validation decomposes the data into k disjoint subsets known as folds of approximately the same size. Each fold is considered unseen data, whereas the rest of k−1 folds are considered seen data. One trains the model on the seen data and validates against the unseen data to find the model's goodness-of-fit. This method reveals the degree to which a model is robust by preventing overfitting and how accurate the model is at classifying new or unique data.

The final step of the causal inference framework is to interpret the results. In a structural causal model (SCM), we can look at the direction of network arcs that flow from one factor to another to interpret the causal relationship between them. Furthermore, the causal effect estimate will indicate the strength of the causal relationships between the factors. Since the causal effects are computed using regression, they can be interpreted the same way the regression coefficients are interpreted. For example, an SCM was developed to study substance abuse, and with the help of regression coefficients, the authors found that taking marijuana causes the odds of taking heroin to be 23 times more than that of people who didn't take marijuana (Lewis and Kuerbis, 2016). In another study, an SCM was developed to identify causes that affect the results of coronally advanced flap procedure (Nieri et al., 2009). Factors such as patients' age and complete root coverage had a negative causal effect on their root canal sensitivity. In contrast, factors such as gender had a positive causal effect on deeper average recessions. Two alternative SCMs were developed to estimate the causal effect of exclusive breastfeeding (EBF) on early childhood caries (ECC) for Ugandan kids aged 5 (Birungi et al., 2017). Based on the regression coefficients, EBF showed considerable impact on ECC in the first SCM, and EBF with marital status showed significant impact on ECC in the second SCM.

The study was approved by the university's Institutional Review Board (IRB). Seventy-two participants were enlisted using posters and oral communication. Monetary incentives ($10/h) were offered for their participation in the study. The participants' demographic distributions are as follows:

The experiments (both IVE and in-situ sessions) were performed in a climate chamber using the experimental procedure described in previous studies (Rentala et al., 2021; Saeidi et al., 2021). The details of the climate chamber and its IVE model have also been explained in those earlier studies. Following the submission of their signed consent forms, the immersive tendencies and TPB questionnaires were sent to the participants. After completing these questionnaires, the participants were asked to arrive at the climate chamber dressed in a long-sleeved shirt or T-shirt, trousers, socks, and shoes. Then they were asked to sit in the resting area of the chamber (temperature and humidity set at 75°F/23.8°C and 55% respectively), where they were screened for cigarette or alcohol consumption and vigorous physical activity before the experimental sessions. This session took about 10 mins, and the participants were excluded from the experiments if any of the above were discovered. During this session, the participants' individual factors were also recorded. After the screening, they entered the control area of the chamber, sat in a chair, and got the temperature sensors taped to the surface of their skin.

The experiments were conducted from December 2019 to March 2020 following the experiment procedure shown in Figure 2. Each IVE and in-situ experiment had a cooling sequence or a heating sequence. The IVE and in-situ experiments were conducted on the same day, but the heating and cooling sequence experiments were conducted at least 2 weeks apart. To offset and reduce the order effect, participants were arbitrarily allocated to each of the experimental conditions. The IVE experiments were conducted using an HTC Vive HMD device where the 3D model of the chamber was displayed using Unreal Engine 4 (for details and pictures, refer to previous studies). During an IVE experiment, participants only view the chamber's virtual environment through an HMD (no interaction), while the climate chamber controls the environmental conditions. The heating sequence started from 65°F/18.3°C and then increased to 75°F/23.8°C and finally to 85°F/29.4°C. The cooling sequence started from 85°F/29.4°C and then decreased to 75°F/23.8°C and finally to 65°F/18.3°C. The psychrometric chart of ASHRAE (Bhattacharya and Milne, 2009) was used to select these step temperatures. In this chart, for both summer and winter seasons, 75°F/23.8°C falls in the comfort zone range, while 65°F/18.3°C and 85°F/29.4°C fall outside the comfort zone range. The humidity and CO₂ levels were maintained at 55% and 1,000 ppm during the experiments, respectively.

The indoor room temperature around the participants and their skin temperatures were continuously measured at the one-second interval. After each temperature step change, the participants' thermal state responses (i.e., thermal sensation, comfort, and acceptability) were recorded after the temperature stabilized. Throughout the experiments, the participants' were constrained to sedentary physical activity levels. Furthermore, the participants completed the presence and cybersickness questionnaires inside the chamber after an IVE experimental session. Finally, the participants' energy behavioral intentions (i.e., their three most preferred behavioral choices) were collected after exposure to the cool discomfort condition, i.e., 65°F/18.3°C in cooling sequence in both IVE and in-situ (referred to as Cool/IVE and Cool/in-situ conditions). The behavioral intentions were also collected after exposure to the warm discomfort condition, i.e., 85°F/29.4°C in heating sequence in both IVE and in-situ (referred to as Warm/IVE and Warm/in-situ conditions). The energy behavioral intentions are collected in discomfort conditions because the occupants usually react to thermal discomfort situations to maintain thermal comfort.

The raw data (i.e., indoor room and skin temperatures) collected from the sensors were averaged by taking into account data from the precise start time when the indoor room temperature stabilized at the target step temperature through the end of the thermal state questionnaires. Furthermore, to ensure that the actual range of the indoor room temperature represents the corresponding step temperature zone, the data was cleaned by comparing the averaged indoor room temperature between IVE and in-situ settings and removing those data points whose differences were more than ±3°F/1.7°C. This threshold was selected based on the fact that an individual may be unaware of a temperature increase of 4–5°C (7.2–9°F) when the rate of change of temperature is <0.5°C (0.9°F) per min (Darian-Smith and Johnson, 1977). In the climate chamber used in this study, a temperature change of about 0.5°C/1°F frequently takes 1 min or more. At this rate, it was assumed that the participants were unaware of at least a 3°F/1.7°C change in about 10 mins (10 mins is the estimated time to rise/fall for 5.6°C/10°F).

This study selected the key treatment factors that may affect the behavioral intentions and thermal states (Table 1) based on the current domain knowledge and the models reviewed from various peer-reviewed occupant energy behavior and thermal state literature (Background section). It must be noted that the participants' other individual factors such as age, BMI, ethnicity, education levels, and employment status were not included in the model because the samples in those factor groups were not evenly distributed.

The outcome factors of this study are the participants' energy behavioral intentions (EBI) and their thermal states, i.e., general thermal sensation (GTS), general thermal acceptability (GTA), and general thermal comfort (GTC).

As described in Experiment Procedure section, the energy behavioral intention data were collected outside the comfort zone range, i.e., at 65°F/18.3°C (cool discomfort condition) and 85°F/29.4°C (warm discomfort condition). In contrast, the thermal state data was collected under all three-step temperatures (both inside and outside comfort zone range). Thus, for analyzing energy behavioral intentions (EBI), the data is segmented based on cool and warm discomfort conditions in IVE and in-situ, and appropriate treatment factors such as gender, TPB, and virtual experience factors are assigned (Figure 4). For analyzing thermal states, the data is segmented based on IVE and in-situ conditions only, and factors such as indoor room temperature (all three levels), overall skin temperature, and virtual experience factors are assigned as treatment factors (Figure 5).

After factor assignment, the continuous data are discretized into categories, as shown in Table 2. Data discretization was performed because the algorithms used for structure learning and approximate causal inference work well with the categorical data. The entire discretization process is explained in detail in the following sections.

Based on the discussions presented Casual Model Specification section, the elicitation method, in conjunction with a PC-Stable algorithm, was used to learn the SCMs from the experimental data. First, the data were randomly divided into two datasets assuming a 50:50 ratio. The first 50% of the data is the seen data on which the SCMs were learnt, and the other 50% of the data is the unseen data that was used for validating the causal models learnt on the seen data. The following sub-sections describe the structure learning process in detail.

Given the causal relationships encoded in the SCMs, the average causal effects (ACE) of only those treatment factors that have a direct path to the outcome factors, i.e., direct average causal effects, were estimated using logistic regression on the seen data. The causal risk odds ratio obtained from the logistic regression coefficients is used to represent the average causal effect. If the causal risk odds ratio is not equal to 1, the factors have a causal effect. The logistic regression equation is given below:

Where, P(O_i = o|T_i) is the probability of the outcome Oi given the treatment factor Ti, α denotes the intercept, and β is the coefficient of the treatment factor.

Since there are no confounding factors between the treatment and outcomes, the average causal effects of treatment factors with only direct paths to the outcome were estimated by regressing the former directly on the latter. Since the outcome factors of EBI, GTA, and GTC had two levels, binomial logistic regression was used. In contrast, multinomial logistic regression (Menard, 2002) was used for estimating the average causal effect of treatment factors on GTS, which has three levels.

The SCM shown in Figure 6 suggests that in Cool/in-situ and Cool/IVE, there is a causal relationship from Control to EBI and Attitude to EBI, respectively. In contrast, SCM shown in Figure 7 suggests that in Warm/in-situ and Warm/IVE, there are no causal relationships from any of the treatment factors to EBI. In Figure 8, the SCM suggests that in in-situ and IVE conditions, there is a causal relationship from RoomTemp to GTS and from GTC to GTA. Tables 3, 4 show the average causal effects of those treatment factors on the outcomes. Since there are no paths from any of the treatment factors to EBI in Warm/in-situ and Warm/IVE, the average causal effects on EBI were not estimated in those conditions.

Validation is essential for testing the strength of causal models. In this study, a two-fold cross-validation process was applied where the data was randomly divided based on a 50:50 ratio into seen data (50%) and unseen data (50%). The choice of a two-fold with a 50:50 ratio was based on the need for simplicity and consistency in the results of the hypothesis tests, which depend on the sample size. The relationships identified by SCMs learnt on the seen data are cross-validated against the previously unseen data. The SCMs provide relationships between the factors in the form of testable implications, and these implications can be verified by performing hypothesis tests. The following steps describe the validation process:

Depending on the direction of the arcs, the implications test the null hypothesis that a factor does not affect another factor. The alternate hypothesis states that a factor affects another factor. The null hypothesis is accepted if the p-value is larger than α = 0.05; otherwise, it is rejected. Since all the factors are categorical and depend on the number of levels in a factor, binomial and multinomial logistic regression with Type II tests were used to test the hypotheses. The “Type II tests” test for the significance of each independent factor (treatment) on the outcome, under the assumption that there are no interactions between the factors (Langsrud, 2003). The Type II tests will provide the Wald chi-square statistic and its associated p-value.

A total of 22 testable implications were identified from the causal models, and Type II tests were performed for each of those implications. The validation results (p-values) are shown in Tables 5, 6. Out of 22 Type II tests, 19 were significant, i.e., the p-values were below 0.05. Thus the null hypothesis was rejected in all those cases. However, there were only three cases, one in Cool/In-situ, one in Warm/IVE, and one in IVE condition, respectively, where the p-values exceeded 0.05, indicating that those relationships were inconsistent against the unseen data. Therefore, the 2-fold cross-validation results revealed that the SCMs have a high accuracy of about 86.3% in identifying the causal relationships in the unseen data.