(*:Mathematica Version: 3.0 *)

(*:Package Version: 1.0 *)

(*:Name: RotationLine *)

(*:Author: Karin Isler *)

(*:Keywords: Allometry, Major Axis, Robust Regression, Line-fitting by Rotation *)

(*:Requirements: None. *)

(*:Warnings: None*)

(*:Sources: K. Isler, A.D. Barbour, R. D. Martin: Line-fitting by rotation: A nonparametric method for bivariate allometric analysis. Biometrical Journal 44: 289-304.*)

(*:Summary: This package provides methods for line-fitting in bivariate allometric analyses: least-squares regression, major axis regression and a new nonparametric method, line-fitting by rotation. *)


BeginPackage["RotationLine`"]



ReadFile::usage = "ReadFile[file_,form_,heading_] reads data from a file. The filename must be given in quotation marks. In form, the type of data in the columns must be specified, e.g. {Word,Number,Number}. The number of heading rows that contain no data can also be give."

Data::usage = "Data[list_,name_,xvar_,yvar_] yields a (nx3) data matrix from a list. name, xvar and yvar specify the columns that are to be included. Xvar is the independent variable and yvar the dependent variable."

LinePlot::usage = "LinePlot[data_] yields a plot of the data with the least-squraes regression (dotted line) and the major axis (broken line). Optionally, in LinePlot[data_,{x_,y_}] the axes can be labeled in quotation marks."

RotationLinePlot::usage = "RotationLinePlot[data_,theta_] yields a lot of the data with the least-squares regression (dotted line), the major axis (broken line) and the line of the rotation method (solid line) with the rotation angle theta. Theta is required in radians. Arctan(theta) is the slope of the rotation line. Optionally, in RotationLinePlot[data_,{x_,y_}] the axes can be labeled in quotation marks."

RotationValue::usage = "RotationValue[data_] gives the value of the measure of dependence D, i.e. the sum of the absolute deviations of the empirical distribution of the data points from the distribution that would be expected if the margin distributions were independent, in the quantile squares."

RotationGraph::usage = "RotationGraph[data_,{theta1_:0,theta2_:Pi/2}] yields a graph of the measure of dependence D(theta) of the data which are rotated from theta1 to theta2."

RotationMinimum::usage = "RotationMinimum[data_] yields a list of ten minima of the measure of dependence D(theta), which are found with ten initial values of theta between 0 and Pi/2. "


Begin["RotationLine`"]


(* Read Data *)

ReadFile[file_,form_,heading_] := Module[{tfile,data},
If[heading =!= 0,
(tfile = OpenRead[file];
Do[Skip[tfile,Record],{i,1,heading}];
data = ReadList[tfile,form];
Close[tfile]),
data = ReadList[file,form]];
Return[data]]

