HLM多层线性模型讲义

R-practice session 7 CS&SS 560

Marijtje van Duijn

Winter 2006 The commands used in this session are available as R syntax file (Session7.R) at the website.

Data input and preparation

We continue with the data used in Snijders & Bosker. For a description see Example 4.1 (p. 46).

Download the data file SBbookR2.csv from the class website. Also get the file session7.r and execute the commands under data preparation (also used in practice session 4).

Modeling heteroscedasticity at level 1

We will first estimate the models as specified in table 8.1.

>model8.11ML<-lme(testscore~ IQc+sesc+gender+meanIQc, random=~1+IQc| +schoolNR, data= datasb7,method="ML")

>model8.12ML<-lme(testscore~ IQc+sesc+gender+meanIQc, random = ~1+IQc| +schoolNR, data=datasb7, weights=varIdent(form=~1|gender), method="ML") Although it is immediately clear that the second model in table 8.1 is not an improvement of the first model:

>anova(model8.12ML,model8.11ML)

we will check that it is the exact same model, although differently parameterized. For almost all of the variance components this is clear from

>VarCorr(model8.12ML)

and remember that the covariance between the intercept and the IQ slope variance can be computed by multiplying the reported correlation by the standard deviations.

We now have to check that the -1.21 covariance of the gender effect is retrieved in the lme model calculation. The estimated variance for girls is 36.30 (38.72-

2*1.21) from table 8.1

This number is obtained in our analysis by multiplying 38.72 by 0.968^2. Check that this indeed correct.

Note that the weights command does not work with lmer(). I‟m pretty sure it must be possible to specify an equivalent model, but I have not found out how (yet).

Exploring the level 1 error structure

In chapter 8 it may have seemed somewhat surprising to assume different residual variances for boys and girls, even though we could immediately test whether the assumption was useful and drop it. In chapter 9 where the assumptions of the hierarchical linear model are investigated some methods for inspecting the level 1 residuals and deciding whether they could be related to level 1 explanatory variables (if available) are given.

In section 9.4.1 a test for (level 1) heteroscedasticity is given, based on the OLS residuals within each group, after estimating a certain model (with only level 1 explanatory variables). We can use the lmList() method to produce the required input for this test statistic H. See the commands in session7.r. What is done is that the var iances of …eligible‟ groups (with enough degrees of freedom, defined as the number of observations – the number of estimated parameters, here at least 10) are …normalized‟. The sum of the squares of these normalized variances are chi-square distributed under the null hypothesis of equal variances with df equal to the number of contributing groups minus 1.

The resulting H test can be performed and interpreted easily. Not surprisingly, it is sensitive to the normality assumption, so it can be wise to check the normality of the dj-elements, by a simple qq plot. It is also advisable to check the normality of the within-group residuals (as we did before) to make sure that the variances are an adequate summary of the within-group variability to start with.

Further inspection of level 1 residuals

Although on the basis of the test statistic H there is actually no good reason for further inspection of the level 1 residuals, the procedure to do this is illustrated in section 9.5. The idea is to see how the “smoothed” residuals are associated with the explanatory variables. Smoothing here means summarizing the residuals for certain (categorical) values of the explanatory variable. I will follow example 9.2. See session7.r for how the plots are obtained.

The curvature in the obtained graph is weaker than what is shown in figure 9.1. This may be due to the definition of the IQ groups. My conclusion here would be that there is not enough evidence for exploring a nonlinear relation between testscore and IQ.

Obviously the same thing can be done for SES, and for modified models.