A Note on Parameter Estimation for Misspecified Regression Models with Heteroskedastic Errors

Finding the best fitting sinusoid is closely related to fitting a linear model. Using the sine angle addition formula we can rewrite the maximum likelihood estimator from Equation (2.2) as

The sum over $i$ can be simplified by noting the linearity of the model and reparameterizing. Let $Y=(y_{1},\ldots,y_{n})^{T}$. Let $\beta_{k1}=a_{k}\cos(\phi_{k})$ and $\beta_{k2}=a_{k}\sin(\phi_{k})$. Define $\beta=(\beta_{0},\beta_{11},\beta_{12},\ldots,\beta_{K1},\beta_{K2})^{T}\in\mathbb{R}^{2K+1}$. Let $\Sigma$ be a $n\times n$ diagonal matrix where $\Sigma_{ii}=\sigma_{i}^{2}$. Define

We rewrite the ML estimator as

Every frequency $\omega$ in the grid of frequencies $\Omega$ determines a design matrix $X(\omega)$. At a particular $\omega$, the $\beta$ which minimizes the objective function is the weighted least squares estimator

The frequency estimate may then be written as

Thus estimating frequency involves performing a weighted least squares regression (Equation (2.3)) at every frequency in the grid $\Omega$. The motivation for the procedure is maximum likelihood. As discussed earlier, in cases where the model is misspecified, there is no theoretical support for using $\Sigma^{-1}$ as the weight matrix in either Equation (2.3) or (2.4).

Let $X\in\mathbb{R}^{n\times p}$ be the matrix with row $i$ equal to $x_{i}^{T}$. Let $Y=(y_{1},\ldots,y_{n})^{T}$. Let $\Sigma$ be the diagonal matrix of observation variances such that $\Sigma_{ii}=\sigma_{i}^{2}$. Let $\widehat{W}$ be a diagonal positive definite matrix. The weighted least squares estimator is

In this work we seek $\widehat{W}$ which minimize error in estimating $\beta=\mathbb{E}[xx^{T}]^{-1}\mathbb{E}[xf(x)]$.

There is a long history of studying estimators for misspecified models, often in the context of sandwich estimators for asymptotic variances. In [10], it was shown that when the true data generating distribution $\theta_{t}$ is not in the model, the MLE converges to the distribution $\theta_{0}$ in the model $\Theta$ which minimizes Kullback–Liebler divergence, ie

The asymptotic variance has a “sandwich” form which is not the inverse of the information matrix. [30] and [31] studied this behavior in the context of the linear regression model and the OLS estimator, proposing consistent estimators of the asymptotic variance. See [18] and [15] for sandwich estimators with improved finite sample performance and [27] for recent work on sandwich estimators in a Bayesian context. [2] provides a summary of sandwich estimators and proposes a bootstrap estimator for the asymptotic variance. By specializing our asymptotic theory from the weighted to the unweighted case, we rederive some of these results. However our focus is different in that we find weightings for least squares which minimize asymptotic variance, rather than estimating the asymptotic variance of unweighted procedures.

Other work has focused on correcting model misspecification, often by modeling deviations from a parametric regression function with some non–parametric model. [1] studied model misspecification when response variances are known up to a constant due to repeated measurements, ie $Var(y_{i})=\sigma^{2}/m_{i}$ where $m_{i}$ is known. A Gaussian prior was placed on $\beta$ and the non–linear component $g$ was modeled as being drawn from a Gaussian process. See [13] for an example with homoskedastic errors in the context of computer simulations. See [6] for an example in astronomy with known heteroskedastic errors. Our focus here is different in that instead of correcting model misspecification we consider how weighting observations affects estimation of the linear component of $f$.

Heteroskedasticity in the partial linear model

is studied in [17] and [16]. Here $\text{Var}(\epsilon_{i})=\xi(x_{i},z_{i})$ for some function $\xi$. The parameter $h$ is some unknown function. The response $y$ depends on the $x$ covariates linearly and the $z$ covariates nonlinearly. When $h$ is estimated poorly, weighting by the inverse of the observation variances causes parameter estimates of $\beta$ to be inconsistent. In contrast, ignoring observation variances leads to consistent estimates of $\beta$. Qualitatively these conclusion are similar to our own in that they caution against using weights in the standard way.

