Robust estimation of non-linear aspects

Christos P. Kitsos Athens University of Business and Economics, Athens, Greece and Christine H. Müller Free University of Berlin, Berlin, Germany

Summary

As a first step for dealing with efficient robust estimation in non-linear models, we regard the problem of efficient robust estimation of non-linear aspects (functions) φ(β) of the unknown parameter β of a linear model. For robust estimation of a general non-linear aspect we propose estimators which are based on one-step-M-estimators and derive their asymptotic behaviour at the contaminated linear model, where the errors have contaminated normal distributions. The asymptotic behaviour provides criteria for robustness and optimality of the estimators and the corresponding designs. Because it is impossible to find globally optimal robust estimators and designs locally optimal solutions are used for efficiency comparisons. Simple formulas for the efficiency rates are given for the general case. Using these results the efficiency rates for estimating robustly the relative variation of a circadian rhythm are calculated. These efficiency rates are very similar to those for non-robust estimation although on principle there is an important difference.

Key words: Contaminated linear model, non-linear aspect, locally optimal robust estimators, locally optimal designs, efficiency rate, circadian rhythm.

1 Introduction

A general linear model

Y_{nN} = x(t_{nN})′β + Z_{nN}, n = 1, ..., N

is considered where Y_{nN} are observations, t_{nN} ∈ T are experimental conditions, x : T → IR^{p} is a known ”regression” function, β ∈ B ⊂ IR^{p} is an unknown parameter vector, Z_{nN} are error variables. In this model a general non-linear aspect φ(β) of β shall be estimated, where φ : B → IR^{r} is a given non-linear function.

In classical linear models it is assumed that the error variables Z_{1N}, ..., Z_{NN} are independent and identically distributed and often it is assumed that they are normally distributed with mean 0 and known or unknown variance σ^{2}, i.e.

\frac{1}{σ} Z_{nN} ∼ P := n_{(0,1)},

where n_{(0,1)} denotes the standard normal distribution. But if some outlying observations (gross errors) may appear the normal distribution is not correct and the classical estimators based on least squares estimators may be baised very much.

For estimating the whole unknown parameter vector several approaches were developed to derive outlier robust estimators (see Hampel et al. (1986) and Rousseeuw and Leroy (1987)). One often used approach is the approach of Hampel (1978) and Krasker (1980) who derived the robustness criterion by calculating the influence function of the estimator, where mainly M-estimators were regarded. A similar approach which leads to the same robustness criterion is the more model oriented and more general approach of Bickel (1981, 1984) and Rieder (1985, 1987). They studied the asymptotic behaviour of general asymptotically linear estimators at contaminated linear models and in particular at conditionally contaminated linear models where the errors Z_{nN} can have different contaminated normal distributions for different experimental conditions. The class of asymptotically linear estimators in particular includes one-step-M-estimators which generalize M-estimators. This approach is especially appropiate for designed experiments and for estimation general linear aspects Cβ, C ∈ IR^{r×p}, of β as Kurotschka and Müller (1992) showed. Therefore we here also use this model oriented approach to study the robustness properties and efficiency properties of estimators for general non-linear aspects φ(β) of β. Because some problems in non-linear models can be transformed in a linear model with an interesting non-linear aspect (see also the example in Section 4 of this paper) the approach of this paper may also a first step for regarding the efficiency of robust estimators in non-linear models.

In Section 2 we describe the conditionally contaminated linear model and introduce the class of regarded estimators for the non-linear aspect which are based on one-step-M-estimators. In this section we also derive the asymptotic behaviour of the estimators at the contaminated linear model which provides the robustness and optimality criterion. As for classical non-robust estimation (see for example Silvey (1980), Buonaccorsi et al. (1986), Kitsos (1986, 1992), Kitsos et al. (1988), McDonald and Studden (1990)) it is impossible to find globally optimal designs for robust estimation. Moreover, in opposite to non-robust estimation, it is also impossible to find globally optimal robust estimators. Only locally optimal designs and robust estimators can be derived and used for efficiency comparisons as is shown in Section 3. In this section it is also shown that the locally optimal solutions have a simple form so that general formulas for the efficiency rates can be derived. In Section 4 the efficiency rates are investigated for the problem of estimating the relative variation of a circadian rhythm. The efficiency comparisons show that on principle there is an important difference to those for classical non-robust estimation which were given by Kitsos et al. (1988). But nevertheless the calculated efficiency rates are very similar to those for classical non-robust estimation. This may imply, as we discuss in Section 5, that strategies which were developed for deriving optimal designs for non-robust non-linear estimation problems can be transferred to robust non-linear estimation problems.

