unit WLSUnit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
  StdCtrls, Buttons, ExtCtrls, ComCtrls,
  Globals, MainUnit, DictionaryUnit, Matrixlib, DataProcs,
  MathUnit, ReportFrameUnit, ChartFrameUnit, BasicStatsParamsFormUnit;

type

  { TWLSFrm }

  TWLSFrm = class(TBasicStatsParamsForm)
    OptionsBevel: TBevel;
    DepInBtn: TBitBtn;
    DepOutBtn: TBitBtn;
    IndInBtn: TBitBtn;
    IndOutBtn: TBitBtn;
    PageControl: TPageControl;
    OLSPage: TTabSheet;
    ResRegPageControl: TPageControl;
    ResidualsRegPage: TTabSheet;
    ResRegReportPage: TTabSheet;
    UserWeightsChk: TRadioButton;
    WeightChk: TRadioButton;
    WLSPage: TTabSheet;
    WeightInBtn: TBitBtn;
    WeightOutBtn: TBitBtn;
    SaveChk: TCheckBox;
    OriginChk: TCheckBox;
    Origin2Chk: TCheckBox;
    DepVarEdit: TEdit;
    WeightVarEdit: TEdit;
    OptionsGroup: TGroupBox;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    IndVarList: TListBox;
    VarList: TListBox;
    procedure DepInBtnClick(Sender: TObject);
    procedure DepOutBtnClick(Sender: TObject);
    procedure IndInBtnClick(Sender: TObject);
    procedure IndOutBtnClick(Sender: TObject);
    procedure IndVarListDblClick(Sender: TObject);
    procedure UserWeightsChkChange(Sender: TObject);
    procedure VarListDblClick(Sender: TObject);
    procedure VarListSelectionChange(Sender: TObject; {%H-}User: boolean);
    procedure WeightInBtnClick(Sender: TObject);
    procedure WeightOutBtnClick(Sender: TObject);

  private
    OLSReportFrame: TReportFrame;
    ResidualsRegReportFrame: TReportFrame;
    WLSReportFrame: TReportFrame;

    procedure AddPredictedStuffToGrid(AIndepCols: IntDyneVec; ANumIndepCols: Integer;
      BWeights: DblDyneVec);

    procedure CreateOrGetChartFrame(AColIndex: Integer; AVarName: String;
      out AMemo: TMemo; out AChartFrame: TChartFrame);

    function GetPageCaption(AVarName: String): String;

    procedure PlotSquaredResiduals(AIndepCols: IntDyneVec;
      ANumIndepCols, ADepCol: Integer; AConfLevel: Double);

    procedure PlotXY(AChartFrame: TChartFrame; const XPoints, YPoints: DblDyneVec;
      const ARegressionResults: TBivariateRegressionResults; const XLabel, YLabel: String);

    procedure PredictIt(ColNoSelected: IntDyneVec; NoVars: integer;
      Means, StdDevs, BetaWeights: DblDyneVec;
      StdErrEst: double; NoIndepVars: integer);

    function PrepareData(out ADepCol, ANumIndepCols: Integer;
      out AIndepCols: IntDyneVec; out AWeightCol: Integer;
      out ARowLabels: StrDyneVec): Boolean;

    function Process_OLSRegression(AIndepCols: IntDyneVec; ANumIndepCols, ADepCol: Integer;
      ARowLabels: StrDyneVec; ANumCases: Integer; PrintAll: Boolean): Boolean;

    function Process_SquaredResidualsRegression(AIndepCols: IntDyneVec;
      ANumIndepCols, ADepCol: Integer; ARowLabels: StrDyneVec;
      BWeights: DblDyneVec; ANumCases: Integer; PrintAll: Boolean): Boolean;

    procedure WriteDescriptiveReport(AMemo: TMemo;
      const ARegressionResults: TBivariateRegressionResults;
      const XLabel, YLabel: String);

  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
  WLSFrm: TWLSFrm;


implementation

{$R *.lfm}

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

