diff --git a/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas b/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas
index 4d7bda5b7..d262a807f 100644
--- a/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas
+++ b/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas
@@ -57,9 +57,11 @@ type
     ResidualsRegReportFrame: TReportFrame;
     WLSReportFrame: TReportFrame;
 
-    procedure AddVariable(AVarName: String; AData: DblDyneVec; ANumFormat: String);
+    procedure AddVariable(AVarName: String; AData: DblDyneVec;
+      ANumFormat: String; const ABadRows: IntDyneVec);
 
-    procedure AddWeightsToGrid(const ASqrPredictedResiduals, AWeights: DblDyneVec);
+    procedure AddWeightsToGrid(const ASqrPredictedResiduals, AWeights: DblDyneVec;
+      const ABadRows: IntDyneVec);
 
     procedure CalcWeights(xValues: DblDyneMat; ACoeffs: DblDyneVec;
       out ASquaredPredictedResiduals: DblDyneVec; out AWeights: DblDyneVec);
@@ -76,19 +78,20 @@ type
       const ARegressionResults: TBivariateRegressionResults; const XLabel, YLabel: String);
 
     procedure Predict(const xData: DblDyneMat; const yData: DblDyneVec;
-      ARegressionResults: TMultipleRegressionResults);
+      const ABadRows: IntDyneVec; ARegressionResults: TMultipleRegressionResults);
 
     function PrepareData(out AIndepCols: IntDyneVec; out ADepCol: Integer;
       out AWeightCol: Integer; out ARowLabels: StrDyneVec;
-      out xValues: DblDyneMat; out yValues: DblDyneVec): Boolean;
+      out xValues: DblDyneMat; out yValues: DblDyneVec; out AWeights: DblDyneVec;
+      out ABadRows: IntDyneVec): Boolean;
 
     function Process_OLSRegression(AIndepCols: IntDyneVec; ADepCol: Integer;
       const ARowLabels: StrDyneVec; const xValues: DblDyneMat;
-      const yValues: DblDyneVec): Boolean;
+      const yValues: DblDyneVec; const ABadRows: IntDyneVec): Boolean;
 
     function Process_SquaredResidualsRegression(AIndepCols: IntDyneVec;
       const ARowLabels: StrDyneVec; const xValues: DblDyneMat;
-      out AWeights: DblDyneVec): Boolean;
+      out AWeights: DblDyneVec; const ABadRows: IntDyneVec): Boolean;
 
     function Process_WeightedRegression(AIndepCols: IntDyneVec;
       const ARowLabels: StrDyneVec; const xValues: DblDyneMat;
@@ -173,10 +176,27 @@ end;
 
 
 { Adds a new variable names AColTitle after the last grid column,
-  and writes the specified data to the grid (in the specified number format). }
-procedure TWLSFrm.AddVariable(AVarName: String; AData: DblDyneVec; ANumFormat: String);
+  and writes the specified data to the grid (in the specified number format).
+  Rows mentioned in ABadRows must be omitted because they are not contained in
+  AData. }
+procedure TWLSFrm.AddVariable(AVarName: String; AData: DblDyneVec;
+  ANumFormat: String; const ABadRows: IntDyneVec);
+
+  function IsBadRow(ARow: Integer): Boolean;
+  var
+    j: Integer;
+  begin
+    for j := 0 to High(ABadRows) do
+      if ARow = ABadRows[j] then
+      begin
+        Result := true;
+        exit;
+      end;
+    Result := false;
+  end;
+
 var
-  i, colIndex: Integer;
+  i, j, colIndex, row: Integer;
 begin
   colIndex := GetVariableIndex(OS3MainFrm.DataGrid, AVarname);
   if colIndex = -1 then
@@ -188,17 +208,25 @@ begin
     OS3MainFrm.DataGrid.Cells[colIndex, 0] := AVarName;
     OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
   end;
+  row := 1;
   for i := 0 to High(AData) do
-    OS3MainFrm.DataGrid.Cells[colIndex, i+1] := Format(ANumFormat, [AData[i]]);
+  begin
+    while IsBadRow(row) do inc(row);
+    if row >= OS3MainFrm.DataGrid.RowCount then
+      raise Exception.Create('Bad row error.');
+    OS3MainFrm.DataGrid.Cells[colIndex, row] := Format(ANumFormat, [AData[i]]);
+    inc(row);
+  end;
 end;
 
 
 { Calculate predicted squared residuals and save recipricols to grid as weights }
