// Test file: srh.tab
// Created with data from https://rcompanion.org/handbook/F_14.html
// - Dependent variable: Yield
// - Factor 1: Crop
// - Factor 2: Fertilizer

unit SRHTestUnit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
  StdCtrls, Buttons, ExtCtrls, ComCtrls,
  TASources, TACustomSeries, TAStyles,
  MainUnit, Globals, FunctionsLib, BasicStatsReportAndChartFormUnit;

type

  TSRHResults = record
    TotalCount: Integer;  // N

    CellCounts: IntDyneMat;
    CellSums: DblDyneMat;
    CellVars: DblDyneMat;

    RowSums, ColSums: DblDyneVec;
    RowCount, ColCount: IntDyneVec;

    MeanDep, MeanF1, MeanF2: double;

    SSDep, SSErr, SSF1, SSF2, SSF1F2: double;
    MSDep, MSErr, MSF1, MSF2, MSF1F2: double;
    DFTot, DFErr, DFF1, DFF2, DFF1F2: double;
    Omega, OmegaF1, OmegaF2, OmegaF1F2: double;
    FF1, FF2, FF1F2: Double;
    ProbF1, ProbF2, ProbF1F2: double;
  end;


  { TSRHTestForm }

  TSRHTestForm = class(TBasicStatsReportAndChartForm)
    Bevel1: TBevel;
    ChartStyles: TChartStyles;
    DepIn: TBitBtn;
    DepOut: TBitBtn;
    DepVar: TEdit;
    Fact1In: TBitBtn;
    Fact1Out: TBitBtn;
    Fact2In: TBitBtn;
    Fact2Out: TBitBtn;
    Factor1: TEdit;
    Factor2: TEdit;
    Label1: TLabel;
    Label3: TLabel;
    OverallAlpha: TEdit;
    StaticText1: TStaticText;
    StaticText2: TStaticText;
    StaticText3: TStaticText;
    VarList: TListBox;
    procedure DepInClick(Sender: TObject);
    procedure DepOutClick(Sender: TObject);
    procedure Fact1InClick(Sender: TObject);
    procedure Fact1OutClick(Sender: TObject);
    procedure Fact2InClick(Sender: TObject);
    procedure Fact2OutClick(Sender: TObject);
    procedure VarListDblClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);

  private
    F1Labels: StrDyneVec;
    F2Labels: StrDyneVec;

    procedure GetColsSelected(out AColsSelected: IntDyneVec);

    function GetLabels(AColsSelected: IntDyneVec): Boolean;

    function Calc2Way(const AColsSelected: IntDyneVec;
      var AResults: TSRHResults) : Boolean;

    procedure TwoWayTable(const AResults: TSRHResults);

    procedure TwoWayPlot(const AResults: TSRHResults);

  private
    FChartSeries: TChartSeries;
    FChartSources: array[0..3] of TListChartSource;
    FSeriesButtons: array[0..3] of TToolButton;
    procedure PrepareStyles(const ALabels: StrDyneVec);
    procedure SelectPlot(Sender: TObject);

  protected
    procedure AdjustConstraints; override;
    procedure Compute; override;
    procedure UpdateBtnStates; override;
    function Validate(out AMsg: String; out AControl: TWinControl): Boolean; override;

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

var
  SRHTestForm: TSRHTestForm;


implementation

{$R *.lfm}

uses
  Math,
  TAChartUtils, TACustomSource, TALegend, TASeries,
  Utils, MathUnit, GridProcs, ChartFrameUnit;


{ TSRHTestForm }

constructor TSRHTestForm.Create(AOwner: TComponent);
var
  tb: TToolButton;
