lazarus-ccr/applications/lazstats/source/forms/analysis/cross-classification/loglinscreenunit.pas

// File for testing: ABCLogLinData.laz
// Use all variables in original order
// Click the button and define Min/max values:
//   variable "Row":   1 .. 2
//   variable "Col":   1 .. 2
//   variable "Slice": 1 .. 3
//   variable  "X":    1 .. 9
//
// NOTE: the calculation crashes

unit LogLinScreenUnit;

{$mode objfpc}{$H+}

interface

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

type

  { TLogLinScreenForm }

  TLogLinScreenForm = class(TBasicStatsReportForm)
    Bevel2: TBevel;
    InBtn: TBitBtn;
    Label1: TLabel;
    OutBtn: TBitBtn;
    AllBtn: TBitBtn;
    MarginsChk: TCheckBox;
    GenlModelChk: TCheckBox;
    OptionsGroup: TGroupBox;
    Label2: TLabel;
    Label3: TLabel;
    SelectList: TListBox;
    MinMaxGrid: TStringGrid;
    VarList: TListBox;
    CountVarChk: TCheckBox;
    procedure AllBtnClick(Sender: TObject);
    procedure CountVarChkChange(Sender: TObject);
    procedure InBtnClick(Sender: TObject);
    procedure OutBtnClick(Sender: TObject);
    procedure SelectListSelectionChange(Sender: TObject; {%H-}User: boolean);
    procedure VarListDblClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);

  private
    function ArrayPosition(ANumDims: integer;
      const Subscripts, DimSize: IntDyneVec): integer;

    procedure UpdateMinMaxGrid;

    procedure Screen(VAR NVAR : integer;
                 VAR MP : integer; VAR MM : integer;
                 VAR NTAB : integer; VAR TABLE : DblDyneVec;
                 VAR DIM : IntDyneVec; VAR GSQ : DblDyneVec;
                 VAR DGFR : IntDyneVec; VAR PART : DblDyneMat;
                 VAR MARG : DblDyneMat; VAR DFS : IntDyneMat;
                 VAR IP : IntDyneMat; VAR IM : IntDyneMat;
                 VAR ISET : IntDyneVec; VAR JSET : IntDyneVec;
                 VAR CONFIG : IntDyneMat; VAR FIT : DblDyneVec;
                 VAR SIZE : IntDyneVec; VAR COORD : IntDyneVec;
                 VAR X : DblDyneVec; VAR Y : DblDyneVec;
                 VAR IFAULT : integer);

    procedure CONF(VAR N : integer; VAR M : integer;
               VAR MP : integer;
               VAR MM : integer;
               VAR ISET : IntDyneVec; VAR JSET : IntDyneVec;
               VAR IP : IntDyneMat; VAR IM : IntDyneMat; VAR NP : integer);

    procedure COMBO(VAR ISET : IntDyneVec;
                   N, M : Integer;
                   VAR LAST : boolean);

    procedure EVAL(VAR IAR : IntDyneMat;
                 NC, NV, IBEG, NVAR, MAX : integer;
                 VAR CONFIG : IntDyneMat;
                 VAR DIM : IntDyneVec; VAR DF : integer);

    procedure RESET(VAR FIT : DblDyneVec; NTAB : Integer;
                 AVG : Double);

    procedure LIKE(var GSQ: Double; const FIT, TABLE: DblDyneVec; NTAB: integer);

    procedure LOGFIT(NVAR, NTAB, NCON : integer;
                 VAR DIM : IntDyneVec;
                 VAR CONFIG : IntDyneMat; VAR TABLE : DblDyneVec;
                 VAR FIT : DblDyneVec; VAR SIZE : IntDyneVec;
                 VAR COORD : IntDyneVec; VAR X : DblDyneVec;
                 VAR Y : DblDyneVec);

    procedure MaxCombos(NumDims: integer; out MM, MP: integer);

    procedure Marginals(NumDims, ArraySize: integer; const Indexes: IntDyneMat;
      const Data: DblDyneVec; const Margins: IntDyneMat);

  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
  LogLinScreenForm: TLogLinScreenForm;

implementation

{$R *.lfm}

uses
  Math, LCLIntf, LCLType, 
  Utils, DataProcs, GridProcs;


{ TLogLinScreenForm }

procedure TLogLinScreenForm.AdjustConstraints;
begin
  inherited;
  ParamsPanel.Constraints.MinWidth := MaxValue([
    4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left,
    OptionsGroup.Width,
    Label2.Width * 2 + AllBtn.Width + 2*varList.BorderSpacing.Right +
      MinMaxGrid.Width + MinMaxGrid.BorderSpacing.Left
    ]);
  ParamsPanel.Constraints.MinHeight := AllBtn.Top + AllBtn.Height +
    CountVarChk.BorderSpacing.Top + CountVarChk.Height +
    VarList.BorderSpacing.Bottom + OptionsGroup.Height + ButtonBevel.Height +
    CloseBtn.BorderSpacing.Top + CloseBtn.Height;
end;
 

procedure TLogLinScreenForm.AllBtnClick(Sender: TObject);
var
  i: integer;
begin
  for i := 0 to VarList.Items.Count-1 do
    SelectList.Items.Add(VarList.Items[i]);
  VarList.Clear;
  UpdateBtnStates;
  UpdateMinMaxGrid;
end;


function TLogLinScreenForm.ArrayPosition(ANumDims: integer;
  const Subscripts, DimSize: IntDyneVec): integer;
var
  Pos         : integer;
  i, j        : integer;
  PriorSizes  : IntDyneVec = nil;
begin
  // allocate space for PriorSizes
  SetLength(PriorSizes, ANumDims);

  // calculate PriorSizes values
  for i := 0 to ANumDims - 2 do
    PriorSizes[i] := 1; // initialize

  for i := ANumDims - 2 downto 0 do
    for j := 0 to i do PriorSizes[i] := PriorSizes[i] * DimSize[j];

  Pos := Subscripts[0] - 1;
  for i := 0 to ANumDims - 2 do
    Pos := Pos + (PriorSizes[i] * (Subscripts[i+1]-1));

  Result := Pos;
  PriorSizes := nil;
end;



{  SUBROUTINE COMBO(ISET, N, M, LAST)

   ALGORITHM AS 160.2  APPL. STATIST. (1981) VOL.30, NO.1

   Subroutine to generate all possible combinations of M of the
   integers from 1 to N in a stepwise fashion.   Prior to the first
   call, LAST should be set to .FALSE.   Thereafter, as long as LAST
   is returned .FALSE., a new valid combination has been generated.
   When LAST goes .TRUE., there are no more combinations.

   LOGICAL LAST
   INTEGER N, M, ISET(M)
}
procedure TLogLinScreenForm.COMBO(var ISET: IntDyneVec; N, M: Integer;
  var LAST: boolean);
label
  100, 110, 130, 150;
var
  I, K, L : integer;
begin
      IF (LAST) then GOTO 110;
//
//     Get next element to increment
//
      K := M;
100:  L := ISET[K] + 1;
      IF (L + M - K <= N) then GOTO 150;
      K := K - 1;
//
//     See if we are done
//
      IF (K <= 0) then GOTO 130;
      GOTO 100;
//
//     Initialize first combination
//
110:  for I := 1 to M do ISET[I] := I;
130:  LAST := NOT LAST;
      exit;
//
//     Fill in remainder of combination.
//
150:  for I := K to M do //DO 160 I = K, M
      begin
	     ISET[I] := L;
	     L := L + 1;
      end; //160   CONTINUE
end;


procedure TLogLinScreenForm.Compute;
var
  ArraySize: integer;
  N: integer;
  index, index2, i, j, k, l, nVars: integer;
  count: integer;
  Data: DblDyneVec = nil;
  Subscripts: IntDyneVec = nil;
  DimSize: IntDyneVec = nil;
  GridPos: IntDyneVec = nil;
  Labels: StrDyneVec = nil;
  Margins: IntDyneMat = nil;
  Expected: DblDyneVec = nil;
  WorkVec: IntDyneVec = nil;
  Indexes: IntDyneMat = nil;
  LogM: DblDyneVec = nil;
  M: DblDyneMat = nil;
  astr, HeadStr : string;
  MaxDim, MP, MM  : integer;
  U, Mu : Double;
  Chi2, G2 : double;
  DF : integer;
  ProbChi2, ProbG2 : double;
  GSQ : DblDyneVec = nil;
  DGFR : IntDyneVec = nil;
  PART : DblDyneMat = nil;
  MARG : DblDyneMat = nil;
  DFS : IntDyneMat = nil;
  IP : IntDyneMat = nil;
  IM : IntDyneMat = nil;
  ISET : IntDyneVec = nil;
  JSET : IntDyneVec = nil;
  CONFIG : IntDyneMat = nil;
  FIT : DblDyneVec = nil;
  SIZE : IntDyneVec = nil;
  COORD : IntDyneVec = nil;
  X: DblDyneVec = nil;
  Y : DblDyneVec = nil;
  IFAULT : integer;
  TABLE : DblDyneVec = nil;
  DIM : IntDyneVec = nil;
  nDims: Integer;
  Minimums: IntDyneVec = nil;
  Maximums: IntDyneVec = nil;
//   Response: BoolDyneVec;
//   Interact: BoolDyneVec;
  lReport: TStrings;
