note MLwiN macro for producing results in examples in Chapter 10,
note Sections 10.4 - 10.6, of Snijders & Bosker, 2nd edition.
echo 1
logo 0
logo mlbook1_ch10.log
clear
note Enter data.
obey mlbook2_r.obe
note IQ_verb and SES already have been centered.
note original means were 11.8668 for IQ_verb and 27.7273 for ses.
note Basic model definition
pref 0
post 0
iden 1 c2
iden 2 c1
expl 1 "cons"
setv 1 "cons"
setv 2 "cons"
resp "lang_post"
maxi 100
batch 1
tole 3
expl 1 "IQ_verb"
expl 1 "mean_iqv"
setv 2 "IQ_verb"
calc c14 = 'IQ_verb'*'ses'
name c14 "iqv*ses"
expl 1 c14
dumm "denomina" c15-c19
name c15 'd1' c16 'd2' c17 'd3' c18 'd4' c19 'd5'
calc c20 = c5 * c10
calc c21 = c5 * c9
calc c22 = c4 * c10
calc c23 = c4 * c9
calc c24 = c9 * c10
name c20 "iqv*m_iqv" c21 "iqv*m_ses"
name c22 "ses*m_iqv" c23 "ses*m_ses"
name c24 "m_iqv*m_ses"
calc c20 = c5 * c10
calc c21 = c5 * c9
calc c22 = c4 * c10
calc c23 = c4 * c9
calc c24 = c9 * c10
name c20 "iqv*m_iqv" c21 "iqv*m_ses"
name c22 "ses*m_iqv" c23 "ses*m_ses"
name c24 "m_iqv*m_ses"
expl 1 c4 c6 c9 c24
calc c13 = "IQ_verb" - "mean_iqv"
name c13 "IQ_dev"
calc c37 = 'iq_verb'^2
name c37 'iqv^2'
calc c38 = 'iq_verb'*abso('iq_verb')
calc c38 = (c38 + 'iqv^2')/2
name c38 "iqv^2+"
calc c39 = c37 - c38
name c39 "iqv^2-"
note The current specification is as Model 1 in Table 8.1.
start
next
fixe
rand
like
note Now exclude the level-two variables.
expl 0 c9 c10 c24
next
fixe
rand
like
note Example 10.1, heteroscedasticity test:
obey hetsced1.obe
note Put denominations of schools in C212:
take c1 c8 c211 c212
note Only for schools used in heteroscedasticity test:
merge c211 c212 c201 c213
note Take average of lang_post in used schools:
mlav c1 c3 c214
take c1 c214 c211 c215
merge c211 c215 c201 c216
plot standardized within-school standard deviations against average lang_post:
plot c204 c216
note If you wish to see the values of
note identifier, d.f., s_j^2, d_j, denomination, average language test,
note then drop the "note" in the following command:
note write c201-c204 c213 c216
note Now also a plot of within-school standard deviations against average ses:
take c1 c9 c10 c211 c217 c218
merge c211 c217 c218 c201 c219 c220
plot c204 c219
note Example 10.2.
note Use macro for calculation of residuals:
set b11 10
obey residuals_1.obe
note Groups with too low residual d.f. are omitted from residuals.
note Residuals are now in C226, standardised residuals are in C229,
note for the reduced data set (in which small groups are omitted).
note Explanatory variables of the original data set are in Group 1,
note Explanatory variables of the reduced data set are in Group 4.
note In this case, Group 4 is C301-C304.
note You can find this out as follows:
note In MLwiN open the Data Manipulation - Groups window to see which are the
note variables in groups G1 and G4, and in the Data window
note you may have a look at their values to check the correspondence.
note Now make some plots of the residuals.
name c301 'riqv' c302 'riqv_x_ses' c303 'rses' c304 'rsex'
note First treat IQ and SES as continuous variables for smoothing:
note Smooth residuals as function of iq:
set b1 301
set b2 226
set b21 50
obey smooth_2.obe
note Variables & residuals can be written to file; this now is suspended by
note writing "note" before the dwrite command and filename.
note dwrite c404 c406
note riq.dat
note smooth them as function of SES:
set b1 303
obey smooth_2.obe
note dwrite c404 c406
note rses.dat
note Now treat IQ and SES as categorical variables for smoothing:
note First round them to integer values
calc c301 = 2*(c301 - 0.63)
calc c303 = c303 - .27
set b1 301
set b2 226
set b3 10
obey table_2.obe
note dwrite c411 c413 c415
note triq.dat
set b1 303
obey table_2.obe
note dwrite c411 c413 c415
note trses.dat
set b1 304
obey table_2.obe
note dwrite c411 c413 c415
note trsex.dat
note if you want to save the normal scores for standardized residuals
note as calculated by residuals_1.obe, drop the "note" before the next two lines:
note dwrite c229 c228
note rnorm0.dat
note Now extend fixed part with non-linear functions of IQ:
expl 1 'iqv^2-' 'iqv^2+'
start
fixe
rand
like
note Calculate residuals again for the new model:
obey residuals_1.obe
note The preceding plot is a
note normal probability plot of semi-standardized level-1 residuals,
note similar to Figure 10.3 but for a different model.
note To write the data for this plot to a file use the next two commands:
note dwrite c229 c228
note rnorm.dat
note Now the explanatory variables in the reduced data set
note are in C301 - C306.
note This can be checked by looking at the <> window.
note Give them corresponding names:
name c301 'c301' c302 'c302' c303 'c303' c304 'c304'
name c301 'riqv' c302 'riqv_x_ses' c303 'rses' c304 'rsex'
name c305 'riq^2-' c306 'riq^2+'
note smooth these residuals as function of IQ:
set b1 301
set b2 226
set b21 50
obey smooth_2.obe
note For writing the results to file:
note dwrite c401 c403
note riq1.dat
note smooth them as function of SES:
set b1 303
obey smooth_2.obe
note dwrite c404 c406
note rses1.dat
note Now the same for the standardized residuals:
note w.r.t. iq:
set b1 301
set b2 228
set b21 50
obey smooth_2.obe
note dwrite c404 c406
note riq1s.dat
note w.r.t. ses:
set b1 303
obey smooth_2.obe
note dwrite c404 c406
note rses1s.dat
note Now treat IQ and SES as categorical variables for smoothing:
note First transform them back to integer values
calc c301 = 2*(c301 - 0.63)
calc c303 = c303 - .27
note For IQ:
set b1 301
set b2 226
set b3 12
obey table_2.obe
note For producing the data for Figure 10.1 (left):
note Transform back to earlier IQ scale
note calc c411 = 0.5*c411 + 0.63
note dwrite c411 c413 c416 c417
note triq1.dat
set b1 303
calc c303 = 0.5*c303
obey table_2.obe
calc c303 = 2*c303
note For producing the data for Figure 10.1 (right):
note Transform back to earlier SES scale
note calc c411 = 2*c411 + 0.27
note dwrite c411 c413 c416 c417
note trses1.dat
note Extend the model with non-linear effect of IQ:
expl 1 'iqv^2+' 'iqv^2-'
start
fixe
rand
like
note Calculate residuals again for the new model:
set b11 10
obey residuals_1.obe
note Go through the same as before:
name c301 'c301' c302 'c302' c303 'c303' c304 'c304'
name c301 'riqv' c302 'riqv_x_ses' c303 'rses' c304 'rsex'
name c305 'riq^2-' c306 'riq^2+'
note Treat IQ and SES as categorical variables for smoothing:
note First transform them back to integer values
calc c301 = 2*(c301 - 0.63)
calc c303 = c303 - .27
note For IQ:
set b1 301
set b2 226
set b3 12
obey table_2.obe
note For producing the data for Figure 10.2 (left):
note Transform back to earlier IQ scale
note calc c411 = 0.5*c411 + 0.63
note dwrite c411 c413 c416 c417
note triq2.dat
set b1 303
calc c303 = 0.5*c303
obey table_2.obe
calc c303 = 2*c303
note For producing the data for Figure 10.2 (right):
note Transform back to earlier SES scale
note calc c411 = 2*c411 + 0.27
note dwrite c411 c413 c416 c417
note trses2.dat
note Data for Figure 10.3:
note dwrite c228 c229
note rnorm_2.dat
note Now go on to level-two residuals (Section 10.6)
note First the settings for calculation of posterior means:
expl 1 c9 c10 c24
next
rlev 2
rcov 1
rtyp 0
rout c231-c234
note Actually calculate the residuals:
resi
calc c235 = c231/sqrt(c233)
calc c236 = c232/sqrt(c234)
note c231 and c232 unstandardized residuals
note c233 and c234 diagnostic variances
note c235 and c236 standardized residuals
note Calculate group-level variables
take c1 c10 c237 c238
take c1 c9 c237 c239
name c238 "sch_iqv" c239 "sch_ses"
note plot posterior intercepts as a function of average IQV
plot c231 c238
note smoothed:
set b1 238
set b2 231
set b21 30
obey smooth_2.obe
note plot posterior intercepts as a function of average ses
plot c231 c239
note smoothed:
set b1 239
set b2 231
obey smooth_2.obe
note plot posterior slopes
plot c232 c238
set b1 238
set b2 232
set b21 30
obey smooth_2.obe
plot c232 c239
set b1 239
set b2 232
obey smooth_2.obe
(The lowess approximations in the pictures in the book were obtained using R.)
note Write data for Figure 10.6.
nsco c235 c240
nsco c236 c241
dwrite c240 c235
resi2_int_norm.txt
dwrite c241 c236
resi2_sl_norm.txt
note check that diagnostic variances plus comparative variances
note are equal to random-effect variances:
rtyp 1
rout c242-c245
resi
calc c246 = c244 + c233
calc c247 = c245 + c234
aver c246
aver c247