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="") > # 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