lazarus-ccr/applications/lazstats/source/forms/analysis/nonparametric/kwanovaunit.pas

// File for testing: "kwanova.laz"

unit KWANOVAUnit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
  StdCtrls, Buttons, ExtCtrls,
  MainUnit, FunctionsLib, Globals, DataProcs, BasicStatsReportFormUnit;

type

  { TKWAnovaForm }

  TKWAnovaForm = class(TBasicStatsReportForm)
    AlphaEdit: TEdit;
    OptionsGroup: TGroupBox;
    Label5: TLabel;
    MWUChk: TCheckBox;
    GrpEdit: TEdit;
    DepEdit: TEdit;
    GrpIn: TBitBtn;
    GrpOut: TBitBtn;
    DepIn: TBitBtn;
    DepOut: TBitBtn;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    VarList: TListBox;
    procedure DepInClick(Sender: TObject);
    procedure DepOutClick(Sender: TObject);
    procedure GrpInClick(Sender: TObject);
    procedure GrpOutClick(Sender: TObject);
    procedure VarListDblClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);
  private
	
  protected
    procedure AdjustConstraints; override;
    procedure Compute; override;
    procedure UpdateBtnStates; override;
    function Validate(out AMsg: String; out AControl: TWinControl): Boolean; override;

  public
  procedure Reset; override;
  end; 

var
  KWAnovaForm: TKWAnovaForm;


implementation

{$R *.lfm}

uses
  Math,
  Utils, GridProcs, MatrixUnit;


{ TKWAnovaForm }

procedure TKWAnovaForm.AdjustConstraints;
begin
  inherited;
  ParamsPanel.Constraints.MinWidth := Max(
    4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left,
    2*OptionsGroup.Width);
  ParamsPanel.Constraints.MinHeight := DepOut.Top + DepOut.Height +
    VarList.BorderSpacing.Bottom + OptionsGroup.Height +
    ButtonBevel.Height + CloseBtn.BorderSpacing.Top + CloseBtn.Height;
end;


procedure TKWAnovaForm.Compute;
var
  i, j, k, m, ind_var, dep_var, min_grp, max_grp, group, total_n : integer;
  NoTies, NoTieGroups, nogroups, NoSelected, npairs, n1, n2 : integer;
  largestN : integer;
  ColNoSelected : IntdyneVec = nil;
  group_count : IntDyneVec = nil;
  Ranks: DblDyneMat = nil;
  X : DblDyneMat = nil;
  RankSums: DblDyneVec = nil;
  score, t, SumT, Avg, Probchi, H, CorrectedH, value : double;
  Correction, Temp, TieSum, alpha, U, U2, SD, z, prob : double;
  cellstring, outline: string;
  lReport: TStrings;
