kickstarter_analysis.qmd

---
title: "Kickstarter Dataset Analysis"
author:
    - name: Thibault Senegas
      email: senegas.th@gmail.com
format:
  html:
    code-fold: true
    code-summary: "Show the code"
    toc: true
    toc-location: left
    warning: false
    message: false
---

# Executive Summary

Dans le cadre de cette analyse technique en science des données, j'ai choisi d'utiliser et découvrir le nouvel IDE de Posit (anciennement Rstudio). Cet outil novateur offre une intégration fluide entre différents langages de programmation, notamment R et Python. En tant qu'utilisateur principal de R, cette approche m'a permis de réaliser l'analyse de manièere rapide en utilisant R tout en fournissant une version de Python au sein du même projet.

J'ai voulu faire l'analyse du jeu de donnée Kickstarter de manière simple et rapide. La réalisation de l'analyse en R a pris environ 2h30. Dans l'analyse en Python, curieux de voir si je pouvais améliorer mon premier modèle réalisé en R, j'ai décidé de tester directement un modèle lightGBM et sélectionné des variables un peu différentes - la traduction basique du code R en Python ne me semblant pas très intéressante.


# Packages for Analysis
```{r}
library(tidyverse)
library(dplyr)
library(skimr)
library(lubridate)
library(highcharter)
library(tidytext)
library(textdata)
library(tidymodels)
library(vip)
```


# EDA

## Load and clean raw dataset
```{r}
df_raw = read.csv("data/ks_dataset.csv")

```

Il apparaît que certaines observations dans le jeu de données sont décalées, générant ainsi quatre colonnes presque vides (X à X.3). Nous pourrions prendre le temps de corriger ces observations, mais étant donné qu'elles sont peu nombreuses, j'ai choisi de ne conserver que les observations correctes.

```{r}
existring_category = as.data.frame(table(df_raw$category)) |>
    subset(Freq >= 10 )

df = df_raw |>
    subset(category %in% existring_category$Var1) |>
    select(-X, -X.1, -X.2, -X.3)
```

Nous disposons désormais d'un jeu de données propre, avec une perte limitée à seulement 632 observations.

## Cleaning

```{r}
df = df |>
    mutate(
        category = factor(category),
        main_category = factor(main_category),
        currency = factor(currency),
        state = factor(state),
        country = factor(country),
        deadline = lubridate::ymd_hms(deadline),
        launched = lubridate::ymd_hms(launched),
        goal = as.numeric(goal),
        pledged = as.numeric(pledged),
        backers = as.numeric(backers),
        usd.pledged = as.numeric(usd.pledged)
        )

skimr::skim(df)

```

Examinons la distribution des projets lancés au fil des années :

```{r}
table(lubridate::year(df$launched)) 
```

Nous constatons que 7 projets ont été lancés en 1970 et seulement 1 324 en 2009. Nous allons supprimer ces projets afin de travailler avec des années complètes.

```{r}
df = df |>
    dplyr::filter(!(lubridate::year(launched) %in% c(1970, 2009)))
```

Faisons une dernière vérification pour nous assurer qu'il n'y a pas de valeurs manquantes.

```{r}
sapply(df, function(x) sum(is.na(x)))
```

Nous avons une observation avec un 'name' manquant et 3 790 observations avec un 'usd.pledged' manquant. Par souci de simplicité, nous allons retirer ces observations. Il aurait également été possible de retrouver le taux de change du jour et de réintégrer la valeur correcte pour 'usd.pledged'.

```{r}
df = df[complete.cases(df), ]
```


Pour simplifier l'analyse, nous allons recoder notre variable cible 'state' en une variable binaire. Par la même occasion, nous supprimerons les projets ayant un état 'undefined' ou 'live'.

```{r}
df = df |>
    dplyr::filter(state %in% c("failed", "canceled", "successful", "suspended")) |>
    mutate(
        state_binary = factor(ifelse(state == "successful", 1, 0))
    )

table(df$state_binary)
```

Sur les 313,569 projets dans notre jeu de données, 112,400 ont été financés, ce qui représente environ 36 %. Bien que notre jeu de données ne soit pas parfaitement équilibré, nous allons le garder tel quel.

## Data Visualizations

```{r}

plot_success_per_year <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
     group_by(year_launch, state_binary) |>
     summarise(n_proj = n()) |>
    highcharter::hchart('column', hcaes(x = year_launch, y = n_proj, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Nombre de projets 'successful' et 'failed' par années") |>
        hc_xAxis(title = list(text = "Année de lancement")) |>
        hc_yAxis(title = list(text = "Nombre de projets")) 
 
 
plot_success_per_year
```

