projet4_mlogit

Author

Team 4

Définition des Objectifs :

Section 1

Chargez le jeu de données Yogurt.
Utilisez summarise ou stargazer pour générer le tableau des statistiques descriptives, en vous inspirant de la Table 1 de l’article “Generalized Reverse Discrete Choice Models”

Section 2

Estimez le modèle Logit Multinomial (MNL) sur les données Yogurt avec mlogit.
Estimez aussi le modèle Nested Logit (NL), considérant la corrélation entre alternatives.
Formatez les résultats de vos estimations selon la Table 2 de l’article en utilisant stargazer.

Afin d’arriver à bout de ces objectifs, nous allons part la collecte de données.

Collecte des Données

La source de données provient de Ecdata. Ecdata est un recueil de datasets pour l’economotries et les sciences politiques.

En l’occurence, notre projet portera sur le dataset Yogurt.

Exploration des Données

Nombres d’observation : 2412

Pays : Etats-Unis

Organisation de la donnees :

id : identite des individus
choice : le choix effectué par l’individu pouvant etre un des possibilites suivantes :
- yoplait
- dannon
- weight
- hiland
feat.{z} : nombre d’article de presse sur z ou {z} pouvant être une valeur de la variable choice,
price.{z} : contient le prix de z ou {z} pouvant être une valeur de la variable choice.

Introduction

Le but du projet

Les dependances R du projet :

mlogit
stargazer
tidyverse
Ecdat
ggplot2

library(mlogit)
library(tidyverse)
library(Ecdat)
library(stargazer)
library(ggplot2)
library(kableExtra)

Avant de proceder a une analyse, on a besoin de charge la donnée.

 data("Yogurt")
raw_data = Yogurt
head(raw_data, 10)

   id feat.yoplait feat.dannon feat.hiland feat.weight price.yoplait
1   1            0           0           0           0          10.8
2   1            0           0           0           0          10.8
3   1            0           0           0           0          10.8
4   1            0           0           0           0          10.8
5   1            0           0           0           0          12.5
6   1            0           0           0           0          10.8
7   1            0           0           0           0          10.3
8   1            0           0           0           0          10.8
9   2            0           0           0           0          10.8
10  2            0           0           0           0          10.8
   price.dannon price.hiland price.weight  choice
1           8.1          6.1          7.9  weight
2           9.8          6.4          7.5  dannon
3           9.8          6.1          8.6  dannon
4           9.8          6.1          8.6  dannon
5           9.8          4.9          7.9  dannon
6           9.2          5.0          7.9  dannon
7           8.1          4.9          7.9  dannon
8           8.6          5.4          7.9  weight
9           9.8          5.0          7.9 yoplait
10          9.8          5.0          7.9 yoplait

On constate que la donnée est au format : wide.

Avant de continuer pour une question d’integrité de la donnée, veuillez ajouter un id à chaque choix effectuer.

raw_data <- raw_data %>% mutate(choiceid = row_number())

head(raw_data, 10)

   id feat.yoplait feat.dannon feat.hiland feat.weight price.yoplait
1   1            0           0           0           0          10.8
2   1            0           0           0           0          10.8
3   1            0           0           0           0          10.8
4   1            0           0           0           0          10.8
5   1            0           0           0           0          12.5
6   1            0           0           0           0          10.8
7   1            0           0           0           0          10.3
8   1            0           0           0           0          10.8
9   2            0           0           0           0          10.8
10  2            0           0           0           0          10.8
   price.dannon price.hiland price.weight  choice choiceid
1           8.1          6.1          7.9  weight        1
2           9.8          6.4          7.5  dannon        2
3           9.8          6.1          8.6  dannon        3
4           9.8          6.1          8.6  dannon        4
5           9.8          4.9          7.9  dannon        5
6           9.2          5.0          7.9  dannon        6
7           8.1          4.9          7.9  dannon        7
8           8.6          5.4          7.9  weight        8
9           9.8          5.0          7.9 yoplait        9
10          9.8          5.0          7.9 yoplait       10

Prétraitement des Données :

Analyse de la qualité de la donnée :

Le premier point qui sera analysé est la présence de valeur nulle commodement appélées na dans notre dataset. La présence de ces valeurs non renseignées peut alterer la qualité du modèle qui sera mis en place

library(dplyr)
colonnes <- colnames(raw_data)
head(is.na(raw_data),1)

     id feat.yoplait feat.dannon feat.hiland feat.weight price.yoplait
1 FALSE        FALSE       FALSE       FALSE       FALSE         FALSE
  price.dannon price.hiland price.weight choice choiceid
1        FALSE        FALSE        FALSE  FALSE    FALSE

Après le check a travers la fonction is.na, on s’est rendu compte qu’on n’a pas de valeur na.

Analyse de la présence de donnée abbrérantes :

une donnée aberrante (anglais outlier) est une valeur ou une observation qui est « distante » des autres observations effectuées sur le même phénomène, c’est-à-dire qu’elle contraste grandement avec les valeurs « normalement » mesurées. Une donnée aberrante peut être due à la variabilité inhérente au phénomène observé, ou indiquer une erreur expérimentale.

