statistical analysis.Rmd

---
title: "statistical test"
author: "Jiawen Zhao"
date: "12/5/2022"
output: html_document
---

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

```{r, echo = FALSE, include=FALSE}
library(tidyverse)
library(corrplot)
library(glmnet)
library(modelr)
```

Correlation matrix between predictor variables

```{r}
data <- read.csv("./data/worldcup_final.csv") %>% 
  na.omit() %>% 
  select(-country)
#hist.data.frame(data, nclass = 8)
num = data[2:14] %>% 
  select(-gd, -player)
#summary(num)
corrplot(cor(num),  diag = FALSE)
```

Model selection: stepwise method

Model selection that didn't work.
```{r}
model_step_wise = step(glm(w ~ ., data = data), direction = "both")
model_forward = step(glm(w ~ 1, data = data), direction = "forward", scope = formula(lm(gf ~ ., data = data)), test = "F")
model_backward = step(glm(w ~ ., data = data), direction = "backward", test = "F")
model_step_wise
model_forward
model_backward
```

```{r}
predictor_state = num%>%
  select(-w)

lambda_seq <- 10^seq(-3, 0, by = .1)
set.seed(2022)
cv_object <- cv.glmnet(as.matrix(predictor_state), num$w,
lambda = lambda_seq,
nfolds = 5)
cv_object
tibble(lambda = cv_object$lambda,
mean_cv_error = cv_object$cvm) %>%
ggplot(aes(x = lambda, y = mean_cv_error)) +
geom_point()
lambda = cv_object$lambda.min

fit_bestcv <- glmnet(as.matrix(predictor_state), num$gf, lambda = cv_object$lambda.min)
coef(fit_bestcv)
lambda

```


Cross-validation of models

Cross validation - 10 fold
```{r}
cv_df = 
  crossv_kfold(data, k = 10) %>%
  mutate(
    train = map(train, as_tibble),
    test = map(test, as_tibble)
  ) %>%
  mutate(
    fit = map(train, ~glm(w ~ pld + d + rank+ gf, data = .x)),
    rmse = map2_dbl(fit, test, ~rmse(model = .x, data = .y))
  )

models <- map(cv_df$train, ~ glm(w ~ pld + d + rank+ gf, data = .))
summary(map2_dbl(models, cv_df$test, modelr::rmse))

cv_df %>%
  ggplot(aes(x = rmse)) +
  geom_density()


cv_df2 = 
  crossv_kfold(data, k = 10) %>%
  mutate(
    train = map(train, as_tibble),
    test = map(test, as_tibble)
  ) %>%
  mutate(
    fit = map(train, ~glm(w ~ pld + d + rank + gf + ga, data = .x)),
    rmse = map2_dbl(fit, test, ~rmse(model = .x, data = .y))
  )

models <- map(cv_df2$train, ~ glm(w ~ pld + d + rank + gf+ ga, data = .))
summary(map2_dbl(models, cv_df2$test, modelr::rmse))

cv_df2 %>%
  ggplot(aes(x = rmse)) +
  geom_density()
```


cross validation with violin plot
```{r}
cv_df =
  crossv_mc(data, 100) %>% 
  mutate(
    train = map(train, as_tibble),
    test = map(test, as_tibble))

cv_df =
  cv_df %>% 
  mutate(
    train = map(train, as_tibble),
    test = map(test, as_tibble))

cv_df = 
  cv_df %>% 
  mutate(
    fit1  = map(train, ~lm(w ~ pld + d + rank, data = .x)),
    fit2  = map(train, ~lm(w ~ pld + d + rank + gf, data = .x)),
    fit3  = map(train, ~lm(w ~ pld + d + rank + gf + ga, data = .x))) %>% 
  mutate(
    rmse_fit1 = map2_dbl(fit1, test, ~rmse(model = .x, data = .y)),
    rmse_fit2 = map2_dbl(fit2, test, ~rmse(model = .x, data = .y)),
    rmse_fit3 = map2_dbl(fit3, test, ~rmse(model = .x, data = .y)))

cv_df %>% 
  select(starts_with("rmse")) %>% 
  pivot_longer(
    everything(),
    names_to = "model", 
    values_to = "rmse",
    names_prefix = "rmse_") %>% 
  mutate(model = fct_inorder(model)) %>% 
  ggplot(aes(x = model, y = rmse)) + geom_violin()+
  ylab("rmse Value") +
  xlab("Models") +
  labs(title = "Violin Plot of rmse Values for Three Models")
```


Final regression model output


```{r}
summary(glm(w ~ pld + d + rank + gf+ ga, data = data))
```

Results/output from regression model


This model will be w = -1.554369 + 0.661219pld - 0.622216d + 0.015920rank + 0.153794gf -0.225378ga.
That means if the number of number of goals scored against a country in the World Cup (since 1990) increases by 1 goal, the would increase by 2.7789.


```{r}
filter_data = data %>% 
  select( -player) %>%  
  na.omit()
m=glm(w~., data = filter_data)
m=lm(w~pld+d+rank, data = num)
summary(m)

```