Practice session5HLM多层线性模型讲义
- 格式：doc
- 大小：46.50 KB
- 文档页数：4
R-practice session 5 CS&SS 560
Marijtje van Duijn
Winter 2006
The commands used in this session are available as R syntax file (Session5.R)
at the website.
Data input and preparation
We use the data as provided in example 12.1 of Snijders & Bosker. The
dependent variable is student evaluations of 51 teachers over 4 years. The data
are very unbalanced. Explanatory variables are teacher experience (coincides
with ‘time’) own evaluation and gender.
After getting the necessary libraries, download the data file SBbook12.csv from
the class website. Also get the file session5.r and execute the commands under
data preparation.
Note that a groupedData object ‘longdata’ is created.
Some more commands for investigating grouped data
We can extract some useful information about the grouping with the following
commands:
>table(getGroups(longdata))
>getGroups(longdata)
>unique(getGroups(longdata))
Another useful trick in this case may be to split the data in separate objects for
men and women. This is especially attractive for making pictures. (The same can
be done for the men; see the file session5.R)
>longdataf<-groupedData(prox_stu~occ|teacher,data=subset(data12,gender>0),
+labels = list(x="time", y="proximity score")) Estimating multilevel models for longitudinal data
One of the important questions in modeling multilevel data is how to treat ‘time’.
In this case with fixed occasions we have to choose what kind of ‘contrasts’ we
want to have. In the example in the book, a ‘dummy’ coding is used, giving
separate estimates for the means at each time point.
Here we will also explore other contrasts that are statistically equivalent, but lead
to different parameter estimates.
We will start with the ‘empty’ compound symmetry model (or random intercept
model) with equal means for all time points (model 1 in table 12.1).
>model1.1<-lme(prox_stu~1,data=longdata,random=~1|teacher)
>summary(model1.1)
We then estimate three different specifications for model 2: P for polynomial, F
for factor (estimating differences with respect to the first time point), and D, for
dummy. Look carefully at the model specification, ‘transforming’ occasion (the
time variable).
Have a close look at the differences and similarities of the results! Also estimate
the dummy coded model with maximum likelihoods. Due to the small groups,
there is a larger difference between the results obtained with ML and REML than
we have seen sofar.
>model1.2P<-lme(prox_stu~ordered(occ),data=longdata,random=~1|teacher)
>model1.2F<-lme(prox_stu~factor(occ),data=longdata,random=~1|teacher)
>model1.2D<-lme(prox_stu~factor(occ)-1,data=longdata,random=~1|teacher)
>model1.2DML<-lme(prox_stu~factor(occ)-1, data=longdata, random=~1|
+teacher, method="ML")
Finally we look at what happens if we forget to indicate what kind of contrast we
want for occasion:
>model1.2L<-lme(prox_stu~occ,data=longdata,random=~1|teacher)
The next step is to define a random effect for occasion, which implies that the
variance increases over time. Note the different treatment of occasion in the fixed
and in the random part.
>model2.3D<-lme(prox_stu~factor(occ)-1, data= longdata, random = ~1+occ|
+teacher)
>model2.3DML<-lme(prox_stu~factor(occ)-1, data = longdata, random = ~1 +
+occ| teacher,method="ML")
In the next model (4) gender is added, as well as its interaction with experience
(measured through occasion). Note again the difference in treatment of occasion.
The other peculiar thing is that we have to ‘subtract’ occ (automatically added
when we specify the interaction between occ and gender) because it was
included in the model in a different way – through the dummies belonging to
factor(occ).
>model2.4D<-lme(prox_stu~factor(occ)-1 + gender + occ*gender-occ, data =
+longdata, random = ~1+occ|teacher)
A different but somewhat shorter way to arrive at this model is using the update
function:
>model2.4D<-update(model2.3D,~.+gender+occ*gender-occ)
We now want to estimate the fully multivariate model. In a ‘technical’ way, this is
not a multilevel anymore, and that is why we cannot use the lme() function, which
always estimates at least one variance component at level 1. But we can use the
gls() function (gls = generalized least squares), and specify that we want a
‘general’ correlation structure through the use of corSymm for our data (that are
still organized in a multilevel way). Here we use form=~1|teacher, but we could
also use form=~ordered(occ)|teacher, showing a little more clearly that we want
to have general correlation structure for occasion. The weights part specifies that
the variances are different for each occasion. I agree that this is not simple to
