lazarus-ccr/applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas

// Use file "Sickness.laz" for testing

// to do Not sure about the "red" slope line...
// to do: There are also two code positions where a negative intercept is made positive. Why? Correct?

unit ResistanceLineUnit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
  StdCtrls, ExtCtrls, Buttons, Printers, ComCtrls, Grids,
  MainUnit, Globals, DataProcs, DictionaryUnit,
  MatrixUnit, RegressionUnit, BasicStatsReportAndChartFormUnit;

type

  { TResistanceLineForm }

  TResistanceLineForm = class(TBasicStatsReportAndChartForm)
    OptionsBevel: TBevel;
    GridChk: TCheckBox;
    PlotMediansChk: TRadioButton;
    ConfChk: TCheckBox;
    ConfEdit: TEdit;
    DescChk: TCheckBox;
    OptionsGroup: TGroupBox;
    ConfLabel: TLabel;
    LineChk: TCheckBox;
    PlotPointsChk: TRadioButton;
    VarList: TListBox;
    MeansChk: TCheckBox;
    XInBtn: TBitBtn;
    YInBtn: TBitBtn;
    XOutBtn: TBitBtn;
    YOutBtn: TBitBtn;
    XEdit: TEdit;
    YEdit: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    procedure PlotMediansChkChange(Sender: TObject);
    procedure StdCorChkChange(Sender: TObject);
    procedure VarListDblClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);
    procedure XInBtnClick(Sender: TObject);
    procedure XOutBtnClick(Sender: TObject);
    procedure YInBtnClick(Sender: TObject);
    procedure YOutBtnClick(Sender: TObject);

  private
    procedure PlotMedians(const XMedians, YMedians: DblDyneVec;
      ASlope, AIntercept: Double);
    procedure PlotXY(const XPoints, YPoints: DblDyneVec;
      ARegressionResults: TBivariateRegressionResults);

    function PrepareData(ADataGrid: TStringGrid;
      out xCol, ycol: Integer; out XData, YData: DblDyneVec;
      out ColNoSelected: IntDyneVec): Boolean;

    procedure ResidualsToGrid(xCol, yCol: Integer;
      const AColNoSelected: IntDyneVec; Slope, Intercept: Double);

    procedure ResistantLineAnalysis(const XData, YData: DblDyneVec;
      out XMedians, yMedians: DblDyneVec; out GrpSize: IntDyneVec);

    procedure WriteMedianReport(AReport: TStrings;
      const XMedians, YMedians: DblDyneVec; const GrpSize: IntDyneVec;
      Slope1, Slope2, Slope, Intercept: Double);
    procedure WriteXYReport(AReport: TStrings; AFileName: String;
      const ARegressionResults: TBivariateRegressionResults);

  protected
    procedure AdjustConstraints; override;
    procedure Compute; override;		
    procedure UpdateBtnStates; override;

  public
    constructor Create(AOwner: TComponent); override;
    procedure Reset; override;
  end; 

var
  ResistanceLineForm: TResistanceLineForm;

implementation

{$R *.lfm}

uses
  Math,
  TATypes, TAChartUtils, TAChartAxisUtils, TALegend, TASources, TACustomSeries,
  ChartFrameUnit, Utils, GridProcs;

{ TResistanceLineForm }

constructor TResistanceLineForm.Create(AOwner: TComponent);
begin
  inherited;
  if DictionaryFrm = nil then Application.CreateForm(TDictionaryFrm, DictionaryFrm);

  FChartFrame.Chart.Legend.Alignment := laBottomCenter;
  FChartFrame.Chart.Legend.ColumnCount := 3;
  FChartFrame.Chart.Legend.TextFormat := tfHTML;
  with FChartFrame.Chart.AxisList.Add do
  begin
    Alignment := calRight;
    Marks.Source := TListChartSource.Create(self);
    Marks.Style := smsLabel;
    Grid.Visible := false;
    TickColor := clNone;
  end;
  with FChartFrame.Chart.AxisList.Add do
  begin
    Alignment := calTop;
    Marks.Source := TListChartSource.Create(self);
    Marks.Style := smsLabel;
    Grid.Visible := false;
    TickColor := clNone;
  end;

  PageControl.ActivePage := ChartPage;
