# MMreg

MMreg computes MM estimator of regression coefficients

## Syntax

• out =MMreg(y,X)example
• out =MMreg(y,X,Name,Value)example
• [out , varargout]=MMreg(___)example

## Description

 out =MMreg(y, X) MMreg with all default options.

 out =MMreg(y, X, Name, Value) MMreg with optional input arguments.

 [out , varargout] =MMreg(___) MMreg with optional input arguments.

## Examples

expand all

### MMreg with all default options.

Run this code to see the output shown in the help file

    n=200;
p=3;
randn('state', 123456);
X=randn(n,p);
% Uncontaminated data
y=randn(n,1);
% Contaminated data
ycont=y;
ycont(1:5)=ycont(1:5)+6;
[out]=MMreg(ycont,X);


### MMreg with optional input arguments.

MMreg using the hyperbolic rho function.

    % Run this code to see the output shown in the help file
n=200;
p=3;
randn('state', 123456);
X=randn(n,p);
% Uncontaminated data
y=randn(n,1);
% Contaminated data
ycont=y;
ycont(1:5)=ycont(1:5)+6;
[out]=MMreg(ycont,X,'Srhofunc','optimal');


### MMreg with optional input arguments.

MMreg using the OLS estimates ac InitialEst.

    % Run this code to see the output shown in the help file
