---
title: "Homework 03"
author: "Minji Jeong"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(MASS)
library(simcf)           
library(tile)            # For some graphics; used by plot.binPredict()library(MASS)
library(boot)            # For cv.glm()
library(pscl)            # Alternative PCP code
library(verification)    # For ROC area
library(RColorBrewer)    # For nice colors
source("https://faculty.washington.edu/cadolph/mle/binaryGOF.R")    # Percent correctly predicted and concordance indexes
source("https://faculty.washington.edu/cadolph/mle/binPredict.R")   # Code for making predicted vs actual plots
library(tidyverse)
library(data.table)
library(kableExtra)


set.seed(2020)
```


##  Modeling Cy Young winners using logit

```{r setup2}
df <- read_csv("https://faculty.washington.edu/cadolph/mle/cyyoung.csv")

# Get nice colors
col <- brewer.pal(5, "Set1")
blue <- col[2]
orange <- col[5]
```


## Answer a.
```{r warning=FALSE, message=FALSE}
m1 <- cy ~ era + winpct
m1df <- extractdata(m1, df, na.rm=TRUE) 

y <- df$cy
x1 <- m1df[,2:ncol(m1df)] %>% as.matrix()


# Likelihood function for logit
llk.logit <- function(param,y,x) {
  os <- rep(1,length(x[,1]))
  x <- cbind(os,x)
  b <- param[ 1 : ncol(x) ]
  xb <- x%*%b
  sum( y*log(1+exp(-xb)) + (1-y)*log(1+exp(xb)))
  # OR  -sum(y*xb - log(1 + exp(xb))) They are the same (Do the math)
}

# Fit logit model using optim
ls.result <- lm(y~x1)  # use ls estimates as starting values
stval <- ls.result$coefficients  # initial guesses
optim.m1 <- optim(par=stval,fn=llk.logit,method="Nelder-Mead",hessian=T,y=y,x=x1)
                   # call minimizer procedure
pe.m1 <- optim.m1$par   # point estimates
vc.m1 <- solve(optim.m1$hessian)  # var-cov matrix
se.m1 <- sqrt(diag(vc.m1))    # standard errors
ll.m1 <- -optim.m1$value  # likelihood at maximum

# Alternative estimation technique:  GLM
glm.m1 <- glm(m1, family="binomial", data=m1df)

```

```{r warning=FALSE, message=FALSE}

tibble(
  term = c("Intercept","era","winpct" ),
  "glm PE" = coef(glm.m1),
  "optim PE" = pe.m1,
  "glm SE" = summary(glm.m1)$coefficients[,"Std. Error"],
  "optim SE" = se.m1 
) %>% 
  kable() 

tibble(
  "glm LL" = logLik(glm.m1),
  "optim LL" = ll.m1
) %>% 
  kable() 

```

## Answer b.:

```{r warning=FALSE, message=FALSE, out.width='0.7\\textwidth', fig.align='center'}
## simulate parameters from model
sims <- 10000
simbetas <- mvrnorm(n=sims,mu=pe.m1,Sigma =vc.m1)

## create a matrix with scenarios to loop
x <- 
  tibble(
    # intercept
    os = rep(1, 15),
    # hypothetical scenarios
    era = seq(1.5, 5, 0.25),
    # higher-order hypo scens
    winpct = mean(df$winpct) %>% round(3)
  )

## create a matrix for simulated expected values
sim.ev <- tibble(pe = rep(NA, nrow(x)), 
                 lwr = rep(NA, nrow(x)), 
                 upr = rep(NA, nrow(x)))

## function of inverse logit for re-parametrization
inverse.logit <- function(xb){
  1 / (1 + exp (- xb))}

## loop the EV over the hypothetical scenarios

for (i in 1:nrow(x)) {
  # matrix product for each row's scenario
  xb <- simbetas %*% as.matrix(t(x[i,]))
  # re-parametrization
  res <- inverse.logit(xb)
  # expected value
  sim.ev[i,1] <- mean(res)
  # lower confidence interval at 0.95
  sim.ev[i,2] <- quantile(res,probs=0.025)
  # upper confidence interval at 0.95
  sim.ev[i,3] <- quantile(res,probs=0.975)
}