end;


procedure TResistanceLineForm.AdjustConstraints;
begin
  ParamsPanel.Constraints.MinHeight :=  YOutBtn.Top + YOutBtn.Height +
    VarList.BorderSpacing.Bottom + OptionsGroup.Height + ButtonBevel.Height +
    CloseBtn.Height;

  ParamsPanel.Constraints.MinWidth := Max(
    4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left,
    OptionsGroup.Width);
end;


procedure TResistanceLineForm.Compute;
var
  xValues: DblDyneVec = nil;
  yValues: DblDyneVec = nil;
  xMedians: DblDyneVec = nil;
  yMedians: DblDyneVec = nil;
  colNoSelected : IntDyneVec = nil;
  grpSize: IntDyneVec = nil;
  xCol, yCol: Integer;
  regressionRes: TBivariateRegressionResults;
  confBand: Double;
  slope1, slope2, slopeRL: Double;
  c1, c2, c3, c: Double;
  lReport: TStrings;
  dataGrid: TStringGrid;
begin
  dataGrid := OS3MainFrm.DataGrid;

  // Get values from the grid
  if not PrepareData(dataGrid, xCol, yCol, xValues, yvalues, colNoSelected) then
    exit;

  // Sort on the x values
  SortOnX(xValues, yValues);

  // Calculate bivariate regression
  confBand := StrToFloat(ConfEdit.Text) / 100;
  BivariateRegression(xValues, yValues, confBand, regressionRes);

  // Do the resistant line analysis
  ResistantLineAnalysis(xValues, yValues, xMedians, yMedians, grpSize);

  slope1 := (yMedians[1] - yMedians[0]) / (xMedians[1] - xMedians[0]);
  slope2 := (yMedians[2] - yMedians[1]) / (xMedians[2] - xMedians[1]);
  slopeRL  := (yMedians[2] - yMedians[0]) / (xMedians[2] - xMedians[0]);

  // obtain estimate of the constant for the prediction equation
  c1 := slopeRL * xMedians[0] - yMedians[0];
  c2 := slopeRL * xMedians[1] - yMedians[1];
  c3 := slopeRL * xMedians[2] - yMedians[2];
  c := (c1 + c2 + c3) / 3.0;
  if c < 0 then c := -c;         // wp: why is this necessary?

  if GridChk.Checked then
    ResidualsToGrid(xCol, yCol, colNoSelected, slopeRL, c);

  // Write results to report in ReportFrame
  lReport := TStringList.Create;
  try
    WriteMedianReport(lReport,
      xMedians, yMedians, grpSize,
      slope1, slope2, slopeRL, c
    );

    if DescChk.Checked then
       WriteXYReport(lReport,
         OS3MainFrm.FileNameEdit.Text,
         regressionRes
       );

    FReportFrame.DisplayReport(lReport);
  finally
    lReport.Free;
  end;

  // Plot results in ChartFrame
  if PlotPointsChk.Checked then
    PlotXY(xValues, yValues, regressionRes);

  if PlotMediansChk.Checked then
    PlotMedians(xMedians, yMedians, slopeRL, c);
end;


procedure TResistanceLineForm.PlotMedians(const XMedians, YMedians: DblDyneVec;
  Aslope, AIntercept: Double);
var
  ser: TChartSeries;
  rightLabels: TListChartSource;
  topLabels: TListChartSource;
  labels: StrDyneVec = nil;
  xpts: DblDyneVec = nil;
  ypts: DblDyneVec = nil;
  i: Integer;