{ TWLSFrm }

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

  OLSReportFrame := TReportFrame.Create(self);
  OLSReportFrame.Name := '';
  OLSReportFrame.Parent := OLSPage;
  OLSReportFrame.Align := alClient;
  OLSReportFrame.BorderSpacing.Left := 0;
  OLSReportFrame.BorderSpacing.Top := 0;
  OLSReportFrame.BorderSpacing.Bottom := 0;
  OLSReportFrame.BorderSpacing.Right := 0;
  InitToolbar(OLSReportFrame.ReportToolbar, tpRight);

  ResidualsRegReportFrame := TReportFrame.Create(self);
  ResidualsRegReportFrame.Name := '';
  ResidualsRegReportFrame.Parent := ResRegReportPage;
  ResidualsRegReportFrame.Align := alClient;
  ResidualsRegReportFrame.BorderSpacing.Left := 0;
  ResidualsRegReportFrame.BorderSpacing.Top := 0;
  ResidualsRegReportFrame.BorderSpacing.Bottom := 0;
  ResidualsRegReportFrame.BorderSpacing.Right := 0;
  InitToolbar(ResidualsRegReportFrame.ReportToolbar, tpRight);

  WLSReportFrame := TReportFrame.Create(self);
  WLSReportFrame.Name := '';
  WLSReportFrame.Parent := WLSPage;
  WLSReportFrame.Align := alClient;
  WLSReportFrame.BorderSpacing.Left := 0;
  WLSReportFrame.BorderSpacing.Top := 0;
  WLSReportFrame.BorderSpacing.Bottom := 0;
  WLSReportFrame.BorderSpacing.Right := 0;
  InitToolbar(WLSReportFrame.ReportToolbar, tpRight);
end;


{ Get predicted squared residuals and save recipricols to grid as weights }
procedure TWLSFrm.AddPredictedStuffToGrid(AIndepCols: IntDyneVec;
  ANumIndepCols: Integer; BWeights: DblDyneVec);
var
  col: Integer;
  i, j: Integer;
  X, predicted: Double;
begin
  col := NoVariables + 1;
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.NewVar(col);
  DictionaryFrm.DictGrid.Cells[1, col] := 'PredResid2';
  OS3MainFrm.DataGrid.Cells[col, 0] := 'PredResid2';

  col := NoVariables + 1;
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.NewVar(col);
  DictionaryFrm.DictGrid.Cells[1, col] := 'WEIGHT';
  OS3MainFrm.DataGrid.Cells[col, 0] := 'WEIGHT';
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);

  for i := 1 to NoCases do
  begin
    if (DataProcs.ValidValue(i, col-2)) then // do we have a valid squared OLS residual?
    begin
      predicted := 0.0;
      for j := 0 to ANumIndepCols - 1 do
      begin
        X := StrToFloat(OS3MainFrm.DataGrid.Cells[AIndepCols[j], i]);
        predicted := predicted + BWeights[j] * X;
      end;
      predicted := predicted + BWeights[ANumIndepCols];
      predicted := abs(predicted);
      OS3MainFrm.DataGrid.Cells[col-1, i] := Format('%.3f', [predicted]);
      if (predicted > 0.0) then
        predicted := 1.0 / sqrt(predicted)
      else
        predicted := 0.0;
      OS3MainFrm.DataGrid.Cells[col, i] := Format('%.3f', [predicted]);
    end;
  end;
end;


procedure TWLSFrm.AdjustConstraints;
begin
  ParamsPanel.Constraints.MinHeight := DepInBtn.Top + (IndOutBtn.Top - DepInBtn.Top)*2 + DepInBtn.Top +
    OptionsGroup.Height + ButtonBevel.Height + CloseBtn.BorderSpacing.Top + CloseBtn.Height;

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

procedure TWLSFrm.Compute;
var
  i, j, noIndep, depCol, weightCol, oldDepCol, NCases, pos, col: integer;
  IndepCols: IntDyneVec = nil;
  RowLabels: StrDyneVec = nil;
  X, Y: double;
  Means: DblDyneVec = nil;
  Variances: DblDyneVec = nil;
  StdDevs: DblDyneVec = nil;
  BWeights: DblDyneVec = nil;
  BetaWeights: DblDyneVec = nil;
  BStdErrs: DblDyneVec = nil;
  BtTests: DblDyneVec = nil;
  tProbs: DblDyneVec = nil;
  predicted: Double;
  lReport: TStrings;
  StdErrEst: Double = 0.0;
  R2: Double = 0.0;
  errorcode: Boolean = false;
  PrintDesc: boolean = true;
