R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ############################################################################
> # MLwiN User Manual
> #
> # 9 Logistic Models for Binary and Binomial Responses . . . . . . . . .117
> #
> # Rasbash, J., Steele, F., Browne, W. J. and Goldstein, H. (2012).
> # A User's Guide to MLwiN, v2.26. Centre for Multilevel Modelling,
> # University of Bristol.
> ############################################################################
> # R script to replicate all analyses using R2MLwiN
> #
> # Zhang, Z., Charlton, C., Parker, R, Leckie, G., and Browne, W.J.
> # Centre for Multilevel Modelling, 2012
> # http://www.bristol.ac.uk/cmm/software/R2MLwiN/
> ############################################################################
>
> library(R2MLwiN)
R2MLwiN: A package to run models implemented in MLwiN from R
Copyright 2013-2017 Zhengzheng Zhang, Christopher M. J. Charlton, Richard M. A. Parker, William J. Browne and George Leckie
Support provided by the Economic and Social Research Council (ESRC)
(Grants RES-149-25-1084, RES-576-25-0032 and ES/K007246/1)
To cite R2MLwiN in publications use:
Zhengzheng Zhang, Richard M. A. Parker, Christopher M. J. Charlton,
George Leckie, William J. Browne (2016). R2MLwiN: A Package to Run
MLwiN from within R. Journal of Statistical Software, 72(10), 1-43.
doi:10.18637/jss.v072.i10
A BibTeX entry for LaTeX users is
@Article{,
title = {{R2MLwiN}: A Package to Run {MLwiN} from within {R}},
author = {Zhengzheng Zhang and Richard M. A. Parker and Christopher M. J. Charlton and George Leckie and William J. Browne},
journal = {Journal of Statistical Software},
year = {2016},
volume = {72},
number = {10},
pages = {1--43},
doi = {10.18637/jss.v072.i10},
}
The MLwiN_path option is currently set to C:/Program Files/MLwiN v3.05/
To change this use: options(MLwiN_path="<path to MLwiN>")
> # MLwiN folder
> mlwin <- getOption("MLwiN_path")
> while (!file.access(mlwin, mode = 1) == 0) {
+ cat("Please specify the root MLwiN folder or the full path to the MLwiN executable:\n")
+ mlwin <- scan(what = character(0), sep = "\n")
+ mlwin <- gsub("\\", "/", mlwin, fixed = TRUE)
+ }
> options(MLwiN_path = mlwin)
>
>
> # 9.1 Introduction and description of the example data . . . . . . . . . 117
>
> data(bang, package = "R2MLwiN")
> summary(bang)
woman district use
Min. : 1.0 Min. : 1.00 Not_using:1728
1st Qu.: 717.5 1st Qu.:14.00 Using :1139
Median :1434.0 Median :29.00
Mean :1434.0 Mean :29.25
3rd Qu.:2150.5 3rd Qu.:45.00
Max. :2867.0 Max. :61.00
use4 lc age
Sterilization : 302 None : 774 Min. :-14.0000
Modern_reversible_method: 555 One_child : 517 1st Qu.: -8.0000
Traditional_method : 282 Two_children: 461 Median : -2.0000
Not_using_contraception :1728 Three_plus :1115 Mean : -0.3279
3rd Qu.: 6.0000
Max. : 19.0000
urban educ hindu d_lit
Rural:2063 None :1806 Muslim:2480 Min. :0.0000
Urban: 804 Lower_primary : 357 Hindu : 387 1st Qu.:0.0850
Upper_primary : 265 Median :0.1100
Secondary_and_above: 439 Mean :0.1115
3rd Qu.:0.1400
Max. :0.3000
d_pray cons
Min. :0.1000 Min. :1
1st Qu.:0.2900 1st Qu.:1
Median :0.4100 Median :1
Mean :0.4253 Mean :1
3rd Qu.:0.5500 3rd Qu.:1
Max. :0.7800 Max. :1
>
> # 9.2 Single-level logistic regression . . . . . . . . . . . . . . . . . 119
>
> # Link functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
>
> # Interpretation of coeficients . . . . . . . . . . . . . . . . . . . . .120
>
> # Fitting a single-level logit model in MLwiN . . . . . . . . . . . . . .121
>
> addmargins(with(bang, table(lc, use)))
use
lc Not_using Using Sum
None 584 190 774
One_child 283 234 517
Two_children 234 227 461
Three_plus 627 488 1115
Sum 1728 1139 2867
>
> (mymodel1 <- runMLwiN(logit(use) ~ 1 + lc, D = "Binomial", 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
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.05) multilevel model (Binomial)
Estimation algorithm: IGLS MQL1 Elapsed time : 0.05s
Number of obs: 2867 (from total 2867) The model converged after 4 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc
Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -1.12288 0.08348 -13.45 3.05e-41 *** -1.28650 -0.95926
lcOne_child 0.93275 0.12156 7.67 1.676e-14 *** 0.69450 1.17100
lcTwo_children 1.09250 0.12509 8.73 2.466e-18 *** 0.84733 1.33768
lcThree_plus 0.87224 0.10302 8.47 2.523e-17 *** 0.67033 1.07416
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> if (!require(car)) {
+ warning("car package required to use linearHypothesis() function")
+ } else {
+ linearHypothesis(mymodel1, "FP_lcOne_child = FP_lcTwo_children")
+ }
Loading required package: car
Loading required package: carData
Linear hypothesis test
Hypothesis:
FP_lcOne_child - FP_lcTwo_children = 0
Model 1: restricted model
Model 2: logit(use) ~ 1 + lc
Res.Df Df Chisq Pr(>Chisq)
1 2863
2 2862 1 1.5481 0.2134
>
> # A probit model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
>
> (mymodel2 <- runMLwiN(probit(use) ~ 1 + lc, D = "Binomial", 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
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.05) multilevel model (Binomial)
Estimation algorithm: IGLS MQL1 Elapsed time : 0.05s
Number of obs: 2867 (from total 2867) The model converged after 4 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
probit(use) ~ 1 + lc
Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -0.68879 0.04915 -14.01 1.299e-44 *** -0.78513 -0.59245
lcOne_child 0.56972 0.07396 7.70 1.328e-14 *** 0.42476 0.71468
lcTwo_children 0.66976 0.07631 8.78 1.69e-18 *** 0.52018 0.81933
lcThree_plus 0.53191 0.06195 8.59 8.977e-18 *** 0.41049 0.65332
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> (mymodel3 <- runMLwiN(logit(use) ~ 1 + lc + age, D = "Binomial", 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
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.05) multilevel model (Binomial)
Estimation algorithm: IGLS MQL1 Elapsed time : 0.06s
Number of obs: 2867 (from total 2867) The model converged after 4 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc + age
Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -1.25598 0.09776 -12.85 8.821e-38 *** -1.44758 -1.06438
lcOne_child 0.99131 0.12376 8.01 1.15e-15 *** 0.74874 1.23388
lcTwo_children 1.22356 0.13480 9.08 1.12e-19 *** 0.95935 1.48777
lcThree_plus 1.11655 0.13824 8.08 6.649e-16 *** 0.84560 1.38750
age -0.01629 0.00609 -2.67 0.007514 ** -0.02823 -0.00435
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> # 9.3 A two-level random intercept model . . . . . . . . . . . . . . . . 128
>
> # Model specification . . . . . . . . . . . . . . . . . . . . . . . . . .128
>
> # Estimation procedures . . . . . . . . . . . . . . . . . . . . . . . . .128
>
> # Fitting a two-level random intercept model in MLwiN . . . . . . . . . .129
>
> (mymodel4 <- runMLwiN(logit(use) ~ 1 + lc + age + (1 | district), D = "Binomial", 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
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
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 MQL1 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) ~ 1 + lc + age + (1 | district)
Level 2: district Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -1.36711 0.12338 -11.08 1.557e-28 *** -1.60893 -1.12529
lcOne_child 0.98998 0.12643 7.83 4.869e-15 *** 0.74218 1.23777
lcTwo_children 1.27523 0.13816 9.23 2.711e-20 *** 1.00443 1.54603
lcThree_plus 1.21568 0.14245 8.53 1.413e-17 *** 0.93648 1.49487
age -0.01878 0.00625 -3.00 0.002659 ** -0.03102 -0.00653
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.27409 0.07138
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> (mymodel5 <- runMLwiN(logit(use) ~ 1 + lc + age + (1 | district), D = "Binomial", estoptions = list(nonlinear = c(N = 1,
+ M = 2), startval = list(FP.b = mymodel4@FP, FP.v = mymodel4@FP.cov, RP.b = mymodel4@RP, RP.v = mymodel4@RP.cov)),
+ 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 -1.3671107908183 0.989975145198351 1.27523297062681 1.21567634144111 -0.0187756435513459 '_FP_b'
JOIN 0.0152221068761287 -0.00857801335032769 0.0159844374382819 -0.010151723167991 0.00868106697570791 0.0190895131628167 -0.0124154753668934 0.0096459226360998 0.012442006063399 0.0202921736477393 0.000326589121711455 -0.000146849092291606 -0.000325928858650572 -0.000596278835798142 3.9048334047441e-05 '_FP_v'
JOIN 0.274088563383711 1 '_RP_b'
JOIN 0.0050956166989184 0 0 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
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.1s
Number of obs: 2867 (from total 2867) The model converged after 5 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc + age + (1 | district)
Level 2: district Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -1.46602 0.12791 -11.46 2.059e-30 *** -1.71671 -1.21532
lcOne_child 1.06285 0.12882 8.25 1.575e-16 *** 0.81037 1.31533
lcTwo_children 1.37010 0.14167 9.67 4.014e-22 *** 1.09242 1.64778
lcThree_plus 1.30391 0.14595 8.93 4.104e-19 *** 1.01786 1.58997
age -0.02005 0.00640 -3.13 0.001747 ** -0.03260 -0.00749
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.30776 0.07899
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> if (!require(car)) {
+ warning("car package required to use linearHypothesis() function")
+ } else {
+ linearHypothesis(mymodel5, "RP2_var_Intercept = 0")
+ }
Linear hypothesis test
Hypothesis:
RP2_var_Intercept = 0
Model 1: restricted model
Model 2: logit(use) ~ 1 + lc + age + (1 | district)
Res.Df Df Chisq Pr(>Chisq)
1 2861
2 2860 1 15.181 9.769e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> # Variance partition coeficient . . . . . . . . . . . . . . . . . . . . .131
>
> set.seed(1)
>
> invlogit <- function(x) exp(x)/(1 + exp(x))
>
> u <- sqrt(coef(mymodel5)["RP2_var_Intercept"]) * qnorm(runif(5000))
>
> p1 <- invlogit(coef(mymodel5)["FP_Intercept"] + u)
>
> p2 <- invlogit(coef(mymodel5)["FP_Intercept"] + coef(mymodel5)["FP_lc3plus"] + coef(mymodel5)["FP_age"] * -9.7 + u)
>
> p3 <- invlogit(coef(mymodel5)["FP_Intercept"] + coef(mymodel5)["FP_age"] * 15.3 + u)
>
> v1 <- p1 * (1 - p1)
> lev2var1 <- sd(p1)^2
> lev1var1 <- mean(v1)
>
> v2 <- p2 * (1 - p2)
> lev2var2 <- sd(p2)^2
> lev1var2 <- mean(v2)
>
> v3 <- p3 * (1 - p3)
> lev2var3 <- sd(p3)^2
> lev1var3 <- mean(v3)
>
> cat(paste0("VPC = ", lev2var1/(lev2var1 + lev1var1)))
VPC = 0.0491754301807682>
> cat(paste0("VPC for a young women with 3+ children (low probability use) = ", lev2var2/(lev2var2 + lev1var2)))
VPC for a young women with 3+ children (low probability use) = NA>
> cat(paste0("VPC for an old woman with no children (high probability use) = ", lev2var3/(lev2var3 + lev1var3)))
VPC for an old woman with no children (high probability use) = 0.0419272279288916>
> # Adding further explanatory variables . . . . . . . . . . . . . . . . . 134
>
> table(bang$educ)
None Lower_primary Upper_primary Secondary_and_above
1806 357 265 439
>
> (mymodel6 <- runMLwiN(logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 | district), D = "Binomial", estoptions = list(nonlinear = c(N = 1,
+ M = 2), startval = list(FP.b = mymodel5@FP, FP.v = mymodel5@FP.cov, RP.b = mymodel5@RP, RP.v = mymodel5@RP.cov)),
+ 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 -1.46601726374528 1.06284833330899 1.37009902258045 1.30391318250182 -0.0200481959375968 0 0 0 0 0 '_FP_b'
JOIN 0.0163603695419892 -0.00890762958205187 0.0165943547340283 -0.0106419958708061 0.00897007438476784 0.0200716252826363 -0.0130375955794902 0.0100149271005426 0.0130197255583241 0.0213006854805154 0.000347319944589179 -0.000157647285641369 -0.000346447992241024 -0.00063025959287392 4.10233254113274e-05 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 '_FP_v'
JOIN 0.307757058623237 1 '_RP_b'
JOIN 0.00623907692791556 0 0 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
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.12s
Number of obs: 2867 (from total 2867) The model converged after 5 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc + age + 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.05249 0.13819 -14.85 6.704e-50 *** -2.32334 -1.78164
lcOne_child 1.15141 0.13413 8.58 9.12e-18 *** 0.88853 1.41429
lcTwo_children 1.51227 0.14734 10.26 1.026e-24 *** 1.22349 1.80105
lcThree_plus 1.50191 0.15271 9.83 7.972e-23 *** 1.20260 1.80122
age -0.01736 0.00665 -2.61 0.009033 ** -0.03039 -0.00433
urbanUrban 0.53306 0.10482 5.09 3.665e-07 *** 0.32762 0.73850
educLower_primary 0.24654 0.12836 1.92 0.05478 . -0.00505 0.49812
educUpper_primary 0.72433 0.14380 5.04 4.731e-07 *** 0.44248 1.00618
educSecondary_and_above 1.17020 0.12716 9.20 3.506e-20 *** 0.92096 1.41944
hinduHindu 0.43282 0.12765 3.39 0.0006973 *** 0.18263 0.68301
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.23364 0.06534
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> # 9.4 A two-level random coeficient model . . . . . . . . . . . . . . . .135
>
> (mymodel7 <- runMLwiN(logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 + urban | district), D = "Binomial",
+ estoptions = list(nonlinear = c(N = 1, M = 2), startval = list(FP.b = mymodel6@FP, FP.v = mymodel6@FP.cov, RP.b = mymodel6@RP,
+ RP.v = mymodel6@RP.cov)), 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 -2.0524878126349 1.1514091831992 1.51227328126239 1.50190931077976 -0.0173564094857298 0.533063022676533 0.246535751852934 0.724328418970923 1.17019898626514 0.432817014945708 '_FP_b'
JOIN 0.0190969057181895 -0.010249298461099 0.0179897690949509 -0.0123640566867248 0.00988274764015615 0.0217090961479245 -0.0151256474103302 0.0110367966927285 0.0143831177429424 0.0233214607718276 0.000362885872261872 -0.000160504805163484 -0.000363883885511234 -0.000670143247882038 4.41947424646325e-05 -0.00280269546137152 0.0005261144156233 0.000849328751980673 0.000702752807306554 -1.93563397557956e-05 0.0109869565822594 -0.00283420656434017 0.00036229338945785 0.000199158958934153 -0.000127364884113156 4.59999223868248e-05 -0.00059630716642084 0.0164773314339917 -0.00265716065127714 -0.000132304943914623 -1.94267258205956e-05 -0.00012744129452743 7.55732006752395e-05 -0.000986714124785436 0.00312301236770023 0.0206794937660382 -0.00433453673756542 0.00178163571915053 0.00218751328572626 0.00282048124886048 3.59899532831458e-05 -0.00290733241170694 0.00339389445448919 0.00374076332064882 0.0161709104496365 -0.00310427091118153 0.000291548581329297 0.000563339535587039 0.00123544577460758 -5.5691862556009e-05 6.03577341367192e-05 2.13428023572961e-05 -0.000294041097998737 0.000123287082941657 0.0162946578344111 '_FP_v'
JOIN 0.233641534263787 0 0 1 '_RP_b'
JOIN 0.00426981357741411 0 0 0 0 0 0 0 0 8.77180296005631e-36 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 7
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
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.21s
Number of obs: 2867 (from total 2867) The model converged after 8 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 + urban |
district)
Level 2: district Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -2.09365 0.14823 -14.12 2.704e-45 *** -2.38418 -1.80312
lcOne_child 1.16597 0.13489 8.64 5.44e-18 *** 0.90159 1.43036
lcTwo_children 1.52627 0.14841 10.28 8.31e-25 *** 1.23539 1.81715
lcThree_plus 1.52278 0.15408 9.88 4.924e-23 *** 1.22079 1.82477
age -0.01818 0.00670 -2.71 0.00667 ** -0.03131 -0.00504
urbanUrban 0.57420 0.13647 4.21 2.58e-05 *** 0.30674 0.84167
educLower_primary 0.24518 0.12951 1.89 0.05834 . -0.00866 0.49901
educUpper_primary 0.73327 0.14533 5.05 4.523e-07 *** 0.44843 1.01811
educSecondary_and_above 1.17969 0.12839 9.19 4.003e-20 *** 0.92804 1.43134
hinduHindu 0.50956 0.13292 3.83 0.0001263 *** 0.24904 0.77007
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.35974 0.09867
cov_Intercept_urbanUrban -0.25801 0.11151
var_urbanUrban 0.34896 0.17340
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> if (!require(car)) {
+ warning("car package required to use linearHypothesis() function")
+ } else {
+ linearHypothesis(mymodel7, "RP2_cov_Intercept_urbanUrban = 0")
+ }
Linear hypothesis test
Hypothesis:
RP2_cov_Intercept_urbanUrban = 0
Model 1: restricted model
Model 2: logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 + urban |
district)
Res.Df Df Chisq Pr(>Chisq)
1 2854
2 2853 1 5.3534 0.02068 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> if (!require(car)) {
+ warning("car package required to use linearHypothesis() function")
+ } else {
+ linearHypothesis(mymodel7, "RP2_var_urbanUrban = 0")
+ }
Linear hypothesis test
Hypothesis:
RP2_var_urbanUrban = 0
Model 1: restricted model
Model 2: logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 + urban |
district)
Res.Df Df Chisq Pr(>Chisq)
1 2854
2 2853 1 4.0499 0.04417 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> if (!require(car)) {
+ warning("car package required to use linearHypothesis() function")
+ } else {
+ linearHypothesis(mymodel7, c("RP2_cov_Intercept_urbanUrban = 0", "RP2_var_urbanUrban = 0"))
+ }
Linear hypothesis test
Hypothesis:
RP2_cov_Intercept_urbanUrban = 0
RP2_var_urbanUrban = 0
Model 1: restricted model
Model 2: logit(use) ~ 1 + lc + age + urban + educ + hindu + (1 + urban |
district)
Res.Df Df Chisq Pr(>Chisq)
1 2855
2 2853 2 5.4717 0.06484 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> (mymodel8 <- runMLwiN(logit(use) ~ 1 + lc + age + urban + educ + hindu + d_lit + d_pray + (1 + urban | district),
+ D = "Binomial", estoptions = list(nonlinear = c(N = 1, M = 2), startval = list(FP.b = mymodel7@FP, FP.v = mymodel7@FP.cov,
+ RP.b = mymodel7@RP, RP.v = mymodel7@RP.cov)), 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 -2.09365114401487 1.16597447430765 1.52626872549308 1.52277975007742 -0.0181756030954108 0.574203474031062 0.245175315697699 0.733269026333489 1.17968787196353 0.509555790433903 0 0 '_FP_b'
JOIN 0.0219733196802088 -0.0104541888031348 0.0181959598975733 -0.0126078929973787 0.0100424237520441 0.0220256026701775 -0.0154396033053014 0.0112258857845742 0.0146092310342419 0.0237399309100986 0.000370979834478551 -0.000165444373814286 -0.000370526959534443 -0.000682628446709666 4.48870021916425e-05 -0.00767744730675773 0.000480161961342514 0.000870968622887474 0.000659233695422207 -2.36492406938521e-05 0.018622783474728 -0.00291720245689676 0.000407641708114746 0.00022877717133059 -7.54876972434869e-05 4.6508417010466e-05 -0.000705193071727363 0.0167726708681131 -0.0028148671892566 -2.73405942759159e-05 7.53379457881814e-05 -4.51502866793914e-05 7.59301974829547e-05 -0.00108779055559781 0.00324217544128519 0.0211210684327045 -0.00459151677866056 0.00186306059818434 0.00229861328230719 0.00297792844854554 3.349332574253e-05 -0.00294660060949399 0.00351102245403737 0.00388494224260141 0.0164852284249341 -0.00349614929466233 0.000407001163384805 0.000674570662564585 0.00135662192237635 -5.71220594618246e-05 -2.59130801425283e-05 2.33979892126268e-05 -0.000157509096792477 0.000341272982594661 0.0176668751018998 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 '_FP_v'
JOIN 0.359738458340248 -0.258012324937657 0.348955050974729 1 '_RP_b'
JOIN 0.00973512093334196 -0.00773262176575147 0.0124351943508086 0.00611415126486965 -0.0150050721588829 0.0300670987206531 0 0 0 0 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
iteration 5
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
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.16s
Number of obs: 2867 (from total 2867) The model converged after 6 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use) ~ 1 + lc + age + 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% Conf. Interval]
Intercept -1.72295 0.26327 -6.54 5.976e-11 *** -2.23896 -1.20695
lcOne_child 1.17020 0.13503 8.67 4.462e-18 *** 0.90555 1.43486
lcTwo_children 1.53401 0.14861 10.32 5.6e-25 *** 1.24273 1.82529
lcThree_plus 1.52828 0.15425 9.91 3.857e-23 *** 1.22595 1.83061
age -0.01814 0.00670 -2.71 0.006761 ** -0.03127 -0.00501
urbanUrban 0.52822 0.13814 3.82 0.0001314 *** 0.25747 0.79897
educLower_primary 0.23770 0.12995 1.83 0.06738 . -0.01700 0.49239
educUpper_primary 0.74232 0.14559 5.10 3.421e-07 *** 0.45697 1.02768
educSecondary_and_above 1.19596 0.12895 9.27 1.779e-20 *** 0.94323 1.44869
hinduHindu 0.50955 0.13249 3.85 0.0001201 *** 0.24987 0.76922
d_lit 2.07495 1.70615 1.22 0.2239 -1.26905 5.41894
d_pray -1.40837 0.53393 -2.64 0.008346 ** -2.45486 -0.36188
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.30501 0.08818
cov_Intercept_urbanUrban -0.23343 0.10549
var_urbanUrban 0.35162 0.17408
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
>
> # 9.5 Modelling binomial data . . . . . . . . . . . . . . . . . . . . . .139
>
> # Modelling district-level variation with district-level proportions . . 139
>
> # Creating a district-level data set . . . . . . . . . . . . . . . . . . 140
>
> if (!require(doBy)) {
+ warning("package doBy required to run this example")
+ } else {
+ bangshort <- summaryBy(use + cons ~ district + d_lit + d_pray, FUN = c(mean, sum), data = bang)
+ bangshort$use.sum <- NULL
+ colnames(bangshort) <- c("district", "d_lit", "d_pray", "use", "cons", "denom")
+ bangshort$use <- bangshort$use - 1
+
+ # Fitting the model . . . . . . . . . . . . . . . . . . . . . . . . . . .142
+
+ (mymodel9 <- runMLwiN(logit(use, denom) ~ 1 + d_lit + d_pray + (1 | district), D = "Binomial", data = bangshort))
+ print(mymodel9)
+
+ (mymodel10 <- runMLwiN(logit(use, denom) ~ 1 + d_lit + d_pray + (1 | district), D = "Binomial", estoptions = list(nonlinear = c(N = 1,
+ M = 2), startval = list(FP.b = mymodel9@FP, FP.v = mymodel9@FP.cov, RP.b = mymodel9@RP, RP.v = mymodel9@RP.cov)),
+ data = bangshort))
+ }
Loading required package: doBy
/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
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.05) multilevel model (Binomial)
N min mean max N_complete min_complete mean_complete max_complete
district 60 1 1 1 60 1 1 1
Estimation algorithm: IGLS MQL1 Elapsed time : 0.04s
Number of obs: 60 (from total 60) The model converged after 4 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use, denom) ~ 1 + d_lit + d_pray + (1 | district)
Level 2: district Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -0.39797 0.23395 -1.70 0.08892 . -0.85650 0.06056
d_lit 3.75968 1.61790 2.32 0.02014 * 0.58865 6.93071
d_pray -1.19571 0.50106 -2.39 0.01702 * -2.17777 -0.21366
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.21061 0.05890
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
/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 -0.397970286407069 3.75968303276531 -1.19571286717458 '_FP_b'
JOIN 0.054731962557466 -0.164897319586594 2.61760668722093 -0.0718793179017356 -0.308613346250095 0.251058756151548 '_FP_v'
JOIN 0.210612498304647 1 '_RP_b'
JOIN 0.00346930817161946 0 0 '_RP_v'
TOLE 2
MAXI 20
NEXT
iteration 2
iteration 3
iteration 4
Convergence achieved
ECHO 0
Execution completed
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
MLwiN (version: 3.05) multilevel model (Binomial)
N min mean max N_complete min_complete mean_complete max_complete
district 60 1 1 1 60 1 1 1
Estimation algorithm: IGLS PQL2 Elapsed time : 0.06s
Number of obs: 60 (from total 60) The model converged after 5 iterations.
Log likelihood: NA
Deviance statistic: NA
---------------------------------------------------------------------------------------------------
The model formula:
logit(use, denom) ~ 1 + d_lit + d_pray + (1 | district)
Level 2: district Level 1: l1id
---------------------------------------------------------------------------------------------------
The fixed part estimates:
Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval]
Intercept -0.42665 0.24072 -1.77 0.07633 . -0.89845 0.04515
d_lit 3.99760 1.68806 2.37 0.01788 * 0.68906 7.30614
d_pray -1.25083 0.52205 -2.40 0.01658 * -2.27402 -0.22763
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.22523 0.06234
---------------------------------------------------------------------------------------------------
The random part estimates at the l1id level:
Coef. Std. Err.
var_bcons_1 1.00000 0.00000
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
> ############################################################################
>
> proc.time()
user system elapsed
4.03 0.81 7.85