begin
  inherited;

  FChartSources[0] := TListChartSource.Create(FChartFrame);
  FChartSources[1] := TListChartSource.Create(FChartFrame);
  FChartSources[2] := TListChartSource.Create(FChartFrame);
  FChartSources[3] := TListChartSource.Create(FChartFrame);

  FChartFrame.ChartToolbar.ShowCaptions := true;
  FChartFrame.Charttoolbar.ButtonHeight := 40;

  tb := TToolButton.Create(FChartFrame.ChartToolbar);
  tb.Style := tbsDivider;
  AddButtonToToolbar(tb, FChartFrame.ChartToolbar);

  FSeriesButtons[0] := TToolButton.Create(FChartFrame.ChartToolbar);
  FSeriesButtons[0].Caption := 'Factor 1';
  FSeriesButtons[0].ShowCaption := true;
  FSeriesButtons[0].Style := tbsCheck;
  FSeriesButtons[0].Grouped := true;
  FSeriesButtons[0].Tag := 10;
  FSeriesButtons[0].OnClick := @SelectPlot;
  AddButtonToToolbar(FSeriesButtons[0], FChartFrame.ChartToolbar);

  FSeriesButtons[1] := TToolButton.Create(FChartFrame.ChartToolbar);
  FSeriesButtons[1].Caption := 'Factor 2';
  FSeriesButtons[1].ShowCaption := true;
  FSeriesButtons[1].Style := tbsCheck;
  FSeriesButtons[1].Grouped := true;
  FSeriesButtons[1].Tag := 11;
  FSeriesButtons[1].OnClick := @SelectPlot;
  AddButtonToToolbar(FSeriesButtons[1], FChartFrame.ChartToolbar);

  FSeriesButtons[2] := TToolButton.Create(FChartFrame.ChartToolbar);
  FSeriesButtons[2].Caption := 'Factor 1 * Factor 2';
  FSeriesButtons[2].ShowCaption := true;
  FSeriesButtons[2].Style := tbsCheck;
  FSeriesButtons[2].Grouped := true;
  FSeriesButtons[2].Tag := 12;
  FSeriesButtons[2].OnClick := @SelectPlot;
  AddButtonToToolbar(FSeriesButtons[2], FChartFrame.ChartToolbar);

  FSeriesButtons[3] := TToolButton.Create(FChartFrame.ChartToolbar);
  FSeriesButtons[3].Caption := 'Factor 1 * Factor 2 (rotated)';
  FSeriesButtons[3].ShowCaption := true;
  FSeriesButtons[3].Style := tbsCheck;
  FSeriesButtons[3].Grouped := true;
  FSeriesButtons[3].Tag := 13;
  FSeriesButtons[3].OnClick := @SelectPlot;
  AddButtonToToolbar(FSeriesButtons[3], FChartFrame.ChartToolbar);

  PageControl.ActivePageIndex := 0;
end;


procedure TSRHTestForm.AdjustConstraints;
begin
  inherited;

  ParamsPanel.Constraints.MinWidth := Max(
    4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left,
    2*(Label3.Width + OverallAlpha.BorderSpacing.Left + OverallAlpha.Width) - Fact1In.Width);
  ParamsPanel.Constraints.MinHeight := OverallAlpha.Top + OverallAlpha.Height +
    ButtonBevel.Height + CloseBtn.BorderSpacing.Top + CloseBtn.Height;
end;


// Returns false if an error has occured.
function TSRHTestForm.Calc2Way(const AColsSelected: IntDyneVec;
  var AResults: TSRHResults): Boolean;
var
  X, Xsq: Double;
  i, j: integer;
  grpName: String;
  grp1, grp2: integer;
  constant, rowsTotCnt, colsTotCnt, SSCells: double;
  groupSize: Integer;
  numF1, numF2: Integer;