begin
  alpha := StrToFloat(AlphaEdit.Text);

  // Allocate array memory
  SetLength(Ranks, NoCases, 2);
  SetLength(X, NoCases, 2);

  // Get column numbers of the independent and dependent variables
  ind_var := GetVariableIndex(OS3MainFrm.DataGrid, GrpEdit.Text);
  dep_var := GetVariableIndex(OS3MainFrm.DataGrid, DepEdit.Text);
  NoSelected := 2;
  SetLength(ColNoSelected, 2);
  ColNoSelected[0] := ind_var;
  ColNoSelected[1] := dep_var;

  // Get minimum and maximum group codes
  total_n := 0;
  min_grp := 10000; //atoi(MainForm.Grid.Cells[ind_var,1].c_str);
  max_grp := -10000;
  for i := 1 to NoCases do
  begin
    if (not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected)) then continue;
    group := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ind_var,i])));
    if (group < min_grp) then min_grp := group;
    if (group > max_grp) then max_grp := group;
    total_n := total_n + 1;
  end;
  nogroups := max_grp - min_grp + 1;

  NoTieGroups := 0;
  SumT := 0.0;

  // Initialize arrays
  SetLength(RankSums,nogroups);
  SetLength(group_count,nogroups);
  for i := 0 to nogroups-1 do
  begin
    group_count[i] := 0;
    RankSums[i] := 0.0;
  end;

  // Setup for printer output
  lReport := TStringList.Create;
  try
    lReport.Add('KRUSKAL-WALLIS ONE-WAY ANALYSIS OF VARIANCE');
    lReport.Add('See pages 184-194 in S. Siegel: Nonparametric Statistics for the Behavioral Sciences');
    lReport.Add('');

    // Get data
    for i := 1 to NoCases do
    begin
      if (not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected)) then continue;
      score := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[dep_var, i]));
      group := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ind_var, i])));
      group := group - min_grp + 1;
      if (group > nogroups) then
      begin
        ErrorMsg('Group codes must be sequential like 1 and 2!');
        exit;
      end;
      group_count[group-1] := group_count[group-1] + 1;
      X[i-1, 0] := score;
      X[i-1, 1] := group;
    end;

    //Sort all scores in ascending order
    for i := 1 to total_n - 1 do
    begin
      for j := i + 1 to total_n do
      begin
        if (X[i-1,0] > X[j-1,0]) then
        begin
          Exchange(X[i-1, 0], X[j-1, 0]);
          Exchange(X[i-1, 1], X[j-1, 1]);
        end;
      end;
    end;

    // Store ranks
    for i := 0 to total_n-1 do
    begin
      Ranks[i,0] := i+1;
      Ranks[i,1] := X[i,1];
    end;

    //Check for ties in ranks - replace with average rank and calculate
    //T for each tie and sum of the T's
    i := 1;
    while i < total_n do
    begin
      j := i + 1;
      TieSum := 0;
      NoTies := 0;
      while (j < total_n) do
      begin
        if (X[j-1,0] > X[i-1,0]) then
          break;
        if (X[j-1,0] = X[i-1,0]) then // match
        begin
          TieSum := TieSum + round(Ranks[j-1,0]);
          NoTies := NoTies + 1;
        end;
        j := j + 1;
      end;

      if (NoTies > 0) then //At least one tie found
      begin
        TieSum := TieSum + Ranks[i-1,0];
        NoTies := NoTies + 1;
        Avg := TieSum / NoTies;
        for j := i to i + NoTies - 1 do Ranks[j-1,0] := Avg;
        t := Power(NoTies,3) - NoTies;
        SumT := SumT + t;
        NoTieGroups := NoTieGroups + 1;
        i := i + (NoTies - 1);
      end;
      i := i + 1;
    end; // next i

    // Calculate sum of ranks in each group
    for i := 0 to total_n-1 do
    begin
      group := round(Ranks[i, 1]);
      RankSums[group-1] := RankSums[group-1] + Ranks[i, 0];
    end;

    // Calculate statistics
    H := 0.0;
    for j := 0 to nogroups-1 do
      H := H + (RankSums[j] * RankSums[j] / (group_count[j]));
    H := H * (12.0 / ( total_n * (total_n + 1)) );
    H := H - (3.0 * (total_n + 1));
    Correction := 1.0 - ( SumT / (Power(total_n,3) - total_n) );
    CorrectedH := H / Correction;
    k := max_grp - min_grp;
    probChi := 1.0 - ChiSquaredProb(H, k);

    // Report results
    lReport.Add('     Score     Rank      Group');
    lReport.Add('');
    for i := 0 to total_n-1 do
      lReport.Add('%10.2f %10.2f %10.0f', [X[i,0], Ranks[i,0], Ranks[i,1]]);
    lReport.Add('');
    lReport.Add('Sum of Ranks in each Group');
    lReport.Add('Group   Sum    No. in Group');
    for i := 0 to noGroups-1 do
      lReport.Add('%3d  %10.2f %5d', [i+min_grp, RankSums[i], group_count[i]]);
    lReport.Add('');
    lReport.Add('No. of tied rank groups           %8d',   [NoTieGroups]);
    lReport.Add('Statistic H uncorrected for ties: %8.4f', [H]);
    lReport.Add('Correction for Ties:              %8.4f', [Correction]);
    lReport.Add('Statistic H corrected for ties:   %8.4f', [CorrectedH]);
    lReport.Add('Corrected H is approx. chi-square with %d degrees of freedom and probability %.4f', [k, ProbChi]);

      if MWUChk.Checked then
      begin
        lReport.Add('');
        lReport.Add('------------------------------------------------------------------------');
        lReport.Add('');

        // Do Mann-Whitney U tests on group pairs
        alpha := StrToFloat(AlphaEdit.Text);
        npairs := nogroups * (nogroups - 1) div 2;
        alpha := alpha / npairs;
        lReport.Add('New alpha for %d paired comparisons: %.3f', [npairs, alpha]);
        for i := 1 to nogroups - 1 do
        begin
          for j := i + 1 to nogroups do
          begin
            // Setup for printer output
            lReport.Add('');
            lReport.Add('');
            lReport.Add('MANN-WHITNEY U TEST');
            lReport.Add('See pages 116-127 in S. Siegel: Nonparametric Statistics for the Behavioral Sciences');
            lReport.Add('');
            lReport.Add('Comparison of group %d with group %d', [i, j]);

            group_count[0] := 0;
            group_count[1] := 0;
            RankSums[0] := 0;
            RankSums[1] := 0;
            total_n := 0;
            for k := 1 to NoCases do
            begin
                if (not GoodRecord(OS3MainFrm.DataGrid, k, ColNoSelected)) then continue;
                score := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[dep_var, k]));
                value := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[ind_var, k]));
                if round(value) = i then
                begin
                     X[total_n,0] := score;
                     X[total_n,1] := value;
                     group_count[0] := group_count[0] + 1;
                     total_n := total_n + 1;
                end;
                if round(value) = j then
                begin
                     X[total_n,0] := score;
                     X[total_n,1] := value;
                     group_count[1] := group_count[1] + 1;
                     total_n := total_n + 1;
                end;
            end; // next case k

            //Sort all scores in ascending order
            for k := 1 to total_n - 1 do
            begin
                for m := k + 1 to total_n do
                begin
                    if (X[k-1,0] > X[m-1,0]) then
                    begin
                        Temp := X[k-1,0];
                        X[k-1,0] := X[m-1,0];
                        X[m-1,0] := Temp;
                        Temp := X[k-1,1];
                        X[k-1,1] := X[m-1,1];
                        X[m-1,1] := Temp;
                    end;
                end;
            end;

            // get ranks for these two groups
            for k := 0 to total_n-1 do
            begin
                Ranks[k,0] := k+1;
                Ranks[k,1] := X[k,1];
            end;

            //Check for ties in ranks - replace with average rank and calculate
            //T for each tie and sum of the T's
            NoTieGroups := 0;
            k := 1;
            while k < total_n do
            begin
                m := k + 1;
                TieSum := 0;
                NoTies := 0;
                while (m < total_n) do
                begin
                    if (X[m-1,0] > X[k-1,0]) then
                      Break;
                    if (X[m-1,0] = X[k-1,0]) then // match
                    begin
                        TieSum := TieSum + round(Ranks[m-1,0]);
                        NoTies := NoTies + 1;
                    end;
                    m := m + 1;
                end;
                if (NoTies > 0) then //At least one tie found
                begin
                    TieSum := TieSum + Ranks[k-1,0];
                    NoTies := NoTies + 1;
                    Avg := TieSum / NoTies;
                    for m := k to k + NoTies - 1 do Ranks[m-1,0] := Avg;
                    t := Power(NoTies,3) - NoTies;
                    SumT := SumT + t;
                    NoTieGroups := NoTieGroups + 1;
                    k := k + (NoTies - 1);
                end;
                k := k + 1;
            end; // next k

            // Calculate sum of ranks in each group
            for k := 1 to total_n do
            begin
                group := round(Ranks[k-1,1]);
                RankSums[group-1] := RankSums[group-1] + Ranks[k-1,0];
            end;

            //Calculate U for larger and smaller groups
            n1 := group_count[0];
            n2 := group_count[1];
            if (n1 > n2) then
            begin
                group := i-1;
                U := (n1 * n2) + ((n1 * (n1 + 1)) / 2.0) - RankSums[group];
            end
            else
            begin
                group := j - 1;
                U := (n1 * n2) + ((n2 * (n2 + 1)) / 2.0) - RankSums[group];
            end;
            U2 := (n1 * n2) - U;
            SD := (n1 * n2 * (n1 + n2 + 1)) / 12.0;
            SD := sqrt(SD);
            if (U2 > U) then z := (U2 - (n1 * n2 / 2)) / SD
            else z := (U - (n1 * n2 / 2)) / SD;
            prob := 1.0 - probz(z);

            // Report results
            lReport.Add('     Score     Rank      Group');
            lReport.Add('');
            for k := 1 to total_n do
                lReport.Add('%10.2f %10.2f %10.0f', [X[k-1,0], Ranks[k-1,0], Ranks[k-1,1]]);
            lReport.Add('');
            lReport.Add('Sum of Ranks in each Group');
            lReport.Add('Group   Sum    No. in Group');
            group := i - 1;
            lReport.Add('%3d  %10.3f %5d', [i, RankSums[group], group_count[0]]);
            group := j - 1;
            lReport.Add('%3d  %10.3f %5d', [j, RankSums[group], group_count[1]]);
            lReport.Add('');
            lReport.Add(            'No. of tied rank groups:          %8d', [NoTieGroups]);
            if (n1 > n2) then largestN := n1 else largestN := n2;
            if (largestN < 20) then
                outline := Format(  'Statistic U:                      %8.4f',[U])
            else
            begin
                if (U > U2) then
                  outline := Format('Statistic U:                      %8.4f',[U])
                else
                  outline := Format('Statistic U:                      %8.4f',[U2]);
            end;
            lReport.Add(outline);
            lReport.Add(            'z Statistic (corrected for ties): %8.4f', [z]);
            lReport.Add(            'Prob. > z:                        %8.4f', [prob]);
            if (n2 < 20) then
            begin
                lReport.Add('z test is approximate.  Use tables of exact probabilities in Siegel.');
                lReport.Add('(Table J or K, pages 271-277)');
            end;
          end; // next group j
        end; // next group i
      end;

    if lReport.Count > 0 then
      FReportFrame.DisplayReport(lReport);

  finally
    lReport.Free;
  end;