Ces valeurs abérrantes peut detorier la qualité du modèle,

En cas de valeur abberante, on peut amener à la supprimer ou à la mettre remplacée par le moyenne ou le médian si cela ,’impacte pas la distribution initiale de la donnée.

Pour detecter une valeur abbérante, on a deux solutions :

Soit à travers la data visualisation (utilisation d’un boxplot)
On peut utilser la méthode z-score
En utilsant : le filtre de Hampel Soit en calcul l’interval non abbérante obtenus comme suit [moyen-3*ecartype; moyen +3 ecartype]; toute les valeurs en dehors cet interval sont extenpiées comme abérrantes, cela

Dans cette section, la prémiere méthode sera utilisée sur la variable prix.

df_long <- pivot_longer(Yogurt,
 cols = 2:9, # Exclure les colonnes id et choice
 names_to = c(".value", "Brand"), # .value pour les valeurs, alternative pour les noms de marque
 names_sep = "\\.", # Séparer les noms de colonnes sur le point littéral
)

#df_long%>%mutate(choice=ifelse(choice ==Brand,1,0))

#donnees <- data.frame(
#  categorie = "Price",
#  valeurs = df_long$price
#)

#ggplot(donnees, aes(x = categorie, y = valeurs, fill = categorie)) +
#  geom_boxplot() +
#  labs(title = "Boxplot des données par catégorie", x = "Catégorie", y = "Valeurs") +
#  theme_minimal() +
#  scale_fill_brewer(palette = "Set3")

On remarque un certain nombre de valeur exterme cependant, chose qui nous semble normal du au caractère niche des produit danont.

La statistique description :

Il sera question dans cette section de reproduire la table 1.

mean_table = aggregate(df_long[, 4], list(df_long$Brand), mean)
st_table = aggregate(df_long[, 4], list(df_long$Brand), sd)
colnames(mean_table)[colnames(mean_table) == 'feat'] <- 'mean'
colnames(st_table)[colnames(st_table) == 'feat'] <- 'std'
colnames(mean_table)[colnames(mean_table) == 'Group.1'] <- 'Brand'
colnames(st_table)[colnames(st_table) == 'Group.1'] <- 'Brand'

#stargazer(df_long, type="text")

#head(mean_table, 4)

kable(mean_table, format = "pipe")

Brand	mean
dannon	0.0377280
hiland	0.0368988
weight	0.0377280
yoplait	0.0559701

#head(st_table, 4)

kable(st_table, format = "pipe")

Brand	std
dannon	0.1905772
hiland	0.1885525
weight	0.1905772
yoplait	0.2299117

final = merge(mean_table,st_table)

df2 <- data.frame(t(final[-1]))
colnames(df2) <- st_table[, 1]
kable(df2, format = "pipe")

	dannon	hiland	weight	yoplait
mean	0.0377280	0.0368988	0.0377280	0.0559701
std	0.1905772	0.1885525	0.1905772	0.2299117

La data visualisation :

La data visualization permet de repérer des modèles, des tendances et des corrélations qui, autrement, passeraient inaperçus dans les rapports, tableaux ou feuilles de calcul traditionnels.undefined

La premiere visualisation qu’on vous propose est le diagramme circulaire(Pie Chart). Elle est souvent utilisée pour représenter graphiquement des données sous forme de proportions ou de pourcentages.

library("ggplot2")
library("ggthemes")
data_pie_chart <- data.frame(
  choice=raw_data$choice,
  choice_possibles=c("dannon", "yoplait", "weight", "hiland")
)
comptage <- table(raw_data$choice)

df_comptage <- as.data.frame(comptage)
names(df_comptage) <- c("marque", "nombre")
df_comptage$pourcentage <- df_comptage$nombre / sum(df_comptage$nombre) * 100


ggplot(df_comptage, aes(x="", y = pourcentage, fill = marque)) +
  geom_bar(stat = "identity", width = 1) +
  geom_text(aes(label = paste0(round(pourcentage, 1), "%")), position = position_stack(vjust = 0.5))+
  coord_polar("y") +
  labs(title = "Part de marchés des marques") +
  theme_minimal() +
  theme_economist() + 
  scale_fill_brewer(palette = "Set3") +
  theme(legend.position = "right",  # Position de la légende
        axis.title.x = element_blank(),  # Supprimer la légende de l'axe x
        axis.title.y = element_blank())

#appliquer au format large (wide) [*]
#changer choices => marque
#effacer x en ordonnee
#effacer count => nombre
#changer les effectifs en %centage

La seconde visusalisation est un diagramme circulaire qui prends en compte dans le diametre des cercles, la proportion de chaque modalité de la variables.

library(ggplot2)

# Exemple de données factices
donnees <- data.frame(
  marque = raw_data$choice)

# Tracer le graphique avec geom_bar et stat=count
ggplot(donnees, aes(x = marque)) +
  geom_bar(color = "black", aes(y = ..count..)) +
  coord_polar("x") +
  theme_economist() +
  scale_fill_brewer(palette = "Set3")

histogram : pour les valeurs continue // par marque

Interpretation a mettre après :

