Unlike most statistical software packages, R often stores the results of an analysis in an object. The advantage of this is that while not all output is shown in the screen ad once, it is neither necessary to estimate the statistical model again if different output is required.
This paragraph will show the kind of data that is stored in a multilevel model estimated by R-Project and introduce some functions that make use of this data.
Inside the model
Let’s first estimate a simple multilevel model, using the nlme-package. For this paragraph we will use a model we estimated earlier: the education model with a random intercept and a random slope. This time though, we will assign it to an object that we call model.01. It is estimated as follows:
model.01 <- lme(fixed = normexam ~ standLRT, data = Exam,
random = ~ standLRT | school)
Basically, this results in no output at all, although the activation of the packages creates a little output. Basic results can be obtained by simply calling the object:
> model.01 Linear mixed-effects model fit by REML Data: Exam Log-restricted-likelihood: -4663.8 Fixed: normexam ~ standLRT (Intercept) standLRT -0.01164834 0.55653379 Random effects: Formula: ~standLRT | school Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 0.3034980 (Intr) standLRT 0.1223499 0.494 Residual 0.7440699 Number of Observations: 4059 Number of Groups: 65
This gives a first impression of the estimated model. But, there is more. To obtain an idea of the elements that are actually stored inside the model, we use the names() functions, that gives us the names of all the elements of the model.
The output below shows that our model.01 contains seventeen elements. For reasons of space, only some will be described. ‘Contrasts’ contains information on the way categorical variables were handled, ‘coefficients’ contains the model-parameters, in ‘call’ the model formula is stored and in ‘data’ even the original data is stored.
> names(model.01)  "modelStruct" "dims" "contrasts" "coefficients"  "varFix" "sigma" "apVar" "logLik"  "numIter" "groups" "call" "terms"  "method" "fitted" "residuals" "fixDF"  "data" > model.01$method  "REML" > model.01$logLik  -4663.8
In the syntax above two specific elements of the model were requested: the estimation method and the loglikelihood. This is done by sub-setting the model using the $-sign after which the desired element is placed. The output tells us that model.01 was estimated using Restricted Maximum Likelihood and that the loglikelihood is -4663.8 large.
All the information we could possibly want is stored inside the models, as we have seen. In order to receive synoptic results, many functions exist to extract some of the elements from the model and present them clearly. The most basic of these extractor-functions is probably summary():
> summary(model.01) Linear mixed-effects model fit by REML Data: Exam AIC BIC logLik 9339.6 9377.45 -4663.8 Random effects: Formula: ~standLRT | school Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 0.3034980 (Intr) standLRT 0.1223499 0.494 Residual 0.7440699 Fixed effects: normexam ~ standLRT Value Std.Error DF t-value p-value (Intercept) -0.0116483 0.04010986 3993 -0.290411 0.7715 standLRT 0.5565338 0.02011497 3993 27.667639 0.0000 Correlation: (Intr) standLRT 0.365 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.8323045 -0.6316837 0.0339390 0.6834319 3.4562632 Number of Observations: 4059 Number of Groups: 65
The last extractor function that will be shown here is anova. This is a very general function that can be used for a great variety of models. When it is applied to a multilevel model, it results in a basic test for statistical significance of the model-parameters, as is showed below.
model.02 <- lme(fixed = normexam ~ standLRT, data = Exam,
random = ~ 1 | school)
In the syntax above an additional model is estimated that is very similar to our model.01, but does not have a random slope. It is stored in the object model.02. This is done to show that it is possible to test whether the random slope model fits better to the data than the fixed slope model. The output below shows that this is indeed the case.
> anova(model.01) numDF denDF F-value p-value (Intercept) 1 3993 124.3969 <.0001 standLRT 1 3993 765.4983 <.0001 > > model.02 <- lme(fixed = normexam ~ standLRT, data = Exam, + random = ~ 1 | school) > > anova(model.02,model.01) Model df AIC BIC logLik Test L.Ratio model.02 1 4 9376.765 9401.998 -4684.383 model.01 2 6 9339.600 9377.450 -4663.800 1 vs 2 41.16494 p-value model.02 model.01 <.0001
- - -- --- ----- --------
- Discuss this article and pose additional questions in the R-Sessions Forum
- Find the original article embedded in the manual.
- - -- --- ----- --------
R-Sessions is a collection of manual chapters for R-Project, which are maintained on Curving Normality. All posts are linked to the chapters from the R-Project manual on this site. The manual is free to use, for it is paid by the advertisements, but please refer to it in your work inspired by it. Feedback and topic requests are highly appreciated.
-------- ----- --- -- - -