begin
  lReport := TStringList.Create;
  try
    // Allocate space for labels, DimSize and SubScripts
    nDims := MinMaxGrid.RowCount - 1;
    nVars := SelectList.Items.Count;

    SetLength(Labels, nVars);
    SetLength(GridPos, nVars);
    SetLength(Minimums, nDims);
    SetLength(Maximums, nDims);
    SetLength(DimSize, nDims);
    SetLength(Subscripts, nDims);

    for i := 1 to nDims do
    begin
      if not TryStrToInt(MinMaxGrid.Cells[1, i], Minimums[i-1]) then
      begin
        ErrorMsg('Valid Integer number > 0 expected.');
        exit;
      end;
      if not TryStrToInt(MinMaxGrid.Cells[2, i], Maximums[i-1]) then
      begin
        ErrorMsg('Integer number > 0 expected.');
        exit;
      end;
    end;

    // get variable labels and column positions
    for i := 0 to nVars-1 do
    begin
      Labels[i] := SelectList.Items[i];
      GridPos[i] := GetVariableIndex(OS3MainFrm.DataGrid, Labels[i]);
    end;

    // Get no. of categories for each dimension (DimSize)
    MaxDim := 0;
    ArraySize := 1;
    for i := 0 to nDims - 1 do
    begin
      DimSize[i] := Maximums[i] - Minimums[i] + 1;
      if DimSize[i] > MaxDim then MaxDim := DimSize[i];
      ArraySize := ArraySize * DimSize[i];
    end;

    // Allocate space for Data and marginals
    SetLength(WorkVec, MaxDim);
    SetLength(Data, ArraySize);
    SetLength(Margins, nDims, MaxDim);
    SetLength(Expected, ArraySize);
    SetLength(Indexes, ArraySize+1, nDims);
    SetLength(LogM, ArraySize);
    SetLength(M, ArraySize, nDims);

    // Read and store frequencies in Data
    for i := 1 to NoCases do
    begin
      if GoodRecord(OS3MainFrm.DataGrid, i, GridPos) then // casewise check
      begin
        // get cell subscripts
        for j := 0 to nDims-1 do
        begin
          index := StrToInt(OS3MainFrm.DataGrid.Cells[GridPos[j], i]);
          index := index - Minimums[j] + 1;
          Subscripts[j] := index;
        end;

        index := ArrayPosition(nDims, Subscripts, DimSize);

        // save subscripts for later use
        for j := 0 to nDims-1 do
          Indexes[index, j] := Subscripts[j];

        if CountVarChk.Checked then
        begin
          k := GridPos[nVars-1];
          Data[index] := Data[index] + StrToInt(OS3MainFrm.DataGrid.Cells[k, i]);
        end else
          Data[index] := Data[index] + 1;
      end;
    end;

     // get total N
    N := 0;
    for i := 0 to ArraySize-1 do
      N := N + Round(Data[i]);

    // Get marginal frequencies
    Marginals(nDims, ArraySize, Indexes, Data, Margins);

    // Print Marginal totals if requested
    if MarginsChk.Checked then
    begin
      lReport.Add('FILE: '+ OS3MainFrm.FileNameEdit.Text);
      lReport.Add('');
      for i := 0 to nDims-1 do
      begin
        HeadStr := 'Marginal Totals for ' + Labels[i];
        k := DimSize[i];
        for j := 0 to k-1 do WorkVec[j] := Margins[i, j];
        VecPrint(WorkVec, k, HeadStr, lReport);
      end;
    end;
    lReport.Add('');
    lReport.Add('Total Frequencies: %d', [N]);

    lReport.Add('');
    lReport.Add(DIVIDER);
    lReport.Add('');

    // Get Expected cell values
    U := 0.0; // overall mean (mu) of log linear model
    for i := 0 to ArraySize-1 do // indexes point to each cell
    begin
      Expected[i] := 1.0;
      for j := 0 to nDims-1 do
      begin
        k := Indexes[i, j];
        Expected[i] := Expected[i] * (Margins[j, k-1] / N);
      end;
      Expected[i] := Expected[i] * N;
      LogM[i] := ln(Expected[i]);
    end;
    for i := 0 to ArraySize-1 do U := U + LogM[i];
    U := U / ArraySize;

    // print expected values
    lReport.Add('FILE: '+ OS3MainFrm.FileNameEdit.Text);
    lReport.Add('');
    lReport.Add('EXPECTED CELL VALUES FOR MODEL OF COMPLETE INDEPENDENCE');
    lReport.Add('');

    astr := 'Cell';
    for j := 2 to nDims do astr := astr + '    ';
    lReport.Add(astr + 'Observed  Expected  Log Expected');

    astr := '';
    for j := 1 to nDims do astr := astr + '--- ';
    astr := astr + '---------- ---------- ----------';
    lReport.Add(astr);

    for i := 1 to ArraySize do
    begin
      astr := '';
      for j := 1 to nDims do
        astr := astr + Format('%3d ', [Indexes[i-1, j-1]]);
      astr := astr + Format('%10.0f %10.2f %10.3f',[Data[i-1], Expected[i-1], LogM[i-1]]);
      lReport.Add(astr);
    end;

    // Calculate chi-squared and G squared statistics
    chi2 := 0.0;
    G2 := 0.0;
    for i := 0 to ArraySize-1 do
    begin
      chi2 := chi2 + Sqr(Data[i] - Expected[i]) / Expected[i];
      G2 := G2 + Data[i] * ln(Data[i] / Expected[i]);
    end;
    G2 := 2.0 * G2;
    DF := 1;
    for i := 0 to nDims-1 do
      DF := DF * (DimSize[i] - 1);
    ProbChi2 := 1.0 - ChiSquaredProb(chi2,DF);
    ProbG2 := 1.0 - ChiSquaredProb(G2,DF);
    lReport.Add('Chisquare: %10.3f with probability %10.3f (DF %d)', [chi2, ProbChi2, DF]);
    lReport.Add('G squared: %10.3f with probability %10.3f (DF %d)', [G2, ProbG2, DF]);
    lReport.Add('');
    lReport.Add('U (mu) for general loglinear model: %10.2f', [U]);

    lReport.Add('');
    lReport.Add(DIVIDER);
    lReport.Add('');

    // Get log linear model values for each cell
    // get M's for each cell
    lReport.Add('First Order LogLinear Model Factors and N of Cells in Each');
    astr := 'CELL              ';
    for i := 1 to nDims do
      astr := astr + Format(' U%d  N Cells   ',[i]);
    lReport.Add(astr);
    lReport.Add('');
    for i := 1 to ArraySize do // cell
    begin
      astr := '';
      for j := 1 to nDims do
        astr := astr + Format('%3d ', [Indexes[i-1,j-1]]);
      for j := 1 to nDims do  // jth mu
      begin
        index := Indexes[i-1,j-1]; // sum for this mu
        count := 0;
        Mu := 0.0;
        for k := 1 to ArraySize do
        begin
          if index = Indexes[k-1,j-1] then
          begin
            count := count + 1;
            Mu := Mu + LogM[k-1];
          end;
        end;
        Mu := Mu / count - U;
        astr := astr + format('%10.3f %3d ',[Mu,count]);
      end;
      lReport.Add(astr);
    end;

    lReport.Add('');
    lReport.Add(DIVIDER);
    lReport.Add('');

    // get second order interactions
    lReport.Add('Second Order Loglinear Model Terms and N of Cells in Each');
    astr := 'CELL              ';
    for i := 1 to nDims-1 do
      for j := i + 1 to nDims do
        astr := astr + format('U%d%d  N Cells  ',[i,j]);
    lReport.Add(astr);
    lReport.Add('');
    for i := 1 to ArraySize do // cell
    begin
      astr := '';
      for j := 0 to nDims-1 do
        astr := astr + Format('%3d ', [Indexes[i, j]]);
      for j := 1 to nDims-1 do  // jth
      begin
        index := Indexes[i-1,j-1]; // sum for this mu using j and k
        for k := j+1 to nDims do // with kth
        begin
          index2 := Indexes[i-1,k-1];
          Mu := 0.0;
          count := 0;
          for l := 1 to ArraySize do
          begin
            if ((index = Indexes[l-1,j-1]) and (index2 = Indexes[l-1,k-1])) then
            begin
              Mu := Mu + LogM[l-1];
              count := count + 1;
            end;
          end; // next l
          Mu := Mu / count - U;
          astr := astr + Format('%10.3f %3d', [Mu, count]);
        end; // next k (second term subscript)
      end; // next j (first term subscript)
      lReport.Add(astr);
    end; // next i

    lReport.Add('');
    lReport.Add(DIVIDER);
    lReport.Add('');

    // get maximum no. of interactions in saturated model
    MaxCombos(nDims, MM, MP);

    SetLength(GSQ, nDims+1);
    SetLength(DGFR, nDims+1);
    SetLength(PART, nDims+1, MP+1);
    SetLength(MARG, nDims+1, MP+1);
    SetLength(DFS, nDims+1, MP+1);
    SetLength(IP, nDims+1, MP+1);
    SetLength(IM, nDims+1, MM+1);
    SetLength(ISET, nDims+1);
    SetLength(JSET, nDims+1);
    SetLength(CONFIG, nDims+1, MP+1);
    SetLength(FIT, ArraySize+1);
    SetLength(SIZE, nDims+1);
    SetLength(COORD, nDims+1);
    SetLength(X, ArraySize+1);
    SetLength(Y, ArraySize+1);
    SetLength(TABLE, ArraySize+1);
    SetLength(DIM, nDims+1);

     // Load TABLE and DIM one up from Data
     for i := 1 to ArraySize do Table[i] := Data[i-1];
     for i := 1 to nDims do DIM[i] := DimSize[i-1];

     Screen(
       nDims, MP, MM, ArraySize, TABLE, DIM,
       GSQ, DGFR, PART, MARG, DFS, IP, IM, ISET, JSET, CONFIG, FIT, SIZE,
       COORD, X, Y, IFAULT
     );

    // show results
    lReport.Add('SCREEN FOR INTERACTIONS AMONG THE VARIABLES');
    lReport.Add('Adapted from the Fortran program by Lustbader and Stodola printed in');
    lReport.Add('Applied Statistics, Volume 30, Issue 1, 1981, pages 97-105 as Algorithm');
    lReport.Add('AS 160 Partial and Marginal Association in Multidimensional Contingency Tables');
    lReport.Add('');
    lReport.Add('Statistics for tests that the interactions of a given order are zero');
    lReport.Add('ORDER     STATISTIC      D.F.       PROB.');
    lReport.Add('-----     ----------     ----     ----------');
    for i := 1 to nDims do
    begin
      ProbChi2 := 1.0 - ChiSquaredProb(GSQ[i], DGFR[i]);
      lReport.Add('%5d     %10.3f     %4d     %10.3f',[i,GSQ[i], DGFR[i], ProbChi2]);
    end;
    lReport.Add('');
    lReport.Add('Statistics for Marginal Association Tests');
    lReport.Add('VARIABLE     ASSOC.    PART ASSOC.   MARGINAL ASSOC.   D.F.   PROB');
    lReport.Add('--------     ------    -----------   ---------------   ----   ----------');
    for i := 1 to nDims-1 do
    begin
      for j := 1 to MP do
      begin
        ProbChi2 := 1.0 - ChiSquaredProb(MARG[i,j],DFS[i,j]);
        lReport.Add('%6d        %5d     %10.3f    %12.3f      %3d    %10.3f',
          [i, j, PART[i,j], MARG[i,j], DFS[i,j], ProbChi2]);
      end;
    end;

    FReportFrame.DisplayReport(lReport);

  finally
    lReport.Free;
  end;