```{r}

plot_success_percent_per_year <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1))
    ) |>
     group_by(year_launch) |>
     summarise(succ_rate = round((mean(as.numeric(state_binary) - 1)) * 100)) |>
    highcharter::hchart('column', hcaes(x = year_launch, y = succ_rate),
    stacking = 'normal'
    ) |>
        hc_title(text = "Taux de réussite par années") |>
        hc_xAxis(title = list(text = "Année de lancement")) |>
        hc_yAxis(title = list(text = "Taux réussite (%)")) 
 
 
plot_success_percent_per_year
```

```{r}

plot_usd_invest_per_year <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
     group_by(year_launch, state_binary) |>
     summarise(pledged_per_year = round(sum(usd.pledged) / 1000000)) |> 
    highcharter::hchart('column', hcaes(x = year_launch, y = pledged_per_year, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Total d'argent investis par années") |>
        hc_xAxis(title = list(text = "Année de lancement")) |>
        hc_yAxis(title = list(text = "USD en million")) 
 
plot_usd_invest_per_year
```

```{r}

plot_avg_goal_year_usd <- df |>
    subset(currency == 'USD') |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
     group_by(year_launch, state_binary) |>
     summarise(goal_per_year = round(mean(goal))) |> 
    highcharter::hchart('column', hcaes(x = year_launch, y = goal_per_year, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Moyenne des Goal des projets (Projets en USD)") |>
        hc_xAxis(title = list(text = "Année de lancement")) |>
        hc_yAxis(title = list(text = "USD")) 
 
plot_usd_invest_per_year
```

Avec ces 4 visualisations simples nous pouvons voir les choses suivantes :
- Le nombre de projets déposés sur kickstarter a augmenté pour atteindre son piqu en 2015
- Cette augmentation ne se traduit pas par un taux de succés stable, au contraire. 
- Si l'on regarde la somme d'argent investis par années (en millions de USD), en séparant les projets financés et ceux qui échouent, nous pouvons penser que c'est quitte ou double pour un projet. Soit celui-ci est financé totalement soit il recoit très peu d'investissement et il échoue.
- En nous basant uniquement sur les projets en USD, la moyenne pour les projets qui réussisse de l'objectif (en USD) sont nettement inférieur à la moyenne des objectifs (en USD) des projets qui échouent.

Explorons maintenant les category / main_category des projets, leur taux de réussite, ainsi que la durée des projets déposés et le pays dans lequel sont déposés les projets.

Avec ces quatre visualisations simples, nous pouvons observer les points suivants :

- Le nombre de projets déposés sur Kickstarter a augmenté, atteignant un pic en 2015.
- Cette augmentation n'a pas entraîné un taux de succès stable ; au contraire, le taux de succès semble diminuer.
- En examinant les sommes investies par année (en millions de USD), en distinguant les projets financés de ceux qui échouent, on constate que les projets semblent avoir une issue binaire : soit ils sont complètement financés, soit ils reçoivent très peu d'investissements et échouent.
- En nous basant uniquement sur les projets en USD, la moyenne des objectifs de financement (en USD) pour les projets ayant réussi est nettement inférieure à la moyenne des objectifs de financement (en USD) des projets ayant échoué.

Explorons maintenant les catégories et sous-catégories des projets, leur taux de réussite, ainsi que la durée des projets et les pays d'origine.

```{r}

df <- df |>
    mutate(
        duration = round(deadline - launched)
    )
```

```{r}

plot_main_category_success <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
        group_by(main_category, state_binary) |>
        summarise(n_proj = n()) |>
        arrange(desc(n_proj))|> 
        highcharter::hchart('column', hcaes(x = main_category, y = n_proj, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Nombre de projets Failed & Succed par main_catergory") |>
        hc_xAxis(title = list(text = "Main Category")) |>
        hc_yAxis(title = list(text = "nombre projet")) 
 
plot_main_category_success
```

```{r}

plot_duration_success <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
        group_by(duration, state_binary) |>
        summarise(n_proj = n()) |>
        arrange(desc(n_proj))|> 
        highcharter::hchart('column', hcaes(x = duration, y = n_proj, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Nombre de projets Failed & Succed par durée") |>
        hc_xAxis(title = list(text = "Durée (jours)")) |>
        hc_yAxis(title = list(text = "nombre projet")) 
 
plot_duration_success
```

```{r}

plot_country_success <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
        group_by(country, state_binary) |>
        summarise(n_proj = n()) |>
        arrange(desc(n_proj))|> 
        highcharter::hchart('column', hcaes(x = country, y = n_proj, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Nombre de projets Failed & Succed par pays") |>
        hc_xAxis(title = list(text = "Pays")) |>
        hc_yAxis(title = list(text = "nombre projet")) 
 
plot_country_success
```