Our asymptotic theory makes assumptions on the form of the weight matrix.

These assumptions include both the ordinary least squares (OLS) estimator where $\widehat{W}_{ii}=w(\sigma_{i})=1$ and the standard weighted least squares estimator where $\widehat{W}_{ii}=w(\sigma_{i})=\sigma_{i}^{-2}$ (assuming $\mathbb{E}[\sigma^{-8}]<\infty$). In both these cases $\delta_{nm_{i}}=0$ and $d=0$ for all $n,m$. These additional terms are used in Sections 3.5 and 3.6 to construct adaptive estimators for the known and unknown variance cases.

The major assumption here is independence between $x$ and $\sigma$. We address dependence in Section 3.7.

See Section A.1 for a proof. If the response is linear ($g\equiv 0$) then the variance is

Setting $w(\sigma)=\sigma^{-2}$ we have $\frac{\mathbb{E}[\sigma^{2}w^{2}]}{\mathbb{E}[w]^{2}}=(\mathbb{E}[\sigma^{-2}])^{-1}$. This is the standard weighted least squares estimator. This $w$ can be shown to minimize the variance using the Cauchy Schwartz inequality. With $w(\sigma)=1$, the asymptotic variance can be rewritten

This is the sandwich form of the covariance for OLS derived in [30] and [31] (see [2], specifically Equations 1-3), valid even when $\sigma$ and $x$ are not independent.

For notational simplicity define

The asymptotic variances for OLS ($\widehat{\beta}(I)$) and standard WLS ($\widehat{\beta}(\Sigma^{-1})$) are

Each of these asymptotic variances is composed of the same two terms. The $A$ term is caused by model misspecification while the $B$ term is the standard asymptotic variance in the case of no model misspecification. The coefficient on $A$ is larger for $W=\Sigma^{-1}$ because $\frac{\mathbb{E}[\sigma^{-4}]}{\mathbb{E}[\sigma^{-2}]^{2}}\geq 1$ by Jensen’s Inequality. The coefficient on $B$ is larger for $W=I$ because $\mathbb{E}[\sigma^{2}]\geq\frac{1}{\mathbb{E}[\sigma^{-2}]}$. The relative merits of OLS and standard WLS depend on the size of the coefficients and the precise values of $A$ and $B$. However, qualitatively, OLS and standard WLS suffer from high asymptotic variance in opposite situations which depend on the distribution of the errors. To make matters concrete, consider error distributions of the form

where $\delta_{1},\delta_{2}$ are small non–negative numbers and $c>1$ is large. Note that $A$ and $B$ do not depend on $F_{\sigma}$.

• •

  $\delta_{1}=0,\delta_{2}>0$: In this situation the error standard deviation is usually $1$ and occasionally some large value $c$. The result is large asymptotic variance for OLS. Since $\mathbb{E}[\sigma^{2}]>c^{2}\delta_{2}$,

  $$\nu(I)\succeq A+c^{2}\delta_{2}B$$

  For large $c$ this will be large. In contrast the coefficients on $A$ and $B$ for standard WLS can be bounded. For the coefficient on $B$ we have $\mathbb{E}[\sigma^{-2}]^{-1}\leq(1-\delta_{2})^{-1}$. The coefficient on $A$ with $c>1$ is

  $$\frac{\mathbb{E}[\sigma^{-4}]}{\mathbb{E}[\sigma^{-2}]^{2}}=\frac{\delta_{2}c^{-4}+(1-\delta_{2})}{\delta_{2}^{2}c^{-4}+2\delta_{2}c^{-2}(1-\delta_{2})+(1-\delta_{2})^{2}}<\frac{1}{1-\delta_{2}}.$$

  Therefore

  $$\nu(\Sigma^{-1})\preceq(1-\delta_{2})^{-1}(A+B).$$

  In summary, standard WLS performs better than OLS when there are a small number of observations with large variance.

