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

Misspecified models often provide useful information about the true data generating distribution. For example, if $y$ is a non-linear function of $x$ the least squares estimator $\hat{\beta}$ is an estimate of $\beta$, the slope of the best linear approximation to the non-linear function. Motivated by problems in astronomy, we study how to incorporate observation measurement error variances into fitting parameters of misspecified models. Our asymptotic theory focuses on the particular case of linear regression where often weighted least squares procedures are used to account for heteroskedasticity. We find that when the response is a non-linear function of the independent variable, the standard procedure of weighting by the inverse of the observation variances can be counter-productive. In particular, ordinary least squares may have lower asymptotic variance. We construct an adaptive estimator which has lower asymptotic variance than either OLS or standard WLS. We demonstrate our theory in a small simulation and apply these ideas to the problem of estimating the period of a periodic function using a sinusoidal model.


Introduction
Misspecified models are common. In prediction problems, simple, misspecified models may be used instead of complex models with many parameters in order to avoid overfitting. In big data problems, true models may be computationally intractable, leading to model simplifications which induce some level of misspecification. In many scientific domains there exist sets of well established models with fast computer implementations. A practitioner with a particular data set may have to choose between using one of these models (even when none are exactly appropriate) and devising, testing and implementing a new model. Pressed for time, the practitioner may use an existing misspecified model. In this work we study how to fit a misspecified linear regression model with heteroskedastic measurement error. Problems involving heteroskedastic measurement error and misspecified models are common in astronomy. We discuss an example in Section 2.
Suppose x i ∈ R p ∼ F X independent across i and σ i ∈ R ∼ F σ independent across i for 1 ≤ i ≤ n. Suppose and Var( i ) = 1 ∀i, independent across i and independent of x i and σ i . Define The parameter β is the slope of the best fitting least squares line. The parameter β may be of interest in several situations. For example, β minimizes mean squared error in predicting y from x among all linear functions, ie β = argmin β E[(y−x T β) 2 ]. Define g(x) = f (x)−x T β. The function g is the non-linear component of f .
When the model is correctly specified (ie g(x) ≡ 0), weighted least squares (WLS) using the inverse of the observation variances as weights is asymptotically normal and has minimum asymptotic variance among all WLS estimators. In the case with model misspecification and x i , σ i independent, we show that WLS estimators remain asymptotically normal. However weighting by the inverse of the observation variances can result in a larger asymptotic variance than other weightings, including ordinary least squares. Using the asymptotic variance formula we determine an optimal weighting which has lower asymptotic variance than standard WLS (using the inverse of the observation variances as weights) and OLS. The optimal weighting function has the form w(σ) = (σ 2 + ∆) −1 where ∆ ≥ 0 is a function of the degree of model misspecification and the design. We find adaptive estimators for w in the cases where the error variances are assumed known and where the error variances belong to one of M groups with group membership known. We also briefly consider the case where x i and σ i are dependent. In this setting the OLS estimator is consistent but weighted estimators are generally not consistent.
This work is organized as follows. In Section 2 we introduce a motivating problem from astronomy and offer some heuristic thinking about misspecified models and heteroskedasticity. For those readers primarily interested in the statistical theory, Section 2 can be skipped. In Section 3 we review some relevant literature and develop asymptotic results for the linear model. We present results for simulated data and the astronomy application in Section 4. We conclude in Section 5.