sim.ev %>% kable() %>% kable_styling()# check out

## data wrangling
dat <-
  bind_cols(sim.ev,x) %>% 
  mutate(winpct = factor(winpct)) %>% 
  mutate_if(is.numeric,round,3)

## visualization
dat %>% 
  ggplot(aes(y=pe,
             x=era)) +
  geom_line(color=col[4]) +
  geom_ribbon(aes(ymin=lwr,
                  ymax=upr),
              fill=col[4], alpha = 0.1, linetype = 0)+
  theme_classic() +
  theme(legend.position = "none") 
  
```

## Answer c.:

```{r warning=FALSE, message=FALSE}
sims <- 10000
simbetas <- mvrnorm(sims, pe.m1, vc.m1) 

xhyp <- seq(1.5, 5, 0.25) # era
nscen <- length(xhyp)

# simcf object
winpct.mean <- winpct.0.5  <- winpct.0.9  <- cfMake(m1, m1df, nscen)  


for (i in 1:nscen) {
  # different ways to arrange  
  winpct.mean <- cfChange(winpct.mean, "era", x = xhyp[i], scen = i)
  winpct.0.5  <- cfChange(winpct.0.5,  "era", x = xhyp[i], scen = i)
  winpct.0.9  <- cfChange(winpct.0.9,  "era", x = xhyp[i], scen = i)
  # mean is set as default, so we can leave it as it is for winpct.mean 
  winpct.0.5  <- cfChange(winpct.0.5, "winpct", x=0.5, scen=i)
  winpct.0.9  <- cfChange(winpct.0.9, "winpct", x=0.9, scen=i)
}

# Simulate expected probabilities 
sims.mean <- logitsimev(winpct.mean, simbetas, ci=0.95) 
sims.0.5  <- logitsimev(winpct.0.5,  simbetas, ci=0.95)
sims.0.9  <- logitsimev(winpct.0.9,  simbetas, ci=0.95)

(all.scen <- bind_rows(winpct.mean=data.frame(sims.mean),
                       winpct.0.5=data.frame(sims.0.5),
                       winpct.0.9=data.frame(sims.0.9),
                       .id="scenario") %>% 
    mutate_if(is.numeric,round,3)%>% kable() %>% kable_styling())
```


Plot by `tile`

```{r warning=FALSE, message=FALSE, results='hide'}
trace.fun <- function(sim.result, coln, cfmake, plt){
 lineplot(x=xhyp,
          y=sim.result$pe, 
          lower=sim.result$lower, 
          upper=sim.result$upper, 
          col=col[coln], 
          extrapolate=list
            (data=m1df[,2:ncol(m1df)],
             cfact=cfmake$x[,2:ncol(cfmake$x)],
             omit.extrapolated=FALSE),
          plot=plt) 
}


trace.mean <- trace.fun(sim.result = sims.mean,
                        coln = 1,
                        cfmake = winpct.mean,
                        plt = 1)

trace.0.5 <- trace.fun(sim.result = sims.0.5,
                       coln = 2,
                       cfmake = winpct.0.9,
                       plt = 2)

trace.0.9 <- trace.fun(sim.result = sims.0.9,
                       coln = 3,
                       cfmake = winpct.0.9,
                       plt = 3)


# Plot traces using tile
tile(trace.mean,
             trace.0.5,
             trace.0.9,
             xC = c(1,3),
             limits=matrix(c(1.50,5.00,0,1,1.50,5.00,0,1,1.50,5.00,0,1), 
                           nrow=3, 
                           ncol=4, byrow=TRUE),
             xaxistitle=list(labels=c("Earned Run Average", 
                                      "Earned Run Average", 
                                      "Earned Run Average")),
             columntitle=list(labels=c("Mean Winning %",".50 Winning %", ".90 Winning %")),
             maintitle=list(labels="Probability of Winning Cy Young Award at Different winning rate"))
             
```

Plot by `ggplot`
```{r warning=FALSE, message=FALSE, out.width='0.7\\textwidth', fig.align='center',}

## visualization