• •

  $\delta_{1}>0,\delta_{2}=0$: In this situation the error standard deviation is usually $1$ and occasionally some small value $c^{-1}$. For standard WLS with $c$ large and $\delta_{1}$ small, the coefficient for $A$ is

  $$\frac{\mathbb{E}[\sigma^{-4}]}{\mathbb{E}[\sigma^{-2}]^{2}}=\frac{\delta_{1}c^{4}+(1-\delta_{1})}{\delta_{1}^{2}c^{4}+2\delta_{1}c^{2}(1-\delta_{1})+(1-\delta_{1})^{2}}\approx\frac{1}{\delta_{1}}.$$

  Thus the asymptotic variance induced by model misspecification will be large for standard WLS. In contrast, we can bound the asymptotic variance above for OLS, independently of $c$ and $\delta_{1}$. Since $c>1$, $\mathbb{E}[\sigma^{2}]<1$ and

  $$\nu(I)\preceq A+B.$$

The case where both $\delta_{1}$ and $\delta_{2}$ are non–zero presents problems for both OLS and standard WLS. For example if $\delta=\delta_{1}=\delta_{2}$, both OLS and standard WLS can be made to have large asymptotic variance by setting $\delta$ small and $c$ large. In the following section we construct an adaptive weighting which improves upon both OLS and standard WLS.

Let $\Gamma$ be a linear function from the set of $p\times p$ matrices to $\mathbb{R}$ such that $\Gamma(C)>0$ whenever $C$ is positive definite. We seek some weighting $w=w(\sigma)$ for which $\Gamma(\nu(w))$ (recall that $\nu$ is the asymptotic variance) is lower than OLS and standard WLS. Natural choices for $\Gamma$ include the trace (minimize the sum of variances of the parameter estimates) and the $\Gamma(C)=C_{jj}$ (minimize the variance of one of the parameter estimates).

Section A.2 contains a proof. The proportionality is due to the fact that the estimator is invariant to multiplicative scaling of the weights.

A proof is contained in Section A.3. Thus if we can construct a weight matrix $\widehat{W}$ which satisfies Assumptions 1 with $w(\sigma)=w_{min}(\sigma)$ , then by the preceding theorem the associated weighted estimator will have lower asymptotic variance then either OLS or standard WLS. We now construct such a weighting in the case of known and unknown error variances.

With the $\sigma_{i}$ known we only need to estimate $A$ and $B$ in $w_{min}$ in Equation (3.3). Let $\Delta=\Gamma(A)\Gamma(B)^{-1}$. Let

Let $\widehat{\beta}(\widehat{W})$ be a root $n$ consistent estimator of $\beta$ (eg $\widehat{W}=I$ is root $n$ consistent by Theorem 3.1) and let

Let

Then we have

The estimated optimal weighting matrix is the diagonal matrix $\widehat{W}_{min}$ with diagonal elements

A few notes on this estimator:

• •

  The term $x_{i}x_{i}^{T}\widehat{g}(x_{i})^{2}$ is an estimate of $x_{i}x_{i}^{T}g(x_{i})^{2}$. These estimates are weighted by $\sigma_{i}^{-4}$. The term $(\sum\sigma_{i}^{-4})^{-1}$ normalizes the weights. This weighting is motivated by the fact that

  	$\displaystyle x_{i}x_{i}^{T}\widehat{g}(x_{i})^{2}$	$\displaystyle=x_{i}x_{i}^{T}((y_{i}-x_{i}^{T}\widehat{\beta}(\widehat{W}))^{2}-\sigma_{i}^{2})$
  		$\displaystyle=x_{i}x_{i}^{T}((y_{i}-x_{i}^{T}\beta)^{2}-\sigma_{i}^{2})+O(n^{-1/2})$
  		$\displaystyle=x_{i}x_{i}^{T}((g(x_{i})+\sigma_{i}\epsilon_{i})^{2}-\sigma_{i}^{2})+O(n^{-1/2}).$

  Analysis of the first order term shows

  $$\mathbb{E}[x_{i}x_{i}^{T}((g(x_{i})+\sigma_{i}\epsilon_{i})^{2}-\sigma_{i}^{2})|x_{i},\sigma_{i}]=x_{i}x_{i}^{T}g^{2}(x_{i})$$

  and

  	$\displaystyle\text{Var}(x_{i}x_{i}^{T}((g(x_{i})+\sigma_{i}\epsilon_{i})^{2}-\sigma_{i}^{2})|x_{i},\sigma_{i})_{jk}$
  	$\displaystyle=x_{ij}^{2}x_{ik}^{2}(\sigma_{i}^{4}(\mathbb{E}[\epsilon^{4}]-1)+4g(x_{i})^{2}\sigma_{i}^{2}+4g(x_{i})\sigma_{i}^{3}\mathbb{E}[\epsilon^{3}])$

  Thus by weighting the estimates by $\sigma_{i}^{-4}$, we can somewhat account for the different variances. Unfortunately since the variance depends on $g$, $\mathbb{E}[\epsilon^{3}]$, and $\mathbb{E}[\epsilon^{4}]$ which are unknown, it is not possible to weight by exactly the inverse of the variances. Other weightings are possible and in general adaptivity will hold.