begin
  Result := true;

  numF1 := Length(F1Labels);
  numF2 := Length(F2Labels);

  // Allocate memory for arrays
  AResults.CellCounts := nil;
  AResults.CellSums := nil;
  AResults.CellVars := nil;
  AResults.RowSums := nil;
  AResults.ColSums := nil;
  AResults.RowCount := nil;
  AResults.ColCount := nil;
  SetLength(AResults.CellCounts, numF1, numF2);
  SetLength(AResults.CellSums, numF1, numF2);
  SetLength(AResults.CellVars, numF1, numF2);
  SetLength(AResults.RowSums, numF1);
  SetLength(AResults.ColSums, numF2);
  SetLength(AResults.RowCount, numF1);
  SetLength(AResults.ColCount, numF2);

  AResults.TotalCount := 0;
  AResults.MeanDep := 0.0;
  AResults.MeanF1 := 0;
  AResults.MeanF2 := 0;
  AResults.SSDep := 0.0;
  AResults.SSF1 := 0;
  AResults.SSF2 := 0;
  AResults.SSF1F2 := 0;
  SSCells := 0.0;
  rowsTotCnt := 0.0;
  colsTotCnt := 0.0;

  // Get working totals
  for i := 1 to NoCases do
  begin
    if not GoodRecord(OS3MainFrm.DataGrid, i, AColsSelected) then continue;

    // Note: AColsSelected contains: 0: GrpVar, 1: Factor1, 2: Factor2
    X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[AColsSelected[0], i]));
    Xsq := X*X;
    grpName := Trim(OS3MainFrm.DataGrid.Cells[AColsSelected[1], i]);
    grp1 := IndexOfString(F1Labels, grpName);
    grpName := Trim(OS3MainFrm.DataGrid.Cells[AColsSelected[2], i]);
    grp2 := IndexOfString(F2Labels, grpName);
    AResults.CellCounts[grp1, grp2] := AResults.CellCounts[grp1, grp2] + 1;
    AResults.CellSums[grp1, grp2] := AResults.CellSums[grp1, grp2] + X;
    AResults.CellVars[grp1, grp2] := AResults.CellVars[grp1, grp2] + Xsq;
    AResults.RowSums[grp1] := AResults.RowSums[grp1] + X;
    AResults.ColSums[grp2] := AResults.ColSums[grp2] + X;
    AResults.RowCount[grp1] := AResults.RowCount[grp1] + 1;
    AResults.ColCount[grp2] := AResults.ColCount[grp2] + 1;
    AResults.MeanDep := AResults.MeanDep + X;
    AResults.SSDep := AResults.SSDep + Xsq;
    AResults.TotalCount := AResults.TotalCount + 1;
  end;

  // Check for equal groups
  groupSize := AResults.CellCounts[0, 0];
  for i := 0 to numF1-1 do
    for j := 0 to numF2-1 do
      if AResults.CellCounts[i, j] <> groupSize then
      begin
        ErrorMsg('All groups must have the same size.');
        Result := false;
        exit;
      end;

  // Calculate results
  for i := 0 to numF1-1 do
  begin
    AResults.SSF1 := AResults.SSF1 + (sqr(AResults.RowSums[i]) / AResults.RowCount[i]);
    RowsTotCnt := RowsTotCnt + AResults.RowCount[i];
  end;

  for j := 0 to numF2-1 do
  begin
    AResults.SSF2 := AResults.SSF2 + (sqr(AResults.ColSums[j]) / AResults.ColCount[j]);
    ColsTotCnt := ColsTotCnt + AResults.ColCount[j];
  end;

  for i := 0 to numF1-1 do
    for j := 0 to numF2-1 do
      if AResults.CellCounts[i,j] > 0 then
        SSCells := SSCells + (sqr(AResults.CellSums[i,j]) / AResults.CellCounts[i,j]);

  if AResults.TotalCount > 0 then
    constant := sqr(AResults.MeanDep) / AResults.TotalCount
  else
    constant := 0.0;

  AResults.SSF1 := AResults.SSF1 - constant;
  AResults.SSF2 := AResults.SSF2 - constant;
  AResults.SSF1F2 := SSCells - AResults.SSF1 - AResults.SSF2 - constant;
  AResults.SSErr := AResults.SSDep - SSCells;
  AResults.SSDep := AResults.SSDep - constant;

  if (AResults.SSF1F2 < 0) or (AResults.SSF1 < 0) or (AResults.SSF2 < 0) then
  begin
    ErrorMsg('A negative sum of squares has been found. Unbalanced design? Ending analysis.');
    Result := false;
    exit;
  end;

  // Degrees of freedom
  with AResults do
  begin
    DFTot := TotalCount - 1;
    DFF1 := numF1 - 1;
    DFF2 := numF2 - 1;
    DFF1F2 := DFF1 * DFF2;
    DFErr := DFTot - DFF1 - DFF2 - DFF1F2;
  end;

  // Mean standard error
  with AResults do
  begin
    MSF1 := SSF1 / DFF1;
    MSF2 := SSF2 / DFF2;
    MSF1F2 := SSF1F2 / DFF1F2;
    MSErr := SSErr / DFErr;
    MSDep := SSDep / DFTot;
  end;

  // Omega
  with AResults do
  begin
    OmegaF1 := (SSF1 - DFF1 * MSErr) / (SSDep + MSErr);
    OmegaF2 := (SSF2 - DFF2 * MSErr) / (SSDep + MSErr);
    OmegaF1F2 := (SSF1F2 - DFF1F2 * MSErr) / (SSDep + MSErr);
    Omega := OmegaF1 + OmegaF2 + OmegaF1F2;
  end;

  // Mean of dependent variable
  with AResults do
    MeanDep := MeanDep / TotalCount;

  // F tests for fixed effects
  with AResults do
  begin
    FF1 := abs(MSF1 / MSErr);
    FF2 := abs(MSF2 / MSErr);
    FF1F2 := abs(MSF1F2 / MSErr);

    ProbF1 := ProbF(FF1, DFF1, DFErr);
    ProbF2 := ProbF(FF2, DFF2, DFErr);
    ProbF1F2 := ProbF(FF1F2, DFF1F2, DFErr);
    if (ProbF1 > 1.0) then ProbF1 := 1.0;
    if (ProbF2 > 1.0) then ProbF2 := 1.0;
    if (ProbF1F2 > 1.0) then ProbF1F2 := 1.0;

    // Obtain omega squared (proportion of dependent variable explained)
    if (OmegaF1 < 0.0) then OmegaF1 := 0.0;
    if (OmegaF2 < 0.0) then OmegaF2 := 0.0;
    if (OmegaF1F2 < 0.0) then OmegaF1F2 := 0.0;
    if (Omega < 0.0) then Omega := 0.0;
  end;