2 Robust estimators for non-linear aspects and their asymptotic behaviour

For robust estimation of a non-linear aspect φ(β) we consider estimators of the form φ̂_{N} = φ(β̂_{N}) where β̂_{N} is an estimator for β. An appropriate class of estimators for β which includes the least squares estimator as well as many robust estimators is the class of one-step-M-estimators as were regarded in Bickel (1975) and Müller (1992). These one-step-M-estimators are given by

β̂_{N}(y_{N}, d_{N}) = β̂_{N}^{0}(y_{N}, d_{N}) + \frac{1}{N} ∑_{n=1}^{N} ψ ( \frac{y_{nN} −x(t_{nN})′_{β_{N}^{0}} (y_{N} ,d_{N})}{_{σ_{N}}}, t_{nN} ) σ̂_{N}

where d_{N} := (t_{1N}, ..., t_{NN})′, y_{N} := (y_{1N}, ..., y_{NN})′ and ψ : IR × T → IR^{p} is called the score function of the estimator. The initial estimator β̂_{N}^{0} and the variance estimator σ̂_{N} can be robust or non-robust estimators, for example the initial estimator can be the least squares estimator or some M-estimator, and the variance estimator can be the mean squared residuals or Huber’s Proposal 2 (Huber (1964)).

For these estimators we now derive their asymptotic behaviour at a conditionally contaminated linear model which was also regarded by Bickel (1981, 1984), Rieder (1987) and Müller (1992, 1994). In such a model it is assumed that the error variables Z_{1N}, ..., Z_{NN} of the linear model are also independent but distributed according to a contaminated normal distribution where the contamination may be different for different experimental conditions. I.e.

\frac{1}{σ} Z_{nN} ∼ Q_{nN}(dz) := (1 − N^{−1/2}R ϵ(t_{nN})) n_{(0,1)}(dz) (2.1)

+ N^{−1/2}R ϵ(t_{nN}) g(z, t_{nN}) n_{(0,1)}(dz)

with R > 0, \frac{1}{N} ∑_{n=1}^{N} ϵ(t_{nN}) ≤ 1, ∫ g(z, t) n_{(0,1)}(dz) = 1, g(z, t) ≥ 0 for all z ∈ IR, t ∈ T. Thereby the markov kernel g(·, t)n_{(0,1)} models the form and ϵ(t) ≥ 0 the proportion of contamination. The set P_{R} of all sequences (Q^{N} := ⨂_{n=1}^{N} Q_{nN})_{N∈IN} defines a conditionally contamination neighbourhood around the classical model (P^{N})_{N∈IN}. If (Z_{1N}, ..., Z_{NN}) is distributed according to Q^{N} and β is the true parameter vector then Q_{β}^{N} will denote the distribution of (Y_{1N}, ..., Y_{NN}), in particular P_{β}^{N} the distribution of (Y_{1N}, ..., Y_{NN}) for the classical model.

To derive the asymptotic behaviour of the estimators we assume that a finite set T_{0} and an asymptotic design measure δ exist with δ(T_{0}) = 1, {t_{1N}, ..., t_{NN}} ⊂ T_{0} for all N ∈ IN, lim_{N→∞} \frac{1}{N} ∑_{n=1}^{N} e_{t_{nN}} (\{t\}) = δ(\{t\}) for all t ∈ T_{0} and

∫ ψ(z, t)x(t)′zP (dz)δ(dt) = E_{p×p}, (2.2)

where e_{t} denotes the Dirac measure on t and E_{p×p} the identity matrix. Then under some regularity conditions on the initial and the variance estimator an one-step-M-estimator β̂_{N} for β is asymptotically normally distributed in the conditionally contaminated linear model (see Müller (1992)). If φ(β) is continuously differentiable on an open set B ⊂ IR^{p}, where

φ̇_{β} := \frac{∂}{∂β} φ(β) ∈ IR^{r×p}

denotes the matrix of partial derivatives of φ with respect to the elements of β, then also φ(β̂_{N}) is asymptotically normally distributed in the conditionally contaminated linear model. I.e. for all β ∈ B and all (Q^{N})_{N∈IN} ∈ P_{R} we have