n=200;
p=3;
randn('state', 123456);
X=randn(n,p);
% Uncontaminated data
y=randn(n,1);
% Contaminated data
ycont=y;
ycont(1:5)=ycont(1:5)+6;
% OLS estimates
bols=[ones(n,1) X]\y;
res=y-[ones(n,1) X]*bols;
sols=sqrt((res'*res)/(n-p-1));
InitialEst.beta=bols;
InitialEst.scale=sols;
[out]=MMreg(ycont,X,'InitialEst',InitialEst);


### Comparing the output of different MMreg runs.

    state=100;
randn('state', state);
n=100;
X=randn(n,3);
bet=[3;4;5];
y=3*randn(n,1)+X*bet;
y(1:20)=y(1:20)+13;

%For outlier detection we consider both the nominal individual 1%
%significance level and the simultaneous Bonferroni confidence level.

% Define nominal confidence level
conflev=[0.99,1-0.01/length(y)];
% Define number of subsets
nsamp=3000;
% Define the main title of the plots
titl='';

% MM  estimators
[outMM]=MMreg(y,X,'conflev',conflev(1));
laby='Scaled MM residuals';
resindexplot(outMM.residuals,'title',titl,'laby',laby,'numlab','','conflev',conflev)
% In this example MM estimator seems to detect half of the outlier with a Bonferroni significance level.
% By simply changing the seed to 543 (state=543), using a Bonferroni size
%of 1%, no unit is declared as outlier and just half of them using the 99%
%band.

Total estimated time to complete S estimate:  3.61 seconds


### Comparison between direct call to MMreg and call to Sreg and MMregcore.

    n=30;
p=3;
randn('state', 123456);
X=randn(n,p);
% Uncontaminated data
y=randn(n,1);
% Contaminated data
ycont=y;
ycont(1:5)=ycont(1:5)+6;
% Two different rho functions are used for S and MM
rhofuncS='hyperbolic';
rhofuncMM='hampel';
% Direct call to MMreg
[out]=MMreg(ycont,X,'Srhofunc',rhofuncS,'rhofunc',rhofuncMM,'Snsamp',0);

% Call to Sreg and then to MMregcore
[outS]=Sreg(ycont,X,'rhofunc',rhofuncS,'nsamp',0);
outMM=MMregcore(ycont,X,outS.beta,outS.scale,'rhofunc',rhofuncMM);
disp('Difference between direct call to S and the calls to Sreg and MMregcore')
max(abs([out.beta-outMM.beta]))


## Input Arguments

### y — Response variable. Vector.

A vector with n elements that contains the response variable. y can be either a row or a column vector.

Data Types: single| double

### X — Data matrix of explanatory variables (also called 'regressors') of dimension (n x p-1). Rows of X represent observations, and columns represent variables.

Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) with missing or infinite values will automatically be excluded from the computations.

Data Types: single| double

### Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as  Name1,Value1,...,NameN,ValueN.

Example:  'intercept',1 , 'InitialEst',[] , 'Snsamp',1000 , 'eff',0.99 , 'effshape',1 , 'rhofunc','optimal' , 'rhofuncparam',5 , 'refsteps',10 , 'tol',1e-10 , 'conflev',0.99 , 'nocheck',1 , 'plots',0 , 'yxsave',1 , 'msg',1 

### intercept —Indicator for constant term.scalar.

If 1, a model with constant term will be fitted (default), else no constant term will be included.

Example:  'intercept',1 

Data Types: double

### InitialEst —starting values of the MM-estimator.[] (default) | structure.

InitialEst must contain the following fields

Value Description
beta

v x 1 vector (estimate of the centroid)

scale

scalar (estimate of the scale parameter).

If InitialEst is empty (default) program uses S estimators. In this last case it is possible to specify the options given in function Sreg.

Example:  'InitialEst',[] 

Data Types: struct

### Soptions —options to pass to Sreg.name value pairs.

Options if initial estimator is S and InitialEst is empty.

The options are: Srhofunc,Snsamp,Srefsteps, Sreftol, Srefstepsbestr, Sreftolbestr, Sminsctol, Sbestr.

See function Sreg for more details on these options.

It is necessary to add to the S options the letter S at the beginning. For example, if you want to use the optimal rho function the supplied option is 'Srhofunc','optimal'. For example, if you want to use 3000 subsets, the supplied option is 'Snsamp',3000

Example:  'Snsamp',1000 

Data Types: single | double MM options

### eff —nominal efficiency.scalar.

Scalar defining nominal efficiency (i.e. a number between 0.5 and 0.99). The default value is 0.95 Asymptotic nominal efficiency is:

$(\int \psi' d\Phi)^2 / (\psi^2 d\Phi)$

Example:  'eff',0.99 

Data Types: double

### effshape —location or scale efficiency.dummy scalar.

If effshape=1 efficiency refers to shape efficiency else (default) efficiency refers to location

Example:  'effshape',1 

Data Types: double

### rhofunc —rho function.string.

String which specifies the rho function which must be used to weight the residuals. Possible values are 'bisquare' 'optimal' 'hyperbolic' 'hampel' 'bisquare' uses Tukey's $\rho$ and $\psi$ functions.

See TBrho and TBpsi 'optimal' uses optimal $\rho$ and $\psi$ functions.

See OPTrho and OPTpsi 'hyperbolic' uses hyperbolic $\rho$ and $\psi$ functions.

See HYPrho and HYPpsi 'hampel' uses Hampel $\rho$ and $\psi$ functions.

See HArho and HApsi The default is bisquare

Example:  'rhofunc','optimal' 

Data Types: char

### rhofuncparam —Additional parameters for the specified rho function.scalar | vector.

For hyperbolic rho function it is possible to set up the value of k = sup CVC (the default value of k is 4.5).

For Hampel rho function it is possible to define parameters a, b and c (the default values are a=2, b=4, c=8)

Example:  'rhofuncparam',5 

Data Types: single | double

### refsteps —Maximum iterations.scalar.

Scalar defining maximum number of iterations in the MM loop. Default value is 100.

Example:  'refsteps',10 

Data Types: double

### tol —Tolerance.scalar.

Scalar controlling tolerance in the MM loop.

Default value is 1e-7

Example:  'tol',1e-10 

Data Types: double

### conflev —Confidence level which is used to declare units as outliers.scalar.

Usually conflev=0.95, 0.975 0.99 (individual alpha) or 1-0.05/n, 1-0.025/n, 1-0.01/n (simultaneous alpha).

Default value is 0.975

Example:  'conflev',0.99 

Data Types: double

### nocheck —Check input arguments.scalar.

If nocheck is equal to 1 no check is performed on matrix y and matrix X. Notice that y and X are left unchanged. In other words the additional column of ones for the intercept is not added. As default nocheck=0.

Example:  'nocheck',1 

Data Types: double

### plots —Plot on the screen.scalar | structure.

If plots = 1, generates a plot with the robust residuals against index number. The confidence level used to draw the confidence bands for the residuals is given by the input option conflev. If conflev is not specified a nominal 0.975 confidence interval will be used.

Example:  'plots',0 

Data Types: single | double

### yxsave —the response vector y and data matrix X are saved into the output structure out.scalar.

Default is 0, i.e. no saving is done.

Example:  'yxsave',1 

Data Types: double

### msg —Level of output to display.scalar.

It controls whether to display or not messages on the screen If msg==1 (default) messages are displayed on the screen about step in which signal took place else no message is displayed on the screen.

Example:  'msg',1 

Data Types: double

## Output Arguments

### out — description Structure

A structure containing the following fields:

Value Description
beta

p x 1 vector containing MM estimate of regression coefficients.

auxscale

scalar, S estimate of the scale (or supplied external estimate of scale, if option InitialEst is not empty).

residuals

n x 1 vector containing standardized MM residuals.

fittedvalues

n x 1 vector containing the fitted values.

out.residuals=(y-X*out.beta)/out.auxscale

weights

n x 1 vector. Weights assigned to each observation out.Sbeta : p x 1 vector containing S estimate of regression coefficients (or supplied initial external estimate of regression coefficients, if option InitialEst is not empty)

Ssingsub

Number of subsets without full rank in the S preliminary part. Notice that out.singsub > 0.1*(number of subsamples) produces a warning out.outliers : 1 x k vectors containing the outliers which have been found

conflev

Confidence level that was used to declare outliers

rhofunc

string identifying the rho function which has been used. If a different rho function is specified for S and MM loop then insted of out.rhofunc we will have out.rhofuncS and out.rhofuncMM.

rhofuncparam

vector which contains the additional parameters for the specified rho function which have been used. For hyperbolic rho function the value of k =sup CVC. For Hampel rho function the parameters a, b and c. If a different rho function is specified for S and MM loop then insted of out.rhofuncparam we will have out.rhofuncparamS and out.rhofuncparamMM.

y

response vector Y. The field is present if option yxsave is set to 1.

X

data matrix X. The field is present if option yxsave is set to 1.

class

'MMreg'

## References

Maronna, R.A., Martin D. and Yohai V.J. (2006), "Robust Statistics, Theory and Methods", Wiley, New York.

## Acknowledgements

This function follows the lines of MATLAB/R code developed during the years by many authors.

For more details see http://www.econ.kuleuven.be/public/NDBAE06/programs/ and the R library robustbase http://robustbase.r-forge.r-project.org/ The core of these routines, e.g. the resampling approach, however, has been completely redesigned, with considerable increase of the computational performance.