end;


procedure TSRHTestForm.Compute;
var
  colNoSelected: IntDyneVec = nil;
  res: TSRHResults;
begin
  //  Get column numbers of dependent variable and factors
  GetColsSelected(colNoSelected);

  // get labels for each factor level
  if not GetLabels(colNoSelected) then
    exit;

  // Calculate two-way ANOVA
  res := Default(TSRHResults);
  if not Calc2Way(colNoSelected, res) then
    exit;

  // Display the ANOVA etc in a report
  TwoWayTable(res);

  // Plot the means.
  TwoWayPlot(res);
end;


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


procedure TSRHTestForm.DepOutClick(Sender: TObject);
begin
  if DepVar.Text <> '' then
  begin
    VarList.Items.Add(DepVar.Text);
    DepVar.Text := '';
    UpdateBtnStates;
  end;
end;


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


procedure TSRHTestForm.Fact1OutClick(Sender: TObject);
begin
  if Factor1.Text <> '' then
  begin
    VarList.Items.Add(Factor1.Text);
    Factor1.Text := '';
    UpdateBtnStates;
  end;
end;


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


procedure TSRHTestForm.Fact2OutClick(Sender: TObject);
begin
  if Factor2.Text <> '' then
  begin
    VarList.Items.Add(Factor2.Text);
    Factor2.Text := '';
    UpdateBtnStates;
  end;
end;


procedure TSRHTestForm.GetColsSelected(out AColsSelected: IntDyneVec);
begin
  AColsSelected := nil;
  SetLength(AColsSelected, 3);
  AColsSelected[0] := GetVariableIndex(OS3MainFrm.DataGrid, DepVar.Text);
  AColsSelected[1] := GetVariableIndex(OS3MainFrm.DataGrid, Factor1.Text);
  AColsSelected[2] := GetVariableIndex(OS3MainFrm.DataGrid, Factor2.Text);
end;


function TSRHTestForm.GetLabels(AColsSelected: IntDyneVec): Boolean;
var
  i, nF1, nF2: integer;
  floatVal: Double;
  s: String;