Misspecified Models and Heteroskedastic Error in Astronomy
Periodic variables are stars that vary in brightness periodically over time. Figure  1a shows the brightness of a single periodic variable star over time. This is known as the light curve of the star. Two sigma uncertainties are plotted as vertical bars around each point. Magnitude is inversely proportional to brightness, so lower magnitudes are plotted higher on the y-axis. This is a periodic variable so the changes in brightness over time are periodic. Using this data one may estimate a period for the star. When we plot the brightness measurements as time modulo period (Figure 1b), the pattern in brightness variation becomes clear. Periodic variables play an important role in several areas of astronomy including extra-galactic distance determination and estimation of the Hubble constant [26,21]. Modern surveys, such as OGLE-III, have collected hundreds of thousands of periodic variable star light curves [28].
Accurate period estimation algorithms are necessary for creating the folded light curve (Figure 1b). A common procedure for determining the period is to perform maximum likelihood estimation using some parametric model for light curve variation. One popular model choice is a sinusoid with K harmonics. Let the data for a single periodic variable be D where y i is the brightness at time t i , measured with known uncertainty σ i . Magnitude variation  is modeled as Here ω is the frequency, a k is the amplitude of the k th harmonic, and φ k is the phase of the k th harmonic. Let a = (a 1 , . . . , a K ) and φ = (φ 1 , . . . , φ K ). Let Ω be a grid of possible frequencies.
The maximum likelihood estimate for frequency is Generalized Lomb-Scargle (GLS) is equivalent to this estimator with K = 1 [32]. The analysis of variance periodogram in [23] uses this model with a fast algorithm for computing ω. We used estimator (2.2) with K = 1, 2 to determine the period of the light curve in Figure 1a. The estimates for period were essentially the same for both K = 1 and K = 2 so in Figure 1b we folded the light curve using the K = 1 estimate. The solid orange line is the maximum likelihood fit for the K = 1 model (notice the sinusoidal shape). The blue dashed line is for the K = 2 model.
While the period estimates are accurate, both models are misspecified. In particular, note that the vertical lines around the brightness measurements are four standard deviations (4σ i ) in width. If the model is correct, we would expect about 95% of these intervals to contain the maximum likelihood fitted curves. For the K = 1 model, 10% of the intervals contain the fitted curve. For K = 2, 37% of observations contain the ML fitted curve. The source of model misspecification is the light curve shape which cannot be perfectly represented by a sinusoid with K = 1, 2 harmonics. The light curve has a long, slow decline and a sudden, sharp increase in brightness.
The parameter fits of misspecified models are estimates of an approximation. In the K = 1 case, the parameter fits are the orange line in Figure 1b and the approximation is the sinusoid which is closest to the true light curve shape. In many cases this approximation may be useful. For example the period of the approximation may match the period of the light curve.
When fitting a misspecified model with heteroskedastic measurement error, one should choose a weighting which ensures the estimator has small variance and thus is likely close to the approximation. The use of the inverse of the observation variances as weights (in Equation (2.3)) is motivated by maximum likelihood theory under the assumption that the model is correct. However as we show in Section 3 for the linear model, these weights are generally not optimal when there is model misspecification.
As a thought experiment, consider the case where one observation has extremely small variance and other observations have much larger variance. The maximum likelihood fitted curve for this data will be very close to the observation with small variance. However the best sinusoidal approximation to the true function at this point may not be particularly close to the true function. Thus using the inverse of observation variances as weights may overweight observations with small variance in the case of model misspecification. We make these ideas precise in Section 3.3.
The choice of weights is not critical for the light curve in Figure 1a because it is well sampled (n > 50), so the period is easy to determine. However in many other cases light curves are more poorly sampled (n ≈ 20), in which case weighting may affect period estimation accuracy.

Sinusoidal Fit and Linear Models
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 argmin ω∈Ω min a,φ,β0 The sum over i can be simplified by noting the linearity of the model and reparameterizing. Let Y = (y 1 , . . . , y n ) T . Let β k1 = a k cos(φ k ) and β k2 = a k sin(φ k ). Define β = (β 0 , β 11 , β 12 , . . . , β K1 , β K2 ) T ∈ R 2K+1 . Let Σ be a n × n diagonal matrix where Σ ii = σ 2 i . Define We rewrite the ML estimator as Every frequency ω in the grid of frequencies Ω determines a design matrix X(ω). At a particular ω, the β which minimizes the objective function is the weighted least squares estimator 3) 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 Ω. 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 Σ −1 as the weight matrix in either Equation (2.3) or (2.4).

