Estimation and Inference of Special Types of the Coefficients in Geographically and Temporally Weighted Regression Models

Geographically and temporally weighted regression (GTWR) models have been widely used to explore spatiotemporal nonstationarity where all the regression coefficients are assumed to be varying over both space and time. In reality, however, constant, only temporally varying, and only spatially varying coefficients might also be possible depending on the underlying effects of the explanatory variables on the response variable. Therefore, the development of inference and estimation methods for such special types of the coefficients is essential to the deep understanding of spatiotemporal characteristics of the regression relationship. In this article, an average-based approach, relying on a modified estimation of the conventional GTWR models, is proposed to calibrate the GTWR models with the special types of the coefficients, on which a statistical test is formulated to simultaneously infer constant, temporally varying, and spatially varying coefficients. The simulation study shows that the test method is of valid Type I error and satisfactory power and the average-based estimation method yields more accurate estimators for the special types of the coefficients. A real-life example based on Beijing house prices is given to demonstrate the applicability of the test and estimation methods as well as the extensibility of the test in model selection.

W ith the increasing availability of spatiotemporal data in reality, much attention has been paid to modeling spatiotemporal data for the analysis of spatial and temporal characteristics, among which geographically and temporally weighted regression (GTWR) models are perhaps one of the most popular tools for exploring spatiotemporal nonstationarity, a fundamental property of spatiotemporal data, of a regression relationship.
The GTWR models, originally proposed by Huang, Wu, and Barry (2010), are an extension of the geographically weighted regression (GWR) models (Brunsdon, Fotheringham, and Charlton 1996) created by incorporating the time dimension into the spatially varying coefficients and defining a spatiotemporal distance to calibrate the models the same way as in the GWR technique. Fotheringham, Crespo, and Yao (2015) proposed a different calibration procedure of the GTWR models where only the current and previous observations are used to estimate the spatiotemporally varying coefficients at the current time stamp with a more flexible way of selecting the spatial and temporal bandwidths. This estimation method is in line with the claim that "yesterday can exert influence on today's pattern, but the inverse may not be possible" (Dub e and Legros 2013, 4), and therefore the analysis result is more interpretable. Due to their wide application backgrounds, the GTWR models with their estimation methods have been extended in diverse ways. For example, B. Wu, Li, and Huang (2014) incorporated the spatial autoregressive structure of the response variable into the GTWR models and proposed geographically and temporally weighted autoregressive models with the two-stage least squares estimation procedure. Considering some explanatory variables could influence the response variable globally, J. P. Liu et al. (2017) proposed the so-called mixed GTWR models in which some coefficients are assumed to be constant and the others are retained to be spatiotemporally varying and employed a twostep method to estimate the constant and the spatiotemporally varying coefficients. Extending the multiscale GWR model in Fotheringham, Yang, and Kang (2017), C. Wu et al. (2019) formulated multiscale estimation for the GTWR models, in which the spatial and temporal bandwidths are separately selected for each spatiotemporally varying coefficient to uncover different spatial and temporal scales at which each explanatory variable operates. By optimizing a pair of spatial and temporal bandwidths at each time stamp, Z. Zhang et al. (2021) suggested a multiscale approach with a unilateral weighting scheme to calibrate the GTWR models. The other extensions of the GTWR models include geographically and circle-temporally weighted regression for capturing the periodicity in time dimension (Du et al. 2018) and the geographically weighted temporally correlated logistic regression to identify spatial and temporal correlations of binomial outcome data (Y. Liu et al. 2018). Since their inception, the GTWR models with their variants have been widely applied to a variety of practical areas for spatiotemporal data analysis, including geography (Chu, Kong, and Chang 2018), econometrics (Huang, Wu, and Barry 2010;B. Wu, Li, and Huang 2014;Fotheringham, Crespo, and Yao 2015;C. Wu et al. 2019;Q. Wang et al. 2021;Z. Zhang et al. 2021), environmental science (Guo et al. 2017;Peng et al. 2019), epidemiology (Y. Chen et al. 2021;Hong et al. 2021), and so on.
In the conventional GTWR model, all of the regression coefficients are assumed to be related to both space and time, implying that the effects of the explanatory variables on the response variable are all spatiotemporally varying. In reality, however, in addition to the case considered by J. P. Liu et al. (2017) where the effects of some explanatory variables are global, there is also the possibility that the effects of some explanatory variables are either spatially stationary or temporally stationary. Such multiple effects will result in constant, only temporally varying, and only spatially varying coefficients in the GTWR model, respectively, which we call the multieffect GTWR (MEGTWR) model. Therefore, the inference and estimation of such special types of coefficients provides a deeper understanding of spatiotemporal characteristics of the regression relationship, and also is conducive to the choice of a more appropriate model for a given spatiotemporal data set.
In the GTWR and the mixed GTWR models with their existing calibration methods, however, both spatially varying coefficients and temporally varying ones are all assumed as spatiotemporally varying coefficients to be estimated by the local fitting technique, resulting in such estimators that always vary over both space and time, no matter what the types of the underlying coefficients really are. This can lead not only to accuracy loss on the estimators of these special-type coefficients, but also inexact or even misleading interpretation on the spatiotemporal patterns of the regression relationship. Therefore, it is necessary to develop some new calibration methods for the MEGTWR model.
In practice, the inference of the special-type coefficients in the GTWR model is also indispensable. On this aspect, Huang, Wu, and Barry (2010) as well as Xuan, Li, and Amin (2015) proposed tests for global spatiotemporal stationarity, or global spatial or global temporal stationarity of the whole regression relationship, respectively, by comparing the goodness of fit between the respective null models and the GTWR model, where the F distribution approximation was used to derive the p values of the tests. Similar tests were also formulated by Xiao, Tian, and Guo (2014) where the bootstrap method was used to approximate the p values of the tests. Nevertheless, these tests only focus on inferring whether all of the regression coefficients in a GTWR model are constant, only spatially varying, or only temporally varying, which are all extremely special cases of the GTWR regression relationship and might rarely appear in reality. Hong et al. (2021) extended the bootstrap test for constant coefficients in the GWR models (Mei, Xu, and Wang 2016) to detect possible constant coefficients in a GTWR model, on which a mixed GTWR model was chosen and then estimated by the two-step procedure (J. P. Liu et al. 2017). In fact, due to their specific estimation approaches of the null models, neither of the foregoing tests can directly be extended to simultaneously infer the various special types of the coefficients in a GTWR model. Moreover, the recently developed multiscale estimation methods (C. Wu et al. 2019;Z. Zhang et al. 2021) for the GTWR models yield coefficient-specific spatial and temporal bandwidths that could provide some information for the types of the coefficients. That is, larger spatial or temporal optimal bandwidth sizes of a coefficient imply that the coefficient is possibly constant, spatially or temporally stationary. This judgment is rather intuitive, however, because it is difficult to set proper threshold values of the bandwidths for the judgment of the coefficient types. In summary, it is still necessary to develop a formal statistical test to simultaneously infer these special types of coefficients in a GTWR model.
In this article, an average-based estimation framework is first formulated to fit an MEGTWR model, on which a statistical test is proposed to simultaneously infer the special types of coefficients in a GTWR model. Specifically, we first combine the simplicity in selecting the bandwidths in Huang, Wu, and Barry (2010) and the interpretability in Fotheringham, Crespo, and Yao (2015) that the coefficients at the current time stamp are estimated by using the current and previous observations and propose a modified method to calibrate the conventional GTWR model. Based on the resulting coefficient estimators, the average-based methodology, suggested by W. Y. Zhang, Lee, and Song (2002) to estimate the constant coefficients in a semivarying coefficient model, is extended to derive the estimators of the special types of coefficients including constant coefficients and those only relating to space or time dimension in the MEGTWR model. Then, a test statistic based on the residual sum of squares is constructed for inferring the special types of coefficients and a bootstrap procedure is suggested to compute the p value of the test. A simulation study is consequently conducted to assess the performance of the test and the proposed average-based estimation of the MEGTWR model. A real-life example based on Beijing secondhand house prices from 2014 to 2018 is given to demonstrate the applicability of the proposed estimation and inference methods. The article concludes with discussions on the proposed methods and some related issues that deserve further study.