-procedure TWLSFrm.AddWeightsToGrid(const ASqrPredictedResiduals, AWeights: DblDyneVec);
+procedure TWLSFrm.AddWeightsToGrid(const ASqrPredictedResiduals, AWeights: DblDyneVec;
+  const ABadRows: IntDyneVec);
 begin
   // Create new variables and add to grid
-  AddVariable('Pred SqrResid', ASqrPredictedResiduals, '%.3f');
-  AddVariable('WEIGHTS', AWeights, '%.3f');
+  AddVariable('Pred SqrResid', ASqrPredictedResiduals, '%.3f', ABadRows);
+  AddVariable('WEIGHTS', AWeights, '%.3f', ABadRows);
 end;
 
 
@@ -219,7 +247,7 @@ begin
   sum := 0;
   for i := 0 to n-1 do
   begin
-    ASquaredPredictedResiduals[i] := ACoeffs[m];  // intercept value
+    ASquaredPredictedResiduals[i] := ACoeffs[m];  // intercept coefficient
     for j := 0 to m-1 do
       ASquaredPredictedResiduals[i] += abs(xValues[i, j] * ACoeffs[j]);
     if ASquaredPredictedResiduals[i] <> 0 then
@@ -228,9 +256,6 @@ begin
       AWeights[i] := 0;
     sum := sum + AWeights[i];
   end;
-
-  // Normalize weights to 1.0
-  AWeights := AWeights * (1.0 / sum);
 end;
 
 
@@ -245,26 +270,12 @@ begin
   );
 end;
 
+
 procedure TWLSFrm.Compute;
 var
-  i, j, noIndep, NCases, pos, col: integer;
-  X, weight: 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;
-  lReport: TStrings;
-  StdErrEst: Double = 0.0;
-  R2: Double = 0.0;
-  errorcode: Boolean = false;
-  PrintDesc: boolean = true;
-
   indepCols: IntDyneVec = nil;
   rowLabels: StrDyneVec = nil;
+  badRows: IntDyneVec = nil;
   weights: DblDyneVec = nil;
   xValues: DblDyneMat = nil;
   yValues: DblDyneVec = nil;
@@ -272,164 +283,25 @@ var
   weightCol: Integer = -1;
   useOrigin: Boolean;
 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.
-    if not PrepareData(indepCols, depCol, weightCol, RowLabels, xValues, yValues) then
-      exit;
-
-    // Do the OLS regression
-    if not Process_OLSRegression(indepCols, depCol, RowLabels, xValues, yValues) then
-      exit;
-
-    // Regress the squared residuals on the predictors
-    if WeightChk.Checked then
-    begin
-      if not Process_SquaredResidualsRegression(indepCols, RowLabels, xValues, weights) then
-        exit;
-      useOrigin := OriginChk.Checked;
-    end else
-    begin
-      // Read the weights from the user column
-      weights := CollectVecValues(OS3MainFrm.DataGrid, weightCol, indepCols);
-      useOrigin := Origin2Chk.Checked;
-    end;
-
-    // Do the weighted regression, finally
-    Process_WeightedRegression(indepCols, RowLabels, xValues, yValues, weights, useOrigin);
-
-
+  if not PrepareData(indepCols, depCol, weightCol, RowLabels, xValues, yValues, weights, badRows) then
     exit;
 
+  // Do the OLS regression
+  if not Process_OLSRegression(indepCols, depCol, RowLabels, xValues, yValues, badRows) then
+    exit;
 
+  // Regress the squared residuals on the predictors
+  ResidualsRegPage.TabVisible := WeightChk.Checked;
+  if WeightChk.Checked then
+  begin
+    if not Process_SquaredResidualsRegression(indepCols, RowLabels, xValues, weights, badRows) then
+      exit;
+    useOrigin := OriginChk.Checked;
+  end else
+    useOrigin := Origin2Chk.Checked;
 
-    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
-        weight := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[NoVariables,i]));
-        for j := 0 to Noindep do
-        begin
-          pos := IndepCols[j];
-          X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
-          X := X * weight;
-          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
-        weight := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[weightCol, i]));
-        for j := 0 to Noindep do
-        begin
-          pos := indepCols[j];
-          X := StrToFloat(OS3MainFrm.DataGrid.Cells[pos,i]);
-          X := X * weight;
-          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;
+  // Do the weighted regression, finally
+  Process_WeightedRegression(indepCols, RowLabels, xValues, yValues, weights, useOrigin);
 end;
 
 
@@ -584,6 +456,8 @@ begin
 end;
 
 