```{r}

plot_country_success_perc <- df |>
    mutate(
        year_launch = lubridate::year(launched),
        state_binary = factor(state_binary, levels = c(0, 1), labels = c("Failed", "Success"))
    ) |>
        group_by(country, state_binary) |>
        summarise(n_proj = n()) |>
        arrange(desc(n_proj))|> 
        ungroup() |>
        group_by(country) |>
        mutate(
            total_proj = sum(n_proj),
            perc = n_proj / total_proj * 100
        ) |>
            highcharter::hchart('column', hcaes(x = country, y = perc, group = state_binary),
    stacking = 'normal'
    ) |>
        hc_colors(c("rgb(155, 20, 20)", "rgba(8, 160, 31, 0.5)")) |>
        hc_title(text = "Nombre de projets Failed & Succed par Pays (en %)") |>
        hc_xAxis(title = list(text = "Pays")) |>
        hc_yAxis(title = list(text = "")) 
 
plot_country_success_perc
```

En examinant ces visualisations supplémentaires, nous constatons que la grande majorité des projets ont une durée de 30 jours et sont lancés aux États-Unis. L'objectif de financement (goal) en USD semble avoir un impact significatif sur le critère de réussite ou d'échec du financement du projet.

## Features selections pour le modèle

- main_category
- goal
- duration
- Nous allons filtrer nos données pour faire un modèle uniquement pour les projets aux États-Unis

Enfin, il semble important de prendre en compte le nom du projet. Nous allons calculer un score d'analyse de sentiment à partir des noms des projets et l'inclure dans notre modèle.

# Sentiment Analysis - Name
```{r}
df = df |>
    subset(country == "US")
```

## AFINN Lexicon

Nous aurions pu utiliser différents lexiques et comparer les résultats. Pour l'instant, je vais continuer avec le lexique AFINN.

```{r}
df = df |>
    subset(country == "US")
```

```{r}
df_afinn = df |>
    tidytext::unnest_tokens(word, name) |>
    inner_join(tidytext::get_sentiments("afinn"), by = "word") |>
    group_by(ID) |>
    mutate(
        overall_sent = sum(value)
    ) |>
    select(ID, overall_sent) |>
    unique()

df_afinn_clean <- left_join(df, df_afinn, by = "ID")

df_afinn_clean$overall_sent |> is.na() |> sum()
        
```

Nous constatons que de nombreux projets ont un score de sentiment manquant (NA). Nous ferons l'hypothèse que ces projets ont un sentiment neutre et nous imputerons une valeur de 0 pour le sentiment à ces projets.

```{r}

df_afinn_clean <- df_afinn_clean %>% 
  replace_na(list(overall_sent = 0))
```

Regarons avec un modèle très simple si le score de sentiment du nom d'un projet peut avoir un effet sur la réussite de celui-ci :

```{r}

summary(glm(state_binary ~ overall_sent, data = df_afinn_clean, family = "binomial"))
```

Le score de sentiment a un effet significatif, bien que son poids soit relativement faible. Nous conserverons tout de même cette variable dans notre modèle, notamment pour d'éventuelles analyses plus approfondies par catégorie ou sous-catégorie.

# Data Modeling

## Data prep

Voici les variables que nous allons selectionnés pour notre modèle :
- main_category
- goal
- duration
- overall_sent

```{r}
df_model = df_afinn_clean |>
    select(state_binary, overall_sent, main_category, goal, duration) |>
    mutate(
        duration = as.numeric(duration)
    )
```

## Data Modelling & Prep

Nous allons commencer par tester un modèle de régression logistique pénalisée. Ce modèle présente plusieurs avantages, notamment :

- Intégration des techniques Lasso (régularisation L1) et Ridge (régularisation L2)
- Prévention de l'overfitting
- Sélection des variables (avec Lasso)
- Réduction de la multicolinéarité (avec Ridge)

### Data Splitting & resampling

Ici, plutôt que d'utiliser plusieurs itérations de rééchantillonnage, nous avons créé un seul rééchantillon, appelé val_set. Le graphique ci-dessous illustre notre partitionnement des données et notre rééchantillonnage.

![](validation-split.svg)


```{r}
set.seed(123)

splits <- initial_split(df_model, strata = state_binary)

df_other <- training(splits)
df_test <- testing(splits)

val_set <- validation_split(df_other, 
                            strata = state_binary, 
                            prop = 0.80)
```

##  penalized logistic regression

```{r}
lr_mod <- logistic_reg(penalty = tune(), mixture = 1) |>
    set_engine("glmnet")
```

Ici, nous réglons le paramètre de mixture à 1 pour obtenir le modèle le plus simple possible, en privilégiant ainsi la régularisation Lasso.

```{r}
lr_recipe <- 
  recipe(state_binary ~ ., data = df_other) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_zv(all_predictors()) %>% 
  step_normalize(all_predictors())

lr_workflow <- 
  workflow() %>% 
  add_model(lr_mod) %>% 
  add_recipe(lr_recipe)
```

