Updating models and Likelihood Ratio tests

  Model selection is always a difficult procedure, made more difficult in the case of longitudinal data analysis with the additional complication of a model for the variance structure as well as for the mean structure. A common approach to the problem is to estimate the variance parameters $\nu$, $\tau$ and $\phi$ by specifying the most complex mean model you are willing to entertain (say, the group*as.factor(time) model, with one mean parameter per time point in each group), and then keep the variance parameters fixed at those estimates when selecting simpler mean models.

The S model selection function update works with remlfit objects and encourages this kind of analysis. You can use update to specify a new mean model relative to another fit (by dropping one or more terms, for example) but keeping the variance parameters fixed. (Incidentally, such a submodel with fixed variance parameters is comparatively cheap to calculate, as it requires only a single likelihood calculation rather that an iterative procedure to maximise the likelihood.) The update function, like the chitest function, takes two arguments: a remlfit object specifying the ``base'' model, and an update formula specifying the new model relative to the previous model.

Of course, our ultimate aim when updating models in this fashion is inferential comparison of the two models. We can do this quite easily with a Likelihood Ratio test (as described below), but this is only valid when ML likelihoods (not REML likelihoods) are maximised. For this reason we use the ml.fit function for the analyses in the example below.

Example 5.7: plausible mean model for the lograts data might be for a common polynomial within each time series, with a constant offset for each treatment group. We begin by fitting the saturated model to estimate the variance parameters. We then fit a common quadratic model by specifying a new formula to create rats.quad. To create rats.lin we drop the quadratic term (time^2) from rats.quad.

To compare the models thus created, we can use the anova function to test a submodel within a more complex model using a likelihood ratio test. Note that this is only appropriate if ML (not REML) likelihoods were calculated. (If REML likelihoods are used, the anova function gives an error message.) The function to calculate the likelihood ratio statistic is anova, and takes two arguments which are the model objects to be compared.

To test whether the quadratic term is necessary in the lograts example, we use a likelihood ratio test to compare rats.quad with rats.lin: