Implementing GLS Recall the assumptions of Approach 9: E( Y | x ) = - - PowerPoint PPT Presentation

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Implementing GLS

Recall the assumptions of Approach 9: E(Y |x) = f (x, β), var(Y |x) = σ2g(β, θ, x)2. Here: β is, as before, the vector of parameters in the mean function; θ is a possible additional parameter in the structure of the variance function; σ2 is an additional dispersion parameter.

1 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Ad hoc estimation scheme:

i

Get initial estimate of β using OLS;

ii

Get initial estimate of θ, if needed, and construct initial estimated weights ˆ wj = 1 g

• ˆ

βOLS, ˆ θ, x 2

iii Re-estimate β using WLS, treating ˆ

wj as fixed: solve the estimating equation

n

• j=1

ˆ wj {Yj − f (xj, β)} fβ(xj, β) = 0.

2 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Digression

Gauss-Newton method for WLS (including OLS): The equation

n

• j=1

wj {Yj − f (xj, β)} fβ(xj, β) = 0. generally cannot be solved in closed form. But if β∗ is close to the solution β, f (xj, β) ≈ f (xj, β∗) + fβ(xj, β∗)T (β − β∗) , fβ(xj, β) ≈ fβ(xj, β∗) + fββ(xj, β∗) (β − β∗)

3 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Omitting terms likely to be small, the WLS estimating equation becomes

• n
• j=1

wjfβ(xj, β∗) fβ(xj, β∗)T

• (β − β∗)

≈

n

• j=1

wj {Yj − f (xj, β∗)} fβ(xj, β∗)

4 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Write X(β)

(n×p)

=    fβ(x1, β)T . . . fβ(xn, β)T    f(β)

(n×1)

=    f (x1, β) . . . f (xn, β)    W

(n×n) =

     w1 . . . w2 . . . . . . . . . ... . . . . . . wn     

5 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Then the approximate equation may be written

• X(β∗)TWX(β∗)
• (β − β∗) ≈ X(β∗)TW {Y − f(β∗)}
• r (if the inverse exists)

β ≈ β∗ +

• X(β∗)TWX(β∗)

−1 X(β∗)TW {Y − f(β∗)} Iterative solution: β(a+1) = β(a)+

• X
• β(a)

T WX

• β(a)

−1 X

• β(a)

T W

• Y − f
• β(a)
• Note that W is fixed throughout the inner circle of iteration.

6 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Note that if the iteration converges, β(a) → β(∞), and X

• β(a)

T WX

• β(a)
• converges to a non-singular matrix, then

X

• β(∞)

T W

• Y − f
• β(∞)
• = 0

which is the original estimating equation written in matrix form. That is, the limit β(∞) does solve the original WLS problem. This algorithm, suitably refined, is the default method in both R’s nls() and SAS’s proc nlin.

7 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Instead of the nested iterations, we could solve the last equation directly. Note This is not equivalent to minimizing Sg(β) =

n

• j=1

1 g(β, θ, xj)2 {Yj − f (xj, β)}2 , because differentiating Sg(β) w.r.t. β also brings in the derivative

• f g(·).

8 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

The equation can be solved by a Gauss-Newton method; a similar approximation leads to the iteration β(a+1) = β(a)+

• X
• β(a)

T W(a)X

• β(a)

−1 X

• β(a)

T W(a)

• Y − f
• β(a)
• where the weight matrix W is now updated using the current value of

β(a). That is, the weights are updated within the iteration, instead of being held constant, as in the WLS iteration; a nonlinear instance of Iteratively Reweighted Least Squares, IRWLS (note the redundancy!)

• r IWLS.

9 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Estimating σ2

By analogy with WLS, the natural estimator of σ2 is either 1 n

n

• j=1
• Yj − f
• xj, ˆ

βGLS 2 g

• ˆ

βGLS, θ, xj 2

• r

1 n − p

n

• j=1
• Yj − f
• xj, ˆ

βGLS 2 g

• ˆ

βGLS, θ, xj 2 .

10 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

For fixed weights, the latter is unbiased. For fixed weights and Gaussian Y , the former is ML and the latter is REML. The latter is reported by most software.

11 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Summary

Model E(Y |x) = f (x, β), var(Y |x) = σ2g(β, θ, x)2. Generalized least squares (GLS)

n

• j=1

wj {Yj − f (xj, β)} fβ(xj, β) = 0.

12 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Motivation Loss function Sg(β) =

n

• j=1

wj {Yj − f (xj, β)}2 , with plugged-in weights wj = g(β, θ, xj)−2. Gaussian distribution with mean f (xj, β) and plugged-in weights wj.

13 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Implementations

3-step GLS

1

Get initial estimate ˆ β

(0)

2

Repeat until convergence, or for a fixed number of C steps:

1

Update the weights ˆ wj = g(ˆ β

(t), x)−2

2

Estimate β using WLS: Repeat until convergence β(a+1) = β(a)+

• X
• β(a)

T WX

• β(a)

−1 X

• β(a)

T W

• Y − f
• β(a)
• 3

Update the estimate ˆ β

(t+1)

14 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Comment There are 2 loops In the inner loop (step 2.2), the weight matrix W is fixed

15 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Implementations (continued)

Iteratively reweighted least squares (IRWLS)

1

Get initial estimate ˆ β

(0)

2

Repeat until convergence:

1

Update the weights ˆ wj = g(ˆ β

(t), x)−2

2

Estimate β using WLS β(a+1) = β(a)+

• X
• β(a)

T W(a)X

• β(a)

−1 X

• β(a)

T W(a)

• Y − f
• β(a)
• 3

Update the estimate ˆ β

(t+1)

16 / 17

ST 762 Nonlinear Statistical Models for Univariate and Multivariate Response

Comment There is only 1 loop The weight matrix W is iteratively updated

17 / 17