• •

  Since $A$ and $B$ are positive semi–definite, $\Gamma(A)\Gamma(B)^{-1}\geq 0$. Thus for estimating $\Delta$, we use the maximum of a plug–in estimator and $0$ (Equation (3.4)).

See Section A.4 for a proof. Theorem 3.3 shows it is possible to construct better estimators than both OLS and standard WLS. In practice, it may be best to iteratively update estimates of $\widehat{W}_{min}$ starting with a known root $n$ consistent estimator such as $\widehat{W}=I$. We take this approach in our numerical simulations in Section 4.1.

For the purposes of making confidence regions we need estimators of the asymptotic variance. Above we developed consistent estimators for $A$ and $B$. We take a plug–in approach to estimating the asymptotic variance for a particular weighting $W$. Specifically

We also define the oracle $\widehat{\nu}_{OR}(\widehat{W})$ which is the same as $\widehat{\nu}_{1}$ but uses $A$ and $B$ rather than $\widehat{A}$ and $\widehat{B}$. While $\widehat{\nu}_{OR}$ cannot be used in practice, it is useful for evaluating the performance of $\widehat{\nu}_{1}$ in simulations.

Finally suppose the error variance is known up to a constant, i.e. $\sigma_{i}^{2}=k\tau_{i}^{2}$ where $\tau_{i}^{2}$ is known but $k$ and $\sigma_{i}^{2}$ are unknown. In the case without model misspecification, one can simply use weights $\tau_{i}^{-2}$ since the weighted estimator is invariant up to rescaling of the weights. The situation is more complicated when model misspecification is present. Simulations and informal mathematical derivations (not included in this work) suggest that replacing the $\sigma_{i}$ with $\tau_{i}$ in Equation (3.5) results in weights that are suboptimal. In particular, when $k>1$ (underestimated errors), the resulting weights are closer to OLS than optimal while if $k<1$ (overestimated errors), the resulting weights are closer to standard WLS than optimal.

Suppose for observation $i$ we observe $m_{i}\in\{1,\ldots,M\}$, the group membership of observation $i$. Observations in group $m$ have the same (unknown) variance $\sigma_{m}^{2}>0$. See [8], [5], and [9] for work on grouped error models in the case where the response is linear.

The $m_{i}$ are assumed independent of $x_{i}$ and $\epsilon_{i}$, with probability mass function $f_{m}$ (supported on $1,\ldots,M$). While the $\sigma_{m}$ for $m=1,\ldots,M$ are fixed unknown parameters, the probability mass function $f_{m}$ induces the probability distribution function $F_{\sigma}$ on $\sigma$. So we can define

for any function $h$.

Theorem 3.1 shows that even if the $\sigma_{m}$ were known, standard weighted least squares is not generally optimal for estimating $\beta$ in this model. It is not possible to estimate $w_{min}$ as proposed in Section 3.5 because that method requires knowledge of $\sigma_{m}$. However we can re–express the optimal weight function as