Data[list_,name_,xvar_,yvar_] :=Map[{#[[name]], #[[xvar]],#[[yvar]]}&,list]


(* Descriptive Statistics *)

Mean[list_] := Apply[Plus, list] / Length[list]

Variance[list_] :=
Apply[Plus, (list - Mean[list])^2]/(Length[list] - 1)

Covariance[list1_,list2_] := Apply[Plus,(list1-Mean[list1])*(list2 - Mean[list2])] / (Length[list1]-1)

Correlation[list1_,list2_] := Covariance[list1,list2]/Sqrt[Variance[list1]*Variance[list2]]

RSquared[list1_,list2_] := Correlation[list1,list2]^2

RSquared[list_] := Correlation[Transpose[list][[2]],Transpose[list][[3]]]^2

RSquaredMax[data_] :=
Module[{datax,datay,slope,datanew,v1,v2},
datax = Transpose[data][[2]];
datay = Transpose[data][[3]];
v1 = Variance[datay] - Variance[datax];
v2 = Covariance[datax,datay];
slope = (v1 + Sqrt[v1^2+ 4* v2^2]) /(2*v2);
slope = ArcTan[slope];
datanew = Rot2D[data,N[slope-Pi/4]];
Return[RSquared[datanew]]]


(* Major Axis and LS-Regression *)

Regression[data_] :=
Module[{slope,intercept,datax,datay},
datax = Transpose[data][[2]];
datay = Transpose[data][[3]];
slope = Covariance[datax,datay]/Variance[datax];
intercept = Mean[datay] - slope * Mean[datax];
Print["LS-Regression: ",N[intercept]+N[slope]*"x"," (Angle: ", N[ArcTan[slope]], " rad)"];
Return[intercept + slope * x]]

MajorAxis[data_] :=
Module[{datax,datay,slope,intercept,v1,v2},
datax = Transpose[data][[2]];
datay = Transpose[data][[3]];
v1 = Variance[datay] - Variance[datax];
v2 = Covariance[datax,datay];
slope = (v1 + Sqrt[v1^2+ 4* v2^2]) /(2*v2);
intercept = (Mean[datay] - slope * Mean[datax]);
Print["Major Axis: ", N[intercept] + N[slope] * "x"," (Angle: ", N[ArcTan[slope]], " rad)"];
Return[intercept + slope * x]]


(* Line-fitting by Rotation *)


Rot2D[data_,alpha_:0] := Map[{#[[1]],#[[2]]* Cos[alpha] + #[[3]] * Sin[alpha],- #[[2]]* Sin[alpha] + #[[3]] * Cos[alpha]}&, data]

logfct[{x_,y_},{a1_,a2_},{b1_,b2_}] :=
((1-1/(1+Exp[(20/(a2-a1))*(-x+a2)]))*(1-1/(1+Exp[(20/(b2-b1))*(-y+b2)])))/((1+Exp[(20/(a2-a1))*(-x+a1)])*(1+Exp[(20/(b2-b1))*(-y+b1)]))

Quantile[list_, q_] :=
Part[Sort[list],-Floor[-q*Length[list]]] /;
VectorQ[list] && 0 < q < 1 && Length[list] > 2

qua3[t_] := {Min[t],Quantile[t,0.333],Quantile[t,0.666],Max[t]}

qua4[t_] := {Min[t], Quantile[t,0.25], Quantile[t,0.5],Quantile[t,0.75], Max[t]}

qua5[t_] := {Min[t],Quantile[t,0.2],Quantile[t,0.4],Quantile[t,0.6],Quantile[t,0.8],Max[t]}

RotationValue3[data_] := Module[{quax,quay,data1},
quax = qua3[Transpose[data][[2]]];
quay = qua3[Transpose[data][[3]]];
data1 = Map[{#[[2]],#[[3]]}&,data];
Apply[Plus,Abs[Table[Table[N[(Apply[Plus,Map[logfct[#,{quax[[i]],quax[[i+1]]}, {quay[[j]],quay[[j+1]]}]&,data1]]/ Length[data])-0.111111],{i,3}],{j,3}]],{0,1}]]

RotationMinimum3[data_,theta1_] :=
FindMinimum[RotationValue3[Rot2D[data,theta]],{theta,theta1}]

RotationValue4[data_] := Module[{quax,quay,data1},
quax = qua4[Transpose[data][[2]]];
quay = qua4[Transpose[data][[3]]];
data1 = Map[{#[[2]],#[[3]]}&,data];
Apply[Plus,Abs[Table[Table[N[(Apply[Plus,Map[logfct[#,{quax[[i]],quax[[i+1]]}, {quay[[j]],quay[[j+1]]}]&,data1]]/ Length[data])-0.0625],{i,4}],{j,4}]],{0,1}]]

RotationMinimum4[data_,theta1_] :=
FindMinimum[RotationValue4[Rot2D[data,theta]],{theta,theta1}]

RotationValue5[data_] := Module[{quax,quay,data1},
quax = qua5[Transpose[data][[2]]];
quay = qua5[Transpose[data][[3]]];
data1 = Map[{#[[2]],#[[3]]}&,data];
Apply[Plus,Abs[Table[Table[N[(Apply[Plus,Map[logfct[#,{quax[[i]],quax[[i+1]]}, {quay[[j]],quay[[j+1]]}]&,data1]]/ Length[data])-0.04],{i,5}],{j,5}]],{0,1}]]

RotationMinimum5[data_,theta1_] :=
FindMinimum[RotationValue5[Rot2D[data,theta]],{theta,theta1}]

RotationValue[data_] :=
If[Length[data]>1024,RotationValue5[data],If[Length[data]>243,RotationValue4[data],RotationValue3[data]]]

RotationMinimum[data_] :=
If[Length[data]>1024,Sort[Table[RotationMinimum5[data,{theta,theta+0.12}],{theta,0.1,Pi/2,0.15}]], If[Length[data]>243,Sort[Table[RotationMinimum4[data,{theta,theta+0.12}],{theta,0.1,Pi/2,0.15}]], Sort[Table[RotationMinimum3[data,{theta,theta+0.12}],{theta,0.1,Pi/2,0.15}]]]]

RotationGraph[data_,{theta1_:0,theta2_:Pi/2}] :=
Plot[RotationValue[Rot2D[data,theta]],{theta,theta1,theta2}]

RotationLine[data_,theta_] := Module[{med,rotmed,rotationline},
med = {Quantile[Transpose[Rot2D[data,theta]][[2]],0.5],
Quantile[Transpose[Rot2D[data,theta]][[3]],0.5]};
rotmed = med. {{Cos[-theta], -Sin[-theta]},{Sin[-theta],Cos[-theta]}};
rotationline = (rotmed[[2]] - Tan[theta]* rotmed[[1]])+ Tan[theta]*x;
Print["RotationLine: ",N[(rotmed[[2]] - Tan[theta]* rotmed[[1]])] +N[Tan[theta]] * "x"," (Angle: ", N[theta], " rad)"];
Return[rotationline]]


(* Plots *)


LinePlot[data_,axeslabel_:None] :=
Module[{datax,datay,line1,line2,plot1,plot2},
datax = Transpose[data][[2]];
datay = Transpose[data][[3]];
Print["RSquared: ",RSquared[datax,datay]];
Print["RSquaredMax: ",RSquaredMax[data]];
line1 = Regression[data];
line2 = MajorAxis[data];
plot1 = ListPlot[Transpose[{datax,datay}],
DisplayFunction->Identity];
plot2 = Plot[{line1,line2},{x,Min[datax],Max[datax]}, PlotStyle->{{AbsoluteThickness[1], AbsoluteDashing[{1,2}]},{AbsoluteThickness[1],AbsoluteDashing[{8,4}]}}, DisplayFunction->Identity];
Show[plot1,plot2,AxesLabel->axeslabel,AxesOrigin->{Min[datax],Min[datay]},PlotRange->All, DisplayFunction->$DisplayFunction]]

RotationLinePlot[data_,theta_,axeslabel_:None] :=
Module[{datax,datay,line1,line2,rotationline,plot1,plot2},
datax = Transpose[data][[2]];
datay = Transpose[data][[3]];
Print["RSquared: ",RSquared[datax,datay]];
Print["RSquaredMax: ",RSquaredMax[data]];
line1 = Regression[data];
line2 = MajorAxis[data];
rotationline = RotationLine[data,theta];
plot1 = ListPlot[Transpose[{datax,datay}],DisplayFunction->Identity];
plot2 = Plot[{line1,line2,rotationline},{x,Min[datax],Max[datax]},PlotStyle->{{AbsoluteThickness[1], AbsoluteDashing[{1,2}]},{AbsoluteThickness[1],AbsoluteDashing[{8,4}]},{AbsoluteThickness[1]}}, DisplayFunction->Identity];
Show[plot1,plot2,AxesLabel->axeslabel,AxesOrigin->{Min[datax],Min[datay]},PlotRange->All, DisplayFunction->$DisplayFunction]]