Why Reproducibility Matters in GLM Comparison
A comparison of generalized linear models (GLMs) is meaningful only when the fitted models are being compared on a suitable common basis. R can return log-likelihoods, AIC values, BIC values, and analysis-of-deviance tables with little code, but those outputs do not establish that the underlying models used the same response or observations. A reproducible comparison makes those checks explicit.
This tutorial focuses on a practical R workflow: fitting candidate models on shared data, checking their fitted samples, and choosing the comparison function that matches the question. The examples use hypothetical claim-incidence data and ordinary, unweighted GLMs; they do not introduce exposure offsets or weighted likelihoods.
The R Functions and What They Compare
For a fitted GLM, logLik extracts the maximized log-likelihood and reports a degrees-of-freedom attribute that counts estimated model parameters for information-criterion calculations. AIC and BIC compare fit with a penalty for model complexity. In the usual maximum-likelihood formulation,
where \(\ell\) is the maximized log-likelihood, \(k\) is the number of estimated parameters, and \(n\) is the number of observations used by the fit. Smaller values are preferred among models for which the criterion is appropriately comparable. R obtains the relevant values from each fitted object; inspecting logLik and the fitted sample helps verify what went into the result.
The anova function serves a different purpose. For nested GLMs fit to the same response observations with a compatible likelihood, it can report changes in deviance and, when requested, a likelihood-ratio test. Nesting means that the smaller model can be obtained from the larger one by setting one or more of the larger model’s parameters to specified values, commonly zero. An analysis-of-deviance comparison is not a general test for any two models.
A Reproducible Comparison Workflow
A reliable workflow separates data preparation, model fitting, sample checks, and interpretation. The following sequence is useful when candidates differ in their predictor terms.
Decide whether the goal is relative information-criterion ranking, a test of added terms in nested models, or another stated objective. Do not select a function before deciding what question it answers.
Include the response and every variable used by any candidate. Resolve missing and invalid values consistently before fitting.
State the response, formula, family and link, and missing-value behavior. For a common prefiltered data set, setting na.action to na.fail can prevent R from silently dropping additional rows.
Compare nobs, response definitions, and the row names in the model frames. Inspect logLik and its degrees-of-freedom attribute.
State the comparison method, the values compared, and relevant limitations. A lower criterion is a relative result among the candidates, not proof that the model is correct.
A practical R check is to compare the row names of each model frame, not only the number of rows: two models can both use 150 observations but use different 150 observations. Also inspect the response recorded in each model frame. This is particularly important if one candidate formula contains a different response, a transformation that creates missing values, or a factor whose treatment has changed.
Worked Comparisons
Worked Example: Comparing AIC and BIC on Shared Rows
A hypothetical insurer models whether a small commercial policy has a claim during the coming year. Two logistic GLMs use the same 240 policy records and the same binary claim response. Model A has four estimated parameters and maximized log-likelihood \(-118.4\); Model B has six estimated parameters and maximized log-likelihood \(-113.9\). Compare their AIC and BIC values. The reported log-likelihoods are stipulated summaries of these hypothetical fits.
Method. Since both fits use the same response observations and likelihood specification, calculate the two information criteria for each model. The common sample size is \(n=240\).
For Model A, \(-2\ell=-2(-118.4)=236.8\). Thus \[ \mathrm{AIC}_A=236.8+2(4)=244.8. \] For Model B, \(-2\ell=-2(-113.9)=227.8\), so \[ \mathrm{AIC}_B=227.8+2(6)=239.8. \] The AIC difference, Model A minus Model B, is \(244.8-239.8=5.0\); AIC favors Model B.
For BIC, \(\log(240)\approx5.480639\). Therefore \[ \mathrm{BIC}_A=236.8+4(5.480639)=258.7226, \] and \[ \mathrm{BIC}_B=227.8+6(5.480639)=260.6838. \] Rounded to one decimal place, the BIC values are 258.7 and 260.7. BIC favors Model A. This disagreement is possible because BIC applies a larger complexity penalty here; neither criterion says that its preferred model is certainly the true data-generating model.
R implementation. After fitting the models on the same data, commands such as logLik(fitA), AIC(fitA, fitB), and BIC(fitA, fitB) provide the corresponding summaries. Check that nobs(fitA) and nobs(fitB) both equal 240 and inspect the model frames before interpreting the table. The common mistake avoided is choosing the smaller printed criterion without first checking that the fits are comparable.
Worked Example: Correcting a Missing-Data Comparison
A hypothetical service plan study models whether a customer makes a claim. One logistic model uses age and sales channel; another uses age and region. Initially, the first fit uses 180 rows and has log-likelihood \(-72.0\) with four parameters. The second uses 164 rows and has log-likelihood \(-65.0\) with five parameters. Explain why comparing their AIC values directly is not a sound model comparison, then describe a common-row refit.
The initial calculations are \(\mathrm{AIC}_1=-2(-72.0)+2(4)=144+8=152.0\), and \(\mathrm{AIC}_2=-2(-65.0)+2(5)=130+10=140.0\). The second value is smaller, but the models were fitted to different customer records. The difference reflects both model fit and a change in the data contributing to the likelihood; it cannot be attributed cleanly to the predictor choice.
A defensible workflow first forms a shared set of complete rows for the response, age, channel, and region. In R, a data frame called dat_cmp can be created using dat[complete.cases(dat[, c("claim", "age", "channel", "region")]), ]. Fit both candidate formulas to this data frame with the same binomial family and, for example, na.action = na.fail. Then check that both model frames contain the same row names and response.
Suppose the common-row refits each use 158 records. For illustration, let their summaries be log-likelihood \(-64.8\) with four parameters and \(-63.9\) with five parameters. Their AIC values are \(129.6+8=137.6\) and \(127.8+10=137.8\), respectively. On these common rows, the first model has the lower AIC by 0.2. The small difference should be reported as weak separation by this criterion, not as a decisive practical advantage. These refit summaries are hypothetical and are not inferred from the initial fits.
The R checks are nobs(fit1), nobs(fit2), and identical(rownames(model.frame(fit1)), rownames(model.frame(fit2))). Inspecting the response column in each model frame provides a further check. The trap is to believe that a smaller AIC validates a comparison when missingness has changed the analysis population.
Worked Example: Testing Added Terms with anova
A hypothetical portfolio has 200 policy records and a binary claim response. A logistic model includes age and a three-level region factor; a second model adds the age-by-region interaction. Suppose the first fit has four estimated parameters and log-likelihood \(-90.0\), while the second has six parameters and log-likelihood \(-86.7\). Both use exactly the same 200 rows, response, binomial family, and logit link. Test whether the added interaction improves fit at the 5% level using the likelihood-ratio approach.
Check nesting. The model with age and region is nested in the model with age, region, and their interaction: setting the two interaction coefficients to zero recovers the smaller model. The parameter-count difference is \(6-4=2\).
For a GLM, the residual deviance is twice the difference between the saturated-model log-likelihood and the fitted-model log-likelihood. If the same response data and likelihood are used, the saturated-model term cancels when subtracting the two deviances:
Here the deviance reduction is \(2[-86.7-(-90.0)]=2(3.3)=6.6\). Under the null hypothesis that both interaction coefficients are zero, and the usual regularity conditions for the likelihood-ratio approximation, this statistic is approximately chi-squared with 2 degrees of freedom. For a chi-squared variable with 2 degrees of freedom, the upper-tail probability at 6.6 is \(\exp(-6.6/2)=\exp(-3.3)\approx0.0369\). Since \(0.0369<0.05\), reject the null at the 5% level; the data provide evidence of improved likelihood from adding the interaction.
In R, place the smaller model first in anova(fit_small, fit_large, test = "Chisq"). Verify the common sample and nesting before reading the test row. This approximate test depends on the specified model and regularity conditions; it is not a guarantee that the larger model will predict better on new data.
What to Verify Before Reporting
R makes it easy to obtain output, but model comparison remains an analytical decision. Use these checks before stating which candidate is preferred:
- Same outcome: Confirm that the response has the same meaning, coding, and observations in both fits. For a binary outcome, check which category is treated as the event.
- Same analysis rows: Compare row identities, not just sample sizes. Explicitly create a common analysis data set when candidate predictors have different missingness.
- Compatible likelihood: Check that models use an appropriate common distributional specification and link for the comparison being made. Do not compare likelihood criteria across incompatible response definitions.
- Correct method: Use anova with a likelihood-ratio test only for nested candidates under suitable conditions. For non-nested candidates, it does not become a valid nested-model test simply because both are GLMs.
- Honest interpretation: AIC and BIC rank candidates under their respective penalties. They do not establish causation, remove model misspecification, or replace assessment of the analysis objective.
For a compact reproducibility record, report the formulas, response and coding, family and link, number of observations, how missing data were handled, and the comparison statistic. For a nested test, also state the models’ nesting relationship and the test reference distribution. These details make it possible for another analyst to reconstruct the comparison rather than merely repeat a conclusion.
Check Your Understanding
For each problem, identify the appropriate check or calculation before interpreting the result.
- Two binomial GLMs use the same 120 records. Their log-likelihoods are \(-48.0\) and \(-45.5\), with 3 and 5 estimated parameters. Calculate both AIC values and state which is lower. (Answer: AICs are 102.0 and 101.0; the second model has the lower AIC.)
- Two GLMs each report 90 observations, but the row names of their model frames differ. Can you conclude that their AIC values compare fit on the same data? (Answer: No; equal sample sizes do not establish identical observations.)
- A smaller model has log-likelihood \(-70.2\) and 4 parameters. A nested larger model has log-likelihood \(-67.2\) and 6 parameters, fitted to the same 150 observations. Calculate the likelihood-ratio statistic and its degrees of freedom. (Answer: The statistic is \(2[-67.2-(-70.2)]=6.0\), with 2 degrees of freedom.)
- A Poisson GLM and a binomial GLM use different response variables and have AIC values of 300 and 295. Is the second model necessarily preferred? (Answer: No; the AIC values are not a valid direct comparison of models for different response variables.)
- Two nested logistic models use the same response and 210 observations. Adding three parameters reduces deviance by 2.1. What chi-squared reference degrees of freedom apply, and is the result significant at 5%? (Answer: 3 degrees of freedom; the upper-tail probability is approximately 0.552, so it is not significant at 5%.)