MEGTWR Model
Let Y be the response variable and X 1 , X 2 , Á Á Á , X r be r associated explanatory variables. Suppose that the data are collected at n spatial locations labeled as fðu i , v i Þg n i¼1 and T time points denoted by ft j g T j¼1 : Then, the observations of Y and X 1 , X 2 , :::, X r at the spatiotemporal sampling points ðu i , v i , t j Þ can be denoted by ðy ij ; x ij1 , x ij2 , :::, x ijr Þ where i ¼ 1, 2, :::, n and j ¼ 1, 2, :::, T: With these notations, the conventional GTWR model can be expressed as i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T, where b k ðu i , v i , t j Þ ðk ¼ 1, 2, :::, rÞ are r unknown regression coefficients at ðu i , v i , t j Þ and e ij ði ¼ 1, 2, :::, n; j ¼ 1, 2, :::, TÞ are model errors that are assumed to be independent and identically distributed random variables with Eðe ij Þ ¼ 0 and Varðe ij Þ ¼ r 2 : Generally, we take X 1 ¼ 1 (i.e., x ij1 1) to allow the model to include an intercept term. As mentioned in the introduction, the special types of coefficients including constant, only temporally varying, and only spatially varying coefficients are possible in the conventional GTWR model. Taking these possibilities into account yields the following MEGTWR model: þ e ij , i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T, where we assume, without loss of generality, that the first r 1 coefficients vary over both space and time, the next r 2 Àr 1 coefficients are constant, the following r 3 Àr 2 coefficients vary with time only, and the last rÀr 3 coefficients vary over space only. Here, we use different subscripts i and j for the spatial coordinates (u i , v i ) and the time point t j , respectively, which is slightly different from the notation in the conventional GTWR model, to make the forthcoming average-based estimators of the special types of coefficients in the model in Equation 2 easier to express and more explicit.