begin
  SetLength(Means, NoVariables + 2);
  SetLength(Variances, NoVariables + 2);
  SetLength(StdDevs, NoVariables + 2);
  SetLength(BWeights, NoVariables + 2);      // do not remove!
  SetLength(BetaWeights, NoVariables + 2);
  SetLength(BStdErrs, NoVariables + 2);
  SetLength(Bttests, NoVariables + 2);
  SetLength(tprobs, NoVariables + 2);

  lReport := TStringList.Create;
  try
    NCases := NoCases;

    // Get column indexes and do some validation checks.
    // NOTE that the Length(indepCols) is different from NoIndep.
    if not PrepareData(depCol, noIndep, indepCols, weightCol, RowLabels) then
      exit;

    // Save dependent column so we can re-use DepCol
    oldDepCol := depCol;

    // *** Get OLS regression ***
    Process_OLSRegression(indepCols, noIndep, depCol, RowLabels, nCases, printDesc);

    // *** Regress the squared residuals on the predictors ***
    depCol := NoVariables;
    Process_SquaredResidualsRegression(indepCols, noIndep, depCol, RowLabels,
      BWeights, nCases, printDesc);
    if WeightChk.Checked then
      AddPredictedStuffToGrid(indepCols, noIndep, BWeights);

    // *** Display squared residuals for each independent variable ***
    // NOTE: depCol points to the squared residuals column here
    PlotSquaredResiduals(IndepCols, NoIndep, depCol, 0.95);

    if WeightChk.Checked then
    begin
      // Weight variables and do OLS regression on weighted variables
      DepCol := olddepcol;
      IndepCols[Noindep] := DepCol;
      for i := 1 to NoCases do
      begin
        Y := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[NoVariables,i])); // weight
        for j := 0 to Noindep do
        begin
          pos := IndepCols[j];
          X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
          X := X * Y;
          OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X);        // wp: DON'T OVERWRITE GRID CELLS
        end;
      end;

      // get means of variables and subtract from the values
      if OriginChk.Checked then
      begin
        for j := 0 to NoIndep do
        begin
          Means[j] := 0.0;
          NCases := 0;
          pos := IndepCols[j];
          for i := 1 to NoCases do
          begin
            if (DataProcs.ValidValue(i,pos)) then
            begin
              Means[j] := Means[j] + StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
              NCases := NCases + 1;
            end;
          end;
          Means[j] := Means[j] / NCases;
          for i := 1 to NoCases do
          begin
            if (DataProcs.ValidValue(i,pos)) then
            begin
              X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
              X := X - Means[j];
              OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X);            // wp: DON'T OVERWRITE GRID DATA!
            end;
          end;  // next i
        end;  // next j
      end;  // if origin checked

      lReport.Clear;
      lReport.Add('WEIGHTED LEAST SQUARES (WLS) REGRESSION RESULTS');
      lReport.Add('');
      MReg(Noindep, IndepCols, DepCol, RowLabels, Means, Variances, StdDevs,
         BWeights, BetaWeights, BStdErrs, Bttests, tprobs, R2, stdErrEst,
         NCases, errorcode, PrintDesc, lReport);

      WLSReportFrame.DisplayReport(lReport);
      lReport.Clear;
    end; // if weightschk checked

    // use the weights entered by the user
    if UserWeightsChk.Checked then
    begin
      // Weight variables and do OLS regression on weighted variables
      depCol := olddepcol;
      indepCols[Noindep] := depCol;                                       // wp:  CALCULATION SHOULD NORMALIZE USER WEIGHTS HERE !!!
      for i := 1 to NoCases do
      begin
        Y := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[weightCol,i])); // weight
        for j := 0 to Noindep do
        begin
          pos := indepCols[j];
          X := StrToFloat(OS3MainFrm.DataGrid.Cells[pos,i]);
          X := X * Y;
          OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X);
        end;
      end;
      if (Origin2Chk.Checked) then // get means of variables and subtract from the values
      begin
        for j := 0 to Noindep do
        begin
          Means[j] := 0.0;
          NCases := 0;
          pos := IndepCols[j];
          for i := 1 to NoCases do
          begin
            if (DataProcs.ValidValue(i,pos)) then
            begin
              Means[j] := Means[j] + StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
              NCases := NCases + 1;
            end;
          end;
          Means[j] := Means[j] / NCases;
          for i := 1 to NoCases do
          begin
            if (DataProcs.ValidValue(i,pos)) then
            begin
              X := StrToFloat(OS3MainFrm.DataGrid.Cells[pos,i]);
              X := X - Means[j];
              OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X);                // wp: DON'T OVERWRITE GRID DATA!
            end;
          end;  // next i
        end;  // next j
      end;  // if origin checked

      lReport.Clear;
      lReport.Add('WEIGHTED LEAST SQUARES (WLS) REGRESSION RESULTS');
      lReport.Add('');
      MReg(Noindep, IndepCols, DepCol, RowLabels, Means, Variances, StdDevs,
         BWeights, BetaWeights, BStdErrs, Bttests, tprobs, R2, stdErrEst,
         NCases, errorcode, PrintDesc, lReport);

      WLSReportFrame.DisplayReport(lReport);
      lReport.Clear;
    end;

  finally
    lReport.Free;
  end;
