R version 4.5.0 (2025-04-11 ucrt) -- "How About a Twenty-Six"
Copyright (C) 2025 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-2024 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.14/
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.14)  multilevel model (Multinomial) 
Estimation algorithm:  IGLS MQL1        Elapsed time : 2.46s 
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.14)  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.15s 
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.98547379251303 -1.58839309553552 -2.577769235837 2.15092628532756 0.706366979409297 0.726337603658118 2.68998713122082 0.686672926785605 1.06137891077163 2.65831722355661 0.244769415261588 1.12548085973047 '_FP_b'
JOIN 0.0984579919386151 0.000350943628532463 0.0153143182522688 -0.000948235712980023 -0.00246431519677701 0.0287645963378762 -0.0916258236583623 0.00155143638416763 0.00139056280785581 0.114998743303968 0.00155211785662297 -0.00977370567900765 0.00163478788398754 -0.00434100744782333 0.0206047783029191 0.00139072774343675 0.00163474004823584 -0.0235020874874746 -0.00374332766543246 -0.00447292047185664 0.047206625823428 -0.0917621223323997 0.00154931865277688 0.00139678879244172 0.0918200329600388 -0.00154489097828828 -0.0013979087412892 0.109732113389632 0.00155100920921958 -0.0098762673933995 0.00165749319512585 -0.00154465899016307 0.00991205788866741 -0.00166277797978883 -0.00475956291635227 0.023064249697185 0.00139699744774828 0.0016573392727311 -0.0236356573535803 -0.0013978945043083 -0.00166280201994343 0.0236804918556003 -0.0042764052801552 -0.00483595653275818 0.0452120628251991 -0.0917786094480218 0.00155338936713238 0.00140086875110016 0.0917312240593345 -0.00155059781882863 -0.00139592732735958 0.0918419328627491 -0.0015489019128833 -0.00140098659139411 0.0989404429078702 0.00155524086242794 -0.00988214809182161 0.00165791110167583 -0.0015502137549544 0.00984472647774263 -0.00164887710273951 -0.0015486319367578 0.00992750461550272 -0.00166660993414728 -0.00278517724710793 0.0170898767037855 0.00140096371554592 0.00165773998079984 -0.0236489334145857 -0.00139592734944551 -0.00164892231132228 0.0235978071821533 -0.00140100148206734 -0.00166665074620086 0.0237054905752505 -0.00262110989551881 -0.00289410411986476 0.0315664380764775 '_FP_v'
JOIN 0.349472162721053 0.110516126774998 0.28874244010779 0.0278150292934678 -0.0408544850867873 0.26025879538153 1 1 '_RP_b'
JOIN 0.0126188008421579 0.00184799131408074 0.00483967811472236 0.000228039706805086 0.00136256318239832 0.00705613326423031 -2.67042572956823e-05 -0.000775654395265692 -0.000224304553086714 0.0052775132028894 -2.42629289445046e-05 -0.000126776505265505 -0.00114207720610888 0.000785977397507936 0.0040508846127766 -9.8326899737613e-06 2.80703694932822e-06 0.000184367089551809 -1.1478053032832e-05 -0.00127678793067642 0.00882934915327948 -1.20370621524202e-35 3.00926553810506e-36 -1.20370621524202e-35 7.52316384526264e-37 3.59260460970062e-37 -1.20370621524202e-35 1.38759218019312e-34 0 -1.88079096131566e-36 1.50463276905253e-36 -1.03443502872361e-36 -6.75266403249709e-37 0 -3.58236949089818e-35 1.04991294909306e-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.14)  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.69s 
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.22917     0.31898   -13.26   4.029e-40   ***     -4.85435    -3.60398 
Intercept_Modern_reversible_method        -1.74833     0.13233   -13.21   7.538e-40   ***     -2.00770    -1.48896 
Intercept_Traditional_method              -2.72372     0.17859   -15.25   1.606e-52   ***     -3.07374    -2.37370 
lcOne_child_Sterilization                  2.22552     0.33633     6.62   3.662e-11   ***      1.56634     2.88471 
lcOne_child_Modern_reversible_method       0.77615     0.14248     5.45   5.109e-08   ***      0.49690     1.05541 
lcOne_child_Traditional_method             0.74795     0.22597     3.31   0.0009332   ***      0.30505     1.19084 
lcTwo_children_Sterilization               2.82372     0.32943     8.57    1.02e-17   ***      2.17805     3.46938 
lcTwo_children_Modern_reversible_method    0.80305     0.15011     5.35   8.812e-08   ***      0.50884     1.09727 
lcTwo_children_Traditional_method          1.14890     0.21948     5.23   1.654e-07   ***      0.71873     1.57908 
lcThree_plus_Sterilization                 2.79461     0.31288     8.93   4.185e-19   ***      2.18138     3.40784 
lcThree_plus_Modern_reversible_method      0.33648     0.12963     2.60    0.009439   **       0.08242     0.59055 
lcThree_plus_Traditional_method            1.19088     0.18383     6.48   9.294e-11   ***      0.83057     1.55118 
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.54132     0.15367 
cov_Intercept_Sterilization_Intercept_Modern_reversible_method        0.31485     0.09887 
var_Intercept_Modern_reversible_method                                0.39253     0.10635 
cov_Intercept_Sterilization_Intercept_Traditional_method              0.24924     0.09762 
cov_Intercept_Modern_reversible_method_Intercept_Traditional_method   0.13947     0.07864 
var_Intercept_Traditional_method                                      0.32885     0.11115 
--------------------------------------------------------------------------------------------------- 
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.6830372 
> 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.5907354 
> 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.388203 
> 
> 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.67    1.28   12.93