end;


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


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


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


procedure TKWAnovaForm.GrpOutClick(Sender: TObject);
begin
  if GrpEdit.Text <> '' then
  begin
    VarList.Items.Add(GrpEdit.Text);
    GrpEdit.Text := '';
  end;
  UpdateBtnStates;
end;


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

  GrpEdit.Clear;
  DepEdit.Clear;
  AlphaEdit.Text := FormatFloat('0.00', DEFAULT_ALPHA_LEVEL);
  MWUChk.Checked := false;
  VarList.Items.Clear;
  for i := 1 to NoVariables do
    VarList.Items.Add(OS3MainFrm.DataGrid.Cells[i,0]);
  UpdateBtnStates;
end;


procedure TKWAnovaForm.UpdateBtnStates;
begin
  inherited;
  GrpIn.Enabled := (VarList.Items.Count > 0) and (GrpEdit.Text = '');
  DepIn.Enabled := (VarList.Items.Count > 0) and (DepEdit.Text = '');
  GrpOut.Enabled := (GrpEdit.Text <> '');
  DepOut.Enabled := (DepEdit.Text <> '');
end;


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

  if (NoVariables < 1) then
  begin
    AMsg := 'You must have grid data!';
    AControl := VarList;
    exit;
  end;

  if GrpEdit.Text = '' then
  begin
    AMsg := 'Group variable not specified.';
    AControl := GrpEdit;
    exit;
  end;

  if DepEdit.Text = '' then
  begin
    AMsg := 'Dependent variable not selected.';
    AControl := DepEdit;
    exit;
  end;

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

  Result := true;
end;


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


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


end.