end;


procedure TWLSFrm.CreateOrGetChartFrame(AColIndex: Integer; AVarName: String;
  out AMemo: TMemo; out AChartFrame: TChartFrame);
var
  sheetTitle: String;
  tabSheet: TTabSheet;
  i: Integer;
begin
  sheetTitle := GetPageCaption(AVarName);

  // Find existing sheet first.
  for i := 1 to ResRegPageControl.PageCount-1 do
    if ResRegPageControl.Pages[i].Caption = sheetTitle then begin
      tabSheet := ResRegPageControl.Pages[i];
      AChartFrame := tabSheet.Controls[0] as TChartFrame;
      AMemo := tabSheet.Controls[1] as TMemo;
      exit;
    end;

  // Not found: create new sheet ...
  tabSheet := ResRegPageControl.AddTabSheet;
  tabSheet.Caption := sheetTitle;
  tabSheet.Tag := AColIndex;

  // ... and add ChartFrame
  AChartFrame := TChartFrame.Create(tabSheet);
  AChartFrame.Parent := tabSheet;
  AChartFrame.Align := alClient;
  AChartFrame.Chart.Legend.Alignment := laBottomCenter;
  AChartFrame.Chart.Legend.ColumnCount := 3;
  AChartFrame.Chart.Legend.TextFormat := tfHTML;
  AChartFrame.Chart.BottomAxis.Intervals.MaxLength := 80;
  AChartFrame.Chart.BottomAxis.Intervals.MinLength := 30;
  with AChartFrame.Chart.AxisList.Add do
  begin
    Alignment := calRight;
    Marks.Source := TListChartSource.Create(self);
    Marks.Style := smsLabel;
    Grid.Visible := false;
    TickColor := clNone;
  end;
  with AChartFrame.Chart.AxisList.Add do
  begin
    Alignment := calTop;
    Marks.Source := TListChartSource.Create(self);
    Marks.Style := smsLabel;
    Grid.Visible := false;
    TickColor := clNone;
  end;

  // ... and add memo
  AMemo := TMemo.Create(tabSheet);
  AMemo.Parent := tabSheet;
  AMemo.Align := alBottom;
  AMemo.BorderStyle := bsNone;
  AMemo.Font.Name := 'Courier New';
  AMemo.Font.Size := 8;
  AMemo.ReadOnly := true;
  AMemo.Scrollbars := ssAutoBoth;
  AMemo.WordWrap := false;

  // ... and splitter
  with TSplitter.Create(tabSheet) do
  begin
    Parent := tabSheet;
    Align := alBottom;
  end;
end;


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


procedure TWLSFrm.DepOutBtnClick(Sender: TObject);
begin
  if (DepVarEdit.Text <> '') then
  begin
    VarList.Items.Add(DepVarEdit.Text);
    DepVarEdit.Text := '';
    UpdateBtnStates;
  end;
end;


function TWLSFrm.GetPageCaption(AVarName: String): String;
begin
  Result := 'Plot of ' + AVarName;
end;


procedure TWLSFrm.IndInBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while (i < VarList.Items.Count) do
  begin
    if (VarList.Selected[i]) then
    begin
      IndVarList.Items.Add(VarList.Items[i]);
      VarList.Items.Delete(i);
      i := 0;
    end else
      inc(i);
  end;
  UpdateBtnStates;
end;

procedure TWLSFrm.IndOutBtnClick(Sender: TObject);
var
  i: integer;
begin
  i := 0;
  while (i < IndVarList.Items.Count) do
  begin
    if IndVarlist.Selected[i] then
    begin
      VarList.Items.Add(IndVarList.Items[i]);
      IndVarlist.Items.Delete(i);
      i := 0;
    end else
      inc(i);
  end;
  UpdateBtnStates;
end;


procedure TWLSFrm.IndVarListDblClick(Sender: TObject);
var
  index: Integer;
begin
  index := IndVarList.ItemIndex;
  if index > -1 then
  begin
    VarList.Items.Add(IndVarList.Items[index]);
    IndVarlist.Items.Delete(index);
    UpdateBtnStates;
  end;
end;


procedure TWLSFrm.PlotSquaredResiduals(AIndepCols: IntDyneVec;
  ANumIndepCols, ADepCol: Integer; AConfLevel: Double);