Modified Estimation of the Conventional GTWR Model with the Resulting Residual Sum of Squares
We first propose a modified estimation approach of the conventional GTWR model in Equation 1 by combining the simplicity in selecting the optimal spatial and temporal bandwidths in Huang, Wu, and Barry (2010) and the interpretability of the coefficient estimators in Fotheringham, Crespo, and Yao (2015).
are the design matrix consisting of the observations of the r explanatory variables and the observation vector of the response variable, respectively, and ¼Diagðw 11 ðu 0 ,v 0 ,t 0 Þ,:::,w n1 ðu 0 ,v 0 ,t 0 Þ,w 12 ðu 0 ,v 0 ,t 0 Þ,:::, w n2 ðu 0 ,v 0 ,t 0 Þ,:::,w 1T ðu 0 ,v 0 ,t 0 Þ,:::,w nT ðu 0 ,v 0 ,t 0 ÞÞ is the nT Â nT diagonal weight matrix with its elements in the diagonal being the modified weights in Equation 4. Here, we attach the subscripts b S and b T to W b S ,b T ðu 0 ,v 0 ,t 0 Þ to indicate that the weight matrix is related to the spatial and temporal bandwidths.
In regard to the selection of the optimal bandwidth sizes, rather than optimizing the ratio of the scale factors in Equation 3 to indirectly obtain the optimal size of the spatiotemporal bandwidth h in Huang, Wu, and Barry (2010), we directly optimize the sizes of the spatial bandwidth b S and temporal bandwidth b T by the cross-validation (CV) or AIC c criterion. Specifically, the CV score at b S and b T is computed by where tr(Á) stands for the trace of a matrix. The optimal sizes of the bandwidths b S and b T are jointly selected by minimizing CV The direct selection of b S and b T has the advantage that one can refer the scope of the whole studied spatiotemporal region to designate appropriate ranges of b S and b T for searching their optimal sizes.

Average-Based Estimation of the MEGTWR Model
For the MEGTWR model in Equation 2, we denote the vectors of the four types of coefficients, respectively, by Using the estimators in Equation 6-that is, the modified estimators of the coefficients in the conventional GTWR model at all of the spatiotemporal sampling points-and extending the idea of averaging the local estimators to obtain the estimators of the constant coefficients for semivarying coefficient models (W. Y. Zhang, Lee, and Song 2002), we propose an average-based estimation approach for the MEGTWR model in what follows.
1. The estimators of the spatiotemporally varying coefficients in b ST ðu i , v i , t j Þ are retained to be the corresponding ones in Equation 6. 2. Each constant coefficient in b C is estimated by averaging the local estimators of the corresponding coefficient over all of the spatiotemporal sampling points: 3. Each temporally varying coefficient in b T ðt j Þ is estimated by averaging the local estimators of the corresponding coefficient over the spatial sampling locations. That is, for k ¼ r 2 þ 1, r 2 þ 2, :::, :::, n: (15) With the foregoing coefficient estimators, the fitted values of the response variable Y at all of the spatiotemporal sampling points are computed bŷ i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T, and the residuals arê e ijðMÞ ¼ y ij Àŷ ijðMÞ , i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T, where the subscript M stands for the MEGTWR model. Let y M ¼ ðŷ 11ðMÞ , :::,ŷ n1ðMÞ ,ŷ 12ðMÞ , :::,ŷ n2ðMÞ , :::, y 1TðMÞ , :::,ŷ nTðMÞ Þ T be the fitted vector of the response variable Y. Then, the residual vector iŝ e M ¼ ðê 11ðMÞ , :::,ê n1ðMÞ ,ê 12ðMÞ , :::,ê n2ðMÞ , :::, e 1TðMÞ , :::,ê nTðMÞ Þ T ¼ yÀŷ M , and the residual sum of squares is A notable advantage of this estimation method for the MEGTWR is its simplicity. That is, the special types of coefficients are estimated directly based on the modified coefficient estimators in Equation 6 of the conventional GTWR model without the need to reselect the spatial and temporal bandwidth sizes, which is the main computation cost in a local smoothing estimation procedure.
Furthermore, after some algebraic operations, the hat matrix, which we denote by H M , of the averagebased estimation for the MEGTWR model can be derived, with which the fitted vectorŷ M in Equation 17 can be computed bŷ and the residual sum of squares RSS M can be expressed as Because of the verbose mathematical operations, the derivation of H M is shown in the Appendix.