tb.mean <- 
  sims.mean |>
  bind_rows() |>
  mutate(
    xhyp = xhyp,  # add "xhyp" as a covariate
    winpct = "mean"  # add a column to identify which edu. scenario
  )

tb.0.5 <- sims.0.5 |>
  bind_rows() |>
  mutate(xhyp = xhyp, winpct = "0.5")

tb.0.9 <- sims.0.9 |>
  bind_rows() |>
  mutate(xhyp = xhyp, winpct = "0.9")

# bind them all together, in long format

allSims_tb <- bind_rows(tb.mean, tb.0.5, tb.0.9)

allSims_tb %>% 
  ggplot(aes(y=pe,
             x=xhyp,
             color=winpct)) +
  geom_line() +
  geom_ribbon(aes(ymin=lower,
                  ymax=upper,
                  fill=winpct), 
              alpha = 0.1, linetype = 0) +
  theme_classic() +
  labs(x="Era", y="Probabiilty of Winning")

```


The plots shows that pitchers with lower `era`, earned run averages and higher winning rates are more likely to win Cy young. Because of sample size, `winpct = 0.5 or 0.9` have large 95% CI.(Shadow line may suggest lack of data).


## Answer d.:
I added ‘natleag’ variable to the second model. Given that it’s hard to achieve low ERA in American league due to different rules, simply including ERA without pitcher’s league may obfuscate our understanding. American league pitchers who still achieve low ERA albeit hard rules may have a better chance of winning the Cy Young award, all else equal. Additionally, I included ‘strikeout’ variable because higher number of strikeout can show higher capacity of pitcher, which may increase the chance of winning the award, all else equal.

```{r warning=FALSE, message=FALSE}

m2 <- cy ~ era + winpct + strikeout + natleag
m2df <- extractdata(m2, df, na.rm=TRUE) 


x2 <- m2df[,2:ncol(m2df)] %>% as.matrix()

# Fit logit model using optim
ls.result <- lm(y~x2) 
stval2 <- ls.result$coefficients
logit.m2 <- optim(stval2,llk.logit,method="BFGS",hessian=T,y=y,x=x2)
                   
pe.m2 <- logit.m2$par   # point estimates
vc.m2 <- solve(logit.m2$hessian)  # var-cov matrix
se.m2 <- sqrt(diag(vc.m2))    # standard errors
ll.m2 <- -logit.m2$value  # likelihood at maximum


glm.m2 <- glm(m2, data=df, family="binomial")


```

## (i): A likelihood ratio test

```{r warning=FALSE, message=FALSE}
k.m1 <- length(pe.m1)
k.m2 <- length(pe.m2)

lr.test <- 2 * (ll.m2 - ll.m1)
lr.test.p <- pchisq(lr.test, df = (k.m2 - k.m1), lower.tail=FALSE)
```

The likelihood ratio is `r round(lr.test.p, 2) `. The p-value of the likelihood ratio test is smaller than conventional p=0.05 level. It implies that model 2 is a better fit to the data.

## (ii): BIC and AIC

```{r warning=FALSE, message=FALSE}
bic.m1 <- log(nrow(m1df))*k.m1 - 2*ll.m1
bic.m2 <- log(nrow(m2df))*k.m2 - 2*ll.m2
bic.test <- bic.m2 - bic.m1

aic.m1 <- 2*k.m1 - 2*ll.m1
aic.m2 <- 2*k.m2 - 2*ll.m2
aic.test <- aic.m2 - aic.m1
```

The BIC test result is `r round(bic.test,2)` and the AIC test result is `r round(aic.test,2)`. The BIC test result suggest the original model is better, while the AIC test result implies the new model is better. In instances where the two statistics yield conflicting conclusions, it must be highlighted the trade-off between AIC and BIC, considering their distinct penalties. AIC prioritizes prediction and goodness of fit, with a constant penalization on parameters. Conversely, BIC favors more parsimonious model specifications, with penalization on parameters based on sample size instead of a constant. Conflicting test results of this nature require additional checking of the Goodness-of-Fit measures.

## (iii): in-sample ROC curves


```{r warning=FALSE, message=FALSE}
binned.m1 <- binPredict(glm.m1, col = blue, bins = 15, label = "Original Model", quantiles = T)
binned.m2 <- binPredict(glm.m2, col = orange, bins = 15, label = "New Model", quantiles = T)