var
  xCol, yCol: Integer;
  xLabel, yLabel: String;
  i: Integer;
  colNoSelected: IntDyneVec = nil;
  xPoints: DblDyneVec = nil;
  yPoints: DblDyneVec = nil;
  regressionRes: TBivariateRegressionResults;
  memo: TMemo;
  chartFrame: TChartFrame;
begin
  SetLength(colNoSelected, 2);
  xCol := ADepCol;

  for i := 0 to ANumIndepCols-1 do
  begin
    // Get values
    yCol := AIndepCols[i];
    colNoSelected[0] := xCol;
    colNoSelected[1] := yCol;
    xLabel := OS3MainFrm.DataGrid.Cells[xCol, 0];
    yLabel := OS3MainFrm.DataGrid.Cells[yCol, 0];
    xPoints := CollectValues(OS3MainFrm.DataGrid, xCol, colNoSelected);
    yPoints := CollectValues(OS3MainFrm.DataGrid, yCol, colNoSelected);
    SortOnX(xPoints, yPoints);

    // Regression
    Calc_BivariateRegression(xPoints, yPoints, AConfLevel, regressionRes);

    // Create tab with chart and report controls
    CreateOrGetChartFrame(yCol, yLabel, memo, chartFrame);

    // Plot
    PlotXY(chartFrame, xPoints, yPoints, regressionRes, xLabel, yLabel);

    // Print the descriptive statistics
    WriteDescriptiveReport(memo, regressionRes, xLabel, yLabel);
  end;
end;


procedure TWLSFrm.PlotXY(AChartFrame: TChartFrame; const XPoints, YPoints: DblDyneVec;
  const ARegressionResults: TBivariateRegressionResults; const XLabel, YLabel: String);
var
  xpts: DblDyneVec = nil;
  ypts: DblDyneVec = nil;
  conf: DblDyneVec = nil;
  ext: TDoubleRect;
  i: Integer;
  rightLabels, topLabels: TListChartSource;
  ser: TChartSeries;
begin
  rightLabels := AChartFrame.Chart.AxisList[2].Marks.Source as TListChartSource;
  rightLabels.Clear;
  topLabels := AChartFrame.Chart.AxisList[3].Marks.Source as TListChartSource;
  topLabels.Clear;
  AChartFrame.Clear;

  with ARegressionResults do
    AChartFrame.SetFooter(Format('R(X,Y) = %.3f, Slope = %.3f, Intercept = %.3f', [
      R, Slope, Intercept
    ]));

  // Data points
  AChartFrame.SetXTitle(XLabel);
  AChartFrame.SetYTitle(YLabel);
  AChartFrame.PlotXY(ptSymbols, XPoints, YPoints, nil, nil, 'Data', DATA_COLORS[0]);

  // Regression line
  SetLength(xpts, 2);
  SetLengtH(ypts, 2);
  ext := AChartFrame.Chart.GetFullExtent;
  xpts[0] := ext.a.x;
  xpts[1] := ext.b.x;
  with ARegressionResults do
  begin
    ypts[0] := Intercept + Slope * xpts[0];
    ypts[1] := Intercept + Slope * xpts[1];
  end;
  AChartFrame.PlotXY(ptLines, xpts, ypts, nil, nil, 'Predicted', clBlack);
  rightLabels.Add(ypts[1], ypts[1], 'Predicted');

  // Upper ...
  SetLength(conf, ARegressionResults.Count);
  for i := 0 to High(conf) do
    conf[i] := ARegressionResults.ConfidenceLimits(XPoints[i], true);
  ser := AChartFrame.PlotXY(ptLines, XPoints, conf, nil, nil, 'Upper confidence band', clRed);
  rightLabels.Add(ser.yValue[ser.Count-1], ser.yValue[ser.Count-1], 'UCL');

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

  // Mean lines
  with ARegressionResults do
  begin
    AChartFrame.VertLine(XMean, clGreen, psDashDot, 'Mean ' + XLabel);
    topLabels.Add(XMean, XMean, 'Mean ' + XLabel);

    AChartFrame.HorLine(YMean, clGreen, psDash, 'Mean ' + YLabel);
    rightLabels.Add(YMean, YMean, 'Mean ' + YLabel);
  end;
end;


{ Routine obtains predicted raw and standardized scores and their
  residuals.  It is assumed that the dependent variable is last in the
  list of variable column pointers stored in the ColNoSelected vector.
  Get the z predicted score and its residual }
