---
title: "Multiclass logistic regression"
author: "Milica Cudina"
output:
  pdf_document: default
---

For a similar analysis, look at this
[tutorial](https://stats.oarc.ucla.edu/r/dae/multinomial-logistic-regression/) from UCLA. For more about `glmnet`, consult this [tutorial](https://glmnet.stanford.edu/articles/glmnet.html#multinomial-regression-family-multinomial) by Trevor Hastie. 

First, we import the data.

```{r}
data<-read.csv("hsbdemo.csv")
#View(data)
data
```
The data set contains variables on $200$ students. We will focus on a small subset of predictors. The predictor variables will be socio-economic status `ses` (a three-level categorical variable) and writing score `write` (a quantitative variable). The outcome variable is program type `prog` (a three-level categorical variable). Since I am interested in just this subset, I will create a smaller data frame to analyze. 

```{r}
df=data.frame(data$ses, data$write, data$prog)
colnames(df)<-c("ses", "write", "prog")
attach(df)
```

Some exploratory data analysis is called for. The first idea is, probably to try the `plot` command. 

```{r}
plot(df,
     pch=20, col="lightblue")
```

As we can see, this is not too useful. To look at the *association* between `ses` and `prog`, a two-way table might do the trick.

```{r}
tab=table(ses, prog)
tab
mosaicplot(tab)
```

Which test might we use to test the independence hypothesis? 

```{r}
chi2<-chisq.test(tab)
chi2
```
What about the association between the writing score and the program type? 


Side-by-side boxplots might give insight. 

```{r}
write.ac=write[which(prog=="academic")]
write.gen=write[which(prog=="general")]
write.voc=write[which(prog=="vocation")]

boxplot(write.ac,write.gen, write.voc,
  main = "Writing scores by program type",
  ylab = "Writing scores",
  names = c("Academic", "General", "Vocational"),
  col=c("lightblue", "orange", "lightgreen")
)
```


We see that there is an association, but we will learn more about the extent of the effect if we run a multiclass logistic regression. 

## `multinom` 

Out first option is the `multinom` implementation in the `nnet` package. 

```{r}
#install.packages("nnet")
library(nnet)
```

This implementation requires of us to set a *baseline* level for the response variable. 

```{r}
df$prog<-as.factor(df$prog)
df$prog.base <- relevel(df$prog, ref = "academic")
```

Now, we can run the `multinom` function. 

```{r}
mn.fit<-multinom(prog.base ~ ses + write, data=df)
```
Here is what the object `mn.fit` looks like:

```{r}
summary(mn.fit)
```
Importantly, the coefficients above are **different** from the coefficients one sees in our slides. However, let us look at the fitted probabilities for our data set. 

```{r}
fitted.ps <- fitted(mn.fit)
fitted.ps
```

Also importantly, we do not have $p-$values reported above. We can figure them out easily, though. 

```{r}
z.scores<- summary(mn.fit)$coefficients/summary(mn.fit)$standard.errors
p.val=2*(1-pnorm(abs(z.scores)))
p.val
```

How well did we do on the training set?

```{r}
pred.mn.tr=predict(mn.fit, newdata = df, "class")
#pred.mn.tr
ind.tr=(pred.mn.tr==prog)
mean(ind.tr)
```

How would we predict for a test case? First, we need to define what the inputs for the test case would be. 

```{r}
new.data <- data.frame(ses = c("low", "middle", "high"), write = c(quantile(df$write, 0.25), median(df$write), quantile(df$write, 0.75)))
```

Now, we can use the `predict` command to see what the predicted probabilities would be.

```{r}
pred.mn=predict(mn.fit, newdata = new.data, "probs")
pred.mn
```


## `glmnet`

Here, we need to install and load a package created by Hastie et al. 

```{r}
#install.packages("glmnet")
library(glmnet)
```

We notice that we need to transform our categorical predictor using dummy variables. For that, the following library is useful: 

```{r}
library(caret)
```

First, I will put all of my predictors into a separate data frame. 

```{r}
temp.x=data.frame(data$ses, data$write)
```

Now, I will use the command `dummyVars` to create dummy variables for my one categorical predictor `ses` in my array of predictors.

```{r}
dummy.x <- dummyVars(" ~ .", data = temp.x)
x <- as.matrix(predict(dummy.x, newdata = temp.x))
print(x)
```


Now, we would like to fit the model using `glmnet`. 

```{r}
y <- df$prog
#y
glm.fit <- glmnet(x, y, family = "multinomial")
#summary(glm.fit)
#glm.fit
predict(glm.fit, newx=x, s=0, type="response")
```

What if I try to find, by cross-validation, an optimal tuning parameter? 

```{r}
cv.fit <- cv.glmnet(x, y, family = "multinomial")
plot(cv.fit)
```

Now, we can try to predict at the optimal $\lambda$:

```{r}
predict(cv.fit, newx = x, s = "lambda.min", type = "class")
```

We might want to see how well we did on the training data. 

```{r}
preds=predict(cv.fit, newx = x, s = "lambda.min", type = "class")
#preds
ind=(preds==y)
mean(ind)
```