plot(binned.m1, binned.m2, display = "roc",
thresholds = c(0.6, 0.5, 0.4, 0.3),
labx = 0.05)


```

AUC of the new model is higher than original model, suggesting the new model has a better fit. 



## (iv): insample Actual versus Predicted plots

```{r warning=FALSE, message=FALSE}

plot(binned.m1, binned.m2, display=c("avp", "evp"), hide=TRUE, labx=0.35)

```

The actual vs. predicted plot and error vs. predicted plot show new model performs better in some area, while not in other area. It is hard to identify which model is better solely based on these plots.

## (v): cross-validation

```{r warning=FALSE, message=FALSE}
loocv <- function (obj) {
    data <- obj$data
    m <- dim(data)[1]
    form <- formula(obj)
    fam <- obj$family$family
    loo <- rep(NA, m)
    for (i in 1:m) {
        i.glm <- glm(form, data = data[-i, ], family = fam)
        loo[i] <- predict(i.glm, newdata = data[i,], family = fam, type = "response")
    }
    loo
}

predCVm1 <- loocv(glm.m1)
predCVm2 <- loocv(glm.m2)

# Make cross-validated AVP and ROC plots; note use of newpred input in binPredict
binnedM1cv <- binPredict(glm.m1, newpred=predCVm1, col=blue, bins = 15,
                         label="Original: LOO-CV", quantiles=TRUE)


binnedM2cv <- binPredict(glm.m2, newpred=predCVm2, col=orange, bins = 15,
                         label="New: LOO-CV", quantiles=TRUE)

plot(binnedM1cv, binnedM2cv, display=c("avp","roc"), hide=TRUE, thresholds=c(0.9, 0.8, 0.7, 0.6),labx=0.25)

```

It is hard to tell the difference between models based on actual vs. predicted plot. The ROC curve more clearly shows the new model performs better in general.

## Answer e.:

My variable of interest(Z) is whether the pitcher played in the national league. Cy Young is awarded to one players from each league. I am interested in how the leagues differentiate the other factors such as `era`. I will hold other covariates constant at
their mean.

```{r warning=FALSE, message=FALSE}
simbetas.m2 <- mvrnorm(sims, pe.m2, vc.m2) 
xhyp <- seq(1.5, 5, 0.25)
nscen <- length(xhyp)

aleague.sim <- 
nleague.sim <- 
  cfMake(m2, m2df, nscen)

for (i in 1:nscen) {
  
  aleague.sim <- cfChange(aleague.sim, "era", x = xhyp[i], xpre= xhyp[i], scen = i)
  aleague.sim <- cfChange(aleague.sim, "natleag", x = 0,  scen = i)
  # mean is default, but I still show these 
  aleague.sim <- cfChange(aleague.sim, "winpct", x = mean(m2df$winpct),  scen = i)
  aleague.sim <- cfChange(aleague.sim, "strikeout", x = mean(m2df$strikeout),  scen = i)
  
  nleague.sim <- cfChange(nleague.sim, "era", x = xhyp[i], xpre= xhyp[i], scen = i)
  nleague.sim <- cfChange(nleague.sim, "natleag", x = 1,  xpre = 0, scen = i)
  nleague.sim <- cfChange(nleague.sim, "winpct", x = mean(m2df$winpct),  
                          xpre = mean(m2df$winpct), scen = i)
  nleague.sim <- cfChange(nleague.sim, "strikeout", x = mean(m2df$strikeout),  
                          xpre = mean(m2df$strikeout), scen = i)
}

# Simulate expected probabilities for all scenarios
aleague.prob <- logitsimev(aleague.sim, simbetas.m2, ci=0.95)
nleague.prob <- logitsimev(nleague.sim, simbetas.m2, ci=0.95)


leagueFD <- logitsimfd(nleague.sim, simbetas.m2, ci=0.95)
leagueRR <- logitsimrr(nleague.sim, simbetas.m2, ci=0.95)
```



```{r warning=FALSE, message=FALSE, results='hide'}

