R2MLwiN: Running MLwiN from within R
R2MLwiN is an R command interface to the MLwiN multilevel modelling software package, allowing users to fit multilevel models using MLwiN from within the R environment.
R2MLwiN can also be used to run the script version of MLwiN on a non-Windows system; for more information see the MLwiN System Requirements webpage.
R2MLwiN is available from CRAN (see installation instructions below), although to use it you will also need to obtain and install MLwiN (see the MLwiN webpage for details of how to do so).
News
30th January 2024: The latest version of R2MLwiN (0.8-9) has been released on CRAN
The latest version of R2MLwiN includes the following updates and bug fixes.
Additional datasets
The following additional example datasets have been added:
- cvd
- fysio
- lmdp
These example datasets are also available in MLwiN, are are further described in Leyland, A.H., Groenewegen, P.P. (2020). Multilevel Modelling for Public Health and Health Services Research. Springer, Cham. doi:10.1007/978-3-030-34801-4_11
Examples of handling contrasts for ordered factors
We have also added examples of how to currently handle ordered factors in R2MLwiN. E.g. the example analyses in the help file for the fysio dataset (e.g. see ?R2MLwiN::fysio) illustrate the following:
1) Explicitly setting ordered factors to "treatment" contrast to avoid the warning that, by default, specified contrasts for ordered predictors will be ignored by runMLwiN. As follows:
R> my_contrasts <- options("contrasts")$contrasts
R> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.treatment"))
These can then be reverted to the pre-existing contrast specification:
R> options(contrasts = my_contrasts)
2) As an alternative to changing contrasts (as above), using the C() syntax to specify contrasts for ordered predictors in the formula object, e.g.:
R> (F1 <- logit(referral) ~ 1 + C(pagegrp, "contr.treatment") + patsex + diag + C(patedu, "contr.treatment") + patinsur + gpexper + gpworkload + practype + location + gpphysifr + (1 | gpid)
Bug fix
In addition, a bug fix at the end of the MCMC chapter 16 demo now ensures that the correct model is now fitted (switching the weights used to match those in the corresponding MLwiN manual).
14 March 2023: The latest version of R2MLwiN (0.8-8) has been released on CRAN
This includes updates to how the intercept is specified and handled, and also how contrasts are specified for categorical predictors (see below).
Specification of the intercept
In this latest version (0.8-8), if you do not specify a constant (intercept) in the fixed part of the model formula via "1 + ", then R2MLwiN will look for a constant of ones in the predictors specified in the model. If it finds one, then it will change it to "1" in the formula, and it will report this as the "Intercept" in the summarised output, etc.
We have done this to support better handling of contrasts.
How might this affect your code? In the past, if you included a constant of ones from your data when fitting a model (e.g. as "cons"), then referred to "cons" by name when extracting output from a model fit, or put "cons" in an unorthodox position in the formula (e.g. at the end of the fixed part) and used its index position when extracting output from a model fit, then the results will not be as anticipated following these recent changes (see examples below).
R> library("R2MLwiN")
## Change path to MLwiN executable as appropriate:
# options(MLwiN_path="<path to MLwiN>")
R> data(tutorial)
## The variable 'cons' is a constant of ones: R> table(tutorial$cons)
R> (model_1 <- runMLwiN(normexam ~ cons + standlrt + (1 | student), data = tutorial)) ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The model formula: # normexam ~ 1 + standlrt + (1 | student) # Level 1: student # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # Intercept -0.00119 0.01264 -0.09 0.9249 -0.02596 0.02358 # standlrt 0.59506 0.01273 46.76 0 *** 0.57011 0.62000 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## In the earlier version of R2MLwiN (e.g. 0.8-7) the output would have instead been: ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The model formula: # normexam ~ cons + standlrt + (1 | student) # Level 1: student # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # cons -0.00119 0.01264 -0.09 0.9249 -0.02596 0.02358 # standlrt 0.59506 0.01273 46.76 0 *** 0.57011 0.62000 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## In the earlier version of R2MLwiN (e.g. 0.8-7), if you had specified the fixed intercept as ## "cons", then this would have extracted its estimate: R> model_1@FP["FP_cons"] ## ...but now (v.0.8-8) instead need to specify: R> model_1@FP["FP_Intercept"]
## In the earlier version of R2MLwiN (e.g. 0.8-7) when specifying fixed intercept in an ## unorthodox position, e.g. as the second predictor, then this would have extracted its estimate: (model_2 model_2@FP[2] ## ...but now (v.0.8-8) instead need to specify: R> model_1@FP[1]
Contrasts for categorical predictors
The manner in which R2MLwiN specifies contrasts for categorical predictors is now more consistent with other packages in R.
## Specifying predictor as an ordered factor: R> tutorial$schav <- ordered(tutorial$schav)
## In the latest version (0.8-8), unless specify allowcontrast = TRUE, treatment contrasts ## will be used here - this is in keeping with behaviour of earlier versions of R2MLwiN (e.g. ## 0.8-7), but now we get a warning that this is happening: R> runMLwiN(normexam ~ 1 + schav + (1 | student), data = tutorial) # Warning in runMLwiN(normexam ~ 1 + schav + (1 | student), data = tutorial): # specified contrasts for variable schav will be ignored. To enable this set # allowcontrast to TRUE (this will be the default in future package releases) ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # Intercept -0.41218 0.03820 -10.79 3.887e-27 *** -0.48706 -0.33730 # schavmid 0.35048 0.04327 8.10 5.505e-16 *** 0.26568 0.43529 # schavhigh 0.76076 0.04762 15.98 1.89e-57 *** 0.66742 0.85409 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## However, we can now set allowcontrast = TRUE to specify polynomial contrasts for ## ordered predictors: R> runMLwiN(normexam ~ 1 + schav + (1 | student), allowcontrast = TRUE, data = tutorial) ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # Intercept -0.04177 0.01726 -2.42 0.01551 * -0.07559 -0.00794 # schav_L 0.53794 0.03367 15.98 1.89e-57 *** 0.47194 0.60393 # schav_Q 0.02441 0.02556 0.96 0.3396 -0.02568 0.07450 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## This is the same behaviour as e.g. lme4:
R> library("lme4")
R> lm(normexam ~ 1 + schav, data = tutorial)
# Call:
# lm(formula = normexam ~ 1 + schav, data = tutorial)
#
# Coefficients:
# (Intercept) schav.L schav.Q
# -0.04177 0.53794 0.02441
## An alternative method of specifying polynomial contrasts: R> tutorial$schav <- factor(tutorial$schav, ordered = FALSE) R> contrasts(tutorial$schav) <- contr.poly(3) ## R2MLwiN previously missed the first contrast column, so results were unlikely to be useful; ## now corrected: R> runMLwiN(normexam ~ 1 + schav + (1 | student), data = tutorial) ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # Intercept -0.04177 0.01726 -2.42 0.01551 * -0.07559 -0.00794 # schav_L 0.53794 0.03367 15.98 1.89e-57 *** 0.47194 0.60393 # schav_Q 0.02441 0.02556 0.96 0.3396 -0.02568 0.07450 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## This is the same behaviour as e.g. lme4: R> lm(normexam ~ 1 + schav, data = tutorial) # Call: # lm(formula = normexam ~ 1 + schav, data = tutorial) # # Coefficients: # (Intercept) schav.L schav.Q # -0.04177 0.53794 0.02441
## In the absence of an intercept, R2MLwiN previously (e.g. 0.8-7) did not add the first ## factor with treatment contrasts: R> contrasts(tutorial$schav, 3) <- NULL R> runMLwiN(normexam ~ 0 + schav + (1 | student), data = tutorial) ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # schavmid -0.06170 0.02061 -2.99 0.002753 ** -0.10208 -0.02131 # schavhigh 0.34857 0.02883 12.09 1.191e-33 *** 0.29207 0.40508 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## Now, in the absence of an intercept, R2MLwiN adds the first factor with treatment ## contrasts: ## <Output (truncated)> # -------------------------------------------------------------------------------------------- # The fixed part estimates: # Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] # schavlow -0.41218 0.03820 -10.79 3.887e-27 *** -0.48706 -0.33730 # schavmid -0.06170 0.02032 -3.04 0.002392 ** -0.10152 -0.02188 # schavhigh 0.34857 0.02843 12.26 1.444e-34 *** 0.29286 0.40429 # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # --------------------------------------------------------------------------------------------
## This is in keeping with e.g. lme4: R> lm(normexam ~ 0 + schav, data = tutorial) # Call: # lm(formula = normexam ~ 0 + schav, data = tutorial) # # Coefficients: # schavlow schavmid schavhigh # -0.4122 -0.0617 0.3486
28 May 2019: The latest version of R2MLwiN (0.8-6) has been released on CRAN
This includes the following new features (and also note bug fixes):
- The dependency on the rbugs package has been removed and R2MLwiN now uses the R2WinBUGS package instead for interoperability with BUGS.
06 Feb 2019: R2MLwiN is temporarily archived due to a removed dependency
The rbugs package, which R2MLwiN uses to run models via OpenBUGS/WinBUGS, was recently removed from CRAN. This has resulted in R2MLwiN also being removed as this was defined as a dependency. We are currently investigating how to make the package available again, however in the meantime it is possible to install it from the CRAN archive by running the following in R:
# Install the devtools package to provide install_url function
install.packages("devtools")
# Build and install rbugs
devtools::install_url("https://cran.r-project.org/src/contrib/Archive/rbugs/rbugs_0.5-9.1.tar.gz")
# Build and install R2MLwiN
devtools::install_url("https://cran.r-project.org/src/contrib/Archive/R2MLwiN/R2MLwiN_0.8-5.tar.gz")
Once the installation is complete you can load and use the package as usual.
01 Sep 2017: The latest version of R2MLwiN (0.8-5) has been released on CRAN.
This includes the following new features (and also note bug fixes):
- This is solely a bug-fix version and contains no new features
27 June 2017: The latest version of R2MLwiN (0.8-4) has been released on CRAN.
This includes the following new features (and also note bug fixes):
- Add option for setting the random number generator used by MLwiN
08 Sep 2016: The R2MLwiN paper is now available
If you wish to cite the article (available here), please do so as:
- 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.
07 Sep 2016: The latest version of R2MLwiN (0.8-3) has been released on CRAN.
This includes the following new features (and also note bug fixes):
- Correctly use the new data passed into "newdata" parameter of predict function
- Correctly check whether the model is cross-classified when calculating the model hierarchy
- Drop unused factors when calculating the complete-case group information
- Fix incorrect common coefficients with categorical variables
- Correct the parameter ordering to match that in the formula
- Add initial support for model comparison tables via the memisc package
- Model fit objects now contain contrast information for factor variables in the model
- Update model comparisons in demos to use the texreg package
- Add initial support for model comparison tables via the texreg package
- Allow variables that contain a "." character to work with formulae
(For older news items, see the news archive).
Installation
Both MLwiN and R are required to use R2MLwiN. MLwiN is free to UK academics. A fully functional 30-day free version of MLwiN is available to all other users.
As R2MLwiN depends on a number of third party packages we suggest that before installing the package you run the command:
update.packages(repos = "http://cran.r-project.org")
to ensure that any you already have installed are up to date.
You can then install the package from CRAN with the following command:
install.packages("R2MLwiN", repos = "http://cran.r-project.org")
Alternatively a less tested development version, which will contain the most recent changes and bug fixes, can be installed with the following commands:
install_github("rforge/
To then load it, type:
library(R2MLwiN)
Documentation and Examples
Journal article
An article describing R2MLwiN is available here. If you wish to cite the article, please do so as:
- 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.
Help files
As with all packages on CRAN, R2MLwiN comes with a number of help files; to see a list of the documentation available, type the following at the R command line:
help(package = "R2MLwiN")
The main command for calling MLwiN from R is called runMLwiN, and details of the syntax and options for this can be found with:
help(runMLwiN, package="R2MLwiN")
Demos
R2MLwiN also comes with a range of demos; these replicate all the analyses reported in the (IGLS) User's Guide to MLwiN and the MLwiN MCMC Manual. To see a list of these, type:
demo(package = "R2MLwiN")
To view the script of a specific demo (in this example the code replicating the analyses in Chapter 3 of the MLwiN MCMC Manual) type:
file.show(system.file("demo", "MCMCGuide03.R", package = "R2MLwiN"))
Finally, to run that demo, enter the following:
demo(MCMCGuide03)
These demos and their outputs can also be viewed on the Examples page.
Citations
Please cite R2MLwiN as follows:
We are happy to list all publications that cite R2MLwiN on this web site. More >>
Presentations
We gave a presentation on R2MLwiN at the 10th International Multilevel Conference in Utrecht (9th-10th April 2015). There's a pdf file below (converted from a Powerpoint file, with animations split courtesy of PPspliT):
Slides (PDF, 4.0mb)
User forum
Please use the R2MLwiN user forum to post any questions or feedback you have regarding R2MLwiN.
Acknowledgements
We are very grateful to colleagues at the Centre for Multilevel Modelling and the University of Bristol for their useful comments. This work is funded under the e-STAT node (Grant: RES-149-25-1084) as part of the ESRC Digital Social Research programme.