begin
  rightLabels := FChartFrame.Chart.AxisList[2].Marks.Source as TListChartSource;
  rightLabels.Clear;
  topLabels := FChartFrame.Chart.AxisList[3].Marks.Source as TListChartSource;
  topLabels.Clear;
  FChartFrame.Clear;
  FChartFrame.Clear;

  FChartFrame.Chart.Margins.Left := 20;
  FChartFrame.Chart.Margins.Right := 20;

  FChartFrame.SetTitle('Median plot for three groups');
  FChartFrame.SetFooter(Format('Slope = %.3f', [ASlope]));
  FChartFrame.SetXTitle(XEdit.Text);
  FChartFrame.SetYTitle(YEdit.Text);

  SetLength(labels, Length(XMedians));
  for i := 0 to High(XMedians) do
    labels[i] := Format('M%d', [i+1]);
  ser := FChartFrame.PlotXY(ptLinesAndSymbols, XMedians, YMedians, labels, nil, 'Medians', clNavy, psCircle);
  ser.Marks.Style := smsLabel;
  labels := nil;

  // Draw slope line
  SetLength(xpts, 2);
  SetLength(ypts, 2);
  xpts[0] := XMedians[0];
  ypts[0] := AIntercept + ASlope * xpts[0];
  xpts[1] := XMedians[2];
  ypts[1] := AIntercept + ASlope * xpts[1];
  FChartFrame.PlotXY(ptLines, xpts, ypts, nil, nil, 'Slope line', clRed);
end;


procedure TResistanceLineForm.PlotXY(const XPoints, YPoints: DblDyneVec;
  ARegressionResults: TBivariateRegressionResults);
var
  ser: TChartSeries;
  xPts: DblDyneVec = nil;
  yPts: DblDyneVec = nil;
  conf: DblDyneVec = nil;
  rightLabels: TListChartSource;
  topLabels: TListChartSource;
  i: Integer;
begin
  rightLabels := FChartFrame.Chart.AxisList[2].Marks.Source as TListChartSource;
  rightLabels.Clear;
  topLabels := FChartFrame.Chart.AxisList[3].Marks.Source as TListChartSource;
  topLabels.Clear;
  FChartFrame.Clear;

  FChartFrame.Chart.Margins.Left := 4;
  FChartFrame.Chart.Margins.Right := 4;

  FChartFrame.SetTitle('X vs Y plot using file '+ OS3MainFrm.FileNameEdit.Text);
  with ARegressionResults do
    FChartFrame.SetFooter(Format('R(X,Y) = %.3f, Slope = %.2f, Intercept = %.2f', [R, Slope, Intercept]));
  FChartFrame.SetXTitle(XEdit.Text);
  FChartFrame.SetYTitle(YEdit.Text);

  // Draw means
  if MeansChk.Checked then
    with ARegressionResults do
    begin
      FChartFrame.VertLine(XMean, clGreen, psDashDot, 'Mean ' + XEdit.Text);
      topLabels.Add(XMean, XMean, 'Mean ' + XEdit.Text);

      FChartFrame.HorLine(YMean, clGreen, psDash, 'Mean ' + YEdit.Text);
      rightLabels.Add(YMean, YMean, 'Mean ' + YEdit.Text);
    end;

  // Draw slope line
  if LineChk.Checked then
  begin
    SetLength(xPts, 2);
    SetLength(yPts, 2);
    xPts[0] := XPoints[0];
    xPts[1] := XPoints[High(XPoints)];
    with ARegressionResults do
    begin
      yPts[0] := Intercept + XPoints[0] * Slope;
      yPts[1] := Intercept + XPoints[High(XPoints)] * Slope;
    end;
    ser := FChartFrame.PlotXY(ptLines, xpts, ypts, nil, nil, 'Regression', clBlack);
    rightLabels.Add(ser.XValue[ser.Count-1], ser.YValue[ser.Count-1], 'Regression');
  end;

  if ConfChk.Checked then
  begin
    // Draw upper confidence band
    SetLength(conf, ARegressionResults.Count);
    for i := 0 to ARegressionResults.Count-1 do
      conf[i] := ARegressionResults.ConfidenceLimits(XPoints[i], true);
    ser := FChartFrame.PlotXY(ptLines, XPoints, conf, nil, nil, 'Upper confidence band', clRed);
    rightLabels.Add(ser.yValue[ser.Count-1], ser.YValue[ser.Count-1], 'UCL');

    // Draw lower confidence band
    for i := 0 to ARegressionResults.Count-1 do
      conf[i] := ARegressionResults.ConfidenceLimits(XPoints[i], false);
    ser := FChartFrame.PlotXY(ptLines, XPoints, conf, nil, nil, 'Lower confidence band', clRed);
    rightLabels.Add(ser.yValue[ser.Count-1], ser.YValue[ser.Count-1], 'LCL');
  end;

  // Draw points for x and y pairs
  ser := FChartFrame.PlotXY(ptSymbols, XPoints, YPoints, nil, nil, 'Data', clYellow, psCircle);
  ser.Marks.Style := smsNone;