Inference of the Special Types of Coefficients in the GTWR Model Related Hypotheses
Based on the average-based estimators of the coefficients, we then formulate a formal statistical test to infer the aforementioned special types of coefficients in the GTWR model. In the most general case, the null hypothesis is that the true coefficients in the GTWR model in Equation 1 show in fact those in the MEGTWR model in Equation 2; that is, r 1 spatiotemporally varying coefficients, r 2 Àr 1 constant coefficients, r 3 Àr 2 temporally varying coefficients, and rÀr 3 spatially varying coefficients. Therefore, the null model is the MEGTWR model in Equation 2. For ease of presentation, we simply state the null hypothesis as follows: H 0 : The MEGTWR model in Equation 2 is true: The alternative hypothesis might be diverse. To be specific, we consider here the most common case that all of the coefficients in the conventional GTWR model in Equation 1 vary over both space and time. A note on the extension of the test for other possible alternative hypotheses will be given at the end of this section. Similarly, the currently considered alternative hypothesis is simply represented as H 1 : The conventional GTWR in Equation 1 is true:

Test Statistic with the p Value
In the regression literature, the residual sum of squares is a commonly used statistic to measure the goodness of fit between a model and a data set. In this study, the residual sums of squares from both the null model and the alternative model have been derived earlier. Therefore, a statistic based on the residual sum of squares is constructed to implement the test for the special types of coefficients in the conventional GTWR model. The general form of the statistic is where RSSðH 0 Þ and RSSðH 1 Þ are the residual sums of squares from the null model fit and the alternative model fit, respectively. For the preceding null and alternative hypotheses, we have RSSðH 1 Þ ¼ RSS in Equation 10 and RSSðH 0 Þ ¼ RSS M in Equation 18 or Equation 20. In particular, with the hat matrixes of the modified estimation of the GTWR model and the average-based estimation of the MEGTWR model, Z can be expressed as where ( This kind of statistic in the nonparametric regression literature was originally proposed by Fan, Zhang, and Zhang (2001) and was termed the generalized likelihood ratio statistic. One of the most distinctive properties of the statistic is the so-called Wilks phenomenon, meaning that the null distribution of the statistic is in general irrelevant to the nonparametric fitting methods used for calibrating the models. In addition to wide applications in the nonparametric regression literature, such a statistic has also been used in the GWR literature for spatial nonstationarity diagnosis (see, e.g., Brunsdon, Fotheringham, and Charlton 1999;Leung, Mei, and Zhang 2000;Mei, Xu, and Wang 2016). Because a larger value of Z tends to support the alternative hypothesis, the p value of the test is where z is the observed value of Z computed by Equation 21 and P H 0 means that the probability is computed with the null distribution of Z. Given a significance level a, the null hypothesis is rejected if p<a and is not rejected if otherwise.
Bootstrap Procedure for Deriving the p Value of the Test Because the derivation of the exact or asymptotic null distribution of the statistic Z seems very difficult in the current case, we suggest a residual-based bootstrap procedure to compute the p value of the test, which can be routinely implemented by the following steps: 1. Calibrate the alternative model in Equation 1 by the modified estimation method based on the original data fy ij ; x ij1 , x ij2 , :::, x ijr ; ðu i , v i , t j Þg ði ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T). The optimal sizes of b S and b T and the coefficient estimators at all of the spatiotemporal sampling points are obtained, on which the residual vectorê ¼ ðê 11 , :::,ê n1 ,ê 12 , :::,ê n2 , :::,ê 1T , :::,ê nT Þ T in Equation 9 and the residual sum of squares RSS in Equation 10 are computed. Furthermore, the residual vectorê is centralized asê c ¼ ðê 11c , :::,ê n1c ,ê 12c , :::, e n2c , :::,ê 1Tc , :::,ê nTc Þ T whereê ijc ¼ê ij À 1 nT P T j¼1 P n i¼1ê ij ði ¼ 1, 2, :::, n; j ¼ 1, 2, ::: where IðÁÞ is the indicator function. Similar to the explanation in Mei, Xu, and Wang (2016), both RðH 0 Þ and RðH 1 Þ are unchanged for each bootstrap repetition and the bootstrap value z Ã in Step 3 can be simply computed by substituting y Ã into Equation 21 and there is no need to recalibrate the null and alternative models with the bootstrap sample, which can considerably reduce the computation cost.
This bootstrap method can be extended to test the hypotheses that both null and alternative models are all the MEGTWR models, provided that the three special types of constant, temporally varying, and spatially varying coefficients in the alternative model are nested in their respective special types of coefficients in the null model. In such a case, both null and alternative models are calibrated via the average-based procedure based on the modified coefficient estimators in Equation 6 of the conventional GTWR model, and the bootstrap samples of the model errors are drawn from the centralized residuals derived from the alternative model fit.

Simulation Study
In this section, a simulation study is conducted to assess the performance of the proposed test for inferring the special types of the coefficients in the GTWR model and the effectiveness of the averagebased approach for estimating the coefficients.