begin
  Result := false;

  // Check validity of values of dependent variable
  for i := 1 to NoCases do
  begin
    if not GoodRecord(OS3MainFrm.DataGrid, i, AColsSelected) then continue;
    if not TryStrToFloat(OS3MainFrm.DataGrid.Cells[AColsSelected[0], i], floatVal) then
    begin
      ErrorMsg('No valid number in cell at row %d and column %d', [i, AColsSelected[0]]);
      exit;
    end;
  end;

  // Extract labels of the Factor 1 and Factor 2 cases.
  F1Labels := nil;
  F2Labels := nil;
  Setlength(F1Labels, NoCases);   // over-dimensioned, will be trimmed later.
  SetLength(F2Labels, NoCases);
  nF1 := 0;
  nF2 := 0;
  for i := 1 to NoCases do
  begin
    if not GoodRecord(OS3MainFrm.DataGrid, i, AColsSelected) then continue;

    // Index 1 identifies Factor 1
    s := OS3MainFrm.DataGrid.Cells[AColsSelected[1], i];
    if IndexOfString(F1Labels, s) = -1 then
    begin
      F1Labels[nF1] := s;
      inc(nF1);
    end;

    // Index 2 indentifies Factor 2
    s := OS3MainFrm.DataGrid.Cells[AColsSelected[2], i];
    if IndexOfString(F2Labels, s) = -1 then
    begin
      F2Labels[nF2] := s;
      inc(nF2);
    end;
  end;

  // Trim the label arrays to correct length
  SetLength(F1Labels, nF1);
  SetLength(F2Labels, nF2);

  Result := true;
end;


procedure TSRHTestForm.TwoWayTable(const AResults: TSRHResults);
var
  MinVar, MaxVar, sumvars, sumDFrecip: double;
  i, j, groupSize, numF1, numF2: integer;
  XBar, V, S, RowSS, ColSS: double;
  sumFreqLogVar, c, bartlett, cochran, hartley, chiProb: double;
  H, HProb: double;
  lReport: TStrings;
