Skip to content

Commit 8c6b978

Browse files
pr0m1th3aspr0m1th3as
committed
GeneralizedLinearModel: add comprehensive BIST coverage
 The class shipped with a single construction test; its properties and methods were exercised only indirectly through fitglm. Add a %!shared dataset and 20 tests covering the object directly: * Normal-distribution GLM reproduces OLS (coefficients, fitted values, residuals, deviance) and its standard errors match regress; the dispersion is the estimated SSE/DFE and the log-likelihood matches the closed form. * Poisson fit: coefficients and fit statistics (MATLAB-verified), the scalar count/size properties, unit fixed dispersion, coefficient names, and a symmetric covariance with SE^2 on the diagonal. * Methods: predict (fitted response and confidence bounds), feval, coefCI (default and custom alpha), coefTest, devianceTest, random, AIC/BIC identities, disp output, chained subsref, and the diagnostic/effect plots. * Binomial/logistic fit: coefficients agree with glmfit, fitted values lie in [0, 1], and the deviance is checked; plus an interaction model and an offset. * inst/GeneralizedLinearModel.m: Add the shared dataset and the tests above.
1 parent 3c10335 commit 8c6b978

1 file changed

Lines changed: 155 additions & 0 deletions

File tree

‎inst/GeneralizedLinearModel.m‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1240,3 +1240,158 @@ function disp (this)
12401240
%! mdl = GeneralizedLinearModel ([1 2; 2 1; 3 4; 4 3; 5 6; 6 5], ...
12411241
%! [1;0;2;3;2;4], "linear", ...
12421242
%! "Distribution", "poisson"); mdl(1);
1243+
1244+
## Comprehensive property and method coverage
1245+
%!shared X, yp, yb, yn
1246+
%! X = [ 0.37, 0.06, 1.76; -0.76, -1.52, 0.84; 0.76, -0.19, -0.47; ...
1247+
%! -0.80, -2.74, -0.90; 0.08, 0.39, 1.05; -0.41, -0.03, 0.74; ...
1248+
%! 0.23, 1.21, 0.35; 0.66, 0.94, 0.13; 0.66, -0.12, -0.06; ...
1249+
%! 2.09, 1.33, -0.71; 1.50, 0.08, -0.52; 0.59, 0.07, -1.13; ...
1250+
%! -1.17, -0.35, -1.28; 0.68, 0.63, -0.80; -0.69, 0.08, 0.41; ...
1251+
%! 2.04, 0.96, -0.56];
1252+
%! yp = [5 2 0 3 1 1 0 1 2 1 3 0 0 1 1 3]';
1253+
%! yb = [1 1 1 0 0 1 1 1 1 1 1 0 0 0 0 1]';
1254+
%! yn = [2.1 -0.3 1.2 -1.1 0.8 0.4 1.5 1.1 0.6 2.9 2.0 0.7 -1.3 0.9 -0.2 2.5]';
1255+
1256+
%!test # a normal-distribution GLM with identity link reproduces OLS exactly
1257+
%! mdl = fitglm (X, yn, "Distribution", "normal");
1258+
%! b_ols = [ones(16, 1), X] \ yn;
1259+
%! assert_equal (mdl.Coefficients.Estimate, b_ols, 1e-10);
1260+
%! assert_equal (mdl.Fitted.Response, [ones(16, 1), X] * b_ols, 1e-10);
1261+
%! assert_equal (mdl.Residuals.Raw, yn - [ones(16, 1), X] * b_ols, 1e-10);
1262+
%! assert_equal (mdl.Deviance, sum ((yn - [ones(16, 1), X] * b_ols) .^ 2), 1e-10);
1263+
%! assert_equal (mdl.Link.Name, "identity");
1264+
1265+
%!test # the normal-GLM standard errors match those from regress
1266+
%! mdl = fitglm (X, yn, "Distribution", "normal");
1267+
%! [~, bint] = regress (yn, [ones(16, 1), X], 0.05);
1268+
%! se_reg = (bint(:,1) - bint(:,2)) / 2 / tinv (0.025, 12);
1269+
%! assert_equal (mdl.Coefficients.SE, se_reg, 1e-9);
1270+
1271+
%!test # normal-GLM dispersion is estimated as SSE/DFE; log-likelihood closed form
1272+
%! mdl = fitglm (X, yn, "Distribution", "normal");
1273+
%! rss = mdl.Deviance;
1274+
%! assert_equal (mdl.DispersionEstimated, true);
1275+
%! assert_equal (mdl.Dispersion, rss / mdl.DFE, 1e-12);
1276+
%! assert_equal (mdl.LogLikelihood, -8 * (log (2 * pi * rss / 16) + 1), 1e-6);
1277+
1278+
%!test # Poisson coefficients and fit statistics (verified against MATLAB)
1279+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1280+
%! assert_equal (mdl.Coefficients.Estimate, ...
1281+
%! [-0.3420955; 1.2804868; -1.0743272; 0.8395779], 1e-6);
1282+
%! assert_equal (mdl.Deviance, 7.403008, 1e-5);
1283+
%! assert_equal (mdl.LogLikelihood, -18.543280, 1e-5);
1284+
%! assert_equal (mdl.ModelCriterion.AIC, 45.086559, 1e-5);
1285+
%! assert_equal (mdl.Rsquared.Deviance, 0.6627677, 1e-6);
1286+
1287+
%!test # scalar count/size properties of the Poisson fit
1288+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1289+
%! assert_equal (mdl.NumCoefficients, 4);
1290+
%! assert_equal (mdl.NumEstimatedCoefficients, 4);
1291+
%! assert_equal (mdl.NumPredictors, 3);
1292+
%! assert_equal (mdl.NumObservations, 16);
1293+
%! assert_equal (mdl.DFE, 12);
1294+
%! assert_equal (mdl.ResponseName, "y");
1295+
%! assert_equal (mdl.CoefficientNames, {'(Intercept)', 'x1', 'x2', 'x3'});
1296+
1297+
%!test # Poisson has a fixed unit dispersion (not estimated)
1298+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1299+
%! assert_equal (mdl.Dispersion, 1);
1300+
%! assert_equal (mdl.DispersionEstimated, false);
1301+
%! assert_equal (mdl.Distribution.Name, "poisson");
1302+
%! assert_equal (mdl.Link.Name, "log");
1303+
1304+
%!test # the coefficient covariance is symmetric with SE^2 on its diagonal
1305+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1306+
%! C = mdl.CoefficientCovariance;
1307+
%! assert_equal (size (C), [4, 4]);
1308+
%! assert_equal (C, C', 1e-14);
1309+
%! assert_equal (diag (C), mdl.Coefficients.SE .^ 2, 1e-12);
1310+
1311+
%!test # predict at the training data reproduces the fitted response
1312+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1313+
%! [yhat, yci] = predict (mdl, X);
1314+
%! assert_equal (yhat, mdl.Fitted.Response, 1e-10);
1315+
%! assert_equal (size (yci), [16, 2]);
1316+
%! assert_equal (all (yci(:,1) <= yhat & yhat <= yci(:,2)), true);
1317+
1318+
%!test # feval evaluates the model and agrees with predict
1319+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1320+
%! assert_equal (feval (mdl, X(:,1), X(:,2), X(:,3)), predict (mdl, X), 1e-12);
1321+
1322+
%!test # coefCI matches the t-interval and honours a custom alpha
1323+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1324+
%! b = mdl.Coefficients.Estimate; se = mdl.Coefficients.SE;
1325+
%! t95 = tinv (0.975, mdl.DFE);
1326+
%! assert_equal (coefCI (mdl), [b - t95 * se, b + t95 * se], 1e-12);
1327+
%! t90 = tinv (0.95, mdl.DFE);
1328+
%! assert_equal (coefCI (mdl, 0.10), [b - t90 * se, b + t90 * se], 1e-12);
1329+
1330+
%!test # coefTest gives the Wald F statistic against the constant model
1331+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1332+
%! [p, F, df] = coefTest (mdl);
1333+
%! assert_equal (F, 3.685312, 1e-5);
1334+
%! assert_equal (p, 0.04331745, 1e-7);
1335+
%! assert_equal (df, 3);
1336+
1337+
%!test # devianceTest chi-square equals the drop from the null deviance
1338+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1339+
%! dt = devianceTest (mdl);
1340+
%! assert_equal (class (dt), "table");
1341+
%! assert_equal (dt.chi2Stat(2), dt.Deviance(1) - dt.Deviance(2), 1e-10);
1342+
1343+
%!test # information criteria satisfy their defining identities
1344+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1345+
%! k = mdl.NumEstimatedCoefficients; ll = mdl.LogLikelihood;
1346+
%! assert_equal (mdl.ModelCriterion.AIC, -2 * ll + 2 * k, 1e-9);
1347+
%! assert_equal (mdl.ModelCriterion.BIC, -2 * ll + k * log (16), 1e-9);
1348+
1349+
%!test # raw residuals are response minus fit; random draws match the response size
1350+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1351+
%! assert_equal (mdl.Residuals.Raw, yp - mdl.Fitted.Response, 1e-12);
1352+
%! ysim = random (mdl);
1353+
%! assert_equal (size (ysim), [16, 1]);
1354+
%! assert_equal (all (ysim == round (ysim) & ysim >= 0), true);
1355+
1356+
%!test # binomial/logistic fit: coefficients agree with the glmfit engine
1357+
%! mdl = fitglm (X, yb, "Distribution", "binomial");
1358+
%! assert_equal (mdl.Coefficients.Estimate, glmfit (X, yb, "binomial"), 1e-8);
1359+
%! assert_equal (mdl.Link.Name, "logit");
1360+
%! assert_equal (all (mdl.Fitted.Response >= 0 & mdl.Fitted.Response <= 1), true);
1361+
%! assert_equal (mdl.Deviance, 10.997099, 1e-5);
1362+
1363+
%!test # an interaction model adds the cross term and one coefficient
1364+
%! mdl = fitglm (X, yp, "interactions", "Distribution", "poisson");
1365+
%! assert_equal (any (strcmp (mdl.CoefficientNames, "x1:x2")), true);
1366+
%! assert_equal (mdl.NumCoefficients, 7);
1367+
1368+
%!test # an offset is stored and applied
1369+
%! mdl = fitglm (X, yp, "Distribution", "poisson", "Offset", log (2 * ones (16, 1)));
1370+
%! assert_equal (numel (mdl.Offset), 16);
1371+
1372+
%!test # disp prints the model header and the coefficient table
1373+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1374+
%! s = evalc ("disp (mdl)");
1375+
%! assert_equal (isempty (strfind (s, "Generalized linear regression model")), false);
1376+
%! assert_equal (isempty (strfind (s, "Estimate")), false);
1377+
1378+
%!test # chained subsref reaches property -> table column -> element
1379+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1380+
%! assert_equal (numel (mdl.Coefficients.Estimate), 4);
1381+
%! assert_equal (mdl.Coefficients.Estimate(1), -0.3420955, 1e-6);
1382+
%! assert_equal (mdl.Coefficients.Estimate(2), mdl.Coefficients{2, "Estimate"}, 1e-12);
1383+
1384+
%!test # the diagnostic and effect plots run without error
1385+
%! mdl = fitglm (X, yp, "Distribution", "poisson");
1386+
%! hf = figure ("visible", "off");
1387+
%! unwind_protect
1388+
%! plotResiduals (mdl);
1389+
%! plotResiduals (mdl, "fitted", "ResidualType", "Pearson");
1390+
%! plotDiagnostics (mdl);
1391+
%! plotDiagnostics (mdl, "cookd");
1392+
%! plotEffects (mdl);
1393+
%! plotAdjustedResponse (mdl, 1);
1394+
%! plotAdded (mdl, "x2");
1395+
%! unwind_protect_cleanup
1396+
%! close (hf);
1397+
%! end_unwind_protect

0 commit comments

Comments
 (0)