Design of the Experiment
The Spatiotemporal Layout. Under the threedimensional Cartesian coordinate system with the time axis being the upright one, the spatiotemporal region was designed as the unit cube ½0, 1 Â ½0, 1 Â ½0, 1: The n sampling locations fðu i , v i Þg n i¼1 were devised in such a way that each spatial location (u i , v i ) was determined by independently drawing a pair of random numbers from the uniform distribution U(0, 1). Such spatial sampling points might be more practical than the latticed ones. The T time points ft j g T j¼1 were equally spaced on the interval ½0, 1 with each t j ¼ ðjÀ1Þ=ðTÀ1Þ: In total, there are nT spatiotemporal sampling points ðu i , v i , t j Þ ði ¼ 1, 2, :::, n; j ¼ 1, 2, :::, TÞ in or on the cube. In the simulation, we took n ¼ 300 and T ¼ 8, resulting in the sample size of nT ¼ 2,400. Figure 1 depicts the 300 spatial sampling locations that were taken to be fixed throughout the simulation.
The Model with the Hypotheses. The following model with three explanatory variables X 1 , X 2 , and X 3 was considered in the simulation to generate the experimental data: i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T: For the coefficients with the hypotheses, we based the five functions A nonzero value of c corresponds the following common alternative hypothesis:  The larger the absolute value of c is, the more the alternative hypothesis deviates from each of the preceding four null hypotheses.
The Observations of the Explanatory Variables and Distributions of the Model Errors. The explanatory variable X 1 was set to be 1 (i.e., x ij1 ¼ 1 for i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T) to make the model in Equation 24 include the intercept b 1 (u,v,t). The observations ðx ij2 , x ij3 Þ ði ¼ 1, 2, :::, n; j ¼ 1, 2, :::, TÞ were independently drawn from the standard normal distribution N(0, 1) and the uniform distribution UðÀ1, 1Þ, respectively. To demonstrate that the bootstrap approximation to the null distribution of the test statistic is insensitive to the model error distribution, two kinds of distributions N(0, 1) and U À ffiffi ffi 3 p , ffiffi ffi 3 p À Á that have the common zero mean and unit variance were designated as the model error distribution, respectively. Given each set of the values ð1, x ij2 , x ij3 ; e ij Þ drawn from their respective distributions, the corresponding observation y ij is computed according to the model in Equation 24.
The Other Experimental Settings. The alternative model with each group of the coefficients was calibrated by the modified estimation method, in which the bi-square kernel function was used to generate the weights. The optimal sizes of b S and b T were selected by the CV criterion. To reduce the computation cost, the golden section search method (Kiefer 1953) was used to search for the optimal size of the spatial bandwidth b S at each prespecified equally spaced temporal bandwidth size of b T . Then, the pair (b S , b T ) with the lowest CV score was finally selected as the optimal spatial and temporal bandwidth sizes.
The values of c in the regression coefficients were appointed to be equally spaced on an appropriate symmetric interval depending on the increasing speed of the power for each group of the coefficients. To assess the Type I error and power of the test in identifying the special types of coefficients, we repeatedly run in each experimental setting the experiment for N ¼ 500 times and the rejection rate under a given significance level a was computed, in which both the observations of the explanatory variables and those of the model errors were resampled from their respective distributions. As is known in statistics, the rejection rate is an estimator of the Type I error under the null hypothesis and a measure of the test power under the alternative hypothesis. Therefore, if the test method is well formulated, the rejection rate under the null hypothesis should be close to the significance level a and the rate under the alternative hypothesis should increase rapidly with the increase of the absolute value of c. Moreover, in each experimental repetition, B ¼ 500 bootstrap samples were drawn to calculate the bootstrap p value in Equation 23.