end;


{  SUBROUTINE CONF(N, M, MP, MM, ISET, JSET, IP, IM, NP)

    ALGORITHM AS 160.1 APPL. STATIST. (1981) VOL.30, NO.1

    Set up the arrays IP and IM for a given N and M.   Essentially
    IP contains all possible combinations of (N choose M).   For each
    combination found IM contains all combinations of degree M-1.

    INTEGER ISET(N), JSET(N), IP(N,MP), IM(N,MM)
    LOGICAL ILAST, JLAST
}
procedure TLogLinScreenForm.CONF(var N: integer; var M: integer;
  var MP: integer; var MM: integer; var ISET: IntDyneVec; var JSET: IntDyneVec;
  var IP: IntDyneMat; var IM: IntDyneMat; var NP: integer);
label
  100, 120;
var
  ILAST, JLAST: boolean;
  I, L, NM, JS: integer;
begin
      ILAST := TRUE;
      NP := 0;
      NM := 0;
      //
      //     Get IP
      //
  100:
      COMBO(ISET, N, M, ILAST);
      IF (ILAST) then exit;
      NP := NP + 1;
      for I := 1 to M do IP[I,NP] := ISET[I];
      IF (M = 1) then GOTO 100;
//
//     Get IM
//
      JLAST := TRUE;
      L := M - 1;
  120:
      COMBO(JSET, M, L, JLAST);
      IF (JLAST) then GOTO 100;
      NM := NM + 1;
      for I := 1 to L do //      DO 130 I = 1, L
      begin
	     JS := JSET[I];
	     IM[I,NM] := ISET[JS];
      end; // 130 CONTINUE
      GOTO 120;