where the last equality defines $C_{m}$. Note that $\sigma_{m}$ is a fixed unknown parameter, not a random variable. One can estimate $B$ with $\widehat{B}=(n^{-1}X^{T}X)^{-1}$ and $C_{m}$ with

where $\widehat{\beta}(\widehat{W})$ is a root $n$ consistent estimator of $\beta$ (for example $\widehat{W}=I$ suffices by Theorem 3.1). The estimated weight matrix $\widehat{W}_{min}$ is diagonal with

See Section A.5 for a proof. Thus in the case of unknown errors it is possible to construct an estimator which outperforms standard WLS and OLS. As is the case with known errors, one can iteratively update $\widehat{W}_{min}$, starting with some (possibly inefficient) root $n$ consistent estimate of $\beta$.

For estimating the asymptotic variance we cannot use Equation (3.6) because that method required an estimate of $A$, a quantity for which we do not have an estimate in the unknown error variance setting. Instead note that the asymptotic variance of Equation (3.1) may be rewritten

Thus a natural estimator for the asymptotic variance is

Suppose one drops the independence assumption between $x$ and $\sigma$. This will be the case whenever the error variance is a function of $x$, a common assumption in the heteroskedasticity literature [3, 4, 12]. We require the weight matrix $W$ to be diagonal positive definite with diagonal elements $W_{ii}=w(\sigma_{i})$, some function of the error variance. The estimator for $\beta$ is

Recalling we write $w$ for $w(\sigma)$, we have the following result.

See Section A.6 for a proof. If $x$ and $\sigma$ are independent then the r.h.s is $\mathbb{E}[xx^{T}]^{-1}\mathbb{E}[xf(x)]$ and the estimator is consistent (as demonstrated by Theorem 3.1). Interestingly the estimator is also consistent if one lets $w(\sigma)=1$ (OLS), regardless of the dependence structure between $x$ and $\sigma$. However weighted estimators will not generally be consistent (including standard WLS). This observation suggests the OLS estimator may be preferred in the case of dependent errors. We show an example of this situation in the simulations of Section 4.1.

We conduct a small simulation study to demonstrate some of the ideas presented in the last section.¹¹1Code to reproduce the work in this section can be accessed at http://stat.tamu.edu/~jlong/hetero.zip or by contacting the author. Consider modeling the function $f(x)=x^{2}$ using linear regression with an intercept term. Let $x\sim Unif(0,1)$. The best linear approximation to $f$ is $\beta_{1}+\beta_{2}x$ where $\beta_{1}=-1/6$ and $\beta_{2}=1$. We first suppose $\sigma$ is drawn independently from $x$ from a discrete probability distribution such that $P(\sigma=0.01)=P(\sigma=1)=0.05$ and $P(\sigma=0.1)=0.9$. Since $\sigma$ has support on a finite set of values, we can consider the cases where $\sigma_{i}$ is known (Section 3.5) and where only the group $m_{i}$ of observation $i$ is known (Section 3.6). We let $\Gamma$ be the trace of the matrix.

We generate samples of size $n=100$, $N=1000$ times and make scatterplots of the parameter estimates using weights $W=\Sigma^{-1}$ (standard WLS), $W=I$ (OLS), $\widehat{W}_{min}=(\Sigma+\widehat{\Delta})^{-1}$, and $\widehat{W}_{min,ii}=\frac{\Gamma(\widehat{B})}{\Gamma(\widehat{B}^{T}\widehat{C}_{m_{i}}\widehat{B})}$. The OLS estimator does not require any knowledge about the $\sigma_{i}$. The fourth estimator uses only the group $m_{i}$ of observation $i$. For the two adaptive estimators, we use $\widehat{\beta}(I)$ as an initial root $n$ consistent estimator of $\beta$ and iterate twice to obtain the weights.

Results are shown in Figure 2. The red ellipses are the asymptotic variances. The results show that OLS outperforms standard WLS. Estimating the optimal weighting with or without knowledge of the variances outperforms both OLS and standard WLS. Exact knowledge of the weights (c) somewhat outperforms only knowing the group membership of the variances (d).