Problem Setup and Related Literature
Let X ∈ R n×p be the matrix with row i equal to x T i . Let Y = (y 1 , . . . , y n ) T . Let Σ be the diagonal matrix of observation variances such that Σ ii = σ 2 i . Let W be a diagonal positive definite matrix. The weighted least squares estimator is In this work we seek W which minimize error in estimating 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 θ t is not in the model, the MLE converges to the distribution θ 0 in the model Θ 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 V ar(y i ) = σ 2 /m i where m i is known. A Gaussian prior was placed on β 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 Var( i ) = ξ(x i , z i ) for some function ξ. 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 β to be inconsistent. In contrast, ignoring observation variances leads to consistent estimates of β. Qualitatively these conclusion are similar to our own in that they caution against using weights in the standard way.

Asymptotic Results
Our asymptotic theory makes assumptions on the form of the weight matrix.
These assumptions include both the ordinary least squares (OLS) estimator where W ii = w(σ i ) = 1 and the standard weighted least squares estimator where In both these cases δ nmi = 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.
< ∞, and the variances are bounded below by a positive constant σ 2 min ≡ inf{σ 2 : F σ (σ) > 0} > 0. The major assumption here is independence between x and σ. We address dependence in Section 3.7.
See Section A.1 for a proof. If the response is linear (g ≡ 0) then the variance is This is the standard weighted least squares estimator. This w can be shown to minimize the variance using the Cauchy Schwartz inequality. With w(σ) = 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 σ and x are not independent.

OLS and Standard WLS
For notational simplicity define

