{"id":409,"date":"2012-11-23T22:05:13","date_gmt":"2012-11-23T22:05:13","guid":{"rendered":"http:\/\/www.georg-hosoya.de\/wordpress\/?p=409"},"modified":"2012-12-17T15:11:20","modified_gmt":"2012-12-17T15:11:20","slug":"multinomial-logit-regression-as-pgm","status":"publish","type":"post","link":"https:\/\/www.georg-hosoya.de\/wordpress\/?p=409","title":{"rendered":"Multinomial Logit Regression as PGM"},"content":{"rendered":"

I decided to switch to English as I believe that this topic may be interesting for a wider audience. In the last post I have given a simple example of how the t-test for independent samples can be understood from a generalized linear modeling (GLM) perspective, how the model translates into a probabilistic graphcal model (PGM) and how to use Bayesian estimation techniques for parameter estimation. The main motivation behind that post was to show, how generalized linear models, probabilistic graphical models and bayesian estimation techniques can work together for solving statistical inference tasks, particularly in but not limited to Psychology. <\/p>\n

I can hear the one or the other say: “Pah, t-test, peanuts!”. So in this post I would like to turn to something a bit more advanced: Multinomial Logit Regression (MLR). One task of MLR is to model the probability of a categorical response \"c\" out of a set of \"C\" mutually exclusive possible responses based on a set of predictors that may characterize for example the repondents. Typical questions are: Which characteristics of a person predict the choice of a specific party (Political Science) or a specific product (Economics) or a specific category on a multiple choice questionaire (Psychology). This type of model is also quite prominent in the machine learning community. One task is to classify text documents into \"C\" mutually exclusive categories based on a set of text features as predictors. <\/p>\n

Before we can specify the model as a PGM, we should have a look at its formal structure. First we need a probability distribution that is capable of modeling the probabilty of a certain categorical response \"c\" based on a set of parametetrs. One such distribution is:<\/p>\n

(1) <\/span>   <\/span>\"<\/p>\n

This is quite neat. The denominator of that distribution is sometimes called the partition function or in German Zustandssumme<\/em> from its use in statistical physics. It ensures that the probabilites for the possible events, states or responses sum up to one. Please note that the distribiution above has a little identification problem that can be solved by either setting \"\lambda_{1i}\" to zero or by imposing a sum constraint over the parameters.<\/p>\n

Here we set \"\lambda_{0i}\" to zero:<\/p>\n

(2) <\/span>   <\/span>\"<\/p>\n

This also means that we are using category 1 as a reference category for inference. More on that later.<\/p>\n

The trick in MLR is to simply decompose the parameters of the base distribution by simple linear regression equations:<\/p>\n

(3) <\/span>   <\/span>\"<\/p>\n

That’s basically it. The \"\beta\"-weights have to be estimated from data, either by maximum likelihood or MCMC. To clarify the interpretation of the parameters and to elaborate a bit more, let’s assume we have three possible categories \"C=3\" and only one predictor \"x_{1i}\". <\/p>\n

The probabilities of choosing a certain category \"c\" are:<\/p>\n

(4) <\/span>   <\/span>\"<\/p>\n

with partiton function<\/p>\n

\" <\/p>\n

Please note that \"Z\" runs over all three possible states of the system.<\/p>\n

If we write down the log-probability-ratio (or log-risk-ratio) of choosing category 2 over category 1 we get:<\/p>\n

(5) <\/span>   <\/span>\"<\/p>\n

From this equation the meanings of the parameters become clear. \"\beta_{20}\" is nothing but the log-probability-ratio of choosing category 2 over category 1, if \"x_{1i}\" is zero. \"\beta_{21}\" is the expected change in terms of log-probabilty-ratios per unit increase of \"x_{1i}\". <\/p>\n

The log-probability-ratio of choosing category 3 over category 1 is:<\/p>\n

(6) <\/span>   <\/span>\"<\/p>\n

Here we are predicting log-probability-ratios of choosing category 3 over category 1 based on the predictor \"x_{1i}\". \"\beta_{30}\" is the log-probability-ratio of choosing category 3 over category 1, if \"x_{1i}\" is zero and \"\beta_{31}\" is the expected change in terms of log-probability-ratios per unit increase of \"x_{1i}\". Quite easy and very similar to multiple regression. The only thing one has to keep in mind is that we are predicting \"C-1\" log-probability-ratios instead of changes in one continuous variable.<\/p>\n

We are now in the position of having a graphical look at this simple model.
\n\"\"<\/a>
\nFig. 1: Graphical representation of the simple multinomial logit regression model in the example<\/em><\/p>\n

The variable y[i] is categorically distributed with parameters l1[i], l2[i] and l3[i]. l1[i] is set to zero for indentification purposes and l2[i] and l3[i] are linearily decomposed, e.g. due to the underlying design, if we are working experimentally. x0[i] is a dummy variable for the constants (\"\beta_{20}\" and \"\beta_{30}\") of the regression equations and x1[i] is the predictor. The arrows symbolize dependencies between the nodes of the network.<\/p>\n

To show how the model’s parameters of are estimated via MCMC in a Bayesian fashion, let’s introduce a little example. Let’s assume we have sampled \"n=41\" individuals that are falling into two groups (\"x_i \in \left\{0,1 \right\}\"). Each individual had the choice between 3 alternatives (\"y_i \in \left\{1,2,3 \right\}\"). We are interested in the question: Do the groups differ with respect to the chosen alternatives? Or: Does the chosen category depend on group membership? If we crosstabulate the data in R we find:<\/p>\n

\r\n> table(y,x)\r\n   x\r\ny    0  1\r\n  1  5  3\r\n  2 13  4\r\n  3  2 14\r\n<\/pre>\n
By visual inspection we see that group 0 seems to have preference for alternative 2 and group 1 seems to have a preference for alternative 3. To model these data and to estimate the parameters via MCMC we have to translate the model into OpenBUGS code. This is straight forward from the graphical representation or the equations:<\/p>\n
\r\n#uo_multinom_reg.txt\r\nmodel\r\n{\r\n  # Likelihood\r\n  for(i in 1:n)\r\n  {\r\n    \r\n    lambda[1,i]<-0\r\n    lambda[2,i]<-beta_20+beta_21*x[i]\r\n    lambda[3,i]<-beta_30+beta_31*x[i]\r\n     \r\n    # Partition function\r\n    Z[i]<-exp(lambda[1,i])+exp(lambda[2,i])+exp(lambda[3,i])\r\n\r\n    for(c in 1:3)\r\n    {\r\n      p[i,c]<-exp(lambda[c,i])\/Z[i]\r\n    }\r\n    y[i]~dcat(p[i,1:3])\r\n  }\r\n  \r\n  # Priors\r\n  beta_20~dnorm(0,1.0E-6)\r\n  beta_21~dnorm(0,1.0E-6)\r\n  beta_30~dnorm(0,1.0E-6)\r\n  beta_31~dnorm(0,1.0E-6)\r\n}\r\n<\/pre>\n
The choice of normal priors for the parameters is motivated by the theoretical result that maximum likelihood etimates are asymptotically normally distributed. In addition, the precision of the priors is set to a very low value. So these distributions encode the prior information that we really do not know anything about the distributions of the parameters in the population yet and want to let the data speak for themselves as much as possible.<\/p>\n
The R package R2OpenBUGS is used to control the OpenBUGS script:<\/p>\n
\r\n# uo_multinom_reg.R\r\nlibrary(R2OpenBUGS)\r\ndata<-read.table(\"uo_multinom_reg.csv\", sep=\",\", head=TRUE)\r\n\r\n# data preparation\r\ny<-data$wahl\r\nx<-data$geschl\r\nn<-length(y)\r\n\r\ndata<-list(\"y\",\"x\", \"n\")\r\n\r\nparameters<-c(\"beta_20\", \"beta_30\", \"beta_21\", \"beta_31\")\r\n\r\n# Initial values for the markov chains\r\ninits <- function(){\r\nlist(beta_20 = rnorm(1), beta_21=rnorm(1), beta_30=rnorm(1), beta_31=rnorm(1))\r\n}\r\n\r\n# Handing over to OpenBUGS\r\noutput<-bugs(data, inits=inits, parameters, model.file=\"uo_multinom_reg.txt\",\r\n    n.chains=2, n.iter=8000, n.burnin=1000)\r\n\r\n# Printing the results\r\nprint(output, digits=2)\r\n<\/pre>\n
The output of OpenBUGS is<\/p>\n
\r\nInference for Bugs model at \"uo_multinom_reg.txt\", \r\nCurrent: 2 chains, each with 11000 iterations (first 1000 discarded)\r\nCumulative: n.sims = 20000 iterations saved\r\n          mean   sd  2.5%   25%   50%   75% 97.5% Rhat n.eff\r\nbeta_20   1.01 0.55 -0.02  0.64  0.98  1.36  2.13    1  2200\r\nbeta_30  -1.09 0.94 -3.22 -1.67 -1.01 -0.44  0.49    1  2600\r\nbeta_21  -0.68 1.01 -2.66 -1.35 -0.69 -0.01  1.31    1  1400\r\nbeta_31   2.78 1.16  0.76  1.98  2.68  3.48  5.35    1  2100\r\ndeviance 74.84 3.01 71.11 72.61 74.17 76.25 82.69    1 20000\r\n\r\nFor each parameter, n.eff is a crude measure of effective sample size,\r\nand Rhat is the potential scale reduction factor (at convergence, Rhat=1).\r\n\r\nDIC info (using the rule, pD = Dbar-Dhat)\r\npD = 4.1 and DIC = 79.0\r\nDIC is an estimate of expected predictive error (lower deviance is better).\r\n<\/pre>\n
The statistics \"\hat{R}\" indicate that the markov process seems to have reached a stationary distribution. <\/p>\n
The parameter beta_20<\/code> indicates a preference of category 2 over category 1 for group 0. But this difference is statistically not significant, probably due to small sample size. beta_30<\/code> indicates a preference of category 1 over category 3 for group 0, but from the percentiles of the posterior distributions we can see that this difference is also not significant statistically.<\/p>\n
beta_21<\/code> is the contrast of group 1 to beta_20<\/code>. As we can see, group 1 seems to have a slightly lower preference for category 2 over category 1 than group 0 in terms of log-probability-ratios, but this contrast is also not significant statistically. <\/p>\n
beta_31<\/code> is statistically significant. Group 1 has a much stronger preference for category 3 over category 1 than group 0 in terms of log-probabilty-ratios. <\/p>\n
Please note that this setup can be extended to include a larger set of continuous and categorical predictors. Also a multilevel extension of the model seems to be straight forward by \"simply\" including random effects into the linear regression equations. <\/p>\n
The main objective of this post was to show in a simple way, how categorical data analysis techniques translate into a graphical representation that in terms can be used to estimate model parameters with MCMC-methods.<\/p>\n
Grab the code here<\/a>.<\/p>\n","protected":false},"excerpt":{"rendered":"
I decided to switch to English as I believe that this topic may be interesting for a wider audience. In the last post I have given a simple example of how the t-test for independent samples can be understood from … Continue reading →<\/span><\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/409"}],"collection":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=409"}],"version-history":[{"count":182,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/409\/revisions"}],"predecessor-version":[{"id":706,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/409\/revisions\/706"}],"wp:attachment":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=409"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=409"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=409"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}