# Compare expected probabilities (optional)
# Set up lineplot traces of expected probabilities 
aleague.Trace <- lineplot(x=xhyp,
                         y=aleague.prob$pe,
                         lower=aleague.prob$lower,
                         upper=aleague.prob$upper,
                         col=col[1],
                         extrapolate=list(data=m2df[,2:ncol(m2df)],
                           cfact=aleague.sim$x[,2:ncol(aleague.sim$x)],
                           omit.extrapolated=TRUE),
                         plot=1)

nleague.Trace <- lineplot(x=xhyp,
                         y=nleague.prob$pe,
                         lower=nleague.prob$lower,
                         upper=nleague.prob$upper,
                         col=col[2],
                         ci = list(mark="dashed"),
                         extrapolate=list(data=m2df[,2:ncol(m2df)],
                           cfact=nleague.sim$x[,2:ncol(nleague.sim$x)],
                           omit.extrapolated=TRUE),
                         plot=1)

# Set up traces with labels and legend
labelTrace <- textTile(labels=c("American League", "National League"),
                       x=c( 3.75, 2.9),
                       y=c( 0.65, 0.05),
                       col=col,
                       plot=1)

legendTrace <- textTile(labels=c("Logit estimates:", "95% confidence", "interval is shaded"),
                        x=c(4.5, 4.5, 4.5),
                        y=c(0.95, 0.9, 0.85),
                        cex=0.9,
                        plot=1)

# Plot traces using tile
tile(aleague.Trace,
     nleague.Trace,
     labelTrace,
     legendTrace,
     limits=c(2,5,0,1),
     xaxis=list(at=c(2,3,4,5)),
     xaxistitle=list(labels="Era"),
     yaxistitle=list(labels="Probability of Cy Young Award"),
     width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5)
     )

```

Since the American league has DH system which increases the risk of homerun, the result confirmed that pitchers in AL could be awarded Cy young with higher `era`. 

```{r warning=FALSE, message=FALSE, results='hide'}

# Plot First Difference

# Set up lineplot trace of relative risk
leagueFDTrace <- lineplot(x=xhyp,
                         y=leagueFD$pe,
                         lower=leagueFD$lower,
                         upper=leagueFD$upper,
                         col=col[2],
                         extrapolate=list(data=m2df[,2:ncol(m2df)],
                           cfact=nleague.sim$x[,2:ncol(nleague.sim$x)],
                           omit.extrapolated=TRUE),
                         plot=1)


# Set up baseline: for first difference, this is 0
baseline <- linesTile(x=c(2,5),
                      y=c(0,0),
                      plot=1)

# Set up traces with labels and legend
labelFDTrace <- textTile(labels=c("National league compared \n to American league"),
                       x=c(3.5),
                       y=c(0.3),
                       col=col[2],
                       plot=1)

legendfDTrace <- textTile(labels=c("Logit estimates:", "95% confidence", "interval is shaded"),
                        x=c(4.5, 4.5, 4.5),
                        y=c(0.9, 0.8, 0.7),
                        cex=0.9,
                        plot=1)

# Plot traces using tile
tile(leagueFDTrace,
     labelFDTrace,
     legendfDTrace,
     baseline,
     limits=c(2,5,-1,1),
     xaxis=list(at=c(2,3,4,5)),
     yaxis=list(label.loc=-0.5, major=FALSE),
     xaxistitle=list(labels="Era"),
     yaxistitle=list(labels="Difference in Probability of Cy Young Award"),
     width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5)
    
     )


# Plot Relative Risk

# Set up lineplot trace of relative risk
leagueRRTrace <- lineplot(x=xhyp,
                          y=leagueRR$pe,
                          lower=leagueRR$lower,
                          upper=leagueRR$upper,
                          col=col[2],
                          extrapolate=list(data=m2df[,2:ncol(m2df)],
                          cfact=nleague.sim$x[,2:ncol(nleague.sim$x)],
                          omit.extrapolated=TRUE),
                          plot=1)


# Set up baseline: for relative risk, this is 1
baseRRline <- linesTile(x=c(2,5),
                      y=c(1,1),
                      plot=1)

