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.13/
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.13) multilevel model (Multinomial)
Estimation algorithm: IGLS MQL1 Elapsed time : 3.85s
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.13) 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.81s
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.13) 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 : 4.21s
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
9.57 1.23 23.35