+{ We will plot the selected vector of the independent values vertically,
+  and the dependent values horizontally. }
 procedure TWLSFrm.PlotSquaredResiduals(AIndepCols: IntDyneVec;
   ADepCol: Integer; const AIndepValues: DblDyneMat; const ADepValues: DblDyneVec);
 var
@@ -595,9 +469,6 @@ var
   xLabel, yLabel: String;
   numIndepCols: Integer;
 begin
-  // We will plot the selected vector of the independent values vertically,
-  // and the dependent values horizontally.
-
   xCol := ADepCol;
   x := VecCopy(ADepValues);
   xLabel := OS3MainFrm.DataGrid.Cells[xCol, 0];
@@ -623,47 +494,7 @@ begin
     WriteDescriptiveReport(memo, regressionRes, xLabel, yLabel);
   end;
 end;
-(*
-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 := CollectVecValues(OS3MainFrm.DataGrid, xCol, colNoSelected);
-    yPoints := CollectVecValues(OS3MainFrm.DataGrid, yCol, colNoSelected);
-    SortOnX(xPoints, yPoints);
-
-    // Regression
-    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);
@@ -730,12 +561,13 @@ begin
   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.Predict(const xData: DblDyneMat; const yData: DblDyneVec;
-  ARegressionResults: TMultipleRegressionResults);
+  const ABadRows: IntDyneVec; ARegressionResults: TMultipleRegressionResults);
 var
   means, stddevs, variances: DblDyneVec;
   i, j, n, m: Integer;
@@ -783,11 +615,11 @@ begin
     sqrResid[i] := sqr(rawResid[i]);
   end;
 
-  AddVariable('z Pred', zPred, '%.4f');
+  AddVariable('z Pred', zPred, '%.4f', ABadRows);
 //  AddGridColumn('z Resid', zResid, '%.4f');
-  AddVariable('Raw Pred',  rawPred, '%.3f');
-  AddVariable('Raw Resid', rawResid, '%.3f');
-  AddVariable('Sqr Resid', sqrResid, '%.3f');
+  AddVariable('Raw Pred',  rawPred, '%.3f', ABadRows);
+  AddVariable('Raw Resid', rawResid, '%.3f', ABadRows);
+  AddVariable('Sqr Resid', sqrResid, '%.3f', ABadRows);
 end;
 
 
@@ -802,21 +634,28 @@ end;
     - xValues: matrix with all independent values. The columns of the matrix
         correspond to the variables, the row correspond to the cases.
     - yValues: vector with the dependent variable values
+    - ABadRows: indices of rows in which there is at least one invalid value
+        in the colums specified by AIndepCols. A value is "invalid" when it
+        is filtered, numeric but empty, or contains the missing value code.
 }
 function TWLSFrm.PrepareData(out AIndepCols: IntDyneVec; out ADepCol: Integer;
   out AWeightCol: Integer; out ARowLabels: StrDyneVec;
-  out xValues: DblDyneMat; out yValues: DblDyneVec): Boolean;
+  out xValues: DblDyneMat; out yValues: DblDyneVec; out AWeights: DblDyneVec;
+  out ABadRows: IntDyneVec): Boolean;
 var
-  i: Integer;
+  i, n: Integer;
   msg: String;
   C: TWinControl;
   numIndepCols: Integer;
+  cols: IntDyneVec = nil;
 begin
   Result := false;
   AIndepCols := nil;
   ARowLabels := nil;
   xValues := nil;
   yvalues := nil;
+  AWeights := nil;
+  ABadRows := nil;
 
   if not Validate(msg, C) then
   begin
@@ -864,8 +703,42 @@ begin
     exit;
   end;
 
-  xValues := CollectMatValues(OS3MainFrm.DataGrid, AIndepCols);
-  yValues := CollectVecValues(OS3MainFrm.DataGrid, ADepCol);
+  // Prepare list of all column indices to be loaded: x, y, weights
+  // ADepCol will follow the x columns immediatey, WeightCol is last.
+  if AWeightCol > -1 then
+  begin
+    SetLength(cols, numIndepCols + 2);
+    cols[numIndepCols] := ADepCol;
+    cols[numIndepCols+1] := AWeightCol;
+  end else
+  begin
+    SetLength(cols, NumIndepCols + 1);
+    cols[numIndepCols] := ADepCol;
+  end;
+  for i := 0 to numIndepCols-1 do cols[i] := AIndepCols[i];
+
+  // Determine list of indices of rows containing invalid entries.
+  SetLength(ABadRows, OS3MainFrm.DataGrid.RowCount);
+  n := 0;
+  for i := 1 to OS3MainFrm.DataGrid.RowCount-1 do
+    if not GoodRecord(OS3MainFrm.DataGrid, i, cols) then
+    begin
+      ABadRows[n] := i;
+      inc(n);
+    end;
+  SetLength(ABadRows, n);
+
+  // Extract data values; take care to omit invalid values in both x and y
+  xValues := CollectMatValues(OS3MainFrm.DataGrid, cols);
+  // The y column has index numIndepCols, i.e. follows the x columns.
+  yValues := MatColVector(xValues, numIndepCols);
+  MatColDelete(xValues, numIndepCols);
+  if AWeightCol > -1 then
+  begin
+    // The weight column is the last column
+    AWeights := MatColVector(xValues, High(xValues[0]));
+    MatColDelete(xValues, High(xValues[0]));
+  end;
 
   Result := true;
 end;