end;


procedure TLogLinScreenForm.CountVarChkChange(Sender: TObject);
begin
  UpdateMinMaxGrid;
end;


procedure TLogLinScreenForm.EVAL(var IAR: IntDyneMat; NC, NV, IBEG, NVAR,
  MAX: integer; var CONFIG: IntDyneMat; var DIM: IntDyneVec; var DF: integer);
VAR I, J, K, KK, L : integer;
//      SUBROUTINE EVAL(IAR, NC, NV, IBEG, NVAR, MAX, CONFIG, DIM, DF)
//
//     ALGORITHM AS 160.3 APPL. STATIST. (1981) VOL.30, NO.1
//
//     IAR  = array containing the effects to be fitted
//     NC   = number of columns of IAR to be used
//     NV   = number of variables in each effect
//     IBEG = gebinning column
//     DF   = degrees of freedom
//
//     CONFIG is in a format compatible with algorithm AS 51
//
//      INTEGER IAR(NVAR,MAX), CONFIG(NVAR,NC), DIM(NVAR), DF
//
begin
      DF := 0;
      for J := 1 to NC do //DO 110 J = 1, NC
      begin
	     KK := 1;
	     for I := 1 to NV do //DO 100 I = 1, NV
             begin
	          L := IBEG + J - 1;
	          K := IAR[I,L];
	          KK := KK * (DIM[K] - 1);
	          CONFIG[I,J] := K;
             end; // 100   CONTINUE
	     CONFIG[NV+1,J] := 0;
	     DF := DF + KK;
      end; // 110 CONTINUE
end;


procedure TLogLinScreenForm.InBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while i < VarList.Items.Count do
  begin
    if VarList.Selected[i] then
    begin
      SelectList.Items.Add(VarList.Items[i]);
      VarList.Items.Delete(i);
      i := 0;
    end else
      i := i + 1;
  end;
  UpdateMinMaxGrid;
  UpdateBtnStates;
end;


{  SUBROUTINE LIKE(GSQ, FIT, TABLE, NTAB)
   ALGORITHM AS 160.5 APPL. STATIST. (1981) VOL.30, NO.1

   Compute the likelihood-ration chi-square

   REAL FIT(NTAB), TABLE(NTAB), ZERO, TWO
     DATA ZERO /0.0/, TWO /2.0/
}
procedure TLogLinScreenForm.LIKE(var GSQ: Double; const FIT, TABLE: DblDyneVec;
  NTAB: integer);
const
  ZERO = 0.0;
  TWO = 2.0;
var
  I: integer;
