Code covered by the BSD License

# maregress

### Antonio Trujillo-Ortiz (view profile)

16 Jun 2010 (Updated )

Major Axis Regression (Principal Axis Regression).

maregress(x,y,alpha)
```function [b,bint,l,ang] = maregress(x,y,alpha)
%MAREGRESS Major Axis Regression (Principal Axis Regression).
%   Model II regression should be used when the two variables in the
%   regression equation are random and subject to error, i.e. not
%   controlled by the researcher. Model I regression using ordinary least
%   squares underestimates the slope of the linear relationship between the
%   variables when they both contain error. According to Sokal and Rohlf
%   (1995), the subject of Model II regression is one on which research and
%   controversy are continuing and definitive recommendations are difficult
%   to make. In Sokal and Rohlf (1981, 2nd ed.), the numerical result for
%   major axis regression for the example data set is wrong. The mistake
%   has been corrected in the 1995 edition.
%
%   MAREGRESS is a Model II procedure. When both variables are in the same
%   units of measurement the slope of the major axis (principal axis) of
%   the bivariate sample has been sugested. This method requires a
%   knowledge of correlation techniques. A bivariate normal distribution
%   can be represented by means of concentric ellipses describing the
%   topography of a bell-shaped mound, and an ellipse can be described by
%   two principal axis (major and minor) at rigth angles toeach other.
%   Therefore the main task is to find the slope and equation of the major
%   axis from a sample.
%
%   [B,BINT,L,ANG] = MAREGRESS(X,Y,ALPHA) returns the vector B of
%   regression coefficients in the linear Model II, a matrix BINT of the
%   given confidence intervals for B, the L eigenvalues of the axes, and
%   the angle in degrees of the confidence bounds.
%
%   MAREGRESS treats NaNs in X or Y as missing values, and removes them.
%
%   Syntax: [b,bint,l,ang] = maregress(x,y,alpha)
%
%   Example. From the Box 14.12 (California fish cabezon [Scorpaenichthys
%   marmoratus]) of Sokal and Rohlf (1995). The data are:
%
%   x=[14,17,24,25,27,33,34,37,40,41,42];
%   y=[61,37,65,69,54,93,87,89,100,90,97];
%
%   Calling on Matlab the function:
%                [b,bint,l,ang] = maregress(x,y)
%
%
%   b =
%      6.6566    2.3017
%
%   bint =
%    -36.5408   27.8330
%      1.6043    3.7244
%
%   l =
%    494.6340   17.4933
%
%   ang =
%     66.5172
%
%   Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.edu.mx
%
%   Copyright (C)  June 8, 2010.
%
%   To cite this file, this would be an appropriate format:
%   Trujillo-Ortiz, A. and R. Hernandez-Walls. (2010). maregress:
%      Major Axis Regression. A MATLAB file. [WWW document]. URL
%      http://www.mathworks.com/matlabcentral/fileexchange/27916-maregress
%
%   References:
%   Sokal, R. R. and Rohlf, F. J. (1995), Biometry. The principles and
%              practice of the statistics in biologicalreserach. 3rd. ed.
%              New-York:W.H.,Freeman. [Sections 14.13 and 15.7]
%

if  nargin < 2
error('maregress:TooFewInputs', ...
'MAREGRESS requires at least two input arguments.');
elseif nargin == 2
alpha = 0.05;
end

x = x(:); y = y(:);

% Check that matrix (X) and rigth hand side (Y) have compatible dimensions
[n,ncolx] = size(x);
if ~isvector(y)
error('maregress:InvalidData', 'Y must be a vector.');
elseif numel(y) ~= n
error('maregress:InvalidData', ...
'The number of rows in Y must equal the number of rows in X.');
end

% Remove missing values, if any
wasnan = (isnan(y) | any(isnan(x),2));
havenans = any(wasnan);
if havenans
y(wasnan) = [];
x(wasnan,:) = [];
n = length(y);
end

S = cov(x,y);
b1= ((S(2,2)-S(1,1))+sqrt(((S(2,2)-S(1,1))^2)+4*(S(1,2)^2)))/(2*S(1,2)); %slope
a = mean(y)-b1*mean(x); %intercept
b = [a,b1];
D = sqrt((S(1,1)+S(2,2))^2-4*((S(1,1)*S(2,2))-S(1,2)^2));
l1 = (S(1,1)+S(2,2)+D)/2; %eigenvalue of major (principal) axis
l2 = (S(1,1)+S(2,2)-D)/2; %eigenvalue of minor axis
H = (finv(1-alpha,1,n-2))/(((l1/l2)+(l2/l1)-2)*(n-2));
A = sqrt(H/(1-H));
Li = (b1-A)/(1+b1*A); %confidence lower limit of slope
Ls = (b1+A)/(1-b1*A); %confidence upper limit of slope
l = [l1,l2];
ai = mean(y)-Ls*mean(x); %confidence lower limit of intercept
as = mean(y)-Li*mean(x); %confidence upper limit of intercept
bint = [ai,as;Li,Ls];
ang = atand(b1); %angle in degrees

return,```