R version 4.6.1 (2026-06-24 ucrt) -- "Happy Hop"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ############################################################################
> #     MLwiN User Manual
> #
> # 10  Multinomial Logistic Models for Unordered Categorical Responses . .145
> #
> #     Rasbash, J., Steele, F., Browne, W. J. and Goldstein, H. (2012).
> #     A User's Guide to MLwiN, v2.26. Centre for Multilevel Modelling,
> #     University of Bristol.
> ############################################################################
> #     R script to replicate all analyses using R2MLwiN
> #
> #     Zhang, Z., Charlton, C., Parker, R, Leckie, G., and Browne, W.J.
> #     Centre for Multilevel Modelling, 2012
> #     http://www.bristol.ac.uk/cmm/software/R2MLwiN/
> ############################################################################
> 
> library(R2MLwiN)

R2MLwiN: A package to run models implemented in MLwiN from R
Copyright 2013-2025 Zhengzheng Zhang, Christopher M. J. Charlton, Richard M. A. Parker, William J. Browne and George Leckie

Support provided by the Economic and Social Research Council (ESRC)
(Grants RES-149-25-1084, RES-576-25-0032 and ES/K007246/1)

To cite R2MLwiN in publications use:

  Zhengzheng Zhang, Richard M. A. Parker, Christopher M. J. Charlton,
  George Leckie, William J. Browne (2016). R2MLwiN: A Package to Run
  MLwiN from within R. Journal of Statistical Software, 72(10), 1-43.
  doi:10.18637/jss.v072.i10

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{R2MLwiN}: A Package to Run {MLwiN} from within {R}},
    author = {Zhengzheng Zhang and Richard M. A. Parker and Christopher M. J. Charlton and George Leckie and William J. Browne},
    journal = {Journal of Statistical Software},
    year = {2016},
    volume = {72},
    number = {10},
    pages = {1--43},
    doi = {10.18637/jss.v072.i10},
  }

The MLwiN_path option is currently set to C:/Program Files/MLwiN v3.18/
To change this use: options(MLwiN_path="<path to MLwiN>")
> # MLwiN folder
> mlwin <- getOption("MLwiN_path")
> while (!file.access(mlwin, mode = 1) == 0) {
+   cat("Please specify the root MLwiN folder or the full path to the MLwiN executable:\n")
+   mlwin <- scan(what = character(0), sep = "\n")
+   mlwin <- gsub("\\", "/", mlwin, fixed = TRUE)
+ }
> options(MLwiN_path = mlwin)
> 
> # Change contrasts if wish to avoid warning indicating that, by default,
> # specified contrasts for ordered predictors will be ignored by runMLwiN
> # (they will be fitted as "contr.treatment" regardless of this setting). To
> # enable specified contrasts, set allowcontrast to TRUE (this will be the
> # default in future package releases). NB at the end of this script, the
> # specification for contrasts is changed back.
> my_contrasts <- options("contrasts")$contrasts
> options(contrasts = c(unordered = "contr.treatment",
+                       ordered = "contr.treatment"))
> 
> # As an alternative to changing contrasts, can instead use C() to specify
> # contrasts for ordered predictors in formula object, e.g.:
> 
> # (mymodel1 <- runMLwiN(logit(use4) ~ 1 + C(lc, "contr.treatment"),
> #                       D = "Unordered Multinomial",
> #                       data = bang,
> #                       allowcontrast = TRUE))
> 
> # 10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .145
> 
> data(bang, package = "R2MLwiN")
> 
> addmargins(with(bang, table(use4)))
use4
           Sterilization Modern_reversible_method       Traditional_method 
                     302                      555                      282 
 Not_using_contraception                      Sum 
                    1728                     2867 
> 
> # 10.2 Single-level multinomial logistic regression . . . . . . . . . . .146
> 
> 
> # 10.3 Fitting a single-level multinomial logistic model in MLwiN . . . .147
> 
> addmargins(with(bang, table(lc, use4)))
              use4
lc             Sterilization Modern_reversible_method Traditional_method
  None                    12                      134                 44
  One_child               52                      137                 45
  Two_children            69                      107                 51
  Three_plus             169                      177                142
  Sum                    302                      555                282
              use4