begin
  GSQ := ZERO;
  for I := 1 to NTAB do //DO 100 I = 1, NTAB
  begin
    if (FIT[I] = ZERO) or (TABLE[I] = ZERO) then continue; // GO TO 100
    GSQ := GSQ + TABLE[I] * Ln(TABLE[I] / FIT[I]);
  end; // 100 CONTINUE
  GSQ := TWO * GSQ;
end;


{
  SUBROUTINE LOGFIT(NVAR, NTAB, NCON, DIM, CONFIG, TABLE, FIT, SIZE, COORD, X, Y)

  ALGORITHM AS 160.6 APPL. STATIST. (1981) VOL.30, NO.1

  Iterative proportional fitting of the marginals of a contingency table.
  Relevant code from AS 51 is used.

  REAL TABLE(NTAB), FIT(NTAB), MAXDEV, X(NTAB), Y(NTAB), ZERO
  INTEGER CONFIG(NVAR,NCON), DIM(NVAR), SIZE(NVAR), COORD(NVAR)
  LOGICAL OPTION
  DATA MAXDEV /0.25/, MAXIT /25/, ZERO /0.0/
}
procedure TLogLinScreenForm.LOGFIT(NVAR, NTAB, NCON: integer;
  var DIM: IntDyneVec; var CONFIG: IntDyneMat; var TABLE: DblDyneVec;
  var FIT: DblDyneVec; var SIZE: IntDyneVec; var COORD: IntDyneVec;
  var X: DblDyneVec; var Y: DblDyneVec);
label
  110, 130, 150, 170, 180, 200;
const
  ZERO = 0.0;
  MAXDEV = 0.25;
  MAXIT = 25;
var
  II, K, KK, L, N, J, I: integer;
  OPTION: boolean;
  XMAX, E: double;
  NV1, ISZ: integer;
begin
      for KK := 1 to MAXIT do //DO 230 KK = 1, MAXIT
      begin
           //
           //     XMAX is the maximum deviation between fitted and true marginal
           //
	     XMAX := ZERO;
	     for II := 1 to NCON do //DO 220 II = 1, NCON
             begin
	          OPTION := TRUE;
                  //
                  //     Initialize arrays
                  //
	          SIZE[1] := 1;
	          NV1 := NVAR - 1;
	          for K := 1 to NV1 do //DO 100 K = 1, NV1
                  begin
	               L := CONFIG[K,II];
	               IF (L = 0) then GOTO 110;
	               SIZE[K+1] := SIZE[K] * DIM[L];
                  end; // 100     CONTINUE
	          K := NVAR;
  110:            N := K - 1;
	          ISZ := SIZE[K];
	          for J := 1 to ISZ do //DO 120 J = 1, ISZ
                  begin
	               X[J] := ZERO;
	               Y[J] := ZERO;
                  end; // 120     CONTINUE
                  //
                  //     Initialize co-ordinates
                  //
  130:            for K := 1 to NVAR do COORD[K] := 0;
                  //
                  //     Find locations in tables
                  //
                  I := 1;
  150:            J := 1;
                  for K := 1 to N do //DO 160 K = 1, N
                  begin
	               L := CONFIG[K,II];
	               J := J + COORD[L] * SIZE[K];
                  end; //160     CONTINUE
                  IF (NOT OPTION) then GOTO 170;
                  //
                  //     Compute marginals
                  //
	          X[J] := X[J] + TABLE[I];
	          Y[J] := Y[J] + FIT[I];
	          GOTO 180;
                  //
                  //     Make adjustments
                  //
  170:            IF (Y[J] <= ZERO) then FIT[I] := ZERO;
                  IF (Y[J] > ZERO) then FIT[I] := FIT[I] * X[J] / Y[J];
                  //
                  //     Update co-ordinates
                  //
  180:            I := I + 1;
	          for K := 1 to NVAR do //DO 190 K = 1, NVAR
                  begin
	               COORD[K] := COORD[K] + 1;
	               IF (COORD[K] < DIM[K]) then GOTO 150;
	               COORD[K] := 0;
                  end; //190     CONTINUE
	          IF (NOT OPTION) then GOTO 200;
	          OPTION := FALSE;
	          GOTO 130;
                  //
                  //     Find the largest deviation
                  //
  200:            for I := 1 to ISZ do //DO 210 I = 1, ISZ
                  begin
	               E := ABS(X[I] - Y[I]);
	               IF (E > XMAX) then XMAX := E;
                  end; // 210     CONTINUE
             end; // 220   CONTINUE
             //
             //     Test convergence
             //
	     IF (XMAX < MAXDEV) then exit;
      end; // 230 CONTINUE
end;


procedure TLogLinScreenForm.Marginals(NumDims, ArraySize: integer;
  const Indexes: IntDyneMat; const Data: DblDyneVec; const Margins: IntDyneMat);
var
  i, j, category: integer;
begin
  for i := 1 to ArraySize do
  begin
    for j := 1 to NumDims do
    begin
      category := Indexes[i-1,j-1];
      Margins[j-1,category-1] := Margins[j-1,category-1] + Round(Data[i-1]);
    end;
  end;
