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 > # > # 12 Modelling Count Data . . . . . . . . . . . . . . . . . . . . . . . 181 > # > # 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) > > > # 12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .181 > > data(mmmec, package = "R2MLwiN") > summary(mmmec) nation region county obs Italy :95 Min. : 1.00 Min. : 1.00 Min. : 0.00 France :94 1st Qu.:27.00 1st Qu.: 90.25 1st Qu.: 8.00 UK :70 Median :44.00 Median :178.50 Median : 14.50 W_Germany:30 Mean :42.86 Mean :178.25 Mean : 27.83 Ireland :26 3rd Qu.:60.00 3rd Qu.:266.75 3rd Qu.: 31.00 Denmark :14 Max. :79.00 Max. :355.00 Max. :313.00 (Other) :25 exp cons uvbi Min. : 0.69 Min. :1 Min. :-8.900200 1st Qu.: 11.02 1st Qu.:1 1st Qu.:-4.158400 Median : 18.76 Median :1 Median :-0.886400 Mean : 27.80 Mean :1 Mean : 0.000204 3rd Qu.: 34.39 3rd Qu.:1 3rd Qu.: 3.275525 Max. :258.86 Max. :1 Max. :13.359000 > > # 12.2 Fitting a simple Poisson model . . . . . . . . . . . . . . . . . .182 > > (mymodel1 <- runMLwiN(log(obs) ~ 1 + uvbi + offset(log(exp)), D = "Poisson", data = mmmec)) /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 (Poisson) Estimation algorithm: IGLS MQL1 Elapsed time : 0.04s Number of obs: 354 (from total 354) The model converged after 4 iterations. Log likelihood: NA Deviance statistic: NA --------------------------------------------------------------------------------------------------- The model formula: log(obs) ~ 1 + uvbi + offset(log(exp)) Level 1: l1id --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.07011 0.01105 -6.35 2.195e-10 *** -0.09175 -0.04846 uvbi -0.05719 0.00268 -21.37 2.796e-101 *** -0.06244 -0.05194 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 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > # 12.3 A three-level analysis . . . . . . . . . . . . . . . . . . . . . .184 > > (mymodel2 <- runMLwiN(log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), D = "Poisson", estoptions = list(Meth = 0), + data = mmmec)) /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 (Poisson) N min mean max N_complete min_complete mean_complete max_complete nation 9 3 39.333333 95 9 3 39.333333 95 region 78 1 4.538462 13 78 1 4.538462 13 Estimation algorithm: RIGLS MQL1 Elapsed time : 0.05s Number of obs: 354 (from total 354) The model converged after 5 iterations. Log likelihood: NA Deviance statistic: NA --------------------------------------------------------------------------------------------------- The model formula: log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region) Level 3: nation Level 2: region Level 1: l1id --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept 0.11032 0.16045 0.69 0.4917 -0.20415 0.42480 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------------------------------------------------------------- The random part estimates at the nation level: Coef. Std. Err. var_Intercept 0.21450 0.10894 --------------------------------------------------------------------------------------------------- The random part estimates at the region level: Coef. Std. Err. var_Intercept 0.04537 0.00965 --------------------------------------------------------------------------------------------------- The random part estimates at the l1id level: Coef. Std. Err. var_bcons_1 1.00000 0.00000 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > (mymodel3 <- runMLwiN(log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region), D = "Poisson", estoptions = list(Meth = 0, + nonlinear = c(N = 1, M = 2)), data = mmmec)) /nogui option ignored ECHO 0 Echoing is ON BATC 1 Batch mode is ON MAXI 2 STAR iteration 0 iteration 1 Convergence not achieved TOLE 2 MAXI 20 NEXT iteration 2 iteration 3 iteration 4 iteration 5 iteration 6 Convergence achieved ECHO 0 Execution completed -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- MLwiN (version: 3.05) multilevel model (Poisson) N min mean max N_complete min_complete mean_complete max_complete nation 9 3 39.333333 95 9 3 39.333333 95 region 78 1 4.538462 13 78 1 4.538462 13 Estimation algorithm: RIGLS PQL2 Elapsed time : 0.1s Number of obs: 354 (from total 354) The model converged after 7 iterations. Log likelihood: NA Deviance statistic: NA --------------------------------------------------------------------------------------------------- The model formula: log(obs) ~ 1 + offset(log(exp)) + (1 | nation) + (1 | region) Level 3: nation Level 2: region Level 1: l1id --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.02451 0.15169 -0.16 0.8716 -0.32182 0.27280 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------------------------------------------------------------- The random part estimates at the nation level: Coef. Std. Err. var_Intercept 0.18545 0.09705 --------------------------------------------------------------------------------------------------- The random part estimates at the region level: Coef. Std. Err. var_Intercept 0.05763 0.01246 --------------------------------------------------------------------------------------------------- The random part estimates at the l1id level: Coef. Std. Err. var_bcons_1 1.00000 0.00000 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > (mymodel4 <- runMLwiN(log(obs) ~ 1 + uvbi + offset(log(exp)) + (1 | nation) + (1 | region), D = "Poisson", estoptions = list(Meth = 0, + nonlinear = c(N = 1, M = 2)), data = mmmec)) /nogui option ignored ECHO 0 Echoing is ON BATC 1 Batch mode is ON MAXI 2 STAR iteration 0 iteration 1 Convergence not achieved TOLE 2 MAXI 20 NEXT iteration 2 iteration 3 iteration 4 iteration 5 Convergence achieved ECHO 0 Execution completed -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- MLwiN (version: 3.05) multilevel model (Poisson) N min mean max N_complete min_complete mean_complete max_complete nation 9 3 39.333333 95 9 3 39.333333 95 region 78 1 4.538462 13 78 1 4.538462 13 Estimation algorithm: RIGLS PQL2 Elapsed time : 0.09s Number of obs: 354 (from total 354) The model converged after 6 iterations. Log likelihood: NA Deviance statistic: NA --------------------------------------------------------------------------------------------------- The model formula: log(obs) ~ 1 + uvbi + offset(log(exp)) + (1 | nation) + (1 | region) Level 3: nation Level 2: region Level 1: l1id --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] Intercept -0.08220 0.14226 -0.58 0.5634 -0.36101 0.19662 uvbi -0.02772 0.01132 -2.45 0.01435 * -0.04991 -0.00553 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------------------------------------------------------------- The random part estimates at the nation level: Coef. Std. Err. var_Intercept 0.15728 0.08281 --------------------------------------------------------------------------------------------------- The random part estimates at the region level: Coef. Std. Err. var_Intercept 0.04922 0.01095 --------------------------------------------------------------------------------------------------- The random part estimates at the l1id level: Coef. Std. Err. var_bcons_1 1.00000 0.00000 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > # 12.4 A two-level model using separate country terms . . . . . . . . . .186 > > addmargins(with(mmmec, table(nation))) nation Belgium W_Germany Denmark France UK Italy 11 30 14 94 70 95 Ireland Luxembourg Netherlands Sum 26 3 11 354 > > contrasts(mmmec$nation, 9) <- diag(9) > > (mymodel5 <- runMLwiN(log(obs) ~ 0 + nation + nation:uvbi + offset(log(exp)) + (1 | region), D = "Poisson", estoptions = list(Meth = 0, + nonlinear = c(N = 1, M = 2)), data = mmmec)) /nogui option ignored ECHO 0 Echoing is ON BATC 1 Batch mode is ON MAXI 2 STAR iteration 0 iteration 1 Convergence not achieved TOLE 2 MAXI 20 NEXT iteration 2 iteration 3 iteration 4 iteration 5 iteration 6 Convergence achieved ECHO 0 Execution completed -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- MLwiN (version: 3.05) multilevel model (Poisson) N min mean max N_complete min_complete mean_complete max_complete region 78 1 4.538462 13 78 1 4.538462 13 Estimation algorithm: RIGLS PQL2 Elapsed time : 0.12s Number of obs: 354 (from total 354) The model converged after 7 iterations. Log likelihood: NA Deviance statistic: NA --------------------------------------------------------------------------------------------------- The model formula: log(obs) ~ 0 + nation + nation:uvbi + offset(log(exp)) + (1 | region) Level 2: region Level 1: l1id --------------------------------------------------------------------------------------------------- The fixed part estimates: Coef. Std. Err. z Pr(>|z|) [95% Conf. Interval] nationBelgium 0.69790 0.74658 0.93 0.3499 -0.76538 2.16118 nationW_Germany 0.47761 0.12022 3.97 7.106e-05 *** 0.24198 0.71324 nationDenmark 0.31911 0.87941 0.36 0.7167 -1.40450 2.04273 nationFrance -0.59411 0.05433 -10.94 7.808e-28 *** -0.70059 -0.48762 nationUK 0.61400 0.20697 2.97 0.003011 ** 0.20835 1.01966 nationItaly 0.28151 0.10455 2.69 0.00709 ** 0.07660 0.48643 nationIreland -0.52867 1.30284 -0.41 0.6849 -3.08219 2.02485 nationLuxembourg 14.71263 15.53246 0.95 0.3435 -15.73042 45.15569 nationNetherlands -0.34752 0.92301 -0.38 0.7065 -2.15660 1.46155 nationBelgium:uvbi 0.26472 0.25191 1.05 0.2933 -0.22902 0.75845 nationW_Germany:uvbi -0.01336 0.03211 -0.42 0.6773 -0.07630 0.04957 nationDenmark:uvbi -0.08136 0.15521 -0.52 0.6001 -0.38557 0.22284 nationFrance:uvbi 0.01283 0.01797 0.71 0.4753 -0.02239 0.04804 nationUK:uvbi 0.14230 0.04267 3.33 0.0008534 *** 0.05867 0.22592 nationItaly:uvbi -0.08749 0.01579 -5.54 2.987e-08 *** -0.11842 -0.05655 nationIreland:uvbi -0.00096 0.26276 0.00 0.9971 -0.51595 0.51404 nationLuxembourg:uvbi 6.40748 6.77904 0.95 0.3446 -6.87918 19.69415 nationNetherlands:uvbi -0.11202 0.22159 -0.51 0.6132 -0.54633 0.32230 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 --------------------------------------------------------------------------------------------------- The random part estimates at the region level: Coef. Std. Err. var_Intercept 0.03572 0.00794 --------------------------------------------------------------------------------------------------- The random part estimates at the l1id level: Coef. Std. Err. var_bcons_1 1.00000 0.00000 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- > > xb <- predict(mymodel5) > > plot(mymodel5@data$uvbi, xb, xlab = "UV B radiation", ylab = "prediction", type = "n") > lines(mymodel5@data$uvbi[mymodel5@data$nationBelgium == 1], xb[mymodel5@data$nationBelgium == 1], col = 1) > lines(mymodel5@data$uvbi[mymodel5@data$nationW_Germany == 1], xb[mymodel5@data$nationW_Germany == 1], col = 2) > lines(mymodel5@data$uvbi[mymodel5@data$nationDenmark == 1], xb[mymodel5@data$nationDenmark == 1], col = 3) > lines(mymodel5@data$uvbi[mymodel5@data$nationFrance == 1], xb[mymodel5@data$nationFrance == 1], col = 4) > lines(mymodel5@data$uvbi[mymodel5@data$nationUK == 1], xb[mymodel5@data$nationUK == 1], col = 5) > lines(mymodel5@data$uvbi[mymodel5@data$nationItaly == 1], xb[mymodel5@data$nationItaly == 1], col = 6) > lines(mymodel5@data$uvbi[mymodel5@data$nationIreland == 1], xb[mymodel5@data$nationIreland == 1], col = 7) > lines(mymodel5@data$uvbi[mymodel5@data$nationLuxembourg == 1], xb[mymodel5@data$nationLuxembourg == 1], col = 8) > lines(mymodel5@data$uvbi[mymodel5@data$nationNetherlands == 1], xb[mymodel5@data$nationNetherlands == 1], col = 9) > legend(7, 0.7, c("belgium", "wgermany", "denmark", "france", "uk", "italy", "ireland", "luxembourg", "netherlands"), + lty = 1, col = 1:9) > > # 12.5 Some issues and problems for discrete response models . . . . . . 190 > > # Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .190 > > ############################################################################ > > proc.time() user system elapsed 3.65 0.53 5.54