R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
Copyright (C) 2024 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.11/
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.11)  multilevel model (Multinomial) 
Estimation algorithm:  IGLS MQL1        Elapsed time : 1.71s 
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.11)  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 : 1.16s 
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.98547379251299 -1.58839309553552 -2.57776923583701 2.15092628532752 0.706366979409297 0.726337603658121 2.68998713122079 0.686672926785606 1.06137891077163 2.65831722355658 0.24476941526159 1.12548085973048 '_FP_b'
JOIN 0.0984579919386108 0.000350943628532392 0.0153143182522688 -0.000948235712980013 -0.00246431519677697 0.0287645963378764 -0.0916258236583583 0.00155143638416767 0.00139056280785581 0.114998743303964 0.00155211785662298 -0.00977370567900764 0.00163478788398753 -0.00434100744782334 0.0206047783029191 0.00139072774343673 0.00163474004823583 -0.0235020874874748 -0.00374332766543245 -0.00447292047185663 0.0472066258234283 -0.0917621223323957 0.00154931865277692 0.00139678879244172 0.091820032960035 -0.00154489097828829 -0.00139790874128919 0.109732113389628 0.00155100920921958 -0.0098762673933995 0.00165749319512583 -0.00154465899016308 0.0099120578886674 -0.00166277797978881 -0.00475956291635228 0.023064249697185 0.00139699744774826 0.00165733927273109 -0.0236356573535804 -0.00139789450430829 -0.00166280201994342 0.0236804918556005 -0.00427640528015518 -0.00483595653275816 0.0452120628251993 -0.0917786094480177 0.00155338936713243 0.00140086875110016 0.0917312240593307 -0.00155059781882863 -0.00139592732735957 0.0918419328627452 -0.0015489019128833 -0.00140098659139409 0.0989404429078663 0.00155524086242794 -0.00988214809182161 0.00165791110167581 -0.00155021375495441 0.00984472647774262 -0.0016488771027395 -0.00154863193675781 0.00992750461550272 -0.00166660993414727 -0.00278517724710795 0.0170898767037855 0.0014009637155459 0.00165773998079982 -0.0236489334145858 -0.0013959273494455 -0.00164892231132227 0.0235978071821534 -0.00140100148206733 -0.00166665074620085 0.0237054905752506 -0.0026211098955188 -0.00289410411986475 0.0315664380764777 '_FP_v'
JOIN 0.349472162721045 0.110516126774995 0.288742440107789 0.0278150292934682 -0.0408544850867863 0.260258795381531 1 1 '_RP_b'
JOIN 0.0126188008421577 0.0018479913140807 0.00483967811472231 0.00022803970680508 0.0013625631823983 0.00705613326423024 -2.67042572956688e-05 -0.000775654395265686 -0.000224304553086712 0.00527751320288938 -2.42629289445028e-05 -0.000126776505265498 -0.00114207720610887 0.000785977397507928 0.0040508846127766 -9.83268997376156e-06 2.80703694932664e-06 0.000184367089551809 -1.14780530328208e-05 -0.00127678793067642 0.00882934915327953 3.57350282649975e-36 1.01915360216292e-35 -1.31655367292096e-36 -3.17383474722018e-37 8.14088383073235e-36 -5.07813559555228e-36 1.66810551381706e-34 0 4.51389830715758e-36 0 3.19734463423662e-36 3.04122429076804e-36 0 5.68041862312705e-35 4.19965179637222e-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.11)  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 : 1.83s 
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.81    0.34    9.87