Simulation Results with Analysis
Validity of the Type I Error of the Test. To comprehensively assess the validity of the Type I error of the test, we computed the rejection rates under the three commonly used significance levels of a ¼ 0.01, 0.05, and 0.10. The results for the four null hypotheses under the two kinds of the model error distribution are listed in Table 1. Table 1 indicates that the rejection rates in the various experimental settings are all close to their respective significance levels, indicating that the bootstrap procedure approximates the null distribution of the test statistic reasonably well and yields a valid Type I error. In particular, the rejection rates under the normal and nonnormal error distributions are all comparable for the four null hypotheses, indicating that the bootstrap test seems insensitive to the variation of the model error distribution.
Power of the Test. Figure 2 depicts the rejection rates under the significance level a ¼ 0:05 against different values of c for the four groups of Table 1. Rejection rates of the test under the four null hypotheses and the two kinds of model error distribution  Figure 2 shows that the test is quite powerful in identifying the special types of the coefficients in the GTWR model. Moreover, no evident difference between the powers under the normal and nonnormal model error distributions is observed, demonstrating that the power of the test is also insensitive to the variation of the model error distribution.
Computational Efficiency of the Bootstrap Test. The simulation was performed on our personal computer (Intel Core i7-8700K @ 3.70 GHz of CPU and 32 GB of memory) with R. For the sample size of nT ¼ 2,400, each experimental setting with N ¼ 500 replications for the bootstrap test with B ¼ 500 bootstrap samples drawn took about 7.5 hours in total. On average, each replication took about fifty-four seconds. It is noted that we only need to run one replication in practice to compute the p value of the test. Therefore, the computation cost of the bootstrap test seems acceptable at least for the moderate sample size.
Accuracy of the Coefficient Estimators by the Average-Based Estimation Method. We take the  Huang, Wu, and Barry (2010) and Fotheringham, Crespo, and Yao (2015), which we henceforth refer to as HGTWR and FGTWR, respectively. Specifically, the considered MEGTWR model is þ e ij , i ¼ 1, 2, :::, n; j ¼ 1, 2, :::, T, We repeatedly generated two sets of N ¼ 500 samples of the preceding model with the two kinds of the model error distribution. Here, we adopted the averaged root mean square error (ARMSE) over the nT spatiotemporal sampling points to measure the accuracy of each coefficient estimator. Letb ðqÞ ðu i , v i , t j Þ be the estimator of the coefficient bðu, v, tÞ at ðu i , v i , t j Þ with the qth sample by either of the three estimation methods, the ARMSE of bðu, v, tÞ is defined as ARMSEðbðu, v, tÞÞ For the considered model, bðu, v, tÞ is one of the coefficients b 1 ðu, v, tÞ, b 2 ðtÞ, and b 3 ðu, vÞ: It is noted that both HGTWR and FGTWR yield a spatiotemporally varying estimator of each coefficient no matter what the type of the coefficient really is. The average-based estimator, however, shares the same type with its real coefficient. Therefore, for the average-based method, the ARMSEs for the estimators of b 2 ðtÞ and b 3 ðu, vÞ can respectively be simplified as and ARMSEðb 3 ðu, vÞÞ : Table 2 lists ARMSEs of the coefficient estimators for the set of the N ¼ 500 samples with error distribution being U À ffiffi ffi 3 p , ffiffi ffi 3 p À Á , where the coefficient estimators were derived with the same bi-square kernel and CV criterion for bandwidth selection in the three estimation methods. The results for the samples with the error distribution N(0, 1) are very similar, so we omitted them to save space. Table 2 shows that the average-based method yields median accuracy for the spatiotemporally varying coefficient b 1 (u,v,t), which is reasonable because the average-based estimator of b 1 (u,v,t) is in fact the one yielded by the modified GTWR estimation method, which, as previously mentioned, combines the simplicity in bandwidth selection in HGTWR and the interpretability of the coefficient estimators in FGTWR. For the special-type coefficients b 2 ðtÞ and b 3 ðu, vÞ, however, the proposed average-based method consistently produced much more accurate estimators than both HGTWR and FGTWR with the relative ARMSE decreasing by 37 percent for b 2 ðtÞ and by 14 percent and 21 percent for b 3 ðu, vÞ, respectively. The reason is that both HGTWR and FGTWR do not take into account the coefficient type information and, as aforesaid, yield such estimators that vary over both space and time.
Furthermore, for each estimation method, we averaged the estimators of each coefficient over the N ¼ 500 samples as its final estimator to visually compare the accuracy of the coefficient estimators by means of graphics. As the final estimators of b 2 ðtÞ and b 3 ðu, vÞ by HGTWR and FGTWR are both spatially and temporally varying, their domains are different from those of the final estimators by the average-based method, making less sense for the visual comparison. We therefore turn to compare the estimated curve of b 2 ðtÞ and the estimated surface of b 3 ðu, vÞ by the average-based method with their true curve and true surface, respectively.  Figure 3 depicts the true surfaces and the estimated surfaces of b 1 (u,v,t) at the time points t ¼ t 1 , t 4 , and t 8 . It is observed that, at each time point, the estimated surfaces by the three estimation methods all show a very similar spatial pattern to that of the true surface, which provides further evidence that the three estimation methods are of comparable accuracy on the estimator of the spatiotemporally varying coefficient b 1 (u,v,t). Figure 4 shows the curve of the final estimator of b 2 ðtÞ by the average-based method and its true curve; Figure 5 depicts the surface of the final estimator of b 3 ðu, vÞ with its true surface. It is observed that the average-based method also yields quite Figure 3. True surfaces of the coefficient b 1 ðu, v, tÞ (the first row) and those of its final estimators by the average-based method (the second row), Huang, Wu, and Barry (2010) geographically and temporally weighted regression (HGTWR; the third row), and Fotheringham, Crespo, and Yao (2015) geographically and temporally weighted regression (FGTWR; the fourth row) at t ¼ t 1 , t 4 , and t 8 . accurate estimators on the special-type coefficients. Although the surface of the final estimator of b 3 ðu, vÞ shows some deterioration from its true surface, the typical variation pattern is well followed by the estimator.
Moreover, the simulation study for assessing the performance of the test and the accuracy of the coefficient estimators by average-based method was conducted for a larger sample size and a larger variance of the model error. It is found that the test still performs well in both cases. The accuracy of the estimation method shows obvious improvement under the larger sample size, whereas some deterioration is observed under the larger model error variance. The detailed results are provided as supplementary material due to the limited space.