lc             Not_using_contraception  Sum
  None                             584  774
  One_child                        283  517
  Two_children                     234  461
  Three_plus                       627 1115
  Sum                             1728 2867
> 
> bang$use4 <- relevel(bang$use4, 4)
> 
> (mymodel1 <- runMLwiN(logit(use4) ~ 1 + lc, D = "Unordered Multinomial", data = bang))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8
iteration 9

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.18)  multilevel model (Multinomial) 
Estimation algorithm:  IGLS MQL1        Elapsed time : 2.59s 
Number of obs:  2867 (from total 2867)        The model converged after 10 iterations.
Log likelihood:      NA 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use4) ~ 1 + lc
Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                                             Coef.   Std. Err.        z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept_Sterilization                   -3.88498     0.29093   -13.35   1.127e-40   ***     -4.45519    -3.31477 
Intercept_Modern_reversible_method        -1.47205     0.09500   -15.50   3.729e-54   ***     -1.65825    -1.28586 
Intercept_Traditional_method              -2.58570     0.15523   -16.66   2.674e-62   ***     -2.88994    -2.28146 
lcOne_child_Sterilization                  2.19120     0.32559     6.73   1.696e-11   ***      1.55306     2.82934 
lcOne_child_Modern_reversible_method       0.74692     0.13767     5.43    5.78e-08   ***      0.47709     1.01674 
lcOne_child_Traditional_method             0.74735     0.22004     3.40   0.0006829   ***      0.31607     1.17862 
lcTwo_children_Sterilization               2.66465     0.31885     8.36   6.433e-17   ***      2.03971     3.28959 
lcTwo_children_Modern_reversible_method    0.69036     0.14556     4.74   2.108e-06   ***      0.40507     0.97565 
lcTwo_children_Traditional_method          1.06313     0.21475     4.95   7.399e-07   ***      0.64223     1.48403 
lcThree_plus_Sterilization                 2.57436     0.30267     8.51   1.808e-17   ***      1.98114     3.16759 
lcThree_plus_Modern_reversible_method      0.20768     0.12545     1.66     0.09782   .       -0.03819     0.45355 
lcThree_plus_Traditional_method            1.10103     0.17933     6.14   8.277e-10   ***      0.74954     1.45251 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
             Coef.   Std. Err. 