begin
  lReport := TStringList.Create;
  try
    lReport.Add('TWO WAY ANALYSIS OF VARIANCE');
    lReport.Add('');
    lReport.Add('Variable analyzed:           %s', [DepVar.Text]);
    lReport.Add('Factor A (rows) variable:    %s', [Factor1.Text]);
    lReport.Add('Factor B (columns) variable: %s', [Factor2.Text]);
    lReport.Add('');
    lReport.Add('SOURCE         D.F.      SS          MS        F       PROB.> F  Omega Sqr.    H     H Prob.');
    lReport.Add('-------------  ----  ----------  ----------  --------  --------  ----------  ------  -------');

    with AResults do
    begin
      H := SSF1 / MSDep;
      HProb := 1.0 - ChiSquaredProb(H,round(DFF1));
      lReport.Add('Among Rows     %4.0f  %10.3f  %10.3f  %8.3f  %8.3f     %6.3f   %6.3f   %6.3f',
        [DFF1, SSF1, MSF1, FF1, ProbF1, OmegaF1, H, HProb]);

      H := SSF2 / MSDep;
      HProb := 1.0 - ChiSquaredProb(H,round(DFF2));
      lReport.Add('Among Columns  %4.0f  %10.3f  %10.3f  %8.3f  %8.3f     %6.3f   %6.3f   %6.3f',
       [DFF2, SSF2, MSF2, FF2, ProbF2, OmegaF2, H, HProb]);

      H := SSF1F2 / MSDep;
      HProb := 1.0 - ChiSquaredProb(H,round(DFF1F2));
      lReport.Add('Interaction    %4.0f  %10.3f  %10.3f  %8.3f  %8.3f     %6.3f   %6.3f   %6.3f',
        [DFF1F2, SSF1F2, MSF1F2, FF1F2, ProbF1F2, OmegaF1F2, H, HProb]);
      lReport.Add('Within Groups  %4.0f  %10.3f  %10.3f',
        [DFErr, SSErr, MSErr]);
      lReport.Add('Total          %4.0f  %10.3f  %10.3f',
        [DFTot, SSDep, MSDep]);
      lReport.Add('');
      lReport.Add('Omega squared for combined effects: %10.3f', [Omega]);
    end;

    lReport.Add('');
    lReport.Add(DIVIDER_SMALL_AUTO);
    lReport.Add('');

    lReport.Add('DESCRIPTIVE STATISTICS');
    lReport.Add('');
    lReport.Add('GROUP     Row         Col.      N      MEAN     VARIANCE    STD.DEV. ');
    lReport.Add('-----  ----------  ----------  ----  --------  ----------  ----------');

    MaxVar := -1E304;
    MinVar := 1E304;
    sumvars := 0.0;
    sumFreqLogVar := 0.0;
    sumDFrecip := 0.0;

    // Display cell means, variances, standard deviations
    numF1 := Length(F1Labels);
    numF2 := Length(F2Labels);

    V := 0.0;
    XBar := 0.0;
    S := 0.0;
    for i := 0 to numF1-1 do
    begin
      for j := 0 to numF2-1 do
      begin
        if AResults.CellCounts[i,j] > 1 then
        begin
          XBar := AResults.CellSums[i,j] / AResults.CellCounts[i,j];
          V := AResults.CellVars[i,j] - sqr(AResults.CellSums[i,j]) / AResults.CellCounts[i,j];
          V := V / (AResults.CellCounts[i,j] - 1.0);
          S := sqrt(V);
          sumVars  := sumVars + V;
          if V > MaxVar then MaxVar := V;
          if V < MinVar then MinVar := V;
          sumDFrecip := sumDFrecip + (1.0 / (AResults.CellCounts[i,j] - 1));
          sumFreqLogVar := sumFreqLogVar + ((AResults.CellCounts[i,j] - 1) * ln(V));
        end;
        lReport.Add('Cell   %-10s  %-10s  %4d  %8.3f  %10.3f  %10.3f',
          [F1Labels[i], F2Labels[j], AResults.CellCounts[i, j], XBar, V, S]);
      end;
    end;

    lReport.Add('');

    //Display Row means, variances, standard deviations
    for i := 0 to numF1-1 do
    begin
      XBar := AResults.RowSums[i] / AResults.RowCount[i];
      RowSS := 0.0;
      for j := 0 to numF2-1 do
        RowSS := RowSS + AResults.CellVars[i,j];
      V := RowSS - sqr(AResults.RowSums[i]) / AResults.RowCount[i];
      V := V / (AResults.RowCount[i] - 1);
      S := sqrt(V);
      lReport.Add('Row    %-10s              %4d  %8.3f  %10.3f  %10.3f',
        [F1Labels[i], AResults.RowCount[i], XBar, V, S]);
    end;

    lReport.Add('');

    //Display means, variances and standard deviations for columns
    for j := 0 to numF2-1 do
    begin
      XBar := AResults.ColSums[j] / AResults.ColCount[j];
      ColSS := 0.0;
      for i := 0 to numF1-1 do
        ColSS := ColSS + AResults.CellVars[i,j];

      V := NaN;
      S := NaN;
      if AResults.ColCount[j] > 0 then
      begin
        V := ColSS - sqr(AResults.ColSums[j]) / AResults.ColCount[j];
        if AResults.ColCount[j] > 1 then
        begin
          V := V / (AResults.ColCount[j] - 1);
          if V > 0.0 then
            S := sqrt(V);
        end;
      end;
      lReport.Add('Col                %-10s  %4d  %8.3f  %10.3f  %10.3f',
        [F2Labels[j], AResults.ColCount[j], XBar, V, S]);
    end;

    lReport.Add('');

    lReport.Add(  'TOTAL                          %4d  %8.3f  %10.3f  %10.3f',
      [AResults.TotalCount, AResults.MeanDep, AResults.MSDep, sqrt(AResults.MSDep)]);

    lReport.Add('');
    lReport.Add(DIVIDER_SMALL_AUTO);
    lReport.Add('');

    c := 1.0 + (1.0 / (3.0 * numF1 * numF2 - 1.0)) * (sumDFrecip - (1.0 / AResults.DFErr));
    bartlett := (2.303 / c) * ((AResults.DFErr * ln(AResults.MSErr)) - sumfreqlogvar);
    chiProb := 1.0 - chisquaredprob(bartlett,round(numF1 * numF2 - 1));
    cochran := maxvar / sumvars;
    hartley := maxvar / minvar;
    groupSize := AResults.CellCounts[0, 0];

    LReport.Add('TESTS FOR HOMOGENEITY OF VARIANCE');
    lReport.Add('');
    lReport.Add('Hartley FMax test statistic:   %10.3f', [hartley]);
    lReport.Add('  with %d and %d degrees of freedom', [numF1*numF2, groupsize-1]);
    lReport.Add('');
    lReport.Add('Cochran C statistic:           %10.3f', [cochran]);
    lReport.Add('  with %d and %d degrees of freedom', [numF1*numF2, groupsize - 1]);
    lReport.Add('');
    lReport.Add('Bartlett Chi-square statistic: %10.3f', [bartlett]);
    lReport.Add('  with %d degrees and freedom', [numF1*numF2-1]);
    lReport.Add('  and probability > value      %10.3f', [chiprob]);

    FReportFrame.DisplayReport(lReport);

  finally
    lReport.Free;
  end;