8
The asymptotic variances for OLS ( β(I)) and standard WLS ( β(Σ −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 . 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 δ 1 , δ 2 are small non-negative numbers and c > 1 is large. Note that A and B do not depend on F σ .
• δ 1 = 0, δ 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 E[σ 2 ] > c 2 δ 2 , 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 In summary, standard WLS performs better than OLS when there are a small number of observations with large variance. • δ 1 > 0, δ 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 δ 1 small, the coefficient for A is

9
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 δ 1 . Since c > 1, E[σ 2 ] < 1 and The case where both δ 1 and δ 2 are non-zero presents problems for both OLS and standard WLS. For example if δ = δ 1 = δ 2 , both OLS and standard WLS can be made to have large asymptotic variance by setting δ small and c large. In the following section we construct an adaptive weighting which improves upon both OLS and standard WLS.

Improving on OLS and Standard WLS
Let Γ be a linear function from the set of p × p matrices to R such that Γ(C) > 0 whenever C is positive definite. We seek some weighting w = w(σ) for which Γ(ν(w)) (recall that ν is the asymptotic variance) is lower than OLS and standard WLS. Natural choices for Γ include the trace (minimize the sum of variances of the parameter estimates) and the Γ(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. with strict inequality if E[g 2 (x)xx T ] is positive definite and the distribution of σ is not a point mass.
A proof is contained in Section A.3. Thus if we can construct a weight matrix W which satisfies Assumptions 1 with w(σ) = w min (σ) , 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.

Known Error Variances
With the σ i known we only need to estimate A and B in w min in Equation Let β( W ) be a root n consistent estimator of β (eg W = I is root n consistent by Theorem 3.1) and let Then we have The estimated optimal weighting matrix is the diagonal matrix W min with diagonal elements W min,ii = 1 A few notes on this estimator: These estimates are weighted by σ −4 i . The term ( σ −4 i ) −1 normalizes the weights. This weighting is motivated by the fact that Analysis of the first order term shows and Var Thus by weighting the estimates by σ −4 i , we can somewhat account for the different variances. Unfortunately since the variance depends on g, E[ 3 ], and E[ 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, Γ(A)Γ(B) −1 ≥ 0. Thus for estimating ∆, 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 W min starting with a known root n consistent estimator such as 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 ν OR ( W ) which is the same as ν 1 but uses A and B rather than A and B. While ν OR cannot be used in practice, it is useful for evaluating the performance of ν 1 in simulations. Finally suppose the error variance is known up to a constant, i.e. σ 2 i = kτ 2 i where τ 2 i is known but k and σ 2 i are unknown. In the case without model misspecification, one can simply use weights τ −2 i 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 σ i with τ 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.

Unknown Error Variances
Suppose for observation i we observe m i ∈ {1, . . . , M }, the group membership of observation i. Observations in group m have the same (unknown) variance σ 2 m > 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 i , with probability mass function f m (supported on 1, . . . , M ). While the σ m for m = 1, . . . , M are fixed unknown parameters, the probability mass function f m induces the probability distribution function F σ on σ. So we can define for any function h.
Theorem 3.1 shows that even if the σ m were known, standard weighted least squares is not generally optimal for estimating β in this model. It is not possible to estimate w min as proposed in Section 3.5 because that method requires knowledge of σ m . However we can re-express the optimal weight function as w min (m) = 1 where the last equality defines C m . Note that σ m is a fixed unknown parameter, not a random variable. One can estimate B with B = (n −1 X T X) −1 and C m with where β( W ) is a root n consistent estimator of β (for example W = I suffices by Theorem 3.1). The estimated weight matrix 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 W min , starting with some (possibly inefficient) root n consistent estimate of β.
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 (3.8)

Dependent Errors
Suppose one drops the independence assumption between x and σ. 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(σ i ), some function of the error variance. The estimator for β is Recalling we write w for w(σ), we have the following result.
See Section A.6 for a proof. If x and σ are independent then the r.h.s is E[xx T ] −1 E[xf (x)] and the estimator is consistent (as demonstrated by Theorem 3.1). Interestingly the estimator is also consistent if one lets w(σ) = 1 (OLS), regardless of the dependence structure between x and σ. 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.  Table 1 Fraction of times β is in 95% confidence region.

Simulation
We conduct a small simulation study to demonstrate some of the ideas pre-sented in the last section. 1 Consider modeling the function f (x) = x 2 using linear regression with an intercept term. Let x ∼ U nif (0, 1). The best linear approximation to f is β 1 + β 2 x where β 1 = −1/6 and β 2 = 1. We first suppose σ is drawn independently from x from a discrete probability distribution such that P (σ = 0.01) = P (σ = 1) = 0.05 and P (σ = 0.1) = 0.9. Since σ has support on a finite set of values, we can consider the cases where σ i is known (Section 3.5) and where only the group m i of observation i is known (Section 3.6). We let Γ 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 = Σ −1 (standard WLS), W = I 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 ν 1 (Equation (3.6)) as well as the oracle ν OR for estimating the asymptotic variance when the error variances are known. In Section 3.6 we proposed ν 2 (Equation (3.8)) when the error variances are unknown. The estimator ν 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 W = (Σ + ∆) −1 .
we use only ν 2 because ν 1 requires knowledge of Σ. Table 1 contains the results. While for OLS the nominal coverage probability is approximately attained, the other methods are anti-conservative for ν 1 and ν 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 σ, x independence assumption, we now consider the case where σ is a function of x. Specifically, All other parameters in the simulation are the same as before. Note that the marginal distribution of σ 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.

Analysis of Astronomy Data
[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]. 2 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 where t i is the time of the y i brightness measurement made with uncertainty σ i . In Figure 4 we plot magnitude error (σ 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 σ 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 (Σ −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.  Table 2 Fraction of periods estimated correctly using different weightings for models with K = 1, 2, 3 harmonics. Ignoring the observation uncertainties (I) in the fitting is superior to using them (Σ −1 ). The strategy for determining an optimal weight function (∆) does not provide much improvement over ignoring the weights. More complex models (K = 3) perform worse than simple models (K = 1) when there is limited data (n = 10), but better when the functions are better sampled (n = 40). The standard errors on these accuracies is no larger than 0.5(1 − 0.5)/238 ≈ 0.032 .
The fraction of periods estimated correctly are contained in Table 2. In nearly all cases ignoring observation uncertainties (I) outperforms using the inverse of the observation variances as weights (Σ −1 ). The improvement is greatest for the K = 1 model and least for the K = 3 model, possibly due to the decreasing model misspecification as the number of harmonics increases. The very poor performance of the K = 3 models with 10 magnitude measurements is due to overfitting. With K = 3, there are 8 parameters which is too complex a model for 10 observations. Optimizing the observation weights does not appear to improve performance over not using weights. This is potentially due to the fact that the model is highly misspecified (see Figure 1b).

Other Problems in Astronomy
Heteroskedastic measurement error is ubiquitous in astronomy problems. In many cases some degree of model misspecification is present. In this work, we focused on the problem of estimating periods of light curves. Other problems include: • [22] observe the brightness of galaxies through several photometric filters. Variances on the brightness measurements are heteroskedastic. The brightness measurements for each galaxy are matched to a set of templates. Assuming a normal measurement error model, maximum likelihood would suggest weighting the difference between observed brightness and template brightness by the inverse of the observation variance. In personal communication, [22] stated that galaxy templates contained some level of misspecification. [22] addressed this issue by inflating observation variances, using weights of (σ 2 + ∆) −1 instead of σ −2 . The choice of ∆ > 0 was based on qualitative analysis of model fits. Section 3.3 provides a theoretical justification for this practice. • [20] models spectra of galaxies as linear combinations of simple stellar populations (SSP) and non-linear distortions. While parameters which define an SSP are continuous, a discrete set of SSPs are selected as prototypes and the galaxies are modeled as linear combinations of the prototypes. This is done for computational efficiency and to avoid overfitting. However prototype selection introduces some degree of model misspecification as the prototypes may not be able to perfectly reconstruct all galaxy spectra. Galaxy spectra are observed with heteroskedastic measurement error and the inverse of the observation variances are used as weights when fitting the model (see Equation 2.2 in [20]).

Conclusions
We have shown that WLS estimators can perform poorly when the response is not a linear function of the predictors because observations with small variance have too much influence on the fit. In the misspecified model setting, OLS suffers from the usual problem that observations with large variance induce large asymptotic variance in the parameter estimates. For cases in which some observations have very small variance and other observations have very large variance, procedures which optimize the weights may achieve significant performance improvements as shown in the simulation in Section 4.1. This work primarily focused on the case where x and σ are independent. However results from Section 3.7 showed that when independence fails, weighted estimators will typically be biased. This additional complication makes OLS more attractive relative to weighted procedures.
For practitioners we recommend caution in using the inverse of the observation variances as weights when model misspecification is present. As a check, practitioners could fit models twice, with and without weights, and compare performance based on some metric. More sophisticated methods, such as specifically tuning weights for optimal performance may be attempted. Our asymptotic theory provides guidance on how to do this in the case of the linear model.

Appendix A: Technical Notes
A.1. Proof of Theorem 3.1 Let g(X) ∈ R n be the function g applied to the rows of X. We sometimes write w for w(σ). We have In part 1 we show that q In part 2 we show that

Thus by Slutsky's Theorem
Recall that by Assumptions 1 where h is a bounded function, δ nmi are O P (1), and the d is uniformly (in σ) bounded by an O P (1) random variable.
We show that R 1 , R 2 P → 0. Noting that E[|x ij x ik h(σ i )1 mi=m |] < ∞ because h is bounded and the x have second moments we have Using the fact that |d(σ i , δ n )| < δ n where δ n is O P (1) we have where the last equality follows from the facts that σ and x are independent. The desired result follows from the continuous mapping theorem.
x i ] = 0 and i is independent of all other terms and mean 0. We have The desired result now follows from the CLT and showing that R 3 , R 4 P → 0. Note that Thus because the terms inside the i summand are i.i.d. with expectation 0. Finally recalling that the d(σ i , δ n ) is bounded above by δ n which is uniform O P (1), we have

A.3. Proof of Corollary 3.1
We must show Γ(ν(w min )) ≤ min(Γ(ν(I)), Γ(ν(Σ −1 ))) with strict inequality if E[g 2 (x)xx T ] is positive definite and the distribution of σ is not a point mass. The inequality follows from Theorem 3. By assumption β( W ) = β + n −1/2 δ n , thus . Thus by the CLT and Slutsky's Theorem We have where the last equality follow from the facts that the terms inside the sum are i.i.d with expectation C m = E[(g 2 (x) + σ 2 m )xx T ] and nfm(m) n i=1 1m i =m → P 1. Thus we have W min,ii = w min (m i ) + δ nmi n −1/2 which satisfies the form of Assumptions 1.