L(\sqrt{N}(φ(β̂_{N}) − φ(β))|Q_{β}^{N}) ⟶^{N→∞} N (φ̇_{β} b(ψ, (Q^{N})_{N∈IN}), σ^{2} V (ψ, δ, β)) (2.3)

with

V (ψ, δ, β) := φ̇_{β} ∫ ψψ′ d(P ⊗δ) φ̇′_{β}

and maximum asymptotic bias

sup{|φ̇_{β} b(ψ, (Q^{N})_{N∈IN})|; (Q^{N})_{N∈IN} ∈ P_{R}} = R σ ess sup_{P ⊗δ} |φ̇_{β} ψ| .

If the score function of β̂_{N} satisfies ess sup_{P ⊗δ} |ψ| < ∞, i.e. β̂_{N} is robust, then also ess sup_{P ⊗δ} |φ̇_{β} ψ| < ∞ so that φ(β̂_{N}) is robust. Moreover we can set without loss of generality R = 1 = σ.

3 Efficiency of robust estimators and designs

In (2.3) it is shown that the maximum asymptotic bias ess sup_{P ⊗δ} |φ̇_{β} ψ| and the asymptotic covariance matrix V (ψ, δ, β) of an estimator φ(β̂_{N}) based on an one-step-M-estimator β̂_{N} for β with score function ψ depend on β. If there is no bias bound for the maximum asymptotic bias then the matrix ∫ ψ ψ′ d(P ⊗δ) is minimized uniformly by ψ_{∞}(z, t) := I(δ)^{−1}x(t)z within all score functions ψ fulfilling condition (2.2), where we have I(δ) := ∫ x(t)x(t)′ δ(dt). I.e. the covariance matrix V (ψ_{∞}, δ, β) is the minimum for all β so that the estimator based on the least squares estimator is asymptotically optimal. (Note that the least squares estimator behaves asymptotically like an one-step-M-estimator with score function ψ_{∞}.) But if the maximum asymptotic bias should be bounded by some bias bound b, in principle it is impossible to minimize ∫ ψ ψ′ d(P ⊗δ) uniformly so that an optimal robust estimator cannot be found independently of β. Also optimal designs δ cannot be found independently of β, and this problem already occurs for classical nonrobust estimation. All classical non-linear optimal experimental design theory suffer on this β-dependence so that strategies were invented to overcome this problem. See for an overview Ford et al. (1989).

One basic strategy is to compare a given estimator and design with the locally optimal estimators and designs. This provides the efficiency rate of the estimator and the design and this efficiency rate shall be investigated in this section.

For deriving the locally optimal estimators and designs for robust estimation we regard the local behaviour of an estimator φ(β̂_{N}) based on an robust one-step-M-estimator β̂_{N} for β with score function ψ. Locally, i.e. for given β, this estimator behaves asymptotically like

φ̇_{β}β̂_{N} = φ̇_{β}β̂_{N}^{0} + \frac{1}{N} ∑_{n=1}^{N} φ̇_{β}ψ ( \frac{y_{nN} −x(t_{nN})′_{β_{N}^{0}} (y_{N} ,d_{N})}{_{σ_{N}}}, t_{nN} ) σ̂_{N},

i.e. like an one-step-M-estimator for the linear aspect φ̃(β̃) = φ̇_{β}β̃ with score function φ̇_{β}ψ (see Müller (1992)). Now we can try to find locally optimal robust estimators within all these robust one-step-M-estimators for the linear aspect φ̃(β̃) = φ̇_{β}β̃. I.e., given a bias bound b for estimating this linear aspect, we have to find a score function ψ_{β,b,δ} which solves

ψ_{β,b,δ} = arg min{tr V (ψ, δ, β); φ̇_{β}ψ ∈ Ψ_{δ}(φ̇_{β}) with ess sup_{P ⊗δ} |φ̇_{β}ψ| ≤ b},

where Ψ_{δ}(φ̇_{β}) is the class of all score functions for estimating the linear aspect φ̃(β̃) = φ̇_{β}β̃ and in particular satisfying ∫ ψ(z, t)x(t)′zP (dz)δ(dt) = φ̇_{β}.

Similarly locally optimal designs for robust estimation with bias bound b can be defined as solutions δ_{β,b} of

δ_{β,b} = arg min{tr V (ψ_{β,b,δ}, δ, β); δ ∈ Δ}

where Δ is the set of the regarded design measures.