Real-Life Example
In this section, the Beijing secondhand house price data from 2014 to 2018 are used as an example to demonstrate the applicability of the proposed estimation and inference methods. This data set was analyzed by Z. Zhang et al. (2021) using the multiscale GTWR with a unilateral temporal weighting scheme, which was referred to as MUGTWR in that work. After preprocessing the crude data from Z. Zhang et al. (2021), the data set consists of the yearly inflation-adjusted residential-community (RC)-based average house prices (HPRICE) in about 3,500 RCs in each year and the eight explanatory variables with the Cartesian spatial coordinates (u i , v i ) for each RC attached also. The eight explanatory variables or factors are the following:  To be in accord with the notations in this article, we further removed such RCs with their associated observations that were not been fully recorded from 2014 to 2018 and used the remaining 2,576 RCs with their observations to build the following GTWR model to detect the special types of the coefficients: where i ¼ 1, 2, :::, 2576 and j ¼ 1, 2, :::, 5, yielding the sample size of nT ¼ 2, 576 Â 5 ¼ 12, 880: As done in Z. Zhang et al. (2021), all of the explanatory variables were standardized to make them share a common value range. This model is in fact the same as that in Z. Zhang et al. (2021) except for the mathematical notations. Although the current sample size is about 75 percent of that in Z. Zhang et al. (2021), the spatial pattern of the sampling points is essentially retained. Therefore, the coefficient-specific spatial and temporal bandwidths selected in Z. Zhang et al. (2021) could provide some intuitive information on the types of the coefficients in the model in Equation 26 based on the fact that a relatively large spatial or temporal bandwidth size implies less spatial or temporal heterogeneity of the corresponding coefficient. In MUGTWR, a pair of spatial and temporal bandwidths (b S , b T ) was set at each time stamp for each spatiotemporally varying coefficient, where b T at the Tth year is the number of previous year(s) in which the observations are involved in the estimation and b S is the ratio of the number of the sampling points with their observations involved in the estimation of the coefficients at the Tth year to the number of the total sampling points in the Tth and the previous years. The selected optimal bandwidths are shown in Table 5 in Z. Zhang et al. (2021). To summarize, we averaged for each coefficient the optimal sizes of both b S and b T at the five time stamps, which we denoted here by b S and b T , respectively. We then picked out such values of b S and/or b T that are relatively large and listed in the last two columns of Table 3 the corresponding pairs of ð b S , b T Þ that at least one of the values in each pair is picked out.
Based on the picked pairs of ð b S , b T Þ, we formulated the six null hypotheses for the corresponding individual coefficients, shown in the first column of Table 3. It should be pointed out that, although the value of b T for b 5 (u,v,t) (the coefficient of FLOOR) is relatively small, we still formulated the null hypothesis of b 5 ðu, v, tÞ ¼ b 5 : The reason is that the optimal temporal bandwidth sizes in the first four of the total five years are all largest at the corresponding years, showing b T ¼ 0, 1, 2, and 3. We first tested each of the null hypotheses versus the alternative hypothesis that all of the coefficients in the model in Equation 26 are spatiotemporally varying using the proposed estimation and test methods. The CV criterion with the bi-square kernel was once again used to select the optimal spatial and temporal bandwidth sizes, yielding b S ¼ 3.9 km and b T ¼ 1.5 years. Moreover, B ¼ 500 bootstrap repetitions were run to derive the p value of each test. The p  Table 3, where the corresponding explanatory variables are also attached. We know from the p values that we cannot reject the last five null hypotheses, which is consistent with the intuitive judgment based on the optimal spatial and temporal bandwidth sizes selected in Z. Zhang et al. (2021). We should, however, reject the first null hypothesis because of the very small p value, which conflicts with the intuitive judgment. These results demonstrate that the optimal coefficient-specific bandwidth sizes selected in a multiscale estimation method for the conventional GTWR model do provide some prior information for judging the types of the coefficients, but the judgment based on the coefficient-specific bandwidth sizes is rather intuitive because, unlike the statistical hypothesis testing, it is hard to appoint appropriate threshold values of the bandwidths for achieving the judgment.
Regarding the computation cost, to perform each setting of the preceding tests took about twenty-six to twenty-nine minutes on our personal computer, which seems acceptable in view of the large sample size of nT ¼ 12,880.
We then jointed the previously accepted null hypotheses and tested the comprehensive null hypothesis that and b 9 ðu, v, tÞ ¼ b 9 versus the alternative hypothesis that all of the coefficients in the model in Equation 26 vary over both space and time. The resulting p value is 1.000, suggesting that the null hypothesis cannot be rejected, on which we could build the following MEGTWR model for the analysis of the Beijing house price data: i ¼ 1, 2, :::, 2576; j ¼ 1, 2, :::, 5: It should be noted that, although the preceding model reflects more specific spatiotemporally varying patterns of the explanatory variables influencing the house price than the conventional GTWR model in Equation 26, it cannot be guaranteed that the selected model is the best MEGTWR model for the house price data, because the proposed test method has been used only in a heuristic way in inferring the coefficient types with the prior information of the optimal spatial and temporal bandwidth sizes in Z. Zhang et al. (2021). There is still a gap for the test to be used as a complete model selection tool for an MEGTWR model because too many competing null hypotheses need to be considered. This issue is further discussed in the final section.
We then calibrated the selected MEGTWR model in Equation 27 with the proposed average-based estimation method. For comparison, the model was also calibrated by HGTWR, FGTWR, and the multiscale method MUGTWR under the bi-square kernel and CV criterion. First, we compare the goodness of fits among the estimation methods in terms of trðHÞ, AIC c , residual sum of squares (RSS), and R 2 , where trðHÞ, the trace of the hat matrix H, called the effective degrees of freedom in the nonparametric regression, measures the model complexity (Hurvich, Simonoff, and Tsai 1998). The related results are reported in Table 4, where the values of trðHÞ and AIC c for MUGTWR are not available because its hat matrix remains to be investigated.
It is evident from Table 4 that the average-based method yields the smaller value of trðHÞ, indicating that the MEGTWR model that the average-based method calibrates is the simplest at the optimal sizes of the bandwidths. The average-based method, however, produces the largest value of RSS, leading to a larger value of AIC c and a slightly smaller value of R 2 . These results are reasonable because what HGTWR, FGTWR, and MUGTWR calibrate is in fact the conventional GTWR model in Equation 26, We further show in Table 5 the typical summary statistics of the coefficient estimators by the different estimation methods. It can be observed that, although the values of the statistics deviate among the different methods, their signs are almost coincident, which could result in a basically consistent interpretation on the spatiotemporal patterns of the factors influencing the house price.
More specific spatiotemporal patterns can be exposed, though, by the average-based estimators of the special-type coefficients in the MEGTWR model in Equation 27. To make a visual show, the estimators of the spatially varying coefficients b 6 (u,v) (HERATIO) and b 7 (u,v) (PLOTRATIO) are depicted in Figure 6. First, based on the constant coefficient estimators ofb 5 ¼ 0:006,b 8 ¼ À0:007, andb 9 ¼ À0:023 which correspond to the factors FLOOR, DBUS, and DSUBWAY, respectively, we can conclude that the average floor level of the transacted houses (FLOOR) and the distances from each RC to its nearest bus station (DBUS) and  Fotheringham, Crespo, and Yao (2015) geographically and temporally weighted regression; MUGTWR ¼ multiscale GTWR with a unilateral temporal weighting scheme.
subway station (DSUBWAY) have a spatiotemporally global effect on the house price. FLOOR has a positive effect on the house price, while DBUS and DSUBWAY have a negative influence on the house price. From Figure 6, we observe that the effect of the household's ratio of staircases (HERATIO) and the ratio of the building floor area to the land area (PLOTRATIO) on the house price varies from place to place, with the spatial variation pattern kept over the five years. Caution should be exercised in these interpretations, however, in view of the facts that the values of some goodness-of-fit statistics yielded by the average-based estimation method are somewhat worse and that the proposed test is a global one and local significant effects for space, time, or both might also be possible even if the test has shown that the effects are globally insignificant.

Conclusion
In this article, we first modified the existing estimation methods of the conventional GTWR model by combining the simplicity in bandwidth selection in Huang, Wu, and Barry (2010) and the interpretability of the coefficient estimators in Fotheringham, Crespo, and Yao (2015). Then, based on the coefficient estimators obtained by the modified procedure,

Supplemental Material
Supplemental simulation study for this article can be accessed on the publisher's site at: http://dx.doi. org/10.1080/24694452.2022.2092443. The supplemental material is the performance assessment of the proposed test and the average-based estimation method under a larger sample size or a larger model error variance.
The R codes for performing the proposed test and estimation methods in the simulation study are available in figshare.com via the private link http:// doi.org/10.6084/m9.figshare.16850962.
The raw Beijing secondhand house price data used in this study belong to Homelink, a private real estate agency company. The data are not publicly available unless authorized by Homelink.
The above hat matrix seems quite complex in the very general case in the main text. When we are only interested in identifying a specific type of the coefficients or a single coefficient of a special type, however, the hat matrix can be largely simplified.