@@ -874,7 +747,8 @@ end;
 { Runs the ordinary least squares regression on the grid data }
 function TWLSFrm.Process_OLSRegression(AIndepCols: IntDyneVec;
   ADepCol: Integer; const ARowLabels: StrDyneVec;
-  const xValues: DblDyneMat; const yValues: DblDyneVec): Boolean;
+  const xValues: DblDyneMat; const yValues: DblDyneVec;
+  const ABadRows: IntDyneVec): Boolean;
 var
   lReport: TStrings;
   regressionRes: TMultipleRegressionResults;
@@ -899,7 +773,7 @@ begin
 
     if Result then
     begin
-      Predict(xValues, yValues, regressionRes);
+      Predict(xValues, yValues, ABadRows, regressionRes);
       OLSReportFrame.DisplayReport(lReport);
     end;
 
@@ -910,7 +784,8 @@ end;
 
 
 function TWLSFrm.Process_SquaredResidualsRegression(AIndepCols: IntDyneVec;
-  const ARowLabels: StrDyneVec; const xValues: DblDyneMat; out AWeights: DblDyneVec): Boolean;
+  const ARowLabels: StrDyneVec; const xValues: DblDyneMat; out AWeights: DblDyneVec;
+  const ABadRows: IntDyneVec): Boolean;
 var
   lReport: TStrings;
   sqrResiduals: DblDyneVec;
@@ -919,7 +794,6 @@ var
   i, depCol, numIndepCols: Integer;
 begin
   AWeights := nil;
-  ResidualsRegPage.TabVisible := WeightChk.Checked;
 
   if not WeightChk.Checked then
     exit;
@@ -957,7 +831,7 @@ begin
 
       // Store weights to the grid
       if SaveChk.Checked then
-        AddWeightsToGrid(predSqrResiduals, AWeights);
+        AddWeightsToGrid(predSqrResiduals, AWeights, ABadRows);
     end;
 
   finally
@@ -974,19 +848,27 @@ var
   regressionRes: TMultipleRegressionResults;
   lReport: TStrings;
   means: DblDyneVec;
+  yMean: Double;
 begin
   MatSize(xValues, n, m);
 
   for i :=0 to n-1 do
+  begin
     for j := 0 to m-1 do
       xValues[i, j] := xValues[i, j] * AWeights[i];
+    yValues[i] := yValues[i] * AWeights[i];
+  end;
 
   if SubtractMeans then
   begin
-    means := MatRowMeans(xValues);
+    means := MatColMeans(xValues);
+    yMean := VecMean(yValues);
     for i := 0 to n-1 do
+    begin
       for j := 0 to m-1 do
-        xValues[i, j] := xValues[i, j] - means[i];
+        xValues[i, j] := xValues[i, j] - means[j];
+      yValues[i] := yValues[i] - yMean;
+    end;
   end;
 
   lReport := TStringList.Create;
@@ -1019,9 +901,8 @@ var
 begin
   err := MultipleRegression(xValues, yValues, CONF_LEVEL, ARegressionResults);
   case err of
-    regOK: ;
-    regTooFewValues: ErrorMsg('At least two values required for regression.');
-    regStdDevZero: ErrorMsg('Standard deviation is zero.');
+    regOK           : ;
+    regTooFewValues : ErrorMsg('At least two values required for regression.');
   end;
   Result := (err = regOK);
 