procedure TWLSFrm.PredictIt(ColNoSelected: IntDyneVec; NoVars: integer;
  Means, StdDevs, BetaWeights: DblDyneVec; StdErrEst: double; NoIndepVars: integer);
var
  col, i, j, k, Index: integer;
  predicted, zpredicted, z1, z2, resid, residsqr: double;
begin
  col := NoVariables + 1;
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.DictGrid.ColCount := 8;
  DictionaryFrm.NewVar(col);
  OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.z';
  DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.z';

  col := NoVariables + 1;
  DictionaryFrm.NewVar(col);
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  OS3MainFrm.DataGrid.Cells[col,0] := 'z Resid.';
  DictionaryFrm.DictGrid.Cells[1,col] := 'z Resid.';
  OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 2;
  for i := 1 to NoCases do
  begin
    zpredicted := 0.0;
    for j := 0 to NoIndepVars - 1 do
    begin
      k := ColNoSelected[j];
      z1 := (StrToFloat(OS3MainFrm.DataGrid.Cells[k,i]) -  Means[j]) / StdDevs[j];
      zpredicted := zpredicted + (z1 * BetaWeights[j]);
    end;
    OS3MainFrm.DataGrid.Cells[col-1,i] := Format('%.4f',[zpredicted]);

    if StdDevs[NoVars-1] <> 0.0 then
    begin
      Index := ColNoSelected[NoVars-1];
      z2 := StrToFloat(OS3MainFrm.DataGrid.Cells[Index,i]);
      z2 := (z2 - Means[NoVars-1]) / StdDevs[NoVars-1]; // z score
      OS3MainFrm.DataGrid.Cells[col,i] := Format('%.4f',[z2 - zpredicted]);  // z residual
    end;
  end;

  // Get raw predicted and residuals
  col := NoVariables + 1;
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.NewVar(col);
  DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.Raw';
  OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.Raw';

  // calculate raw predicted scores and store in DataGrid at col
  for i := 1 to NoCases do
  begin   // predicted raw obtained from previously predicted z score
    predicted := StrToFloat(OS3MainFrm.DataGrid.Cells[col-2,i]) * StdDevs[NoVars-1] + Means[NoVars-1];
    OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[predicted]);
  end;

  // Calculate residuals of predicted raw scores end;
  col := NoVariables +1;
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.NewVar(col);
  DictionaryFrm.DictGrid.Cells[1,col] := 'Raw Resid.';
  OS3MainFrm.DataGrid.Cells[col,0] := 'Raw Resid.';

  for i := 1 to NoCases do
  begin
    Index := ColNoSelected[NoVars-1];
    resid := StrToFloat(OS3MainFrm.DataGrid.Cells[col-1,i]) -  StrToFloat(OS3MainFrm.DataGrid.Cells[Index,i]);
    OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[resid]);
  end;

  // get square of raw residuals
  col := NoVariables + 1;
  DictionaryFrm.NewVar(col);
  OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
  DictionaryFrm.DictGrid.Cells[1,col] := 'ResidSqr';
  OS3MainFrm.DataGrid.Cells[col,0] := 'ResidSqr';
  for i := 1 to NoCases do
  begin
    residsqr := StrToFloat(OS3MainFrm.DataGrid.Cells[col-1,i]);
    residsqr := residsqr * residsqr;
    OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[residsqr]);
  end;
end;


function TWLSFrm.PrepareData(out ADepCol, ANumIndepCols: Integer;
  out AIndepCols: IntDyneVec; out AWeightCol: Integer;
  out ARowLabels: StrDyneVec): Boolean;
var
  i: Integer;
  msg: String;
  C: TWinControl;