end;


function TResistanceLineForm.PrepareData(ADataGrid: TStringGrid;
  out xCol, ycol: Integer; out XData, YData: DblDyneVec;
  out ColNoSelected: IntDyneVec): Boolean;
var
  N: Integer;
begin
  Result := false;

  ColNoSelected := nil;
  XData := nil;
  YData := nil;

  xCol := GetVariableIndex(ADataGrid, XEdit.Text);
  yCol := GetVariableIndex(ADataGrid, YEdit.Text);
  if xCol = -1 then
  begin
    ErrorMsg('X variable not specified or not found.');
    exit;
  end;

  if yCol = -1 then
  begin
    ErrorMsg('Y variable not specified or not found.');
    exit;
  end;

  SetLength(ColNoSelected, 2);
  ColNoSelected[0] := xCol;
  ColNoSelected[1] := yCol;

  XData := CollectVecValues(ADataGrid, xCol, colNoSelected);
  YData := CollectVecValues(ADataGrid, ycol, colNoSelected);
  N := Length(XData);
  if N < 3 then
  begin
    ErrorMsg('At least three data points required.');
    exit;
  end;
  if N <> Length(YData) then
  begin
    ErrorMsg('Equal count of cases required for x and y variables.');
    exit;
  end;

  Result := true;
end;


procedure TResistanceLineForm.Reset;
var
  i: integer;
begin
  inherited;

  PlotMediansChk.Checked := true;
  PlotMediansChkChange(PlotMediansChk);
  GridChk.Checked := false;

  XEdit.Text := '';
  YEdit.Text := '';
  VarList.Clear;
  ConfEdit.Text := FormatFloat('0.0', DEFAULT_CONFIDENCE_LEVEL_PERCENT);
  for i := 1 to NoVariables do
    VarList.Items.Add(OS3MainFrm.DataGrid.Cells[i,0]);
  UpdateBtnStates;
end;


{ Get the residuals (Y - predicted Y) for each X value and place them in the grid }
procedure TResistanceLineForm.ResidualsToGrid(xCol, yCol: Integer;
  const AColNoSelected: IntDyneVec; Slope, Intercept: Double);
var
  s: String;
  i: Integer;
  x, y, yPred, yRes: Double;
begin
  // New column for predicted value
  s := 'Pred.' + OS3MainFrm.DataGrid.Cells[yCol, 0];
  DictionaryFrm.NewVar(NoVariables+1);
  DictionaryFrm.DictGrid.Cells[1, NoVariables] := s;
  OS3MainFrm.DataGrid.Cells[NoVariables, 0] := s;

  // New column for residual value
  s := 'Residual';
  DictionaryFrm.NewVar(NoVariables+1);
  DictionaryFrm.DictGrid.Cells[1, NoVariables] := s;
  OS3MainFrm.DataGrid.Cells[NoVariables, 0] := s;

  for i := 1 to NoCases do
  begin
    if not DataProcs.GoodRecord(i, Length(AColNoSelected), AColNoSelected) then continue;
    x := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[xCol, i]));
    y := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[yCol, i]));

    // Calculate predicted value
    if Intercept >= 0 then
      yPred := Slope * x + Intercept        // wp: why two equations?
    else
      yPred := Slope * x - Intercept;

    // Calculate residual
    yRes := y - yPred;

    // Write to grid
    OS3MainFrm.DataGrid.Cells[NoVariables-1, i] := Format('%.3f', [yPred]);
    OS3MainFrm.DataGrid.Cells[NoVariables, i]  := Format('%.3f', [yRes]);
  end;
