---
title: 'Edpsy 589:  Discrete Choice Models:  Conditional Multinomial Logistic Regression'
author: "Carolyn J. Anderson"
date: "6/29/2022"
output: 
    html_document:
       toc: yes
---

How to fit discrete choice models, sometimes called conditional multinomial logistic regression models, are illustrated below. Special cases of these models include the baseline multinomial logistic regression model.  

Two methods are illustrated, one of which uses glm and the other is mlogit.  The latter has two depreciated functions (mlogit.data and mFormula), so code that previously worked will no longer work unless you use archived version of mlogit.  However, the commands below use the non-depreciated functions and package that replace these.  The mlogit.data is replace by dfidx package and mFormula is replaced by Formula.

Three data sets are used in this document:


1.  Chocolates (from SAS) 
2.  Brands and prices (from SAS)
3.  Modes of Transportion (from Powers & Xie)


# Set up

The libraries needed
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(mlogit)
library(dfidx)    # for indexing data (replaces mlogit.data)
library(Formula)  # for createing formaula (resplace mFormula)
```

# Chocolate Data Set

This is a hypothetical data set from SAS book on using logistic regression.  A person chooses one of eight types of chocolates,  Attributes of the chocolates are type of chocolate (milk/dark), whether is contain nuts (nuts/no nuts), and the texture (hard/soft).  

```{r choco}
choc <- read.table("chocolates_data.txt",header=TRUE)
head(choc, n=10)
```

## a. glm (AS a Poisson regression)

The model can be fit simply using
```{r choc_glm}
summary( glm.choc <- glm(choice ~ dark + soft + nuts, data=choc, family=poisson) )
```

We want the probabilities for each type

```{r choc_probs}
beta <- glm.choc$coefficients

#For the probabilities:		 
pi0 <- exp(beta[1]+ beta[2]*choc$dark[1:8] +beta[3]*choc$soft[1:8] + beta[4]*choc$nuts[1:8])
(pi.hat <- pi0/sum(pi0))

# For easier of knowing what's what
probs <- cbind(choc[1:8,],pi.hat)
(pi.hat <- probs[order(probs$pi.hat),])
```

In this example, chooses have the same choice probabilities.  

## b. mlogit


Need to create a logical variable where "TRUE" for chocolate picked "FALSE" otherwise
```{r  }
choc$pick <-  choc$choice==1
```

Create strings for attributes of chocolates
```{r}
choc$dark.s <- ifelse(choc$dark==1, "Dk","Mk")
choc$soft.s <- ifelse(choc$soft==1, "Sft","Hrd")
choc$nuts.s <- ifelse (choc$nuts==1, "Nut","Nnts")
```

Concatenate strings to create a variables for each attribute
```{r}
choc$alt<- paste(choc$dark.s,choc$soft.s,choc$nuts.s,sep="_")
head(choc$alt, n=10)
```

This is what replaced mlogit.data
```{r} 
c <- dfidx(choc, shape = "long", idx=list("subj", "alt"), choice="pick")
head(c, n=10)
```

The basic syntax of mlogit formulas is 
  
  
  (logical variable for choice ~   generic alternative specific | individual specific |  alternative specific )
```{r}
f <- Formula(pick ~ dark + soft + nuts | 0 | 0) 

model <- mlogit(f,c)
summary(model)
```

What I am interested in  are 
```{r}
probs <- model$probabilities[1,]
probs <- sort(probs)

as.matrix(probs,nrow=8,ncol=1)
```

Exactly the same as found using glm.

# Brands and Prices

Objects in this example has different brands of the same product and different prices.

```{r brands}
brands <- read.table("brands_price_data.txt",header=TRUE)
head(brands)
```

## a.  glm

Need to do a bit of data manipulation
```{r}
brands$combo <- as.factor(brands$combo)
brands$brand <- as.factor(brands$brand)
```

I changed the reference category
```{r }
brands <- within(brands, brand <- relevel(brand, ref = "5"))
head(brands)
```

The variable combo is a nuisance variable that will be ignored but needed in the glm.

### Simple Model: Additive Effects

```{r}
summary( mdc1 <- glm(choice ~ combo + brand + price, data=brands, family=poisson) )
anova(mdc1,test="LRT")
```

### Interaction between Brand and Price
```{r}
summary( mdc2 <- glm(choice ~ combo + brand + price + brand*price, data=brands, family=poisson) )
anova(mdc2,test="LRT")