We construct 95% confidence regions using estimates of the asymptotic variance and determine the fraction of times (out of the $N$ simulations) that the true parameters are in the confidence regions. Recall that in Section 3.5 we proposed $\widehat{\nu}_{1}$ (Equation (3.6)) as well as the oracle $\widehat{\nu}_{OR}$ for estimating the asymptotic variance when the error variances are known. In Section 3.6 we proposed $\widehat{\nu}_{2}$ (Equation (3.8)) when the error variances are unknown. The estimator $\widehat{\nu}_{2}$ can also be used when the error variances are known. We use all three of these methods for constructing confidence regions for standard WLS, OLS, and $\widehat{W}=(\Sigma+\widehat{\Delta})^{-1}$. For $\widehat{W}_{ii}=\frac{\Gamma(\widehat{B})}{\Gamma(\widehat{B}^{T}\widehat{C}_{m_{i}}\widehat{B})}$ we use only $\widehat{\nu}_{2}$ because $\widehat{\nu}_{1}$ requires knowledge of $\Sigma$. Table 1 contains the results. While for OLS the nominal coverage probability is approximately attained, the other methods are anti–conservative for $\widehat{\nu}_{1}$ and $\widehat{\nu}_{2}$. Estimates for WLS are especially poor. The performance of the oracle is rather good, suggesting that the problem lies in estimating $A$ and $B$.

To illustrate the importance of the $\sigma$, $x$ independence assumption, we now consider the case where $\sigma$ is a function of $x$. Specifically,

All other parameters in the simulation are the same as before. Note that the marginal distribution of $\sigma$ is the same as the first simulation. We know from Section 3.7 that weighted estimators may no longer be consistent. In Figure 3 we show a scatter plot of parameter estimates using standard WLS and OLS. We see that the WLS estimator has low variance but is highly biased. The OLS estimator is strongly preferred.

[25] identified 483 RR Lyrae periodic variable stars in Stripe 82 of the Sloan Digital Sky Survey III. We obtained 450 of these light curves from a publicly available data base [11].²²2We use only the g–band data for determining periods. Figure 1a shows one of these light curves. These light curves are well observed ($n>50$), so it is fairly easy to estimate periods. For example, [25] used a method based on the Supersmoother algorithm of [7]. However there is interest in astronomy in developing period estimation algorithms that work well on poorly sampled light curves [29, 19, 14, 24]. Well sampled light curves offer an opportunity to test period estimation algorithms because ground truth is known and they can be artificially downsampled to create realistic simulations of poorly sampled light curves.

As discussed in Section 2, each light curve can be represented as $\{(t_{i},y_{i},\sigma_{i})\}_{i=1}^{n}$ where $t_{i}$ is the time of the $y_{i}$ brightness measurement made with uncertainty $\sigma_{i}$. In Figure 4 we plot magnitude error ($\sigma_{i}$) against magnitude ($y_{i}$) for all observations of all 450 light curves. For higher magnitudes (less bright observations), the observation uncertainty is larger. In an attempt to ensure independence between $\sigma$ and $x$ assumed by our asymptotic theory, we use only the bright stars in which all magnitudes are below 18 (left of the vertical black line in Figure 4). In this region, magnitude and magnitude error are approximately independent. This reduces the sample to 238 stars. We also ran our methods on the larger set of stars. Qualitatively, the results which follow are similar.

In order to simulate challenging period recovery settings, we downsample each of these light curves to have $n=10,20,30,40$. We estimate periods using sinusoidal models with $K=1,2,3$ harmonics. For each model we consider three methods for incorporating the error variances. In the first two methods, we weight by the the inverse of the observations variances ($\Sigma^{-1}$) as suggested by maximum likelihood for correctly specified models and the identity matrix ($I$). Since this is not a linear model, it is not possible to directly use the weighting idea proposed in Section 3.5. We propose a modification for the light curve scenario. We first fit the model using identity weights and determine a best fit period. We then determine the optimal weighting at this period following the procedure of Section 3.5. Recall from Section 2 that at a fixed period, the sinusoidal models are linear. Using the new weights, we then refit the model and estimate the period. A period estimate is considered correct if it is within 1% of the true value.