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.

> ############################################################################
> # R2MLWIN - A PROGRAM TO RUN THE MLWIN MULTILEVEL MODELLING SOFTWARE FROM
> #            WITHIN R
> #
> #            Zhengzheng Zhang, Chris Charlton, Richard Parker, George Leckie
> #            and William Browne
> #            Centre for Multilevel Modelling
> #            University of Bristol
> #
> #            http://www.bristol.ac.uk/cmm/software/r2mlwin/
> #
> ############################################################################
> 
> ############################################################################
> # 1.  INTRODUCTION
> ############################################################################
> 
> ############################################################################
> # 1.3 The r2mlwin package
> ############################################################################
> 
> # Install the R2MLwiN package from the Comprehensive R Archive Network
> # (CRAN) archive
> #install.packages("R2MLwiN")
> 
> # Load the R2MLwiN package
> 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>")
> 
> # Specify the file address for mlwin.exe on your computer as the global
> # option MLwiN_path
> options(MLwiN_path="C:/Program Files/MLwiN v3.05/")
> 
> 
> 
> 
> 
> ############################################################################
> # 2.  CONTINUOUS RESPONSE MODELS
> ############################################################################
> 
> # Load the tutorial data into memory
> data(tutorial, package = "R2MLwiN")
> 
> 
> 
> ############################################################################
> # 2.1 Variance components model
> ############################################################################
> 
> # Fit a two-level students-within-schools variance components model to
> # students' age 16 scores
> (model1 <- runMLwiN(normexam ~ 1 + (1 | school) + (1 | student), estoptions=list(debugmode=TRUE), data=tutorial, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Normal) 
        N min     mean max N_complete min_complete mean_complete max_complete
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 3.3s 
Number of obs:  4059 (from total 4059)        The model converged after 3 iterations.
Log likelihood:      -5505.3 
Deviance statistic:  11010.6 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + (1 | school) + (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.01317     0.05363   -0.25      0.806           -0.11827     0.09194 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                  Coef.   Std. Err. 
var_Intercept   0.16863     0.03245 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.84776     0.01897 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> # Display the school-level VPC or ICC statistic
> b <- coef(model1)
> VPC <- b["RP2_var_Intercept"] / (b["RP2_var_Intercept"] + b["RP1_var_Intercept"])
> names(VPC) <- "VPC"
> VPC
      VPC 
0.1659064 
> 
> # Re-fit the previous model , but do not pause in MLwiN
> (model1 <- runMLwiN(normexam ~ 1 + (1 | school) + (1 | student), data=tutorial, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Normal) 
        N min     mean max N_complete min_complete mean_complete max_complete
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 0.03s 
Number of obs:  4059 (from total 4059)        The model converged after 3 iterations.
Log likelihood:      -5505.3 
Deviance statistic:  11010.6 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + (1 | school) + (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z   Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.01317     0.05363   -0.25      0.806           -0.11827     0.09194 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                  Coef.   Std. Err. 
var_Intercept   0.16863     0.03245 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.84776     0.01897 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> 
> 
> ############################################################################
> # 2.2 Random intercept model with predictor variables
> ############################################################################
> 
> # Generate a boys' school dummy variable
> tutorial$boysch = as.integer(tutorial$schgend=="boysch")
> 
> # Generate a girls' school dummy variable
> tutorial$girlsch = as.integer(tutorial$schgend=="girlsch")
> 
> tutorial$girl = as.integer(tutorial$sex=="girl")
> 
> # Fit a two-level random intercept model for student age 16 scores with
> # predictor variables and retrieve the Empirical Bayes estimates of the
> # school random effects
> # (Save the estimation results under the name model2)
> (model2 <- runMLwiN(normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 | school) + (1 | student), 
+   estoptions=list(resi.store=TRUE), data=tutorial, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Normal) 
        N min     mean max N_complete min_complete mean_complete max_complete
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 0.08s 
Number of obs:  4059 (from total 4059)        The model converged after 4 iterations.
Log likelihood:      -4662.7 
Deviance statistic:  9325.4 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 | school) + 
    (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.16815     0.05400   -3.11    0.001846   **      -0.27399    -0.06231 
standlrt     0.55996     0.01244   45.00           0   ***      0.53558     0.58435 
girl         0.16723     0.03408    4.91   9.264e-07   ***      0.10043     0.23403 
boysch       0.17762     0.11075    1.60      0.1088           -0.03945     0.39469 
girlsch      0.15896     0.08725    1.82     0.06848   .       -0.01205     0.32997 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                  Coef.   Std. Err. 
var_Intercept   0.08111     0.01618 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.56227     0.01258 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> # Test whether the boys' school effect is significantly different from the
> # girls' school effect
> if (!require(car)) install.packages("car"); library(car)
Loading required package: car
Loading required package: carData
> linearHypothesis(model2, "FP_boysch = FP_girlsch")
Linear hypothesis test

Hypothesis:
FP_boysch - FP_girlsch = 0

Model 1: restricted model
Model 2: normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 | school) + 
    (1 | student)

  Res.Df Df Chisq Pr(>Chisq)
1   4053                    
2   4052  1 0.023     0.8796
> 
> # Plot a quantile-quantile plot of the school residuals
> u0 <- model2@residual$lev_2_resi_est_Intercept
> 
> qqnorm(u0)
> 
> # Rank the school residuals
> u0se <- sqrt(model2@residual$lev_2_resi_var_Intercept)
> u0rank <- rank(u0)
> u0rankhi <- u0+u0se
> u0ranklo <- u0-u0se
> u0rankno <- order(u0rank)
> 
> # Plot a caterpillar plot of the school residuals 
> plot(1:65, u0[u0rankno], ylim=c(-1,1), pch=15, xlab="Rank",ylab="u0 residual estimate")
> points(1:65, u0rankhi[u0rankno], pch=24,bg="grey")
> points(1:65, u0ranklo[u0rankno], pch=25,bg="grey")
> for(i in 1:65) {lines(rep(i,2), c(u0ranklo[u0rankno[i]],u0rankhi[u0rankno[i]]))}
> 
> # Count how many schools' residuals are significantly below zero
> sum(u0 + 1.96*u0se < 0)
[1] 12
> 
> # Count how many schools' residuals are significantly above zero
> sum(u0 + 1.96*u0se > 0)
[1] 53
> 
> # Drop the variables u0 and u0se
> rm(u0,u0se)
> 
> 
> ############################################################################
> # 2.3 Random slope model
> ############################################################################
> 
> # Fit a two-level random slope model for student age 16 scores and retrieve
> # the empirical Bayes estimates of the school random effects
> # (Save the estimation results under the name model3)
> (model3 <- runMLwiN(normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 + standlrt | school) + (1 | student), 
+   estoptions=list(resi.store=TRUE), data=tutorial, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Normal) 
        N min     mean max N_complete min_complete mean_complete max_complete
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  IGLS        Elapsed time : 0.08s 
Number of obs:  4059 (from total 4059)        The model converged after 4 iterations.
Log likelihood:      -4640.6 
Deviance statistic:  9281.1 
--------------------------------------------------------------------------------------------------- 
The model formula:
normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 + standlrt | 
    school) + (1 | student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z     Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.18882     0.05135   -3.68     0.000236   ***     -0.28947    -0.08817 
standlrt     0.55443     0.01994   27.81   3.524e-170   ***      0.51535     0.59350 
girl         0.16826     0.03382    4.98    6.524e-07   ***      0.10197     0.23455 
boysch       0.17987     0.09914    1.81      0.06964   .       -0.01445     0.37418 
girlsch      0.17482     0.07876    2.22      0.02645   *        0.02045     0.32920 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                           Coef.   Std. Err. 
var_Intercept            0.07954     0.01597 
cov_Intercept_standlrt   0.02010     0.00649 
var_standlrt             0.01465     0.00441 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                  Coef.   Std. Err. 
var_Intercept   0.55020     0.01240 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> # Display the correlation between the random intercepts and the random
> # slopes
> b <- coef(model3)
> corr <- b["RP2_cov_Intercept_standlrt"] / sqrt(b["RP2_var_Intercept"] + b["RP2_var_standlrt"])
> names(corr) <- "corr"
> corr
      corr 
0.06549439 
> 
> # Performs a likelihood-ratio test to compare this model to the previous
> # random intercept model
> if (!require(lmtest)) install.packages("lmtest"); library(lmtest)
Loading required package: lmtest
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

> lrtest(model2, model3)
Likelihood ratio test

Model 1: normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 | school) + 
    (1 | student)
Model 2: normexam ~ 1 + standlrt + girl + boysch + girlsch + (1 + standlrt | 
    school) + (1 | student)
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   7 -4662.7                         
2   9 -4640.6  2 44.306  2.393e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> # Display a table of model results for model2 and model3
> if (!require(texreg)) install.packages("texreg"); library(texreg)
Loading required package: texreg
Version:  1.36.23
Date:     2017-03-03
Author:   Philip Leifeld (University of Glasgow)

Please cite the JSS article in your publications -- see citation("texreg").
> screenreg(c(model2, model3))

================================================
                            Model 1    Model 2  
------------------------------------------------
FP_Intercept                   -0.17      -0.19 
                               (0.05)     (0.05)
FP_standlrt                     0.56       0.55 
                               (0.01)     (0.02)
FP_girl                         0.17       0.17 
                               (0.03)     (0.03)
FP_boysch                       0.18       0.18 
                               (0.11)     (0.10)
FP_girlsch                      0.16       0.17 
                               (0.09)     (0.08)
RP2_var_Intercept               0.08       0.08 
                               (0.02)     (0.02)
RP1_var_Intercept               0.56       0.55 
                               (0.01)     (0.01)
RP2_cov_Intercept_standlrt                 0.02 
                                          (0.01)
RP2_var_standlrt                           0.01 
                                          (0.00)
------------------------------------------------
Num. obs.                    4059       4059    
Log Likelihood              -4662.71   -4640.56 
Deviance                     9325.43    9281.12 
================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
> 
> # Generate predicted age 16 scores based only on age 11 scores and school
> # attended
> u0 <- model3@residual$lev_2_resi_est_Intercept
> u1 <- model3@residual$lev_2_resi_est_standlrt
> 
> ypred = b["FP_Intercept"] + b["FP_standlrt"]*tutorial$standlrt + u0[tutorial$school] + u1[tutorial$school]*tutorial$standlrt
> 
> # Plot the predicted age 16 scores against age 11 scores
> plot(model3@data$standlrt, ypred, type="n")
> for (i in 1:65) {
+   lines(model3@data$standlrt[model3@data$school==i], ypred[model3@data$school==i], col="blue")
+ }
> 
> # Generate the predicted between-school variance for each student
> lev2var = b["RP2_var_Intercept"] + 2*b["RP2_cov_Intercept_standlrt"]*tutorial$standlrt + b["RP2_var_standlrt"]*tutorial$standlrt^2
> 
> # Plot the predicted between-school variance  against age 11 scores
> plot(sort(tutorial$standlrt), lev2var[order(tutorial$standlrt)], type="l")
> 
> 
> 
> 
> 
> 
> ############################################################################
> # 3.  BINARY RESPONSE MODELS
> ############################################################################
> 
> # Load the Bangladeshi contraceptive use data into memory
> data(bang, package = "R2MLwiN")
> 
> 
> 
> ############################################################################
> # 3.1 Variance components model
> ############################################################################
> 
> # Fit a two-level women-within-districts variance components logit model by
> # PQL2 for women's contraceptive use
> (model4 <- runMLwiN(logit(use, cons) ~ 1 + (1 | district), D="Binomial", estoptions=list(nonlinear=c(N=1,M=2)), data=bang, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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 : 0.08s 
Number of obs:  2867 (from total 2867)        The model converged after 5 iterations.
Log likelihood:      NA 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + (1 | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.       z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept   -0.51013     0.08084   -6.31   2.785e-10   ***     -0.66857    -0.35168 
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   0.26510     0.07035 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err. 
var_bcons_1   1.00000     0.00000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> # Fit the previous model by MCMC using the PQL2 estimates from the previous
> # model as starting values for MCMC estimation. Model takes 12 seconds to
> # fit on our computer.
> 
> (model5 <- runMLwiN(logit(use, cons) ~ 1 + (1 | district), D="Binomial",
+   estoptions=list(EstM=1, startval=list(FP.b=model4@FP, FP.v=model4@FP.cov, RP.b=model4@RP, RP.v=model4@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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:  MCMC      Elapsed time : 4.2s 
Number of obs:  2867 (from total 2867)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
3677.104   3635.560   41.544     3718.648   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + (1 | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
               Coef.   Std. Err.      z    Pr(>|z|)       [95% Cred.   Interval]   ESS 
Intercept   -0.52157     0.08811  -5.92   3.229e-09  ***    -0.68823    -0.34903   129 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                  Coef.   Std. Err.   [95% Cred.   Interval]   ESS 
var_Intercept   0.28626     0.08306      0.15856     0.48363   625 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_bcons_1   1.00000       1e-05      1.00000     1.00000   5000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> b <- coef(model5)
> 
> # Display the odds of using contraception in the average district
> exp(b["FP_Intercept"])
FP_Intercept 
   0.5935853 
> 
> # Display the probability of using contraception in the average district
> exp(b["FP_Intercept"])/(1 + exp(b["FP_Intercept"]))
FP_Intercept 
   0.3724842 
> 
> # Display the district-level VPC and ICC
> b["RP2_var_Intercept"]/(b["RP2_var_Intercept"] + pi^2/3)
RP2_var_Intercept 
       0.08004685 
> 
> # Plot a trajectories (trace) plot of the MCMC chains
> trajectories(model5)
> 
> # Plot a fiveway MCMC graphical diagnostic plot
> sixway(model5@chains[, "RP2_var_Intercept",drop=FALSE])
> 
> # Display detailed MCMC summary statistics for the between-district variance
> summary(model5@chains[, "RP2_var_Intercept",drop=FALSE])

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
      0.286258       0.083057       0.001175       0.003322 

2. Quantiles for each variable:

  2.5%    25%    50%    75%  97.5% 
0.1586 0.2266 0.2746 0.3316 0.4836 

> 
> 
> 
> ############################################################################
> # 3.2 Random intercept model with predictor variables
> ############################################################################
> 
> # Fit a two-level random intercept model for women's contraceptive use
> # adjusting for predictor variables by PQL2
> (model6 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 | district), 
+   D="Binomial", estoptions=list(nonlinear=c(N=1,M=2)), data=bang, stdout=NULL))


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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 : 0.11s 
Number of obs:  2867 (from total 2867)        The model converged after 5 iterations.
Log likelihood:      NA 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 | 
    district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                             Coef.   Std. Err.        z    Pr(>|z|)         [95% Conf.   Interval] 
Intercept                 -2.05270     0.13821   -14.85   6.712e-50   ***     -2.32358    -1.78182 
age                       -0.01736     0.00665    -2.61    0.009028   **      -0.03039    -0.00433 
lcOne_child                1.15153     0.13414     8.58   9.097e-18   ***      0.88863     1.41443 
lcTwo_children             1.51244     0.14735    10.26   1.022e-24   ***      1.22363     1.80124 
lcThree_plus               1.50207     0.15272     9.84   7.943e-23   ***      1.20274     1.80141 
urbanUrban                 0.53312     0.10482     5.09   3.661e-07   ***      0.32766     0.73857 
educLower_primary          0.24656     0.12837     1.92     0.05476   .       -0.00503     0.49816 
educUpper_primary          0.72440     0.14381     5.04   4.726e-07   ***      0.44253     1.00627 
educSecondary_and_above    1.17031     0.12717     9.20   3.496e-20   ***      0.92106     1.41957 
hinduHindu                 0.43287     0.12766     3.39   0.0006967   ***      0.18267     0.68307 
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   0.23369     0.06564 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err. 
var_bcons_1   1.00000     0.00000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> 
> # Fit the previous model by MCMC for a burnin of 1000 iterations followed
> # by a monitoring period of 10000 iterations. Model takes 72 seconds to fit
> # on our computer.
> 
> (model7 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 | district), 
+   D="Binomial", estoptions=list(EstM=1, resi.store=TRUE, mcmcMeth=list(burnin=1000,iterations=10000),
+   startval=list(FP.b=model6@FP, FP.v=model6@FP.cov, RP.b=model6@RP, RP.v=model6@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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:  MCMC      Elapsed time : 19.92s 
Number of obs:  2867 (from total 2867)          Number of iter.: 10000  Chains: 1  Burn-in: 1000 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
3392.300   3344.203   48.098     3440.398   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 | 
    district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                             Coef.   Std. Err.       z    Pr(>|z|)       [95% Cred.   Interval]    ESS 
Intercept                 -2.06420     0.13762  -15.00   7.412e-51  ***    -2.32476    -1.79734    144 
age                       -0.01742     0.00652   -2.67    0.007603  **     -0.03043    -0.00476    507 
lcOne_child                1.15234     0.13435    8.58   9.747e-18  ***     0.88695     1.41316    462 
lcTwo_children             1.51414     0.14830   10.21    1.79e-24  ***     1.22923     1.80654    302 
lcThree_plus               1.50584     0.14904   10.10   5.336e-24  ***     1.21538     1.79130    216 
urbanUrban                 0.53930     0.10555    5.11   3.233e-07  ***     0.32637     0.73812    991 
educLower_primary          0.24543     0.12875    1.91     0.05662  .      -0.00304     0.49867   1478 
educUpper_primary          0.72597     0.14552    4.99    6.08e-07  ***     0.44159     1.00731   1478 
educSecondary_and_above    1.17342     0.12691    9.25   2.321e-20  ***     0.91611     1.41426   1199 
hinduHindu                 0.44026     0.12616    3.49   0.0004838  ***     0.19458     0.69441   1110 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                  Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_Intercept   0.25544     0.07591      0.13696     0.43224   1168 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err.   [95% Cred.   Interval]     ESS 
var_bcons_1   1.00000       1e-05      1.00000     1.00000   10000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> 
> ############################################################################
> # 3.3 Random slope models
> ############################################################################
> 
> # Quietly fit a two-level random slope model for contraceptive use by PQL2
> model8 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 + urban | district), 
+   D="Binomial", estoptions=list(nonlinear=c(N=1,M=2)), data=bang, stdout=NULL)

> 
> # Fit the previous model by MCMC using an orthogonal fixed effects
> # parameterisation to make the MCMC chains less autocorrelated. Model takes 
> # 46 seconds to fit on our computer.
> (model9 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 + urban | district), 
+   D="Binomial", estoptions=list(EstM=1, resi.store=TRUE, mcmcOptions=list(orth=1),
+   startval=list(FP.b=model8@FP, FP.v=model8@FP.cov, RP.b=model8@RP, RP.v=model8@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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:  MCMC      Elapsed time : 12.63s 
Number of obs:  2867 (from total 2867)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
3359.126   3294.760   64.366     3423.491   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + (1 + 
    urban | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                             Coef.   Std. Err.       z    Pr(>|z|)       [95% Cred.   Interval]    ESS 
Intercept                 -2.10703     0.15813  -13.32   1.667e-40  ***    -2.42588    -1.80096    287 
age                       -0.01825     0.00673   -2.71    0.006687  **     -0.03161    -0.00535   1122 
lcOne_child                1.17630     0.13538    8.69   3.653e-18  ***     0.91907     1.44008   1024 
lcTwo_children             1.54079     0.15515    9.93   3.059e-23  ***     1.23926     1.83876    881 
lcThree_plus               1.53322     0.16082    9.53   1.521e-21  ***     1.22065     1.83412    916 
urbanUrban                 0.57730     0.13710    4.21   2.543e-05  ***     0.31724     0.85776    350 
educLower_primary          0.24682     0.13336    1.85      0.0642  .      -0.00742     0.50870   1084 
educUpper_primary          0.74061     0.14913    4.97   6.833e-07  ***     0.43837     1.03502    976 
educSecondary_and_above    1.19443     0.13462    8.87   7.147e-19  ***     0.93472     1.47661    942 
hinduHindu                 0.50521     0.13591    3.72   0.0002015  ***     0.22436     0.76555    640 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                              Coef.   Std. Err.   [95% Cred.   Interval]   ESS 
var_Intercept               0.40632     0.11732      0.22378     0.66973   303 
cov_Intercept_urbanUrban   -0.29106     0.11978     -0.57007    -0.10079   157 
var_urbanUrban              0.41135     0.19170      0.15794     0.86430    99 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_bcons_1   1.00000       1e-05      1.00000     1.00000   5000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> b <- coef(model9)
> 
> # Display the between-district variance for women living in rural areas
> b["RP2_var_Intercept"]
RP2_var_Intercept 
        0.4063235 
> 
> # Display the between-district variance for women living in urban areas
> b["RP2_var_Intercept"] + 2*b["RP2_cov_Intercept_urbanUrban"] + b["RP2_var_urbanUrban"]
RP2_var_Intercept 
        0.2355492 
> 
> # Quietly fit the previous model by PQL2 where we add in the two additional
> # predictor variables
> model10 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + d_lit + d_pray + (1 + urban | district), 
+   D="Binomial", estoptions=list(nonlinear=c(N=1,M=2)), data=bang, stdout=NULL)

> 
> # Fit the previous model by MCMC. Model takes 52 seconds to fit on our
> # computer.
> (model11 <- runMLwiN(logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + d_lit + d_pray+ (1 + urban | district), 
+   D="Binomial", estoptions=list(EstM=1, mcmcOptions=list(orth=1),
+   startval=list(FP.b=model10@FP, FP.v=model10@FP.cov, RP.b=model10@RP, RP.v=model10@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Binomial) 
          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:  MCMC      Elapsed time : 13.97s 
Number of obs:  2867 (from total 2867)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
3358.815   3293.307   65.509     3424.324   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use, cons) ~ 1 + age + lc + urban + educ + hindu + d_lit + 
    d_pray + (1 + urban | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                             Coef.   Std. Err.      z    Pr(>|z|)       [95% Cred.   Interval]    ESS 
Intercept                 -1.71880     0.27568  -6.23   4.526e-10  ***    -2.26465    -1.18918    193 
age                       -0.01833     0.00653  -2.81     0.00501  **     -0.03114    -0.00561   1090 
lcOne_child                1.17817     0.13832   8.52    1.63e-17  ***     0.89866     1.44439    987 
lcTwo_children             1.54331     0.14938  10.33   5.061e-25  ***     1.24558     1.83463   1036 
lcThree_plus               1.54268     0.15080  10.23   1.455e-24  ***     1.25103     1.83420   1043 
urbanUrban                 0.51981     0.14102   3.69   0.0002279  ***     0.24287     0.79706    272 
educLower_primary          0.23476     0.13541   1.73     0.08297  .      -0.02051     0.49963   1017 
educUpper_primary          0.74849     0.14890   5.03   4.991e-07  ***     0.45588     1.04285    959 
educSecondary_and_above    1.20271     0.12908   9.32   1.194e-20  ***     0.95763     1.46230   1029 
hinduHindu                 0.50843     0.13534   3.76   0.0001722  ***     0.23290     0.77017    619 
d_lit                      2.13893     1.87041   1.14      0.2528         -1.48382     6.05237    191 
d_pray                    -1.44559     0.57849  -2.50     0.01246  *      -2.62787    -0.35354    167 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                              Coef.   Std. Err.   [95% Cred.   Interval]   ESS 
var_Intercept               0.35184     0.11360      0.18817     0.62610   175 
cov_Intercept_urbanUrban   -0.26773     0.12037     -0.54694    -0.08586   154 
var_urbanUrban              0.43989     0.20560      0.15241     0.91940   115 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
                Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_bcons_1   1.00000       1e-05      1.00000     1.00000   5000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> b <- coef(model11)
> 
> # Display the between-district variance for women living in rural areas
> b["RP2_var_Intercept"]
RP2_var_Intercept 
        0.3518351 
> 
> # Display the between-district variance for women living in urban areas
> b["RP2_var_Intercept"] + 2*b["RP2_cov_Intercept_urbanUrban"] + b["RP2_var_urbanUrban"]
RP2_var_Intercept 
        0.2562748 
> 
> 
> ############################################################################
> # 4.  MORE ADVANCED MODELS
> ############################################################################
> 
> ############################################################################
> # 4.1 Multivariate response model
> ############################################################################
> 
> # Load the tutorial data into memory
> data(tutorial, package = "R2MLwiN")
> 
> tutorial$girl <- as.integer(tutorial$sex == "girl")
> 
> # Generate a binary response variable equal to 1 if normexam is positive,
> # otherwise 0
> tutorial$pass <- as.integer(tutorial$normexam > 0)
> 
> # Quietly fit a two-level students-within-schools bivariate mixed response
> # type model to continuous age 11 scores and binary age 16 scores
> model12 <- runMLwiN(c(standlrt,probit(pass, cons)) ~ 1 + girl + (1 | school) + (1[1] | student), 
+   D=c("Mixed", "Normal", "Binomial"), data=tutorial, stdout=NULL)

> 
> 
> 
> # Fit the previous model by MCMC. Model takes 239 seconds to fit on our
> # computer.
> (model13 <- runMLwiN(c(standlrt,probit(pass, cons)) ~ 1 + girl + (1 | school) + (1[1] | student), 
+   D=c("Mixed", "Normal", "Binomial"), estoptions=list(EstM=1, mcmcOptions=list(orth=1),
+   startval=list(FP.b=model12@FP, FP.v=model12@FP.cov, RP.b=model12@RP, RP.v=model12@RP.cov)), data=tutorial, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
MLwiN (version: 3.05)  multilevel model (Mixed) 
        N min     mean max N_complete min_complete mean_complete max_complete
school 65   2 62.44615 198         65            2      62.44615          198
Estimation algorithm:  MCMC      Elapsed time : 11.8s 
Number of obs:  4059 (from total 4059)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Deviance statistic:  NA 
--------------------------------------------------------------------------------------------------- 
The model formula:
c(standlrt, probit(pass, cons)) ~ 1 + girl + (1 | school) + (1[1] | 
    student)
Level 2: school     Level 1: student      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                   Coef.   Std. Err.      z    Pr(>|z|)       [95% Cred.   Interval]   ESS 
girl_standlrt    0.15522     0.04092   3.79   0.0001489  ***     0.07714     0.23364   255 
girl_pass        0.19430     0.05596   3.47   0.0005158  ***     0.08047     0.30189   137 
Intercept_1     -0.11833     0.04563  -2.59    0.009498  **     -0.20624    -0.03520    94 
Intercept_2     -0.09348     0.06196  -1.51      0.1313         -0.21288     0.03461    51 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the school level: 
                                Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_Intercept_1               0.09747     0.02081      0.06395     0.14431   1060 
cov_Intercept_1_Intercept_2   0.10289     0.02586      0.06025     0.16177    834 
var_Intercept_2               0.20116     0.04369      0.13093     0.30287    660 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the student level: 
                            Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
var_bcons_2               1.00000     0.00001      1.00000     1.00000   5000 
cov_bcons_2_Intercept_1   0.55860     0.01688      0.52595     0.59172    225 
var_Intercept_1           0.89926     0.01996      0.86012     0.93971    331 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> b <- coef(model13)
> 
> # Display the boy predicted probability of passing the age 16 exams
> b["FP_Intercept_2"]
FP_Intercept_2 
   -0.09348274 
> 
> # Display the girl predicted probability of passing the age 16 exams
> pnorm(b["FP_Intercept_2"] + b["FP_girl_pass"])
FP_Intercept_2 
      0.540153 
> 
> 
> ############################################################################
> # 4.2 Ordinal response model
> ############################################################################
> 
> # Load the Bangladeshi contraceptive use data into memory
> data(bang, package = "R2MLwiN")
> 
> # Set Not_using_contraception (4) as the base category
> # bang$use4 <- relevel(bang$use4, ref="Not_using_contraception")
> 
> # Fit a two-level women-within-districts ordinal logit model by MQL1 for
> # contraceptive use
> model14 <- runMLwiN(logit(use4, cons, 4) ~ 1 + age[1:3] + lc[1:3] + urban[1:3] + (1[1:3] | district),
+  D="Ordered Multinomial", data=bang, stdout=NULL)

> 
> # Fit the previous model by MCMC. Model takes 62 seconds to fit on our
> # computer.
> (model15 <- runMLwiN(logit(use4, cons, 4) ~ 1 + age[1:3] + lc[1:3] + urban[1:3] + (1[1:3] | district),
+  D="Ordered Multinomial",  estoptions=list(EstM=1, mcmcOptions=list(orth=1),
+  startval=list(FP.b=model14@FP, FP.v=model14@FP.cov, RP.b=model14@RP, RP.v=model14@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
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:  MCMC      Elapsed time : 12.94s 
Number of obs:  2867 (from total 2867)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
5876.279   5828.449   47.829     5924.108   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use4, cons, 4) ~ 1 + age[1:3] + lc[1:3] + urban[1:3] + 
    (1[1:3] | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                                        Coef.   Std. Err.       z     Pr(>|z|)       [95% Cred.   Interval]    ESS 
Intercept_Sterilization              -3.52680     0.15406  -22.89   5.534e-116  ***    -3.83268    -3.21873     64 
Intercept_Modern_reversible_method   -2.13273     0.14415  -14.79    1.582e-49  ***    -2.42068    -1.84519     79 
Intercept_Traditional_method         -1.64025     0.14208  -11.54    7.855e-31  ***    -1.92434    -1.35883     77 
age_123                              -0.01531     0.00592   -2.59     0.009693  **     -0.02701    -0.00360   1171 
lcOne_child_123                       1.04782     0.12795    8.19    2.628e-16  ***     0.80614     1.30904    883 
lcTwo_children_123                    1.38156     0.13806   10.01     1.42e-23  ***     1.10834     1.66263    862 
lcThree_plus_123                      1.29030     0.14388    8.97    3.029e-19  ***     0.99921     1.57235    926 
urbanUrban_123                        0.67353     0.09358    7.20    6.127e-13  ***     0.48241     0.85385    661 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                      Coef.   Std. Err.   [95% Cred.   Interval]   ESS 
var_Intercept_123   0.25254     0.07705      0.13237     0.42859   485 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
             Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
bcons_1    1.00000       1e-05      1.00000     1.00000   5000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> ############################################################################
> # 4.3 Nominal response model
> ############################################################################
> 
> # Load the Bangladeshi contraceptive use data into memory
> data(bang, package = "R2MLwiN")
> 
> # Set Not_using_contraception (4) as the base category
> bang$use4 <- relevel(bang$use4, ref="Not_using_contraception")
> 
> # Fit a two-level women-within-districts nominal logit model by MQL1 for
> # contraceptive use
> model16 <- runMLwiN(logit(use4, cons) ~ 1 + age + lc + urban + (1 | district),
+  D="Unordered Multinomial", data=bang, stdout=NULL)

> 
> # Fit the previous model by MCMC. Model takes 171 seconds to fit on our
> # computer.
> (model17 <- runMLwiN(logit(use4, cons) ~ 1 + age + lc + urban + (1 | district),
+  D="Unordered Multinomial",  estoptions=list(EstM=1, resi.store=TRUE, mcmcOptions=list(orth=1),
+   startval=list(FP.b=model16@FP, FP.v=model16@FP.cov, RP.b=model16@RP, RP.v=model16@RP.cov)), data=bang, stdout=NULL))
MLwiN is running, please wait......


-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
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:  MCMC      Elapsed time : 31.02s 
Number of obs:  2867 (from total 2867)          Number of iter.: 5000  Chains: 1  Burn-in: 500 
Bayesian Deviance Information Criterion (DIC)
Dbar      D(thetabar)    pD      DIC
5582.304   5477.948   104.356    5686.660   
--------------------------------------------------------------------------------------------------- 
The model formula:
logit(use4, cons) ~ 1 + age + lc + urban + (1 | district)
Level 2: district     Level 1: l1id      
--------------------------------------------------------------------------------------------------- 
The fixed part estimates:  
                                             Coef.   Std. Err.       z    Pr(>|z|)       [95% Cred.   Interval]    ESS 
Intercept_Sterilization                   -4.11203     0.32605  -12.61   1.823e-36  ***    -4.76556    -3.50975    145 
Intercept_Modern_reversible_method        -2.69326     0.16800  -16.03   7.734e-58  ***    -3.01830    -2.36269    323 
Intercept_Traditional_method              -2.78137     0.19659  -14.15   1.922e-45  ***    -3.18495    -2.39819    366 
age_Sterilization                          0.02150     0.00984    2.18     0.02889  *       0.00241     0.04073   1167 
age_Modern_reversible_method              -0.06647     0.00957   -6.94   3.814e-12  ***    -0.08546    -0.04850    706 
age_Traditional_method                     0.00068     0.01010    0.07      0.9464         -0.01897     0.02029   1079 
lcOne_child_Sterilization                  2.20115     0.32969    6.68    2.45e-11  ***     1.57190     2.86066    274 
lcOne_child_Modern_reversible_method       1.06660     0.15943    6.69   2.234e-11  ***     0.75744     1.36366    992 
lcOne_child_Traditional_method             0.80228     0.21932    3.66   0.0002541  ***     0.36028     1.23356    894 
lcTwo_children_Sterilization               2.66174     0.33362    7.98   1.483e-15  ***     2.01968     3.30759    276 
lcTwo_children_Modern_reversible_method    1.37241     0.18085    7.59   3.236e-14  ***     1.01635     1.72837    791 
lcTwo_children_Traditional_method          1.14548     0.24264    4.72   2.348e-06  ***     0.66812     1.61080    818 
lcThree_plus_Sterilization                 2.45488     0.33725    7.28   3.361e-13  ***     1.81292     3.11623    280 
lcThree_plus_Modern_reversible_method      1.32771     0.19111    6.95   3.723e-12  ***     0.96608     1.69694    812 
lcThree_plus_Traditional_method            1.19213     0.23674    5.04   4.761e-07  ***     0.72666     1.64815    831 
urbanUrban_Sterilization                   0.27182     0.16032    1.70     0.08998  .      -0.05797     0.57846    613 
urbanUrban_Modern_reversible_method        1.24600     0.12499    9.97    2.09e-23  ***     1.01741     1.50662    552 
urbanUrban_Traditional_method              0.24574     0.16755    1.47      0.1425         -0.08634     0.55863    640 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
--------------------------------------------------------------------------------------------------- 
The random part estimates at the district level: 
                                                                        Coef.   Std. Err.   [95% Cred.   Interval]   ESS 
var_Intercept_Sterilization                                           0.51829     0.14971      0.27882     0.87025   304 
cov_Intercept_Sterilization_Intercept_Modern_reversible_method        0.27080     0.09127      0.12251     0.48092   331 
var_Intercept_Modern_reversible_method                                0.28688     0.09273      0.14691     0.50744   319 
cov_Intercept_Sterilization_Intercept_Traditional_method              0.22880     0.10495      0.05568     0.46105   297 
cov_Intercept_Modern_reversible_method_Intercept_Traditional_method   0.10516     0.07738     -0.03479     0.27347   342 
var_Intercept_Traditional_method                                      0.36751     0.12498      0.18292     0.67055   293 
--------------------------------------------------------------------------------------------------- 
The random part estimates at the l1id level: 
             Coef.   Std. Err.   [95% Cred.   Interval]    ESS 
bcons_1    1.00000       1e-05      1.00000     1.00000   5000 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- 
> 
> b <- coef(model17)
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Sterilization as a function of age and urban, but
> # holding the number of living children constant at 2 and holding district
> # unobservables constant at zero
> n1 <- exp(b["FP_Intercept_Sterilization"]
+  + b["FP_age_Sterilization"]*bang$age
+  + b["FP_lcTwo_children_Sterilization"]
+  + b["FP_urbanUrban_Sterilization"]*as.integer(bang$urban=="Urban"))
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Modern_reversible_method as a function of age and urban, but
> # holding the number of living children constant at 2 and holding district
> # unobservables constant at zero
> n2 <- exp(b["FP_Intercept_Modern_reversible_method"]
+  + b["FP_age_Modern_reversible_method"]*bang$age
+  + b["FP_lcTwo_children_Modern_reversible_method"]
+  + b["FP_urbanUrban_Modern_reversible_method"]*as.integer(bang$urban=="Urban"))
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Traditional_method as a function of age and urban, but
> # holding the number of living children constant at 2 and holding district
> # unobservables constant at zero
> n3 <- exp(b["FP_Intercept_Traditional_method"]
+  + b["FP_age_Traditional_method"]*bang$age
+  + b["FP_lcTwo_children_Traditional_method"]
+  + b["FP_urbanUrban_Traditional_method"]*as.integer(bang$urban=="Urban"))
> 
> # Generate the denominator of the predicted probability expression
> d <- (1 + n1 + n2 + n3)
> 
> # Generate the predicted probability that use4=1
> p1 <- n1 / d
> 
> # Generate the predicted probability that use4=2
> p2 <- n3 / d
> 
> # Generate the predicted probability that use4=3
> p3 <- n2 / d
> 
> # Generate the predicted probability that use4=4
> p4 <- 1 - p1 - p2 - p3
> 
> # Generate a new variable actualage which stores the actual age of each
> # woman
> actualage <- bang$age + 30
> 
> plotdata <- data.frame(p1=p1,p2=p2,p3=p3,p4=p4,age=actualage,urban=bang$urban)
> plotdata <- plotdata[with(plotdata, order(age)), ]
> 
> 
> if (!require(lattice)) install.packages("lattice"); library(lattice)
Loading required package: lattice
> 
> xyplot(p4+p3+p2+p1~age|urban, data=plotdata, ylab="Probability", xlab="Woman's age (years)", type="l", ylim=c(0,1), 
+  auto.key=list(points=FALSE, lines=TRUE, text=c("No contraception", "Traditional method", "Modern method", "Modern method")))
> 
> # Drop the variables plotdata, n1, n2, n3, d, p1, p2, p3 and p4
> rm(plotdata,n1,n2,n3,d,p1,p2,p3,p4)
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Sterilization as a function of number of living
> # children, urban and district, but holding age constant at 30
> n1 <- exp(b["FP_Intercept_Sterilization"]
+  + b["FP_lcOne_child_Sterilization"]*as.integer(bang$lc == "One_child")
+  + b["FP_lcTwo_children_Sterilization"]*as.integer(bang$lc == "Two_children")
+  + b["FP_lcThree_plus_Sterilization"]*as.integer(bang$lc == "Three_plus")
+  + b["FP_urbanUrban_Sterilization"]*as.integer(bang$urban=="Urban")
+  + model17@residual$lev_2_resi_est_Intercept.Sterilization[bang$district])
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Modern_reversible_method as a function of number of living
> # children, urban and district, but holding age constant at 30
> n2 <- exp(b["FP_Intercept_Modern_reversible_method"]
+  + b["FP_lcOne_child_Modern_reversible_method"]*as.integer(bang$lc == "One_child")
+  + b["FP_lcTwo_children_Modern_reversible_method"]*as.integer(bang$lc == "Two_children")
+  + b["FP_lcThree_plus_Modern_reversible_method"]*as.integer(bang$lc == "Three_plus")
+  + b["FP_urbanUrban_Modern_reversible_method"]*as.integer(bang$urban=="Urban")
+  + model17@residual$lev_2_resi_est_Intercept.Modern_reversible_method[bang$district])
> 
> # Generate the numerator of the predicted probability expression for the
> # predicted probability that use4=Traditional_method as a function of number of living
> # children, urban and district, but holding age constant at 30
> n3 <- exp(b["FP_Intercept_Traditional_method"]
+  + b["FP_lcOne_child_Traditional_method"]*as.integer(bang$lc == "One_child")
+  + b["FP_lcTwo_children_Traditional_method"]*as.integer(bang$lc == "Two_children")
+  + b["FP_lcThree_plus_Traditional_method"]*as.integer(bang$lc == "Three_plus")
+  + b["FP_urbanUrban_Traditional_method"]*as.integer(bang$urban=="Urban")
+  + model17@residual$lev_2_resi_est_Intercept.Traditional_method[bang$district])
> 
> # Generate the denominator of the predicted probability expression
> d <- (1 + n1 + n2 + n3)
> 
> # Generate the predicted probability that use4=1
> p1 <- n1 / d
> 
> # Generate the predicted probability that use4=2
> p2 <- n3 / d
> 
> # Generate the predicted probability that use4=3
> p3 <- n2 / d
> 
> # Generate the predicted probability that use4=4
> p4 <- 1 - p1 - p2 - p3
> 
> # Generate a new variable lc_p4 equal to lc minus 0.15
> lc_p4 <- as.integer(bang$lc) - 0.15
> 
> # Generate a new variable lc_p4 equal to lc minus 0.075
> lc_p3 <- as.integer(bang$lc) - 0.075
> 
> # Generate a new variable lc_p4 equal to lc plus 0.075
> lc_p2 <- as.integer(bang$lc) + 0.075
> 
> # Generate a new variable lc_p4 equal to lc plus 0.15
> lc_p1 <- as.integer(bang$lc) + 0.15
> 
> plotdata <- data.frame(p1=p1,p2=p2,p3=p3,p4=p4,lc_p1=lc_p1,lc_p2=lc_p2,lc_p3=lc_p3,lc_p4=lc_p4,urban=bang$urban)
> 
> xyplot(c(p4,p3,p2,p1)~lc_p4+lc_p3+lc_p2+lc_p1|urban, outer=FALSE, data=plotdata, ylab="Probability", xlab="Number of living children", type="p", ylim=c(0,1), 
+  scales=list(x=list(labels=c(NA,"0","1","2","3+"))), auto.key=list(points=TRUE, lines=FALSE, text=c("No contraception", "Traditional method", "Modern method", "Modern method")))
> 
> # Drop the variables plotdata, n1, n2, n3, lc_p1, lc_p2, lc_p3 and lc_p4
> rm(plotdata,n1,n2,n3,lc_p1,lc_p2,lc_p3,lc_p4)
> 
> ############################################################################
> 
> proc.time()
   user  system elapsed 
   8.96    1.46  127.12