Interpretation des frequences et pourquoi c’est important de faire l’histogramme

La troisième visualisation est l’historigramme des prix. un histogramme est un graphique permettant de représenter la répartition des valeurs d’une variable continue. Chaque colonne de l’ histogramme représente un intervalle de valeurs. La hauteur des colonnes indique le nombre d’instances dans cet intervalle. L’examen de l’histogramme permet de se faire une idée claire sur la distribution des valeurs de la feature analysée.

donnees <- df_long$price

# Créer un dataframe avec les données
df <- data.frame(valeurs = donnees)

# Tracer l'histogramme avec ggplot2
ggplot(df, aes(x = valeurs)) +
  geom_histogram(binwidth = 0.5, fill = "skyblue", color = "black", aes(y = after_stat(count))) + theme_economist() + 
  labs(title = "Histogramme des prix", x = "Valeurs", y = "Fréquence")

La construction du modèle:

Afin de pourvoir, entrainer notre modèle, il faut que la donnée soit pivotée. D’un format wide à un format Long.

Cette prépartion de la donnée sera faite grace à mlogit.data pour la mettre dans un format convenanble pour mlogit.

#data_ready = mlogit.data(Yogurt, shape="wide", varying = c(2:8), choice= "choice" )
data_ready = mlogit.data(Yogurt, shape="wide", varying = c("feat.yoplait", "feat.dannon", "feat.hiland", "feat.weight", "price.yoplait", "price.dannon", "price.hiland", "price.weight"), choice= "choice" )

head(data_ready,4)

~~~~~~~
 first 4 observations out of 9648 
~~~~~~~
  id choice     alt feat price chid    idx
1  1  FALSE  dannon    0   8.1    1 1:nnon
2  1  FALSE  hiland    0   6.1    1 1:land
3  1   TRUE  weight    0   7.9    1 1:ight
4  1  FALSE yoplait    0  10.8    1 1:lait

~~~ indexes ~~~~
  chid     alt
1    1  dannon
2    1  hiland
3    1  weight
4    1 yoplait
indexes:  1, 2

Estimons le modèle Logit Multinomial (MNL) sur les données Yogurt avec mlogit.

# a "mixed" model
m <- mlogit(choice ~ price + feat , data = data_ready, reflevel = "hiland")
#kable(summary(m),format = "pipe") 
summary(m)


Call:
mlogit(formula = choice ~ price + feat, data = data_ready, reflevel = "hiland", 
    method = "nr")

Frequencies of alternatives:choice
  hiland   dannon   weight  yoplait 
0.029436 0.402156 0.229270 0.339138 

nr method
6 iterations, 0h:0m:1s 
g'(-H)^-1g = 0.000763 
successive function values within tolerance limits 

Coefficients :
                     Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept):dannon   3.715595   0.145419  25.5510 < 2.2e-16 ***
(Intercept):weight   3.074411   0.145384  21.1468 < 2.2e-16 ***
(Intercept):yoplait  4.450166   0.187118  23.7827 < 2.2e-16 ***
price               -0.366584   0.024366 -15.0449 < 2.2e-16 ***
feat                 0.491433   0.120063   4.0931 4.256e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-Likelihood: -2656.9
McFadden R^2:  0.062142 
Likelihood ratio test : chisq = 352.09 (p.value = < 2.22e-16)

Estimez aussi le modèle Nested Logit (NL), considérant la corrélation entre alternatives.

mixed_logit_model <- mlogit(choice ~ price + feat , data = data_ready,
                            rpar = c(price = "n", feat= "n"),
                            R = 100,
                            halton = NA)

summary(mixed_logit_model)


Call:
mlogit(formula = choice ~ price + feat, data = data_ready, rpar = c(price = "n", 
    feat = "n"), R = 100, halton = NA)

Frequencies of alternatives:choice
  dannon   hiland   weight  yoplait 
0.402156 0.029436 0.229270 0.339138 

bfgs method
12 iterations, 0h:0m:22s 
g'(-H)^-1g = 2.91E-07 
gradient close to zero 

Coefficients :
                      Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept):hiland  -3.9877853  0.1660274 -24.0188 < 2.2e-16 ***
(Intercept):weight  -0.6639164  0.0549685 -12.0781 < 2.2e-16 ***
(Intercept):yoplait  0.7942609  0.0897850   8.8463 < 2.2e-16 ***
price               -0.3928668  0.0273948 -14.3409 < 2.2e-16 ***
feat                 0.3598263  0.2075419   1.7338   0.08296 .  
sd.price             0.0040998  0.6311035   0.0065   0.99482    
sd.feat              2.3517943  0.5693641   4.1306 3.619e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-Likelihood: -2646
McFadden R^2:  0.065983 
Likelihood ratio test : chisq = 373.85 (p.value = < 2.22e-16)

random coefficients
      Min.   1st Qu.     Median       Mean    3rd Qu. Max.
price -Inf -0.395632 -0.3928668 -0.3928668 -0.3901015  Inf
feat  -Inf -1.226435  0.3598263  0.3598263  1.9460875  Inf

References

Kable Documentation

Ecdata