R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

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.

  Natural language support but running in an English locale

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-2017 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.05/
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)
> 
> 
> # 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.05)  multilevel model (Multinomial) 
Estimation algorithm:  IGLS MQL1        Elapsed time : 1.8s 
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.05)  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.17s 
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.05)  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.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.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
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   4.31    0.50   10.39