# Set up traces with labels and legend
labelRRTrace <- textTile(labels=c("National league compared \n to American league"),
                       x=c(3.5),
                       y=c(1.5),
                       col=col[2],
                       plot=1)

legendRRTrace <- textTile(labels=c("Logit estimates:", "95% confidence", "interval is shaded"),
                          x=c(4.5, 4.5, 4.5),
                          y=c(0.5, 0.4, 0.3),
                          cex=0.9,
                          plot=1)

# Plot traces using tile
tile(leagueRRTrace,
     labelRRTrace,
     legendRRTrace,
     baseRRline,
     limits=c(2,5,0,2),
     xaxis=list(at=c(2,3,4,5)),
     yaxis=list(label.loc=-0.5, major=FALSE),
     xaxistitle=list(labels="Era"),
     yaxistitle=list(labels="Relative Risk of Cy Young Award"),
     width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5))

```

Both the first difference in the probability and the relative risk confirm the differences.The difference is larger for mid-range era. The CI of first difference includes zero for most of the values, indicating it’s not statistically significant. The CI or relative risk also overlaps with 1, suggesting a statistically insignificant result.

Below, I repeat the simulation for different value of a covariate, winpct. Before I hold winpct constant at its mean, now I set winpct at 0.5.

```{r warning=FALSE, message=FALSE}
simbetas.m2 <- mvrnorm(sims, pe.m2, vc.m2) 
xhyp <- seq(1.5, 5, 0.25)
nscen <- length(xhyp)

aleague.sim <- 
nleague.sim <- 
  cfMake(m2, m2df, nscen)

for (i in 1:nscen) {
  
  aleague.sim <- cfChange(aleague.sim, "era", x = xhyp[i], xpre= xhyp[i], scen = i)
  aleague.sim <- cfChange(aleague.sim, "natleag", x = 0,  scen = i)
  # mean is default, but I still show these 
  aleague.sim <- cfChange(aleague.sim, "winpct", x = 0.5,  scen = i)
  aleague.sim <- cfChange(aleague.sim, "strikeout", x = mean(m2df$strikeout),  scen = i)
  
  nleague.sim <- cfChange(nleague.sim, "era", x = xhyp[i], xpre= xhyp[i], scen = i)
  nleague.sim <- cfChange(nleague.sim, "natleag", x = 1,  xpre = 0, scen = i)
  nleague.sim <- cfChange(nleague.sim, "winpct", x = 0.5,  
                          xpre = 0.5, scen = i)
  nleague.sim <- cfChange(nleague.sim, "strikeout", x = mean(m2df$strikeout),  
                          xpre = mean(m2df$strikeout), scen = i)
}

# Simulate expected probabilities for all scenarios
aleague.prob <- logitsimev(aleague.sim, simbetas.m2, ci=0.95)
nleague.prob <- logitsimev(nleague.sim, simbetas.m2, ci=0.95)


leagueFD <- logitsimfd(nleague.sim, simbetas.m2, ci=0.95)
leagueRR <- logitsimrr(nleague.sim, simbetas.m2, ci=0.95)
```

You can visualize plots following the same procedure as above. 

## Answer f.:

The logistic regression model assumes that outcome is a binary and an independent. Also, the model assumes that there is a linear relationship between the logit of the outcome and each covariate. However, the distribution of outcome variable, `cy` is not independent since the winners are limited to two players from each league per year. Specifically, the binary logit model assumes independent Bernoulli trials (conditional on covariates), when in reality the DGP follows a rank-type process. In other words, the independence assumption is violated because only one player in the entire league can win the award, and the probability of winning is conditional on the performance of other contestants during that year, on top of the performance of the player who wants to win the award, all else being equal.


## Answer g.:

The question asks if we have a sampling bias or not. The expert may considers that all covariates are normally distributed and choose the top tier of the sample. If so, our results would not be changed drastically. The full dataset of all pitchers shows that the variance of our estimates would be small by increasing sample size.
However, if the expert's selection is skewed to the specific ability of the pitchers, we may have overfitting problem on our estimation. (though we should take into account that the data is zero-inflated).