For estimating a general linear aspect Cβ the solutions concerning optimal robust estimators are characterized in Kurotschka and Müller (1992) and the solutions concerning optimals designs for robust estimation are characterized in Müller (1994). In particular in Müller (1994) it is shown that the classical A-optimal designs, which minimize tr C I(δ)^{−}C′ (where I(δ)^{−} denotes the g-inverse of I(δ)), are also optimal for robust estimation and this holds independently of the bias bound b. Therefore the locally optimal designs for estimation without bias bound are also locally optimal for robust estimation with bias bound b, i.e.

δ_{β,b} = δ_{β} := arg min{tr φ̇_{β} I(δ)^{−}φ̇′_{β}; δ ∈ Δ}.

Moreover, the score functions

ψ_{β,b} := ψ_{β,b,δ_{β}}

of the locally optimal robust estimators with bias bound b at the locally optimal designs δ_{β} have a very simple form, namely (see Müller (1994))

(3.1) ψ_{β,b}(z, t) = \{ \begin{matrix}I(δ_{β})^{−}x(t) sgn(z) \sqrt{\frac{π}{2}} & for b = \sqrt{tr φ̇_{β} I(δ_{β})^{−}φ̇′_{β} \frac{π}{2}} \\ I(δ_{β})^{−}x(t) sgn(z)_{\frac{min\{|z|,b y_{b,β}\}}{y_{b,β} \sqrt{tr φ̇_{β} I(δ_{β})^{−}φ̇′_{β}}}} & for b > \sqrt{tr φ̇_{β} I(δ_{β})^{−}φ̇′_{β} \frac{π}{2}} \end{matrix} .

Thereby, y_{b,β} is the positive solution of y\sqrt{tr φ̇_{β} I(δ_{β})^{−} φ̇′_{β}} = 2 Φ(b y) − 1 while Φ denotes the distribution function of the normal distribution n_{(0,1)}. Hence, locally optimal robust estimators at locally optimal designs are easy to compute and can be applied for efficiency comparisons.

Now, a given design δ and a given estimator φ(β̂_{N}) for φ(β), where β̂_{N} is a robust onestep-M-estimator for β with bias bound b and score function ψ, can be compared with the locally optimal robust estimator at the locally optimal design. Because, according to (2.3), the maximum asymptotic bias of φ(β̂_{N}) is equal to

b(β) := b(β, ψ, δ) := ess sup_{P ⊗δ} |φ̇_{β} ψ|,

i.e. depends on β, the estimator φ(β̂_{N}) at the design δ should be only compared with locally optimal robust estimators with the same local bias bound b(β) at the locally optimal design δ_{β}, to make fair comparisons. Hence the asymptotic efficiency measure is defined as

E_{β}(ψ, δ) := \frac{tr V (ψ_{β,b(β)}, δ_{β}, β)}{tr V (ψ, δ, β)} ∈ [0, 1]

where ψ_{β,b(β)} is the score function of the locally optimal one-step-M-estimator with local bias bound b(β).

If δ is an A-optimal design for estimating β then the score function ψ_{b} of the optimal robust one-step-M-estimator for β at δ with bias bound b has the form (see Müller (1994))

(3.2) ψ_{b}(z, t) = \{ \begin{matrix}I(δ)^{−1}x(t) sgn(z) \sqrt{\frac{π}{2}} & for b = \sqrt{tr I(δ)^{−1} \frac{π}{2}} \\ I(δ)^{−1}x(t) sgn(z)\frac{min\{|z|,b y_{b}\}}{y_{b} \sqrt{tr I(δ)^{−1}}} & for b > \sqrt{tr I(δ)^{−1} \frac{π}{2}} \end{matrix} ,

where y_{b} is the positive solution of y\sqrt{tr I(δ)^{−1}} = 2 Φ(b y) − 1. Hence, in this case with the equivalence theorem for linear optimality (see Federov (1972), p. 125), we get for the local bias bounds

(3.3) b(β) = \frac{b}{\sqrt{trI(δ)^{−1}}} ess sup_{δ} |φ̇_{β} I(δ)^{−1} x(t)| = ess sup_{δ} \frac{|φ̇_{β} I(δ)^{−1} x(t)|}{|I(δ)^{−1} x(t)|} b,

where b(β) ≥ \sqrt{tr φ̇_{β} I(δ_{β})^{−} φ̇′_{β} \frac{π}{2}} for all b ≥ \sqrt{tr I(δ)^{−1} \frac{π}{2}}. Then the efficiency measure is

E_{β}(ψ_{b}, δ) = \{ \begin{matrix}\frac{g(b(β) y_{b,β})}{y_{b,β}^{2} tr φ̇_{β} I(δ)^{−1} φ̇′_{β}} \frac{2}{π} & for b = \sqrt{tr I(δ)^{−1} \frac{π}{2}} \\ \frac{g(b(β) y_{b,β})}{g(b y_{b})} · \frac{y_{b}^{2}}{y_{b,β}^{2}} · \frac{tr I(δ)^{−1}}{tr φ̇_{β} I(δ)^{−1} φ̇′_{β}} & for b > \sqrt{tr I(δ)^{−1} \frac{π}{2}} \end{matrix} .

Because of lim_{v→∞} g(v) = 1, lim_{b→∞} y_{b} = \frac{1}{\sqrt{tr I(δ)^{−1}}} and lim_{b→∞} y_{b,β} = _{\frac{1}{\sqrt{tr φ̇_{β} I(δ_{β})^{−} φ̇′_{β}}}} it

is easy to see that for b → ∞ the efficiency E_{β}(ψ_{b}, δ) tends to

\frac{tr φ̇_{β} I(δ_{β})^{−} φ̇′_{β}}{tr φ̇_{β} I(δ)^{−1} φ̇′_{β}} , (3.4)

i.e. to the design efficiency of δ relative to δ_{β} as appeared in classical non-robust estimation (see Kitsos et al. (1988)).

4 Example: The rhythmometry problem

Consider the cosinor model of a circadian rhythms which was regarded in Kitsos et al. (1988). In such a model the observations are given by the non-linear model

Y_{nN} = θ_{0} + θ_{1} cos(2πt_{nN} + θ_{2}) + Z_{nN}

with t_{nN} ∈ T := [0, 1) and θ_{0} being the mean, θ_{1} the amplitude and θ_{2} the phase. In Kitsos et al. (1988) the interest was focussed on estimating the ratio

g = g(θ_{0}, θ_{1}) = \frac{θ_{1}}{θ_{0}}

as it gives a potentially useful measure of relative variation, falling usually in practice within the range 0 < g < 1.

By reparameterization one gets a linear model

Y_{nN} = β_{0} + β_{1} cos(2πt_{nN}) + β_{2} sin(2πt_{nN}) + Z_{nN}
= x(t_{nN})′β + Z_{nN}

where β_{0} = θ_{0}, β_{1} = θ_{1} cos(θ_{2}), β_{2} = −θ_{1} sin(θ_{2}), x(t) = (1, cos(2πt), sin(2πt))′ and β = (β_{0}, β_{1}, β_{2})′. Therefore the aspect g(θ_{0}, θ_{1}) of the non-linear cosinor model is a non-linear aspect

φ(β) := \frac{\sqrt{β_{1}^{2} + β_{2}^{2}}}{β_{0}} = g(θ_{0}, θ_{1})

of a linear model. Note that φ̇_{β} = \frac{1}{θ_{0}} (−\frac{θ_{1}}{θ_{0}} , \frac{β_{1}}{θ_{1}} , \frac{β_{2}}{θ_{1}} ).

In Kitsos et al. (1988) is was assumed that the errors Z_{nN} are normally distributed with mean 0 and constant variance σ^{2}. But there may exist some outlying observations so that the errors have a conditionally contaminated normal distribution of form (2.1). Without loss of generality we can here assume that the variance σ^{2} as well as the radius R of the contamination neighbourhood P_{R} are equal to 1.

In a conditionally contaminated cosinor model the non-linear aspect φ(β) should be estimated with a robust estimator, for example with an estimator of the form φ(β̂_{N}), where β̂_{N} is an optimal robust one-step-M-estimator for β. If the design

δ = \frac{1}{3}(e_{0} + e_{\frac{1}{3}} + e_{\frac{2}{3}} )

is used then this design is A- and D-optimal for estimating β (see Kitsos et al. (1988)) and the score function of the optimal robust one-step-M-estimator for β with bias bound b > \sqrt{5\frac{π}{2}} results according to (3.2) in

ψ_{b}(z, t) = (1, 2 cos(2πt), 2 sin(2πt))′ \frac{sgn(z) min\{|z|, b y_{b}\}}{\sqrt{5}y_{b}}.

Thereby y_{b} is the positive solution of y\sqrt{5} = 2Φ(b y) − 1, i.e. y_{b} = 0.216 for b = 3 and y_{b} = 0.444 for b = 6, and \sqrt{5\frac{π}{2}} = \sqrt{tr I(δ)^{−1} \frac{π}{2}} is the minimum possible bias bound. Then, according to (3.3) the maximum asymptotic bias of the estimator φ(β̂_{N}) is

b(β) = \frac{b}{\sqrt{5}θ_{0}} max \{|\frac{−θ_{1}}{θ_{0}} + 2 cos(2πt + θ_{2})|; t ∈ \{0, \frac{1}{3}, \frac{2}{3}\}\} .

These estimator should be compared with the locally optimal robust one-step-M-estimators for φ(β) with local bias bounds b(β) at the locally optimal designs.

Locally optimal designs, i.e. A-optimal designs for estimating φ̇_{β}β, are (see Kitsos et al. (1988))

δ_{β} = \frac{1}{2} ((1 − \frac{θ_{1}}{θ_{0}} ) e_{\frac{−θ_{2}}{2π}} + (1 + \frac{θ_{1}}{θ_{0}} ) e_{\frac{π−θ_{2}}{2π}} ) (4.1)

if 0 < \frac{θ_{1}}{θ_{0}} < 1. According to (3.1) the score function of the locally optimal robust onestep-M-estimator for φ(β) with bias bound b(β) at δ_{β} has the form

(4.2) ψ_{β,b(β)}(z, t) = I(δ_{β})^{−} ( \begin{matrix}1 \\ cos(2πt) \\ sin(2πt) \end{matrix} ) \frac{sgn(z) min\{|z|, b(β) y_{b,β}\} θ_{0}}{y_{b,β}} ,

where y_{b,β} is the positive solution of \frac{y}{θ_{0}} = 2Φ(b(β) y) − 1. Always this solution will exist because b(β) > \frac{1}{θ_{0}} \sqrt{\frac{π}{2}} if b > \sqrt{5\frac{π}{2}}, where \frac{1}{θ_{0}} \sqrt{\frac{π}{2}} = \sqrt{tr φ̇_{β} I(δ_{β})^{−}φ̇′_{β} \frac{π}{2}} is the locally minimum possible bias bound.

Table 1 shows the efficiencies E_{β}(ψ_{b}, δ) of the estimator φ(β̂_{N}) for nine different values for θ_{1}, five different values for θ_{2} and three different values for b including the case b = ∞, i.e. the classical non-robust case. For θ_{0} only the value θ_{0} = 1 was used because the efficiencies are independent of θ_{0}. This is due to the fact that

b(β) y_{b,β} = \frac{b}{\sqrt{5}} max \{|\frac{−θ_{1}}{θ_{0}} + 2 cos(2πt + θ_{2})|; t ∈ \{0, \frac{1}{3}, \frac{2}{3}\}\} \frac{y_{b,β}}{θ_{0}}
= \frac{b}{\sqrt{5}} max \{|\frac{−θ_{1}}{θ_{0}} + 2 cos(2πt + θ_{2})|; t ∈ \{0, \frac{1}{3}, \frac{2}{3}\}\} (2Φ(b(β) y_{b,β}) − 1)

Table 1: Relative efficiencies E_{β}(ψ_{b}, δ) and, in brackets, local bias bounds b(β) for the rhythmometry problem with θ_{0} = 1

b = 3  b = 6  b = ∞
θ_{1}  θ_{2} = 0 θ_{2} = \frac{π}{12} θ_{2} = \frac{π}{6} θ_{2} = \frac{π}{4} θ_{2} = \frac{π}{3}  θ_{2} = 0 θ_{2} = \frac{π}{3}
0.1  0.414 0.414 0.414 0.414 0.413  0.497  0.498
(2.549) (2.458) (2.458) (2.726) (2.817) (5.098) (5.635)
0.2 0.408 0.409 0.408 0.407 0.407 0.490 0.490
(2.415) (2.324) (2.592) (2.860) (2.952) (4.830) (5.903)
0.3 0.399 0.399 0.398 0.397 0.397 0.478 0.478
(2.281) (2.300) (2.726) (2.994) (3.086) (4.562) (6.172)
0.4 0.388 0.386 0.385 0.385 0.384 0.462 0.463
(2.147) (2.434) (2.860) (3.129) (3.220) (4.293) (6.440)
0.5 0.374 0.370 0.369 0.369 0.369 0.444 0.444
(2.012) (2.568) (2.995) (3.263) (3.354) (4.025) (6.708)
0.6 0.355 0.352 0.352 0.352 0.352 0.423 0.424
(2.147) (2.702) (3.129) (3.397) (3.488) (4.293) (6.977)
0.7 0.335 0.334 0.334 0.333 0.333 0.401 0.402
(2.281) (2.837) (3.263) (3.531) (3.622) (4.562) (7.245)
0.8 0.316 0.315 0.315 0.315 0.315 0.378 0.379
(2.415) (2.971) (3.397) (3.665) (3.757) (4.830) (7.513)
0.9 0.296 0.296 0.296 0.296 0.296 0.355 0.356
(2.549) (3.105) (3.531) (3.799) (3.891) (5.098) (7.782)

and

y_{b,β}^{2} tr φ̇_{β}I(δ)^{−1}φ̇′_{β} = tr φ̇_{β}I(δ_{β})^{−}φ̇′_{β} y_{b,β}^{2}  \frac{tr φ̇_{β}I(δ)^{−1}φ̇′_{β}}{tr φ̇_{β}I(δ_{β})^{−}φ̇′_{β}}
= \frac{y_{b,β}^{2}}{θ_{0}^{2}} (\frac{θ_{1}^{2}}{θ_{0}^{2}} + 2)
= (2Φ(b(β) y_{b,β}) − 1)^{2} (\frac{θ_{1}^{2}}{θ_{0}^{2}} + 2)

depend only on the ratio \frac{θ_{1}}{θ_{0}} and the phase θ_{2}. But for non-robust estimation, i.e. for b = ∞, the efficiencies (given by (3.4)) depend only on the ratio \frac{θ_{1}}{θ_{0}} and not on the mean θ_{0} and the phase θ_{2} as already Kitsos et al. (1988) have shown. Hence the main difference between robust and non-robust estimation is the additional dependence on the phase θ_{2} for robust estimation. For b < ∞ Table 1 provides also the values for the local bias bounds b(β).

5 Discussion

As Table 1 for the rhythmometry problem shows, the efficiencies for robust estimation (b < ∞) are very similar to the efficiencies for classical non-robust estimation (b = ∞), in particular for the bias bound b = 6 but also for the bias bound b = 3 which is near to the minimum possible bias bound of \sqrt{5\frac{π}{2}} = 2.802. The main difference to nonrobust estimation is that the efficiency rates depend additionally on the phase θ_{2}. But the influence of the phase θ_{2} is very low, although also those values θ_{2} were used which produce the largest (θ_{2} = \frac{π}{3}) and the smallest (θ_{2} = 0, at least for \frac{θ_{1}}{θ_{0}} ≥ 0.5) values for the local bias bounds b(β). Even at b = 6 the efficiencies were identical at θ_{2} = 0 and θ_{2} = \frac{π}{3}, and this coincides with the behaviour of the efficiencies for non-robust estimation.

From the results for the rhythmometry problem one can expect that for general nonlinear problems the efficiencies for robust estimation are very similar to the efficiencies for non-robust estimation. Hence the strategies which were investigated for deriving optimal designs for non-robust non-linear estimation problems as maximin efficient designs, fully sequential designs, batch sequential designs and two-stage designs (see Silvey (1980), Ford et al. (1989), Kitsos (1986, 1989, 1992)) may be also used for robust non-linear estimation problems. These strategies may also provide strategies for deriving particular robust estimators. For example if we have a two-stage design with N_{1} observations at the first stage and N_{2} observations at the second stage, then we can use in the second stage the locally optimal design δ_{β^{1}} and the locally optimal robust estimator with score function ψ_{_{β^{1}},b}, where the estimator β̂^{1} for β is obtained in the first stage.

For the rhythmometry problem a very practicable two-stage design is the following design: In the first stage a four point design of the form d_{4} := (τ, τ + \frac{1}{4}, τ + \frac{1}{2}, τ + \frac{3}{4})′ with τ ∈ [0, \frac{1}{4}) is used. Because of cos(2πt)^{2} + sin(2πt)^{2} = 1 this design is D-optimal for non-robust estimation of β for all τ ∈ [0, \frac{1}{4}) (see Federov (1972), p.75). Moreover this design is approximately locally optimal for all relevant β and for τ = 0 the least squares estimator β̂_{4}^{∗}, say, for estimating β has a very simple form, namely

β̂_{4}^{∗}((y_{1}, y_{2}, y_{3}, y_{4})′, d_{4}) = ( \begin{matrix}\frac{1}{4}(y_{1} + y_{2} + y_{3} + y_{4}) \\ \frac{1}{2}(y_{1} − y_{3}) \\ \frac{1}{2}(y_{2} − y_{4}) \end{matrix} ) .

Basing on this estimator in the second stage the locally optimal design δ_{β_{4}^{∗}} given by (4.1) and the locally optimal robust one-step-M-estimator with score function ψ_{_{β_{4}^{∗}},b} given by (4.2) can be used. Note that at the locally optimal design δ_{β_{4}^{∗}} the whole parameter vector β is not anymore estimable. But this design has besides its greater efficiency the advantage that only at two different experimental conditions the observations are made, which is more practicable than a three or four point design, in particular for medical investigations of a circadian rhythm.

REFERENCES

[1] BICKEL, P.J. (1975). One-step Huber estimates in the linear model. Journal of the American Statistical Association 70, 428-434.

[2] BICKEL, P.J. (1981). Quelque aspects de la statistique robuste. In École d’Été de Probabilités de St. Flour. Springer Lecture Notes in Mathematics 876, 1-72.

[3] BICKEL, P.J. (1984). Robust regression based on infinitesimal neighbourhoods. The Annals of Statistics 12, 1349-1368.

[4] BUONACCORSI, J.P. and IYER, H.K. (1986). Optimal designs for ratios of linear combinations in the general linear model. Journal of Statistical Planning and Inference 13, 345-356.

[5] FEDEROV, V.V. (1972). Theory of Optimal Experiments. Academic Press, New York.

[6] FORD, I., KITSOS, C.P. and TITTERINGTON, D.M. (1989). Recent advances in nonlinear experimental design. Technometrics 13, 49-60.

[7] HAMPEL, F.R. (1978). Optimally bounding the gross-error-sensitivity and the influence of position in factor space. Proceedings of the ASA Statistical Computing Section, ASA, Washington, D.C., 59-64.

[8] HAMPEL, F.R., RONCHETTI, E.M., ROUSSEEUW, P.J. and STAHEL, W.A. (1986). Robust Statistics - The Approach Based on Influence Functions. John Wiley, New York.

[9] HUBER, P.J. (1964). Robust estimation of a location parameter. Annals of Mathematical Statistics 35, 73-101.

[10] KITSOS, C.P. (1986). Design and inference for nonlinear problems. Ph.D. thesis, University of Glasgow, U.K.

[11] KITSOS, C.P. (1989). Fully-sequential procedures in nonlinear design problems. Comp. Statist. & Data Analysis 8, 13-19.

[12] KITSOS, C.P. (1992). Quasi-sequential procedures for the calibration problem. In COMPSTAT 1992, Vol. 2, Y. Dodge and J. Whittaker (eds.), 227-231, Physica-Verlag.

[13] KITSOS, C.P., TITTERINGTON, D.M. and TORSNEY, B. (1988). An optimal design problem in rhythmometry. Biometrics 44, 657-671.

[14] KRASKER, W.S. (1980). Estimation in linear regression models with disparate data points. Econometrica 48, 1333-1346.

[15] KUROTSCHKA, V. and MÜLLER, Ch.H. (1992). Optimum robust estimation of linear aspects in conditionally contaminated linear models. The Annals of Statistics 20, 331-350.

[16] McDONALD, G.C. and STUDDEN, W.J. (1990). Design aspects of regressionbased ratio estimation. Technometrics 32, 417-424.

[17] MÜLLER, Ch.H. (1992). One-step-M-estimators in conditionally contaminated linear models. Preprint No. A-92-11, Freie Universität Berlin, Fachbereich Mathematik. To appear in Statistics & Decisions.

[18] MÜLLER, Ch.H. (1994). Optimal designs for robust estimation in conditionally contaminated linear models. Journal of Statistical Planning and Inference 38, 125-140.

[19] RIEDER, H. (1987). Robust regression estimators and their least favorable contamination curves. Statistics & Decisions 5, 307-336.

[20] ROUSSEEUW, P.J. and LEROY, A.M. (1987). Robust Regression and Outlier Detection. John Wiley, New York.

[21] SILVEY, S.D. (1980). Optimal Design. Chapman and Hall, London.