Nothing
############################################################################
# MLwiN User Manual
#
# 10 Multinomial Logistic Models for Unordered Categorical Responses . .145
#
# 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)
# 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)
# Change contrasts if wish to avoid warning indicating that, by default,
# specified contrasts for ordered predictors will be ignored by runMLwiN
# (they will be fitted as "contr.treatment" regardless of this setting). To
# enable specified contrasts, set allowcontrast to TRUE (this will be the
# default in future package releases). NB at the end of this script, the
# specification for contrasts is changed back.
my_contrasts <- options("contrasts")$contrasts
options(contrasts = c(unordered = "contr.treatment",
ordered = "contr.treatment"))
# As an alternative to changing contrasts, can instead use C() to specify
# contrasts for ordered predictors in formula object, e.g.:
# (mymodel1 <- runMLwiN(logit(use4) ~ 1 + C(lc, "contr.treatment"),
# D = "Unordered Multinomial",
# data = bang,
# allowcontrast = TRUE))
# 10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . .145
data(bang, package = "R2MLwiN")
addmargins(with(bang, table(use4)))
# 10.2 Single-level multinomial logistic regression . . . . . . . . . . .146
# 10.3 Fitting a single-level multinomial logistic model in MLwiN . . . .147
addmargins(with(bang, table(lc, use4)))
bang$use4 <- relevel(bang$use4, 4)
(mymodel1 <- runMLwiN(logit(use4) ~ 1 + lc, D = "Unordered Multinomial", data = bang))
cat(paste("Pr(y = 1) =", round(exp(mymodel1@FP["FP_Intercept_Sterilization"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) +
exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
cat(paste("Pr(y = 2) =", round(exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) +
exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
cat(paste("Pr(y = 3) =", round(exp(mymodel1@FP["FP_Intercept_Traditional_method"])/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) +
exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) + exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
cat(paste("Pr(y = 4) =", round(1/(1 + exp(mymodel1@FP["FP_Intercept_Sterilization"]) + exp(mymodel1@FP["FP_Intercept_Modern_reversible_method"]) +
exp(mymodel1@FP["FP_Intercept_Traditional_method"])), 4), "\n"))
# 10.4 A two-level random intercept multinomial logistic regression model 154
# 10.5 Fitting a two-level random intercept model . . . . . . . . . . . .155
(mymodel2 <- runMLwiN(logit(use4) ~ 1 + lc + (1 | district), D = "Unordered Multinomial", data = bang))
(mymodel3 <- runMLwiN(logit(use4) ~ 1 + lc + (1 | district), D = "Unordered Multinomial", estoptions = list(nonlinear = c(1,
2), startval = list(FP.b = mymodel2@FP, FP.v = mymodel2@FP.cov, RP.b = mymodel2@RP, RP.v = mymodel2@RP.cov), resi.store = TRUE),
data = bang))
mymodel3@RP["RP2_cov_Intercept_Sterilization_Intercept_Modern_reversible_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Sterilization"] * mymodel3@RP["RP2_var_Intercept_Modern_reversible_method"])
mymodel3@RP["RP2_cov_Intercept_Sterilization_Intercept_Traditional_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Sterilization"] * mymodel3@RP["RP2_var_Intercept_Traditional_method"])
mymodel3@RP["RP2_cov_Intercept_Modern_reversible_method_Intercept_Traditional_method"]/sqrt(mymodel3@RP["RP2_var_Intercept_Modern_reversible_method"] * mymodel3@RP["RP2_var_Intercept_Traditional_method"])
hipos <- rep(0, 2)
hipos[1] <- which(levels(as.factor(bang$district)) == 56)
hipos[2] <- which(levels(as.factor(bang$district)) == 11)
u0 <- mymodel3@residual$lev_2_resi_est_Intercept.Sterilization
u0se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Sterilization)
u0CI95 <- 1.96 * u0se
u0rank <- rank(u0)
u0rankhi <- u0 + u0CI95
u0ranklo <- u0 - u0CI95
u0rankno <- order(u0rank)
plot(1:60, u0[u0rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u0 residual estimate")
points(1:60, u0rankhi[u0rankno], pch = 24, bg = "grey")
points(1:60, u0ranklo[u0rankno], pch = 25, bg = "grey")
for (i in 1:60) lines(rep(i, 2), c(u0ranklo[u0rankno[i]], u0rankhi[u0rankno[i]]))
for (i in 1:2) points(x = which(u0rankno == hipos[i]), y = u0[u0rankno[which(u0rankno == hipos[i])]], pch = 22, bg = i +
1)
u1 <- mymodel3@residual$lev_2_resi_est_Intercept.Modern_reversible_method
u1se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Modern_reversible_method)
u1CI95 <- 1.96 * u1se
u1rank <- rank(u1)
u1rankhi <- u1 + u1CI95
u1ranklo <- u1 - u1CI95
u1rankno <- order(u1rank)
plot(1:60, u1[u1rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u1 residual estimate")
points(1:60, u1rankhi[u1rankno], pch = 24, bg = "grey")
points(1:60, u1ranklo[u1rankno], pch = 25, bg = "grey")
for (i in 1:60) lines(rep(i, 2), c(u1ranklo[u1rankno[i]], u1rankhi[u1rankno[i]]))
for (i in 1:2) points(x = which(u1rankno == hipos[i]), y = u1[u1rankno[which(u1rankno == hipos[i])]], pch = 22, bg = i +
1)
u2 <- mymodel3@residual$lev_2_resi_est_Intercept.Traditional_method
u2se <- sqrt(mymodel3@residual$lev_2_resi_var_Intercept.Traditional_method)
u2CI95 <- 1.96 * u2se
u2rank <- rank(u2)
u2rankhi <- u2 + u2CI95
u2ranklo <- u2 - u2CI95
u2rankno <- order(u2rank)
plot(1:60, u2[u2rankno], ylim = c(-2, 2), pch = 15, xlab = "Rank", ylab = "u2 residual estimate")
points(1:60, u2rankhi[u2rankno], pch = 24, bg = "grey")
points(1:60, u2ranklo[u2rankno], pch = 25, bg = "grey")
for (i in 1:60) lines(rep(i, 2), c(u2ranklo[u2rankno[i]], u2rankhi[u2rankno[i]]))
for (i in 1:2) points(x = which(u2rankno == hipos[i]), y = u2[u2rankno[which(u2rankno == hipos[i])]], pch = 22, bg = i +
1)
# Chapter learning outcomes . . . . . . . . . . . . . . . . . . . . . . .159
# Addendum: changing contrasts back to pre-existing . . . . . . . . . . . NA
# Following re-specification of contrast settings towards the start of this
# script, change contrasts back to pre-existing:
options(contrasts = my_contrasts)
############################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.