end;


procedure TSRHTestForm.TwoWayPlot(const AResults: TSRHResults);
var
  i, j, idx, numF1, numF2: Integer;
  item: PChartDataItem;
begin
  FChartFrame.Clear;
  FChartFrame.SetYTitle('Mean');

  numF1 := Length(F1Labels);
  numF2 := Length(F2Labels);

  // Chart source for Factor A
  FChartSources[0].Clear;
  for i := 0 to numF1-1 do
    FChartSources[0].Add(i+1, AResults.RowSums[i] / AResults.RowCount[i], F1Labels[i]);

  // Chart source for Factor B
  FChartSources[1].Clear;
  for j := 0 to numF2-1 do
    FChartSources[1].Add(j+1, AResults.ColSums[j] / AResults.ColCount[j], F2Labels[j]);

  // Chart source for interaction for Factors A and B: x = Factor 2
  FChartSources[2].Clear;
  FChartSources[2].YCount := numF1;
  for j := 0 to numF2-1 do
  begin
    idx := FChartSources[2].Add(j+1, NaN, F2Labels[j]);    // y values will be added separately
    item := FChartSources[2].Item[idx];
    for i := 0 to numF1-1 do
      item^.SetY(i, AResults.CellSums[i,j] / AResults.CellCounts[i,j]);
  end;

  // Chart source for interaction for Factors A and B: x = Factor 1
  FChartSources[3].Clear;
  FChartSources[3].YCount := numF2;
  for i := 0 to numF1-1 do
  begin
    idx := FChartSources[3].Add(i+1, NaN, F1Labels[i]);     // y values will be added separately
    item := FChartSources[3].Item[idx];
    for j := 0 to numF2-1 do
      item^.SetY(j, AResults.CellSums[i,j] / AResults.CellCounts[i,j]);
  end;

  // Create series if required
  if FChartSeries = nil then
    FChartSeries := FChartFrame.PlotXY(ptLinesAndSymbols, nil, nil, nil, nil, '', DATA_COLORS[0]);

  (FChartSeries as TLineSeries).Styles := ChartStyles;

  // by default display Factor A
  SelectPlot(FSeriesButtons[0]);
end;


procedure TSRHTestForm.PrepareStyles(const ALabels: StrDyneVec);
var
  i: Integer;
begin
  ChartStyles.Styles.Clear;
  for i := 0 to High(ALabels) do
    with TChartStyle(ChartStyles.Styles.Add) do
    begin
      Brush.Color := DATA_COLORS[i mod Length(DATA_COLORS)];
      Pen.Color := DATA_COLORS[i mod Length(DATA_COLORS)];
      Text := ALabels[i];
      UsePen := true;
      UseBrush := true;
    end;
end;


procedure TSRHTestForm.SelectPlot(Sender: TObject);
var
  tb: TToolButton;