end;

{ Do the resistant line analysis }
procedure TResistanceLineForm.ResistantLineAnalysis(const XData, YData: DblDyneVec;
  out XMedians, yMedians: DblDyneVec; out GrpSize: IntDyneVec);
var
  i, N, offs: Integer;
  xVector: DblDyneVec = nil;
  yVector: DblDyneVec = nil;
begin
  N := Length(XData);

  GrpSize := nil;
  XMedians := nil;
  YMedians := nil;

  SetLength(GrpSize, 3);
  SetLength(XMedians, 3);
  SetLength(YMedians, 3);

  GrpSize[0] := N div 3;
  GrpSize[2] := GrpSize[0];
  GrpSize[1] := N - GrpSize[0] - GrpSize[2];

  // Sort the y values. The x values already are sorted.
  //SortOnX(xyPoints[1]);

  // Get median for each group of x and y values

  // First group
  SetLength(xVector, GrpSize[0]);
  SetLength(yVector, GrpSize[0]);
  offs := 0;
  for i := 0 to GrpSize[1] - 1 do
  begin
    xVector[i] := XData[i + offs];
    yVector[i] := YData[i + offs];
  end;
  xMedians[0] := VecMedian(xVector);
  yMedians[0] := VecMedian(yVector);

  // Second group
  SetLength(xVector, GrpSize[1]);
  SetLength(yVector, GrpSize[1]);
  offs := GrpSize[0];
  for i := 0 to GrpSize[1] - 1 do
  begin
    xVector[i] := XData[i + offs];
    yVector[i] := YData[i + offs];
  end;
  xMedians[1] := VecMedian(xVector);
  yMedians[1] := VecMedian(yVector);

  // Third group
  SetLength(xVector, GrpSize[2]);
  SetLength(yVector, GrpSize[2]);
  offs := GrpSize[0] + GrpSize[1];
  for i := 0 to GrpSize[2] - 1 do
  begin
    xVector[i] := XData[i + offs];
    yVector[i] := YData[i + offs];
  end;
  xMedians[2] := VecMedian(xVector);
  yMedians[2] := VecMedian(yVector);
end;


procedure TResistanceLineForm.StdCorChkChange(Sender: TObject);
begin
  //if StdCorChk.Checked then OptionsGroup.Enabled := true else OptionsGroup.Enabled := false;
end;

procedure TResistanceLineForm.PlotMediansChkChange(Sender: TObject);
var
  enable: Boolean;
begin
  enable := ((Sender = PlotPointsChk) and PlotPointsChk.Checked) or
            ((Sender = PlotMediansChk) and (not PlotMediansChk.Checked));
  LineChk.Enabled := enable;
  MeansChk.Enabled := enable;
  ConfChk.Enabled := enable;
  ConfLabel.Enabled := enable;
  ConfEdit.Enabled := enable;
  GridChk.Enabled := enable;
  DescChk.Enabled := enable;
end;


procedure TResistanceLineForm.VarListDblClick(Sender: TObject);
var
  index: Integer;
begin
  index := VarList.ItemIndex;
  if index > -1 then
  begin
    if XEdit.Text = '' then
      XEdit.Text := VarList.Items[index]
    else if YEdit.Text = '' then
      YEdit.Text := VarList.Items[index];
    VarList.Items.Delete(index);
    UpdateBtnStates;
  end;
end;


procedure TResistanceLineForm.UpdateBtnStates;
begin
  inherited;
	
  XInBtn.Enabled := (VarList.ItemIndex > -1) and (XEdit.Text = '');
  YInBtn.Enabled := (VarList.ItemIndex > -1) and (YEdit.Text = '');
  XOutBtn.Enabled := (XEdit.Text <> '');
  YOutBtn.Enabled := (YEdit.Text <> '');
end;


procedure TResistanceLineForm.VarListSelectionChange(Sender: TObject;
  User: boolean);
begin
  UpdateBtnStates;
end;


procedure TResistanceLineForm.WriteMedianReport(AReport: TStrings;
  const XMedians, YMedians: DblDyneVec; const GrpSize: IntDyneVec;
  Slope1, Slope2, Slope, Intercept: Double);