end;


procedure TLogLinScreenForm.MaxCombos(NumDims: integer; out MM, MP: integer);
var
  combos: integer;
  i,j: integer;
begin
  MM := 0;
  MP := 0;
  for i := 1 to NumDims do
  begin
    combos := 1;

    // get numerator factorial products down to i
    for j := NumDims downto i + 1 do
      combos := combos * j;

    // divide by factorial of NumDims - i;
    for j := (NumDims - i) downto 2 do
      combos := combos div j;

    if combos > MP then
      MP := combos;
    if i * combos > MM then
      MM := i * combos;
  end;
end;


procedure TLogLinScreenForm.OutBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while i < SelectList.Items.Count do
  begin
    if SelectList.Selected[i] then
    begin
      VarList.Items.Add(SelectList.Items[i]);
      SelectList.Items.Delete(i);
      i := 0;
    end else
      i := i + 1;
  end;
  UpdateBtnStates;
  UpdateMinMaxGrid;
end;


{  SUBROUTINE RESET(FIT, NTAB, AVG)

   ALGORITHM AS 160.4 APPL. STATIST. (1981) VOL.30, NO.1

   Initialize the fitted values to the average entry

    REAL FIT(NTAB)
}
procedure TLogLinScreenForm.RESET(var FIT: DblDyneVec; NTAB: Integer; AVG: Double);
var
  I: integer;
begin
  for I := 1 to NTAB do //DO 100 I = 1, NTAB
  begin
    FIT[I] := AVG;
  end; // 100 CONTINUE
end;


procedure TLogLinScreenForm.Reset;
begin
  inherited;
  CollectVariableNames(OS3MainFrm.DataGrid, VarList.Items);
  SelectList.Clear;
  UpdateBtnStates;
  UpdateMinMaxGrid;
end;


{  SUBROUTINE SCREEN(NVAR, MP, MM, NTAB, TABLE, DIM, GSQ, DGFR,
   *   PART, MARG, DFS, IP, IM, ISET, JSET, CONFIG, FIT, SIZE,
   *   COORD, X, Y, IFAULT)

   ALGORITHM AS 160 APPL. STATIST. (1981) VOL.30, NO.1

   Screen all efects for partial and marginal association.

   INTEGER NVAR, MP, MM, NTAB, IP(NVAR,MP), IM(NVAR,MM), DGFR(NVAR),
    *   DFS(NVAR,MP), ISET(NVAR), JSET(NVAR), CONFIG(NVAR,MP),
    *   DIM(NVAR), DF, SIZE(NVAR), COORD(NVAR)
   REAL GSQ(NVAR), PART(NVAR,MP), MARG(NVAR,MP), TABLE(NTAB),
    *   FIT(NTAB), X(NTAB), Y(NTAB), ZERO
   DATA ZERO /0.0/
}
procedure TLogLinScreenForm.Screen(var NVAR: integer; var MP: integer;
  var MM: integer; var NTAB: integer; var TABLE: DblDyneVec;
  var DIM: IntDyneVec; var GSQ: DblDyneVec; var DGFR: IntDyneVec;
  var PART: DblDyneMat; var MARG: DblDyneMat; var DFS: IntDyneMat;
  var IP: IntDyneMat; var IM: IntDyneMat; var ISET: IntDyneVec;
  var JSET: IntDyneVec; var CONFIG: IntDyneMat; var FIT: DblDyneVec;
  var SIZE: IntDyneVec; var COORD: IntDyneVec; var X: DblDyneVec;
  var Y: DblDyneVec; var IFAULT: integer);
Label
  160, 170;
const
  ZERO = 0.0;
var
  ISZ, MAX, LIM, I, J, NV1, M, M1, ITP, NP, NP1, L3, DF: integer;
  G21, G22, G23, AVG: double;

begin
//
//     Check for input errors
//
      IFAULT := 1;
      IF (NVAR <= 1) then exit;
      ISZ := 1;
      for I := 1 to NVAR do
      begin
           if (DIM[I] <= 1) then IFAULT := 2;
           ISZ := ISZ * DIM[i];
      end;
      IF (ISZ <> NTAB) then IFAULT := 2;
      MAX := 1;
      LIM := NVAR div 2;
      for I := 1 to LIM do MAX := MAX * (NVAR - I + 1) div I;
      IF (MP < MAX) then IFAULT := 3;
      MAX := 1;
      LIM := (NVAR - 1) div 2;
      for I := 1 to LIM do MAX := MAX * (NVAR - I) div I;
      MAX := MAX * NVAR;
      IF (MM < MAX) then IFAULT := 4;
      IF (IFAULT > 1) then exit;
//
//     Fit the no effect model
//
      DGFR[NVAR] := NTAB - 1;
      AVG := ZERO;
      IFAULT := 5;
      for I := 1 to NTAB do
      begin
	       IF (TABLE[I] < ZERO) then exit; //RETURN
	       AVG := AVG + TABLE[I];
      end;
      IFAULT := 0;
      AVG := AVG / NTAB;
      RESET(FIT, NTAB, AVG);
      LIKE(GSQ[1], FIT, TABLE, NTAB);
