Practice session3HLM多層線性模型講義

1、R-practice session 3cs&ss 560Marijtje van Duijn Winter 2006The commands used in this session are available as R syntax file (Session3.R) at the website.Data input and preparationWe continue with the data used in Snijders & Bosker. For a description see Example 4.1 (p. 46).The data are the same as in

2、 the previous session, except for the treatment of the variable IQc. This variable is the grandmean centered IQ. Then the same procedure is followed as before, resulting in two variables, a centered and an uncentered one, measuring IQ in deviation to its school meanGo the file practice session3.r an

3、d execute the commands under data preparation. One other new variable will be produced, called z2, that is the deviation of the (school variable) groupsize with respect to the overall mean groupsize (equal to 23.1, approximately). Note that the groupsize can be larger than the number of students per

4、 school in the dataset due to missing pupils.Comparing (nested) models and estimating a random slope modelMake sure the library nlme is attached.First we will give a command to format the output a little nicer, restricting the number of decimal places.options(digits=5)This seems to be the best choic

5、e, although R does not always obey” this command. If you like, play around a bit with other choices for the number of digits. In the same options command you can decide to turn of the use of stars to signal significance - in the eyes of many statisticians a despicable practiceoptions(show.signif.sta

6、rs=FALSE)We will first estimate the models from Tables 4.2, 4.4, and 4.5 again.Check that we are now reproducing the output in the tables, and look at the funny p-values.model.4.2summary(model.4.2)VarCorr(model.4.2) model.4.4summary(model.4.4)VarCorr(model.4.4)We can do a deviance test or likelihood

7、 ratio test by giving the commandanova(model.4.4,model.4.2)Here, however, it is not necessary to do a deviance test, because the models only differ in their fixed effects: model 4.4 is model.4.2 + meanIQc. It is then preferable (from a statistical point of view) to use a t-test as reported in the su

8、mmary output, or an F-test (remember its relation to the t-test?) that can be obtained with the commandanova(model.4.4)Because of the small p-values we dont see any difference in the test, but the results can be different - especially when testing multiple fixed effects, or categorical variables.Bec

9、ause going from model.4.2. to model.4.4 means the addition of the variable meanIQc, R allows (for lme() but also for lm() the commandmodel.4.4model.4.4model.5.1summary(model.5.1)VarCorr(model.5.1)We could also get model.5.1 by updating model.4.4 (and changing the estimation method), with the command

10、smodel.5.1model.5.1mode/.4.4anova(model.5.1,model.4.4)Here we cannot use a t-test or F-test for testing the significance of the addition of a random effect, and is the likelihood ratio test the easiest.Note that the likelihood ratio test is done with 2 degrees of freedom. Two extra parameters were e

11、stimated in model.5.1 compared to model.4.4: the variance of the slope and the covariance between slope and intercept (%2and tq1).It is the case that the p-values with the chi-square approximation are conservative for nested multilevel models - the reason for this being that the test considers a nul

12、l hypothesis (a variance equal to zero) on the boundary of the parameter space (a variance cannot be smaller than zero). A solution is to use a bootstrap procedure to obtain a better estimate of the p-value. This procedure can be done in the nlme package (see Pinheiro & Bates, 2000, section 2.4) but

13、 for most purposes the conservative chi-square test seems sufficient.To investigate whether the random slope can be explained using the variable z2, the centered group size, the model of table 5.2 is estimatedmodel.5.2summary(model.5.2)An interaction can be specified with the * sign, but also with t

14、he : sign (as used in the output). Lets see whether we can do without the random slopes now.Unfortunately, the update function does not work for the random part, so we will specify model.5.2asummary(model.5.2a)anova(model.5.2,model.5.2a)So the conclusion is.?I leave the model in table 5.3 as an exer

15、cise.The influence of centeringLast week we saw that it matters for the intercept whether or not we center an explanatory variable. Lets investigate what the effect is on the variance structure. From the theoretical derivations of the between- and within-group variances, we know that this should mak

16、e a difference.First we estimate the model in table 5.1 not using the centered IQ. Here we see the expected effect on the (level 2) intercept variance.model.5.1.asummary(model.5.1.a)VarCorr(model.5.1.a)anova(model.5.1.a,model.5.1)Next we estimate the model in table5.1 using IQ in deviation to its gr

17、oup (school) mean. And now we observe not a difference in the variance structure, but in the estimates of the fixed effects.model.5.1.bsummary(model.5.1.b)VarCorr(model.5.1.b)anova(model.5.1.b,model.5.1)Compare this to model.4.9.model.4.9summary(model.4.9)VarCorr(model.4.9)I admit that the differenc

18、es are not so big, but that is because of the large dataset, with many schools and many students!Visualizing random slopesWe can get a visual impression of the size of the variation by making a graph like figure 5.1. We will draw a sample of 15 school interecepts and slopes from the population with the bivariate normal distributism parameters as estimated in model.5.1. We can draw from a multivariate normal distribution, using the function mvrnorm(), available in the MASS library. We do this after having specified its parameters, taken from the output. (We co