# Test between the two models
anova(mdc1,mdc2)
```

### The probabilities from final model (simple)
```{r}
beta <- glm.choc$coefficients
#For the probabilities:		 
pi0 <- exp(beta[1]+ beta[2]*choc$dark[1:8] +beta[3]*choc$soft[1:8] + beta[4]*choc$nuts[1:8])
(pi.hat <- pi0/sum(pi0))

# For easier of knowing what's what
probs <- cbind(choc[1:8,],pi.hat)
(pi.hat <- probs[order(probs$pi.hat),])
```

## b. mlogit

Use the data that's in individual (long) format
```{r long}
longer <- read.table("D:\\Dropbox\\edps 589\\8 Multicategory_logit\\brands_mdc_data.txt",header=TRUE)
```

Choice to be logical
```{r}
longer$choice <- longer$Y==1

longer$brand <- "brand5"
longer$brand[longer$brand1==1] <- "brand1"
longer$brand[longer$brand2==1] <- "brand2"
longer$brand[longer$brand3==1] <- "brand3"
longer$brand[longer$brand4==1] <- "brand4"
```

Create alternatives' descriptions in terms of price and brand
```{r}
longer$alternative <- paste(longer$brand,longer$price,sep="_")

#Set up data for mlogit
```{r}
long.m <- dfidx(longer, shape="long", idx=list("case","alternative"), choice="choice")
```


### Simpel model:  additive effects

```{r}
fm1 <- Formula(choice ~ brand1 + brand2 + brand3 + brand4 + price | 0 | 0 ) 

summary(model.1<- mlogit(fm1,long.m) )
```


Estimated odds ratios
```{r}
(exp.beta <- exp(model.1$coefficients[1:5]) )

### Model with Interaction between brand an price

```{r}
fm2 <- Formula(choice ~ brand1 + brand2 + brand3 + brand4 + price + brand1*price + brand2*price + brand3+price + brand4*price | 0 | 0 ) 

summary(model.2<- mlogit(fm2,long.m) )

# Test interaction:
(lr <- 2*(as.numeric(model.2$logLik) - as.numeric(model.1$logLik)))
1-pchisq(lr, 3)
```

### Using the simpler model (final model)
```{r bandpprobs}
(beta <- model.1$coefficients)
```


# Modes of Transportation

The data taken from Powers & xie and is a classic example of choice of modes of transportation.  This example includes bother attributes of types of transportation as well as chooser's attributes.  The  Variables used are:


1.   mode= Mode of transportation chosen, 1st=train, 2nd=bus, 3rd=car
2.   ttime= Time in terminal   
3.   invc= Time in vehicle 
4.   gc = Generalized cost 
5.   hinc = Household income 
6.   id = person id


Only mlogit is used in this example

```{r}
trans <- read.table("transportation_data.txt",header=TRUE)

n <- length(unique(trans$id))
trans$alt <- rep(c("Bus","Train","Car"),n) 

# set up data for mlogit
tm <- dfidx(trans, shape="long", idx=list("id","alt"), choice="mode")
```


### Model with alternative specific (i.e.,  Attributes of modes of transportation )
 Recall: 
 
(logical variable for choice ~   generic alternative specific | individual specific |  alternative specific )

 Attributes of modes of transportation 
```{r}
fm.a <- Formula(mode ~ ttme + invc + invt + gc | 0  | 0) 

# Fit the model
summary(model.a <- mlogit(fm.a,tm))
```


### Attributes of modes and individual speicfic

This is reported in lecture notes:

```{r}
fm.b <- Formula(mode ~ ttme  + invc + invt + gc | hinc  | 0 ) 
summary(model.b <- mlogit(fm.b,tm))
```

Run through console to get R script

knitr::purl("R_discrete_choice_models.Rmd", 
          output = "R_discrete_choice_models.R", documentation=1)