//
//     Begin fitting effects
//
      NV1 := NVAR - 1;
      for M := 1 to NV1 do
      begin
          //      DO 200 M = 1, NV1
          //
          //     Set up the arrays IP and IM
          //
          M1 := M;
          CONF(NVAR, M1, MP, MM, ISET, JSET, IP, IM, NP);
          //
          //    Fit the saturated model
          //
          RESET(FIT, NTAB, AVG);
          EVAL(IP, NP, M, 1, NVAR, MP, CONFIG, DIM, DGFR[M]);
          LOGFIT(NVAR, NTAB, NP, DIM, CONFIG, TABLE, FIT, SIZE, COORD, X, Y);
          LIKE(GSQ[M+1], FIT, TABLE, NTAB);
          //
          //     Move the first column of IP to the last
          //
          for I := 1 to M do
          begin
               // DO 150 I = 1, M
	       ITP := IP[I,1];
	       NP1 := NP - 1;
           for J := 1 to NP1 do IP[I,J] := IP[I,J+1];
	       IP[I,NP] := ITP;
          end; //  150   CONTINUE
          L3 := -M + 1;
          for J := 1 to NP do
          begin
              //  DO 190 J = 1, NP
              //
              //  Fit the effects in IP ignoring the last column
              //
	      RESET(FIT, NTAB, AVG);
	      EVAL(IP, NP-1, M, 1, NVAR, MP, CONFIG, DIM, DF);
	      LOGFIT(NVAR, NTAB, NP-1, DIM, CONFIG, TABLE, FIT, SIZE, COORD, X, Y);
	      LIKE(G21, FIT, TABLE, NTAB);
	      DFS[M,J] := DGFR[M] - DF;
	      PART[M,J] := G21 - GSQ[M+1];
              //
              //     For M = 1, partials and marginals are equal
              //
  	      IF (M > 1) then GOTO 160;
	      MARG[1,J] := PART[1,J];
	      GOTO 170;
              //
              //     Fit the last column alone
              //
  160:        RESET(FIT, NTAB, AVG);
	      EVAL(IP, 1, M, NP, NVAR, MP, CONFIG, DIM, DF);
	      LOGFIT(NVAR, NTAB, 1, DIM, CONFIG, TABLE, FIT, SIZE,
                     COORD, X, Y);
	      LIKE(G22, FIT, TABLE, NTAB);
              //
              //     Locate the appropriate columns of IM and fit them
              //
	      L3 := L3 + M;
	      RESET(FIT, NTAB, AVG);
	      EVAL(IM, M, M-1, L3, NVAR, MM, CONFIG, DIM, DF);
	      LOGFIT(NVAR, NTAB, M, DIM, CONFIG, TABLE, FIT, SIZE,
                     COORD, X, Y);
	      LIKE(G23, FIT, TABLE, NTAB);
	      MARG[M,J] := G23 - G22;
              //
              //     Move the next effect to be ignored to the last in IP
              //
  170:        for I := 1 to M do // DO 180 I = 1, M
              begin
	           ITP := IP[I,NP];
	           IP[I,NP] := IP[I,J];
	           IP[I,J] := ITP;
              end;
              //  180     CONTINUE
          end; //  190   CONTINUE
          //
          DGFR[NVAR] := DGFR[NVAR] - DGFR[M];
          GSQ[M] := GSQ[M] - GSQ[M+1];
      end; // 200 CONTINUE
end;


procedure TLogLinScreenForm.SelectListSelectionChange(Sender: TObject;
  User: boolean);
begin
  UpdateBtnStates;
end;


procedure TLogLinScreenForm.UpdateBtnStates;
begin
  inherited;
  
  InBtn.Enabled := AnySelected(VarList);
  OutBtn.Enabled := AnySelected(SelectList);
  AllBtn.Enabled := VarList.Items.Count > 0;
end;


procedure TLogLinScreenForm.UpdateMinMaxGrid;
var
  NumDims, j: Integer;
begin
  if CountVarChk.Checked then
    NumDims := SelectList.Count - 1
  else
    NumDims := SelectList.Count;
  if NumDims < 0 then NumDims := 0;
  MinMaxGrid.RowCount := NumDims + 1;
  for j := 1 to MinMaxGrid.RowCount-1 do
  begin
    MinMaxGrid.Cells[0, j] := SelectList.Items[j-1];
    MinMaxGrid.Cells[1, j] := '';
    MinMaxGrid.Cells[2, j] := '';
  end;
end;


function TLogLinScreenForm.Validate(out AMsg: String;
  out AControl: TWinControl): Boolean;
begin
  Result := false;

  if SelectList.Count = 0 then
  begin
    AMsg := 'No variables selected.';
    AControl := VarList;
    exit;
  end;

  Result := true;
end;


procedure TLogLinScreenForm.VarListDblClick(Sender: TObject);
var
  index: Integer;
begin
  index := VarList.ItemIndex;
  if index > -1 then begin
    SelectList.Items.Add(VarList.Items[index]);
    VarList.Items.Delete(index);
    UpdateBtnStates;
  end;
end;


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


end.