begin
  Result := false;
  AIndepCols := nil;
  ARowLabels := nil;

  if not Validate(msg, C) then
  begin
    C.SetFocus;
    ErrorMsg(msg);
    exit;
  end;

  ANumIndepCols := IndVarList.Items.Count;
  ADepCol := GetVariableIndex(OS3MainFrm.DataGrid, DepVarEdit.Text);
  AWeightCol := GetVariableIndex(OS3MainFrm.DataGrid, WeightVarEdit.Text);

  // The IndepCols store also other variables, in addition to the "real"
  // independent variables. Until I know how this works, this array must be
  // over-dimensions.
  // ARowLabels alike.
  SetLength(AIndepCols, NoVariables + 2);
  SetLength(ARowLabels, NoVariables);

  for i := 0 to ANumIndepCols-1 do
  begin
    AIndepCols[i] := GetVariableIndex(OS3MainFrm.DataGrid, IndVarList.Items[i]);
    if AIndepCols[i] = -1 then
    begin
      ErrorMsg('Dependent variable %s not found.', [IndVarList.Items[i]]);
      exit;
    end;
    ARowLabels[i] := IndVarList.Items[i];
  end;

  // Append dependent column index to the independent columns vector.
  AIndepCols[ANumIndepCols] := ADepCol;

  // Check variable types: all of them must be numeric (float or integer)
  if not IsNumericCol(ADepCol) then
  begin
    ErrorMsg('Incorrect data type of dependent variable.');
    exit;
  end;

  for i := 0 to ANumIndepCols-1 do
    if not IsNumericCol(AIndepCols[i]) then
    begin
      ErrorMsg('Incorrect data type of independent variable "%s"', [ARowLabels[i]]);
      exit;
    end;

  if (AWeightCol > -1) and (not IsNumericCol(AWeightCol)) then
  begin
    ErrorMsg('Incorrect data type of weight variable.');
    exit;
  end;

  Result := true;
end;


{ Runs the ordinary least squares regression on the grid data }
function TWLSFrm.Process_OLSRegression(AIndepCols: IntDyneVec; ANumIndepCols, ADepCol: Integer;
  ARowLabels: StrDyneVec; ANumCases: Integer; PrintAll: Boolean): Boolean;
var
  lReport: TStrings;
  means: DblDyneVec = nil;
  variances: DblDyneVec = nil;
  stdDevs: DblDyneVec = nil;
  BWeights: DblDyneVec = nil;
  BetaWeights: DblDyneVec = nil;
  BStdErrs: DblDyneVec = nil;
  BtTests: DblDyneVec = nil;
  tProbs: DblDyneVec = nil;
  R2, stdErrEst: Double;
  error: Boolean;
begin
  Result := false;

  lReport := TStringList.Create;
  try
    lReport.Add('ORDINARY LEAST SQUARES (OLS) REGRESSION RESULTS');
    lReport.Add('');

    SetLength(means, NoVariables + 2);
    SetLength(variances, NoVariables + 2);
    SetLength(stdDevs, NoVariables + 2);
    SetLength(BWeights, NoVariables + 2);
    SetLength(BetaWeights, NoVariables + 2);
    SetLength(BStdErrs, NoVariables + 2);
    SetLength(Bttests, NoVariables + 2);
    SetLength(tprobs, NoVariables + 2);

    MReg(ANumIndepCols, AIndepCols, ADepCol, ARowLabels, Means, Variances, StdDevs,
      BWeights, BetaWeights, BStdErrs, Bttests, tprobs, R2, stdErrEst,
      ANumCases, error, PrintAll, lReport);

   // if error then               // wp: Why does MReg exit with error???
   //   exit;

    // Get predicted z score, residual z score, predicted raw score,
    // residual raw score and squared raw residual score.  Place in the DataGrid
    PredictIt(AIndepCols, ANumIndepCols+1, means, stdDevs, BetaWeights, stdErrEst, ANumIndepCols);

    OLSReportFrame.DisplayReport(lReport);

    Result := true;
  finally
    lReport.Free;
  end;
end;


function TWLSFrm.Process_SquaredResidualsRegression(AIndepCols: IntDyneVec;
  ANumIndepCols, ADepCol: Integer; ARowLabels: StrDyneVec; BWeights: DblDyneVec;
  ANumCases: Integer; PrintAll: Boolean): Boolean;
var
  lReport: TStrings;
  means: DblDyneVec = nil;
  variances: DblDyneVec = nil;
  stdDevs: DblDyneVec = nil;
  BetaWeights: DblDyneVec = nil;
  BStdErrs: DblDyneVec = nil;
  BtTests: DblDyneVec = nil;
  tProbs: DblDyneVec = nil;
  R2, stdErrEst: Double;
  error: Boolean;