begin
  tb := Sender as TToolButton;
  tb.Down := true;

  if FChartSeries = nil then
    exit;

  FChartSeries.Source := FChartSources[tb.Tag - 10];
  case tb.Tag-10 of
    0: begin  // x = Factor 1
         FChartFrame.SetTitle(Factor1.Text);
         FChartFrame.SetXTitle('Variable ' + Factor1.Text);
         FChartFrame.Chart.Legend.Visible := false;
         FChartSeries.Legend.Multiplicity := lmSingle;
       end;
    1: begin  // x = Factor 2
         FChartFrame.SetTitle(Factor2.Text);
         FChartFrame.SetXTitle('Variable ' + Factor2.Text);
         FChartFrame.Chart.Legend.Visible := false;
         FChartSeries.Legend.Multiplicity := lmSingle;
       end;
    2: begin  // x = Factor 2, multiple curves for all factors 1
         FChartFrame.SetTitle(Factor1.Text + ' * ' + Factor2.Text);
         FChartFrame.SetXTitle('Variable ' + Factor2.Text);
         FChartFrame.Chart.Legend.Visible := true;
         FChartSeries.Legend.Multiplicity := lmStyle;
         PrepareStyles(F1Labels);
       end;
    3: begin  // x = Factor 1, multiple curves for all factors 2
         FChartFrame.SetTitle(Factor1.Text + ' * ' + Factor2.Text);
         FChartFrame.SetXTitle('Variable ' + Factor1.Text);
         FChartFrame.Chart.Legend.Visible := true;
         FChartSeries.Legend.Multiplicity := lmStyle;
         PrepareStyles(F2Labels);
       end;
  end;

  FChartFrame.Chart.BottomAxis.Marks{%H-}.Source := FChartSeries.Source;
  FChartFrame.Chart.BottomAxis.Marks{%H-}.Style := smsLabel;
end;


procedure TSRHTestForm.Reset;
var
  i : integer;
begin
  FreeAndNil(FChartSeries);

  inherited;

  if FChartSources[0] <> nil then
  begin
    //SelectPlot(FSeriesButtons[0]);
    for i := 0 to High(FChartSources) do
      FChartSources[i].Clear;
  end;

  DepVar.Clear;
  Factor1.Clear;
  Factor2.Clear;
  OverAllalpha.Text := FormatFloat('0.00', DEFAULT_ALPHA_LEVEL);
  VarList.Clear;
  for i := 1 to NoVariables do
    VarList.Items.Add(OS3MainFrm.DataGrid.Cells[i, 0]);
		
  UpdateBtnStates;
end;


procedure TSRHTestForm.UpdateBtnStates;
begin
  inherited;

  DepIn.Enabled := (VarList.ItemIndex > -1) and (DepVar.Text = '');
  Fact1In.Enabled := (Varlist.ItemIndex > -1) and (Factor1.Text = '');
  Fact2In.Enabled := (VarList.ItemIndex > -1) and (Factor2.Text = '');

  DepOut.Enabled := DepVar.Text <> '';
  Fact1Out.Enabled := Factor1.Text <> '';
  Fact2Out.Enabled := Factor2.Text <> '';
end;


function TSRHTestForm.Validate(out AMsg: String; out AControl: TWinControl): Boolean;
var
  floatVal: Double;
begin
  Result := false;

  if DepVar.Text = '' then
  begin
    AControl := DepVar;
    AMsg := 'No dependent variable selected.';
    exit;
  end;

  if Factor1.Text = '' then
  begin
    AControl := Factor1;
    AMsg := 'Factor 1 variable not selected.';
    exit;
  end;

  if Factor2.Text = '' then
  begin
    AControl := Factor2;
    AMsg := 'Factor 2 variable not selected.';
    exit;
  end;

  if OverallAlpha.Text = '' then
  begin
    AControl := OverallAlpha;
    AMsg := 'Alpha level not selected.';
    exit;
  end;

  if not TryStrToFloat(OverallAlpha.Text, floatVal) or (floatVal <= 0) or (floatVal >= 1) then
  begin
    AControl := OverallAlpha;
    AMsg := 'Alpha level must be a valid number between 0 and 1.';
    exit;
  end;

  Result := true;
end;


procedure TSRHTestForm.VarListDblClick(Sender: TObject);
var
  index: Integer;
  s: String;
begin
  index := VarList.ItemIndex;
  if index > -1 then
  begin
    s := VarList.Items[index];
    if DepVar.Text = '' then
      DepVar.Text := s
    else if Factor1.Text = '' then
      Factor1.Text := s
    else if Factor2.Text = '' then
      Factor2.Text := s;
    VarList.Items.Delete(index);
    UpdateBtnStates;
  end;
end;


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



end.