bcons_1    1.00000     0.00000 
bcons_2    1.00000     0.00000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> cat(paste("Pr(y = 1) =", round(exp(mymodel1@FP["FP_Intercept_Sterilization"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + 
+   exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
Pr(y = 1) = 0.0155 
> cat(paste("Pr(y = 2) =", round(exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + 
+   exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
Pr(y = 2) = 0.1731 
> cat(paste("Pr(y = 3) =", round(exp(mymodel1@FP["FP_Intercept_Traditional_method"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + 
+   exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
Pr(y = 3) = 0.0568 
> cat(paste("Pr(y = 4) =", round(1/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + 
+   exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
Pr(y = 4) = 0.7545 
> 
> # 10.4 A two-level random intercept multinomial logistic regression model 154
> 
> # 10.5 Fitting a two-level random intercept model . . . . . . . . . . . .155
> 
> (mymodel2 <- runMLwiN(logit(use4) ~ 1 + lc + (1 | district), D = "Unordered Multinomial", data = bang))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8
iteration 9

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.18)  multilevel model (Multinomial) 
          N min     mean max N_complete min_complete mean_complete max_complete
district 60   3 47.78333 173         60            3      47.78333          173
Estimation algorithm:  IGLS MQL1        Elapsed time : 2.1s 
Number of obs:  2867 (from total 2867)        The model converged after 10 iterations.
Log likelihood:      NA 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use4) ~ 1 + lc + (1 | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                                             Coef.   Std. Err.        z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept_Sterilization                   -3.98547     0.31378   -12.70   5.802e-37   ***     -4.60047    -3.37048 
Intercept_Modern_reversible_method        -1.58839     0.12375   -12.84   1.039e-37   ***     -1.83094    -1.34585 
Intercept_Traditional_method              -2.57777     0.16960   -15.20   3.591e-52   ***     -2.91018    -2.24536 
lcOne_child_Sterilization                  2.15093     0.33911     6.34   2.257e-10   ***      1.48627     2.81558 
lcOne_child_Modern_reversible_method       0.70637     0.14354     4.92   8.614e-07   ***      0.42503     0.98771 
lcOne_child_Traditional_method             0.72634     0.21727     3.34   0.0008288   ***      0.30049     1.15218 
lcTwo_children_Sterilization               2.68999     0.33126     8.12   4.642e-16   ***      2.04073     3.33924 
lcTwo_children_Modern_reversible_method    0.68667     0.15187     4.52   6.141e-06   ***      0.38901     0.98433 
lcTwo_children_Traditional_method          1.06138     0.21263     4.99   5.987e-07   ***      0.64463     1.47813 
lcThree_plus_Sterilization                 2.65832     0.31455     8.45   2.883e-17   ***      2.04181     3.27482 
lcThree_plus_Modern_reversible_method      0.24477     0.13073     1.87     0.06116   .       -0.01145     0.50099 
lcThree_plus_Traditional_method            1.12548     0.17767     6.33   2.378e-10   ***      0.77726     1.47371 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                                                                         Coef.   Std. Err. 
var_Intercept_Sterilization                                            0.34947     0.11233 
cov_Intercept_Sterilization_Intercept_Modern_reversible_method         0.11052     0.06957 
var_Intercept_Modern_reversible_method                                 0.28874     0.08400 
cov_Intercept_Sterilization_Intercept_Traditional_method               0.02782     0.07265 
cov_Intercept_Modern_reversible_method_Intercept_Traditional_method   -0.04085     0.06365 
var_Intercept_Traditional_method                                       0.26026     0.09396 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
             Coef.   Std. Err. 
bcons_1    1.00000     0.00000 
bcons_2    1.00000     0.00000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> (mymodel3 <- runMLwiN(logit(use4) ~ 1 + lc + (1 | district), D = "Unordered Multinomial", estoptions = list(nonlinear = c(1, 
+   2), startval = list(FP.b = mymodel2@FP, FP.v = mymodel2@FP.cov, RP.b = mymodel2@RP, RP.v = mymodel2@RP.cov), resi.store = TRUE), 
+   data = bang))
/nogui option ignored
ECHO    0


Echoing is ON
BATC 1

Batch mode is ON
MAXI 2
STAR
iteration 0
iteration 1

Convergence not achieved
JOIN -3.985473792513 -1.58839309553552 -2.57776923583701 2.15092628532752 0.706366979409298 0.726337603658121 2.68998713122079 0.686672926785604 1.06137891077163 2.65831722355658 0.244769415261588 1.12548085973048 '_FP_b'
JOIN 0.0984579919386097 0.000350943628532408 0.0153143182522688 -0.00094823571297995 -0.00246431519677699 0.0287645963378763 -0.0916258236583571 0.00155143638416767 0.00139056280785575 0.114998743303963 0.00155211785662299 -0.00977370567900767 0.00163478788398753 -0.00434100744782335 0.0206047783029191 0.00139072774343668 0.00163474004823583 -0.0235020874874747 -0.00374332766543239 -0.00447292047185663 0.0472066258234282 -0.0917621223323945 0.00154931865277692 0.00139678879244166 0.0918200329600337 -0.0015448909782883 -0.00139790874128913 0.109732113389627 0.00155100920921961 -0.00987626739339952 0.00165749319512583 -0.0015446589901631 0.00991205788866742 -0.00166277797978882 -0.0047595629163523 0.023064249697185 0.0013969974477482 0.00165733927273109 -0.0236356573535804 -0.00139789450430823 -0.00166280201994342 0.0236804918556004 -0.00427640528015512 -0.00483595653275815 0.0452120628251992 -0.0917786094480166 0.00155338936713242 0.00140086875110009 0.0917312240593294 -0.00155059781882864 -0.00139592732735951 0.0918419328627439 -0.00154890191288333 -0.00140098659139403 0.0989404429078651 0.00155524086242795 -0.00988214809182163 0.00165791110167582 -0.00155021375495442 0.00984472647774265 -0.0016488771027395 -0.00154863193675781 0.00992750461550274 -0.00166660993414727 -0.00278517724710795 0.0170898767037855 0.00140096371554584 0.00165773998079983 -0.0236489334145857 -0.00139592734944544 -0.00164892231132228 0.0235978071821533 -0.00140100148206727 -0.00166665074620085 0.0237054905752505 -0.00262110989551873 -0.00289410411986475 0.0315664380764776 '_FP_v'
JOIN 0.349472162721046 0.110516126774996 0.288742440107788 0.0278150292934681 -0.0408544850867867 0.260258795381532 0.999999999999999 1 '_RP_b'
JOIN 0.0126188008421577 0.00184799131408066 0.00483967811472228 0.000228039706805069 0.00136256318239826 0.00705613326423017 -2.67042572956589e-05 -0.000775654395265685 -0.000224304553086707 0.00527751320288938 -2.4262928944501e-05 -0.000126776505265492 -0.00114207720610887 0.000785977397507909 0.00405088461277658 -9.83268997376165e-06 2.8070369493255e-06 0.000184367089551809 -1.14780530328128e-05 -0.00127678793067642 0.00882934915327953 -2.40741243048404e-35 1.50463276905253e-35 -1.20370621524202e-35 4.36343503025233e-35 3.05525675458098e-35 -2.40741243048404e-35 7.50651061212891e-34 -3.00926553810506e-36 1.31655367292096e-36 -3.00926553810506e-36 4.70197740328915e-36 4.25804461486142e-36 0 1.06131751533845e-34 1.66493139320476e-35 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
iteration 8

Convergence achieved
ECHO 0
Execution completed


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.18)  multilevel model (Multinomial) 
          N min     mean max N_complete min_complete mean_complete max_complete
district 60   3 47.78333 173         60            3      47.78333          173
Estimation algorithm:  IGLS PQL2        Elapsed time : 2.61s 
Number of obs:  2867 (from total 2867)        The model converged after 9 iterations.
Log likelihood:      NA 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use4) ~ 1 + lc + (1 | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                                             Coef.   Std. Err.        z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept_Sterilization                   -4.22566     0.31875   -13.26   4.101e-40   ***     -4.85039    -3.60093 
Intercept_Modern_reversible_method        -1.74491     0.13201   -13.22   6.876e-40   ***     -2.00364    -1.48618 
Intercept_Traditional_method              -2.72149     0.17835   -15.26   1.432e-52   ***     -3.07106    -2.37193 
lcOne_child_Sterilization                  2.22291     0.33640     6.61   3.895e-11   ***      1.56359     2.88224 
lcOne_child_Modern_reversible_method       0.77319     0.14241     5.43   5.656e-08   ***      0.49408     1.05231 
lcOne_child_Traditional_method             0.74581     0.22596     3.30   0.0009648   ***      0.30293     1.18869 
lcTwo_children_Sterilization               2.81623     0.32948     8.55   1.259e-17   ***      2.17045     3.46200 
lcTwo_children_Modern_reversible_method    0.79624     0.15009     5.31   1.125e-07   ***      0.50208     1.09041 
lcTwo_children_Traditional_method          1.14329     0.21952     5.21   1.907e-07   ***      0.71304     1.57354 
lcThree_plus_Sterilization                 2.78984     0.31289     8.92   4.822e-19   ***      2.17658     3.40309 
lcThree_plus_Modern_reversible_method      0.33165     0.12958     2.56     0.01048   *        0.07768     0.58562 
lcThree_plus_Traditional_method            1.18728     0.18380     6.46    1.05e-10   ***      0.82704     1.54753 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                                                                        Coef.   Std. Err. 
var_Intercept_Sterilization                                           0.53329     0.15221 
cov_Intercept_Sterilization_Intercept_Modern_reversible_method        0.30899     0.09790 
var_Intercept_Modern_reversible_method                                0.38892     0.10566 
cov_Intercept_Sterilization_Intercept_Traditional_method              0.24327     0.09671 
cov_Intercept_Modern_reversible_method_Intercept_Traditional_method   0.13573     0.07809 
var_Intercept_Traditional_method                                      0.32534     0.11051 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
             Coef.   Std. Err. 
bcons_1    1.00000     0.00000 
bcons_2    1.00000     0.00000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> 
> mymodel3@RP["RP2_cov_Intercept_Sterilization_Intercept_Modern_reversible_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Sterilization"] * mymodel3@RP["RP2_var_Intercept_Modern_reversible_method"])
RP2_cov_Intercept_Sterilization_Intercept_Modern_reversible_method 
                                                         0.6784711 
> mymodel3@RP["RP2_cov_Intercept_Sterilization_Intercept_Traditional_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Sterilization"] * mymodel3@RP["RP2_var_Intercept_Traditional_method"])
RP2_cov_Intercept_Sterilization_Intercept_Traditional_method 
                                                   0.5840432 
> mymodel3@RP["RP2_cov_Intercept_Modern_reversible_method_Intercept_Traditional_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Modern_reversible_method"] * mymodel3@RP["RP2_var_Intercept_Traditional_method"])
RP2_cov_Intercept_Modern_reversible_method_Intercept_Traditional_method 
                                                               0.381576 
> 
> hipos <- rep(0, 2)
> hipos[1] <- which(levels(as.factor(bang$district)) == 56)
> hipos[2] <- which(levels(as.factor(bang$district)) == 11)
> 
> u0 <- mymodel3@residual$lev_2_resi_est_Intercept.Sterilization
> u0se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Sterilization)
> u0CI95 <- 1.96 * u0se
> u0rank <- rank(u0)
> u0rankhi <- u0 + u0CI95
> u0ranklo <- u0 - u0CI95
> u0rankno <- order(u0rank)
> plot(1:60, u0[u0rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
> points(1:60, u0rankhi[u0rankno], pch = 24, bg = "grey")
> points(1:60, u0ranklo[u0rankno], pch = 25, bg = "grey")
> for (i in 1:60) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
> for (i in 1:2) points(x = which(u0rankno == hipos[i]), y = u0[u0rankno[which(u0rankno == hipos[i])]], pch = 22, bg = i + 
+   1)
> 
> u1 <- mymodel3@residual$lev_2_resi_est_Intercept.Modern_reversible_method
> u1se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Modern_reversible_method)
> u1CI95 <- 1.96 * u1se
> u1rank <- rank(u1)
> u1rankhi <- u1 + u1CI95
> u1ranklo <- u1 - u1CI95
> u1rankno <- order(u1rank)
> plot(1:60, u1[u1rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u1 residual estimate")
> points(1:60, u1rankhi[u1rankno], pch = 24, bg = "grey")
> points(1:60, u1ranklo[u1rankno], pch = 25, bg = "grey")
> for (i in 1:60) lines(rep(i, 2), c(u1ranklo[u1rankno[i]], u1rankhi[u1rankno[i]]))
> for (i in 1:2) points(x = which(u1rankno == hipos[i]), y = u1[u1rankno[which(u1rankno == hipos[i])]], pch = 22, bg = i + 
+   1)
> 
> u2 <- mymodel3@residual$lev_2_resi_est_Intercept.Traditional_method
> u2se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Traditional_method)
> u2CI95 <- 1.96 * u2se
> u2rank <- rank(u2)
> u2rankhi <- u2 + u2CI95
> u2ranklo <- u2 - u2CI95
> u2rankno <- order(u2rank)
> plot(1:60, u2[u2rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u2 residual estimate")
> points(1:60, u2rankhi[u2rankno], pch = 24, bg = "grey")
> points(1:60, u2ranklo[u2rankno], pch = 25, bg = "grey")
> for (i in 1:60) lines(rep(i, 2), c(u2ranklo[u2rankno[i]], u2rankhi[u2rankno[i]]))
> for (i in 1:2) points(x = which(u2rankno == hipos[i]), y = u2[u2rankno[which(u2rankno == hipos[i])]], pch = 22, bg = i + 
+   1)
> 
> # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .159
> 
> # Addendum: changing contrasts back to pre-existing . . . . . . . . . . . NA
> 
> # Following re-specification of contrast settings towards the start of this
> # script, change contrasts back to pre-existing:
> options(contrasts = my_contrasts)
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   3.85    1.17   13.48