begin
  lReport := TStringList.Create;
  try
    lReport.Add('REGRESSION OF SQUARED RESIDUALS ON INDEPENDENT VARIABLES');
    lReport.Add('');

    SetLength(means, NoVariables + 2);
    SetLength(variances, NoVariables + 2);
    SetLength(stdDevs, NoVariables + 2);
    SetLength(BetaWeights, NoVariables + 2);
    SetLength(BStdErrs, NoVariables + 2);
    SetLength(Bttests, NoVariables + 2);
    SetLength(tprobs, NoVariables + 2);

    MReg(ANumIndepCols, AIndepCols, ADepCol, ARowLabels, Means, Variances, StdDevs,
      BWeights, BetaWeights, BStdErrs, Bttests, tprobs, R2, stdErrEst,
      ANumCases, error, PrintAll, lReport);

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


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

  VarList.Clear;
  for i := 0 to NoVariables - 1 do
    VarList.Items.Add(OS3MainFrm.DataGrid.Cells[i+1,0]);
  IndVarList.Clear;

  DepVarEdit.Clear;
  WeightVarEdit.Clear;

  if OLSReportFrame <> nil then
    OLSReportFrame.Clear;
  if ResidualsRegReportFrame <> nil then
    ResidualsRegReportFrame.clear;
  if WLSReportFrame <> nil then
    WLSReportFrame.Clear;

  WeightChk.Checked := true;
  UserWeightsChk.Checked := false;
  OriginChk.Checked := false;
  Origin2Chk.Checked := false;

  UpdateBtnStates;
end;


procedure TWLSFrm.UpdateBtnStates;
var
  i: Integer;
  lSelected: Boolean;
begin
  inherited;

  lSelected := false;
  for i:=0 to VarList.Items.Count-1 do
    if Varlist.Selected[i] then
    begin
      lSelected := true;
      break;
    end;
  DepInBtn.Enabled := lSelected and (DepVarEdit.Text = '');
  IndInBtn.Enabled := lSelected;
  WeightInBtn.Enabled := lSelected and (WeightVarEdit.Text = '') and UserWeightsChk.Checked;

  lSelected := false;
  for i:=0 to IndVarList.Items.Count-1 do
    if IndVarList.Selected[i] then
    begin
      lSelected := true;
      break;
    end;
  DepOutBtn.Enabled := (DepVarEdit.Text <> '');
  IndOutBtn.Enabled := lSelected;
  WeightOutBtn.Enabled := (WeightVarEdit.Text <> '') and UserWeightsChk.Checked;

  if OLSReportFrame <> nil then
    OLSReportFrame.UpdateBtnStates;
  if ResidualsRegReportFrame <> nil then
    ResidualsRegReportFrame.UpdateBtnStates;
  if WLSReportFrame <> nil then
    WLSReportFrame.UpdateBtnStates;

end;


procedure TWLSFrm.UserWeightsChkChange(Sender: TObject);
begin
  WeightVarEdit.Enabled := UserWeightsChk.Checked;
  Label4.Enabled := WeightVarEdit.Enabled;
  WeightInBtn.Enabled := UserWeightsChk.Checked and (VarList.ItemIndex > -1) and (WeightVarEdit.Text = '');
  WeightOutBtn.Enabled := UserWeightsChk.Checked and (WeightVarEdit.Text <> '');
end;


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

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

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

  Result := true;
end;


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


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


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


procedure TWLSFrm.WeightOutBtnClick(Sender: TObject);
begin
  if (WeightVarEdit.Text <> '') then
  begin
    VarList.Items.Add(WeightVarEdit.Text);
    WeightVarEdit.Text := '';
  end;
  UpdateBtnStates;
end;


procedure TWLSFrm.WriteDescriptiveReport(AMemo: TMemo;
  const ARegressionResults: TBivariateRegressionResults;
  const XLabel, YLabel: String);
var
  lReport: TStrings;
begin
  lReport := TStringList.Create;
  try
    {  not needed - requires too much space
    lReport.Add('Data file: %s', [OS3MainFrm.FileNameEdit.Text]);
    lReport.Add('');
    lReport.Add('Variables:');
    lReport.Add('  X: %s', [xLabel]);
    lReport.Add('  Y: %s', [yLabel]);
    lReport.Add('');
    }
    lReport.Add('Variable      Mean    Variance  Std.Dev.');
    lReport.Add('----------  --------  --------  --------');
    with ARegressionResults do
    begin
      lReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [xLabel, XMean, XVariance, XStdDev]);
      lReport.Add('%-10s  %8.2f  %8.2f  %8.2f', [yLabel, YMean, YVariance, YStdDev]);
      lReport.Add('');
      lReport.Add('Regression:');
      lReport.Add('  Correlation:                %8.3f', [R]);
      lReport.Add('  Slope:                      %8.3f', [Slope]);
      lReport.Add('  Intercept:                  %8.3f', [Intercept]);
      lReport.Add('  Standard Error of Estimate: %8.3f', [StdErrorPredicted]);
      lReport.Add('  Number of good cases:       %8d',   [Count]);
    end;
    AMemo.Lines.Assign(lReport);
  finally
    lReport.Free;
  end;
end;


end.