var
  i: Integer;
begin
  AReport.Add('RESISTANT LINE ANALYSIS');
  AReport.Add('');
  AReport.Add('Group   X Median   Y Median   Size');
  AReport.Add('-----   --------   --------   -----');
  for i := 0 to 2 do
    AReport.Add('%3d     %8.3f   %8.3f   %3d', [i+1, XMedians[i], YMedians[i], GrpSize[i]]);
  AReport.Add('');

  AReport.Add('Half Slopes:          %10.3f and %.3f', [Slope1, Slope2]);
  AReport.Add('Slope:                %10.3f', [Slope]);
  aReport.Add('Ratio of half slopes: %10.3f', [Slope2 / Slope1]);

  // Prediction equation
  AReport.Add('Equation:  y = %.3f * X + %.3f', [Slope, Intercept]);
end;


procedure TResistanceLineForm.WriteXYReport(AReport: TStrings; AFileName: String;
  const ARegressionResults: TBivariateRegressionResults);
{
  xMean, yMean, xVariance, yVariance, xStdDev, yStdDev, R, slope, intercept, SEPred: double;
  N: Integer);
  }
begin
  AReport.Add('');
  AReport.Add(DIVIDER);
  AReport.Add('');
  AReport.Add('ORIGINAL X vs. Y PLOT');
  AReport.Add('');
  AReport.Add('Data file: ' + AFileName);
  AReport.Add('Variables:');
  AReport.Add('  X: ' + XEdit.Text);
  AReport.Add('  Y: ' + YEdit.Text);
  AReport.Add('');
  AReport.Add('Variable      Mean    Variance  Std.Dev.');
  AReport.Add('----------  --------  --------  --------');
  with ARegressionResults do
  begin
    AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [XEdit.Text, XMean, XVariance, XStdDev]);
    AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [YEdit.Text, YMean, YVariance, YStdDev]);
    AReport.Add('');
    AReport.Add('Correlation:                %8.4f', [R]);
    AReport.Add('Slope:                      %8.3f', [Slope]);
    AReport.Add('Intercept:                  %8.3f', [Intercept]);
    AReport.Add('Standard Error of Estimate: %8.3f', [StdErrorPredicted]);
    AReport.Add('Number of good cases:       %8d',   [Count]);
  end;

  {
  AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [XEdit.Text, xMean, xVariance, xStdDev]);
  AReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [YEdit.Text, yMean, yVariance, yStdDev]);
  AReport.Add('');
  AReport.Add('Correlation:                %8.4f', [R]);
  AReport.Add('Slope:                      %8.3f', [Slope]);
  AReport.Add('Intercept:                  %8.3f', [Intercept]);
  AReport.Add('Standard Error of Estimate: %8.3f', [SEPred]);
  AReport.Add('Number of good cases:       %8d', [N]);
  }
end;


procedure TResistanceLineForm.XInBtnClick(Sender: TObject);
var
  index: integer;
begin
  index := VarList.ItemIndex;
  if (index > -1) and (XEdit.Text = '') then
  begin
    XEdit.Text := VarList.items[index];
    VarList.Items.Delete(index);
  end;
  UpdateBtnStates;
end;


procedure TResistanceLineForm.XOutBtnClick(Sender: TObject);
begin
  if XEdit.Text <> '' then
  begin
    VarList.Items.Add(XEdit.Text);
    XEdit.Text := '';
  end;
  UpdateBtnStates;
end;


procedure TResistanceLineForm.YInBtnClick(Sender: TObject);
var
  index: integer;
begin
  index := VarList.ItemIndex;
  if (index > -1) and (YEdit.Text = '') then
  begin
    YEdit.Text := VarList.items[index];
    VarList.Items.Delete(index);
  end;
  UpdateBtnStates;
end;


procedure TResistanceLineForm.YOutBtnClick(Sender: TObject);
begin
  if YEdit.Text <> '' then
  begin
    VarList.Items.Add(YEdit.Text);
    YEdit.Text := '';
  end;
  UpdateBtnStates;
end;


end.