diff --git a/applications/lazstats/source/units/gridprocs.pas b/applications/lazstats/source/units/gridprocs.pas
index 01b079fbc..23f1f40d9 100644
--- a/applications/lazstats/source/units/gridprocs.pas
+++ b/applications/lazstats/source/units/gridprocs.pas
@@ -75,27 +75,29 @@ end;
 { Extracts the grid values from the columns with indices given by AColIndices
   and puts them into the columns of the result matrix.
   This means: The result matrix contains the variables as columns and the
-  cases as rows. }
+  cases as rows.
+  "Bad" records (filtered, empty) are skipped. }
 function CollectMatValues(AGrid: TStringGrid; AColIndices: IntDyneVec): DblDyneMat;
 var
-  nr, r, c, i, j: Integer;
+  r, c, i, j: Integer;
   val: Double;
 begin
   SetLength(Result, AGrid.RowCount, Length(AColIndices));
-  nr := 0;
+  i := 0;
   for r:= 1 to AGrid.RowCount-1 do
   begin
     if not GoodRecord(AGrid, r, AColIndices) then Continue;
-    i := r - 1;
     for j := 0 to High(AColIndices) do
     begin
       c := AColIndices[j];
       if TryStrToFloat(trim(AGrid.Cells[c, r]), val) then
-        Result[i, j] := val;
+        Result[i, j] := val
+      else
+        Result[i, j] := NaN;
     end;
-    inc(nr);  // count the number of rows in the matrix.
+    inc(i);
   end;
-  SetLength(Result, nr);
+  SetLength(Result, i);
 end;
 
 
diff --git a/applications/lazstats/source/units/matrixunit.pas b/applications/lazstats/source/units/matrixunit.pas
index c40885bd5..f126387b5 100644
--- a/applications/lazstats/source/units/matrixunit.pas
+++ b/applications/lazstats/source/units/matrixunit.pas
@@ -32,6 +32,7 @@ procedure VecSize(A: TDblVector; out n: Integer);
 
 procedure VecMaxMin(const AData: TDblVector;
   out AMax, AMin: Double);
+function VecMean(const AData: TDblVector): Double;
 procedure VecMeanStdDev(const AData: TDblVector;
   out AMean, AStdDev: Double);
 procedure VecMeanVarStdDev(const AData: TDblVector;
@@ -57,6 +58,7 @@ operator * (A: TDblMatrix; v: TDblVector): TDblVector;
 function MatAppendColVector(A: TDblMatrix; v: TDblVector): TDblMatrix;
 procedure MatCheck(A: TDblMatrix);
 procedure MatCheckSquare(A: TDblMatrix; out n: Integer);
+procedure MatColDelete(A: TDblMatrix; ACol: Integer);
 procedure MatColMeanVarStdDev(A: TDblMatrix; out AMeans, AVariances, AStdDevs: TDblVector);
 function MatColMeans(A: TDblMatrix): TDblVector;
 function MatColVector(A: TDblMatrix; AColIndex: Integer): TDblVector;
@@ -227,6 +229,22 @@ begin
 end;
 
 
+function VecMean(const AData: TDblVector): Double;
+var
+  i, n: Integer;
+begin
+  Result := 0;
+  n := Length(AData);
+  if n > 0 then
+  begin
+    for i := 0 to n-1 do
+      Result := Result + AData[i];
+    Result := Result / n;
+  end else
+    Result := NaN;
+end;
+
+
 procedure VecMeanStdDev(const AData: TDblVector; out AMean, AStdDev: Double);
 var
   variance: Double;
@@ -432,6 +450,23 @@ begin
 end;
 
 
+procedure MatColDelete(A: TDblMatrix; ACol: Integer);
+var
+  n, m, i, j: Integer;
+begin
+  MatSize(A, n,m);
+  if (ACol < 0) or (ACol >= m) then
+    raise EMatrix.Create('MatColDelete: illegal column index.');
+
+  for i := 0 to n - 1 do begin
+    for j := 0 to m - 2 do
+      if j >= ACol then
+        A[i, j] := A[i, j+1];
+    SetLength(A[i], m-1);
+  end;
+end;
+
+
 procedure MatColMeanVarStdDev(A: TDblMatrix; out AMeans, AVariances, AStdDevs: TDblVector);
 var
   n, m, i, j: Integer;
diff --git a/applications/lazstats/source/units/regressionunit.pas b/applications/lazstats/source/units/regressionunit.pas
index ad4efb9ac..7e9e500f4 100644
--- a/applications/lazstats/source/units/regressionunit.pas
+++ b/applications/lazstats/source/units/regressionunit.pas
@@ -14,7 +14,7 @@ type
   ERegression = class(Exception);
 
 type
-  TRegressionError = (regOK, regTooFewValues, regStdDevZero);
+  TRegressionError = (regOK, regTooFewValues);
 
   TBivariateRegressionResults = record
   public