- step_dummy(): Convertit les variables catégorielles (caractères ou facteurs) en un ou plusieurs termes binaires numériques représentant chaque niveau des données d'origine.

- step_zv(): Supprime les variables indicatrices qui contiennent uniquement une valeur unique (par exemple, des zéros). Cela est important pour les modèles pénalisés, car les prédicteurs doivent être centrés et mis à l'échelle.

- step_normalize(): Centre et met à l'échelle les variables numériques.
```{r}
lr_reg_grid <- tibble(penalty = 10^seq(-4, -1, length.out = 30))
```

```{r}
lr_res <- 
  lr_workflow %>% 
  tune_grid(val_set,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
```

```{r}
lr_plot <- 
  lr_res %>% 
  collect_metrics() %>% 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  ylab("Area under the ROC Curve") +
  scale_x_log10(labels = scales::label_number())

lr_plot 
```

Ce graphique montre que les performances du modèle sont généralement meilleures lorsque les valeurs de pénalité sont plus faibles. Cela suggère que la plupart des prédicteurs sont importants pour le modèle. Étant donné que nous avons sélectionné un nombre limité de variables pour ce modèle, cela semble cohérent.

```{r}
select_best(lr_res, metric = 'roc_auc')
```

```{r}
lr_best <- lr_res %>% 
  collect_metrics() %>% 
  arrange(penalty) %>% 
  slice(1)

lr_auc <- lr_res |>
    collect_predictions(parameters = lr_best) |>
    roc_curve(state_binary, .pred_0) |>
    mutate(model = "Logistic Regression")
```

Le niveau de performance généré par ce modèle de régression logistique est correct mais pas optimal. La nature linéaire de l'équation de prédiction peut être trop restrictive pour cet ensemble de données. Lors de la prochaine étape, nous pourrions envisager d'utiliser un modèle non linéaire plus complexe, tel qu'une méthode d'ensemble basée sur les arbres.

##  tree-based ensemble

Comparons notre premier modèle à un random forest.

```{r}
cores <- parallel::detectCores()


rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) %>% 
  set_engine("ranger", num.threads = cores/2) %>% 
  set_mode("classification")
```

```{r}
rf_recipe <- 
  recipe(state_binary ~ ., data = df_other)

rf_workflow <- 
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_recipe(rf_recipe)

set.seed(345)
rf_res <- 
  rf_workflow %>% 
  tune_grid(val_set,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
```

```{r}
rf_best <- 
  rf_res %>% 
  select_best(metric = "roc_auc")
rf_best
```

## Compare both models using ROC curves 

```{r}
rf_auc <- 
  rf_res %>% 
  collect_predictions(parameters = rf_best) %>% 
  roc_curve(state_binary, .pred_0) %>% 
  mutate(model = "Random Forest")

bind_rows(rf_auc, lr_auc) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6)
```

La forêt aléatoire est systématiquement meilleure pour tous les seuils de probabilité des événements.

## Fit the Final Model

```{r}
# the last model
last_rf_mod <- 
  rand_forest(mtry = 1, min_n = 38, trees = 1000) %>% 
  set_engine("ranger", num.threads = cores / 2, importance = "impurity") %>% 
  set_mode("classification")

# the last workflow
last_rf_workflow <- 
  rf_workflow %>% 
  update_model(last_rf_mod)

# the last fit
set.seed(345)
last_rf_fit <- 
  last_rf_workflow %>% 
  last_fit(splits)

last_rf_fit %>% 
  collect_metrics()
```


```{r}
last_rf_fit %>% 
  extract_fit_parsnip() %>% 
  vip(num_features = 10)
```

Nous obtenons donc un modèle avec une précision de 66,5 % et un ROC AUC de 0,697. Bien que ces résultats ne soient pas extraordinaires, ils sont tout de même encourageants pour poursuivre l'exploration et améliorer l'efficacité de notre modèle.

# Améliorations possibles et rapides

- Examiner le mois de lancement (en particulier) pour vérifier la saisonnalité ou d'autres facteurs similaires.
- Jeu de données déséquilibré (mais pas de manière excessive) : possibilité de faire du sous-échantillonnage (probablement la méthode la plus simple ici) ou du sur-échantillonnage de la classe réussie. Dans ce cas là, il est préférable d'utiliser des métriques telles que le F1-score, le rappel (recall), la précision (precision) ou le AUC-ROC, car l'accuracy peut être biaisée dans ce cas.
- Envisager de créer un modèle différent par catégorie principale (main_category) ?
- Prendre en compte le jour du mois où le projet est lancé ?
- Privilégier Ridge dans le modèle glmnet étant donné que nous avons déjà très peu de variables.
- Tester divers modèles, notamment XGboost / LightGBM qui permettront probablement d'augmenter l'éfficacité du modèle


# Session info 

```{r}
sessionInfo()
```