Wine_Purchase_Prediction
Pin Lyu
2023-11-21
Project Description
A large wine manufacturer is studying the data in order to predict the number of wine cases ordered based upon the wine characteristics. If the wine manufacturer can predict the number of cases, then that manufacturer will be able to adjust their wine offering to maximize sales. the objective is to build a count regression model to predict the number of cases of wine that will be sold given certain properties of the wine
1. Set Data Files/Directory
# Set data direction
wine <- read.csv("/Users/pin.lyu/Desktop/BC_Class_Folder/Econometrics/DIS_&_ASSIGNMENT/HW#3_Wine/wine_data.csv")
# Split data
set.seed(678)
## Generate a vector of indices for the training set
train_index <- sample(x = nrow(wine), size = round(0.8 * nrow(wine) ) )
## Create the training and testing sets
train <- wine[train_index, ]
test <- wine[-train_index, ]2. Data Exploration
Rename Variables
# Make variable names more readable
colnames(train) <- c("Index",
"Cases_Purchased",
"Fixed_Acidity",
"Volatile_Acidity",
"Citric_Acid",
"Residual_Sugar",
"Chloride",
"Free_Sulfur_Dioxide",
"Total_Sulfur_Dioxide",
"Density",
"pH_Level",
"Sulphate",
"Alcohol_Level",
"Label_Appeal",
"Acidity_Index",
"Quality_Rating"
)
# Delete variable "Index"
train <- train |>
dplyr::select(- Index)# Summary stats table
stargazer(train,
type = "text",
title = "Table 1: Basic Data Summary Statistics"
)##
## Table 1: Basic Data Summary Statistics
## ===============================================================
## Statistic N Mean St. Dev. Min Max
## ---------------------------------------------------------------
## Cases_Purchased 10,236 3.031 1.928 0 8
## Fixed_Acidity 10,236 7.071 6.288 -18.000 34.400
## Volatile_Acidity 10,236 0.320 0.788 -2.790 3.590
## Citric_Acid 10,236 0.305 0.855 -3.240 3.860
## Residual_Sugar 9,729 5.714 33.823 -127.800 141.150
## Chloride 9,731 0.056 0.318 -1.171 1.351
## Free_Sulfur_Dioxide 9,717 30.695 148.768 -555.000 623.000
## Total_Sulfur_Dioxide 9,688 121.071 230.401 -823.000 1,054.000
## Density 10,236 0.994 0.026 0.888 1.099
## pH_Level 9,915 3.210 0.680 0.540 6.130
## Sulphate 9,273 0.532 0.927 -3.130 4.240
## Alcohol_Level 9,714 10.516 3.742 -4.700 26.500
## Label_Appeal 10,236 -0.005 0.895 -2 2
## Acidity_Index 10,236 7.777 1.325 4 17
## Quality_Rating 7,553 2.038 0.904 1 4
## ---------------------------------------------------------------Comments: From the summary statistic table above, we can see that we have multiple variable showing us the same thing, such as acidity of a wine. Other useful information we can extract from this table include the mean pH level of wine in this data set is 3.2 and the mean alcohol level in wine data is around 10%.
Also there are variables that can tell us about how industry experts feel about these wines. Each wine is ranked in two ways, their label appeals, in a variable called “LabelAppeal”, and their overall quality captured in “STARS”. Label Appeal has five different ranks from -2 to 2. And overall quality is ranked from 1 to 4, a total of 4 levels.
Not all variables have the same size which means many of them have missing values.
Lastly, it’s worth to point out that some of the chemical levels are negative which can cause problems in our regressions. Therefore, later I will remove these outliers.
Quality Rating vs. pH Level
train |>
group_by(Quality_Rating,pH_Level,Alcohol_Level) |>
summarise(n = n()) |>
ggplot(aes(x = pH_Level, y = Quality_Rating, size = n)) +
geom_point() +
geom_point(color = "black", alpha = 0.5)+
labs(title = "Wine Rating vs. pH Level",
y = "Rating of Wine",
x = "pH Level") +
theme_classic() ## `summarise()` has grouped output by 'Quality_Rating', 'pH_Level'. You can
## override using the `.groups` argument.## Warning: Removed 2605 rows containing missing values (`geom_point()`).
## Removed 2605 rows containing missing values (`geom_point()`).
Comments: Most of rated wine have a pH level more or less around 2 to 4. However, the range of pH level of 4-stars wines decreased compare to other ranking levels. Additionally, among the highest ranking of wines, many wine in this ranking tend to have a pH level around 3.
Alcohol Level vs. pH Level
train |>
group_by(Quality_Rating,pH_Level,Alcohol_Level) |>
summarise(n = n()) |>
ggplot(aes(x = Alcohol_Level, y = pH_Level, size = n)) +
geom_point(position = "dodge") +
geom_point(color = "aquamarine3", alpha = 0.7) +
labs(title = "pH Level vs. Alcohol Level",
y = "pH Level",
x = "Alcohol Level") +
theme_classic() ## `summarise()` has grouped output by 'Quality_Rating', 'pH_Level'. You can
## override using the `.groups` argument.## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`## Warning: Removed 609 rows containing missing values (`geom_point()`).
## Removed 609 rows containing missing values (`geom_point()`).
Comments: No clear relationship between pH level and alcohol level. Based on the graph, we can tell that most of the wine has a alcohol level centered around 10%. This matches with what we found in the summary statistics table above, it shows us that the mean of alcohol level is 10.489.
Cases Purchased vs. Acidity Index
train |>
group_by(Cases_Purchased, Acidity_Index) |>
summarise(n = n()) |>
ggplot(aes(x = Acidity_Index, y = Cases_Purchased, size = n)) +
geom_point(position = "dodge") +
geom_point(color = "aquamarine3", alpha = 0.7)+
labs(title = "Number of Purchases vs. Acidity Index",
y = "Number of Cases Purchased",
x = "Acid Index") +
theme_classic() ## `summarise()` has grouped output by 'Cases_Purchased'. You can override using
## the `.groups` argument.## Warning: Width not defined
## ℹ Set with `position_dodge(width = ...)`
# Acidity Index range
table(train$Acidity_Index)##
## 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## 3 61 933 3910 3327 1145 442 203 104 50 40 7 5 6# Number of purhcase range
table(train$Cases_Purhcased)## < table of extent 0 >Comments: Based on the tables and the graph of Acid index vs. cases purchased above, we can see that consumers tend to favor wine that are not too acidic (4-10). Most of the purchased wine have an acid index level around 6 to 9. Additionally, majority of the customers tend to purchase wine of this quality, anywhere from, 2 to 6 cases together, in one single setting.
train |>
group_by(Cases_Purchased, Alcohol_Level) |>
summarise(n = n()) |>
ggplot(aes(x = Alcohol_Level, y = Cases_Purchased, size = n)) +
geom_point() +
geom_point(color = "black", alpha = 0.5)+
labs(title = "Number of purhcase by alcohol Level",
y = "Number of cases purhcased",
x = "Alcohol Level") +
theme_classic() ## `summarise()` has grouped output by 'Cases_Purchased'. You can override using
## the `.groups` argument.## Warning: Removed 8 rows containing missing values (`geom_point()`).
## Removed 8 rows containing missing values (`geom_point()`).
Comments: Based on this graph, we can infer some information about wine. What I noticed is that the alcohol level of wine in this data set tend to range anywhere from 2% to 25% but most of wine have a alcohol level around 10% (+/- 5%).
3. Data Preparation
Imputing N/As
# Missing data visulization
missmap(train)
# Impute numeric varibles that have N/As to median
for (col in names(train)) {
# Check if the column contains any N/A values
if (any(is.na(train[[col]])) & col != "Quality_Rating") {
# Calculate the median of the non-missing values in the column
median_value <- mean(train[[col]], na.rm = TRUE)
# Replace N/A values with the median value
train[[col]][is.na(train[[col]])] <- median_value
}
}
# Visualize data set
missmap(train)
Comments: After this procedure, all variables, except “Quality_Rating”, are N/As free. The reason behind not imputing the missing values of “Quality_Rating” to its median is because almost 1/3 of its data are N/As. To impute this large amount of missing values to other average values, median in our case, would cause the data would decrease the truthfulness of the data, hence leading to potential biases.
# subset all rows of Quality_Rating has N/A
train_copy <- train[complete.cases(train$Quality_Rating), ]
# Visulize data set
stargazer(train_copy,
type = "text",
title = "Table 2: Data Summary Statistics (Updated)"
)##
## Table 2: Data Summary Statistics (Updated)
## ==============================================================
## Statistic N Mean St. Dev. Min Max
## --------------------------------------------------------------
## Cases_Purchased 7,553 3.685 1.556 0 8
## Fixed_Acidity 7,553 6.932 6.275 -18.000 32.500
## Volatile_Acidity 7,553 0.291 0.784 -2.750 3.550
## Citric_Acid 7,553 0.310 0.856 -3.160 3.770
## Residual_Sugar 7,553 5.834 32.795 -127.800 141.150
## Chloride 7,553 0.050 0.308 -1.171 1.270
## Free_Sulfur_Dioxide 7,553 33.559 145.835 -555.000 622.000
## Total_Sulfur_Dioxide 7,553 125.648 222.903 -793.000 1,048.000
## Density 7,553 0.994 0.027 0.888 1.099
## pH_Level 7,553 3.205 0.666 0.540 6.050
## Sulphate 7,553 0.517 0.874 -3.130 4.110
## Alcohol_Level 7,553 10.578 3.665 -4.500 26.500
## Label_Appeal 7,553 0.049 0.880 -2 2
## Acidity_Index 7,553 7.643 1.174 4 17
## Quality_Rating 7,553 2.038 0.904 1 4
## --------------------------------------------------------------Comments: After this procedure, a total of 2683 rows were deleted. Now, “Quality_Rating” is also N/A free. Next, I will remove variables that have negative values as any chemical level should only have a minimal level of 0.
# Turn negative values to positive
train_copy$Fixed_Acidity <- abs(train_copy$Fixed_Acidity)
train_copy$Volatile_Acidity <- abs(train_copy$Volatile_Acidity)
train_copy$Citric_Acid <- abs(train_copy$Citric_Acid)
train_copy$Residual_Sugar <- abs(train_copy$Residual_Sugar)
train_copy$Chloride <- abs(train_copy$Chloride)
train_copy$Free_Sulfur_Dioxide <- abs(train_copy$Free_Sulfur_Dioxide )
train_copy$Total_Sulfur_Dioxide <- abs(train_copy$Total_Sulfur_Dioxide)
train_copy$Sulphate <- abs(train_copy$Sulphate)
train_copy$Alcohol_Level <- abs(train_copy$Alcohol_Level)
# Rename data
wine_data <- train_copyComments: Chemical values should have negative values. Here, I will be assuming that the negative sign is a typo. Therefore, I transformed all negative values in to their absolute terms.
- Resources Used For Industry Information:
The next few steps of data cleaning process utilized information from the websites above.
# Delete data outside of industry standard
wine_data <- filter(wine_data, Volatile_Acidity < 1.5 &
Citric_Acid < 0.9 )stargazer(wine_data,
type = "text",
title = "Table 3: Data Summary Statistics (Updated)")##
## Table 3: Data Summary Statistics (Updated)
## ===========================================================
## Statistic N Mean St. Dev. Min Max
## -----------------------------------------------------------
## Cases_Purchased 5,002 3.699 1.541 0 8
## Fixed_Acidity 5,002 7.974 4.940 0.000 32.500
## Volatile_Acidity 5,002 0.505 0.363 0.000 1.490
## Citric_Acid 5,002 0.377 0.200 0.000 0.890
## Residual_Sugar 5,002 22.721 24.833 0.000 141.150
## Chloride 5,002 0.216 0.230 0.000 1.270
## Free_Sulfur_Dioxide 5,002 104.872 107.631 0.000 618.000
## Total_Sulfur_Dioxide 5,002 202.005 157.817 0.000 1,048.000
## Density 5,002 0.994 0.027 0.890 1.099
## pH_Level 5,002 3.195 0.677 0.540 6.050
## Sulphate 5,002 0.800 0.621 0.000 4.110
## Alcohol_Level 5,002 10.619 3.525 0.100 26.500
## Label_Appeal 5,002 0.048 0.880 -2 2
## Acidity_Index 5,002 7.648 1.198 4 17
## Quality_Rating 5,002 2.052 0.905 1 4
## -----------------------------------------------------------# Boxplot for all numeric variables
boxplot_list <- lapply(names(wine_data), function(var) {
ggplot(data = wine_data, aes(x = .data[[var]])) +
geom_boxplot() +
labs(title = paste("Boxplot of", var)) +
theme_minimal()
})
# Boxplot arrangement settings
grid.arrange(grobs = boxplot_list, ncol = 3)
Comments: After this procedure, only 1066 observations remained in the data set.
Correlation Plot
# Correlation chart
corr <- cor(wine_data)
# Present the chart in a graph
corrplot(corr, method = 'square', order = 'FPC', type = 'lower', diag = FALSE)
Comments: Much of the variables in this data set do not correlate with any other variables. However, from this correlation plot, we can extract some very useful information. First, we see that variable “Label_Appeal” and “Quality_Rating” have a significant positive correlation with the total number of cases of wine purchased. This might suggest that lots of consumers’ purchases are influenced by the presentation of the wine’s package label and the rating of the wine during their purchases. “Quality_Rating” and “Label_Appeal” have a strong positive correlation. Additionally, the variable “Fixed_Acidity” has a slight positive correlation with the acidity index. This can suggest that the fixed acidity in a wine is what largely determines the acidity index that is assigned to a wine.
4. Models
Multi-linear Logistic Regression Model A
# A full multi-linear regression model
multl_A <- lm(Cases_Purchased ~ . , data = wine_data)
# Use both "forward" & "backward" selection
linear_modelA <- step(multl_A, direction = "both")## Start: AIC=1307.05
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Residual_Sugar + Chloride + Free_Sulfur_Dioxide + Total_Sulfur_Dioxide +
## Density + pH_Level + Sulphate + Alcohol_Level + Label_Appeal +
## Acidity_Index + Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Sulphate 1 0.01 6456.9 1305.0
## - Residual_Sugar 1 0.18 6457.1 1305.2
## - Free_Sulfur_Dioxide 1 0.28 6457.1 1305.3
## - Chloride 1 0.52 6457.4 1305.5
## - Total_Sulfur_Dioxide 1 1.36 6458.2 1306.1
## - Density 1 2.00 6458.9 1306.6
## - Fixed_Acidity 1 2.41 6459.3 1306.9
## <none> 6456.9 1307.0
## - pH_Level 1 4.53 6461.4 1308.5
## - Volatile_Acidity 1 5.87 6462.7 1309.6
## - Citric_Acid 1 11.24 6468.1 1313.8
## - Alcohol_Level 1 30.70 6487.6 1328.8
## - Acidity_Index 1 175.14 6632.0 1438.9
## - Label_Appeal 1 1473.97 7930.8 2333.5
## - Quality_Rating 1 1815.54 8272.4 2544.4
##
## Step: AIC=1305.05
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Residual_Sugar + Chloride + Free_Sulfur_Dioxide + Total_Sulfur_Dioxide +
## Density + pH_Level + Alcohol_Level + Label_Appeal + Acidity_Index +
## Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Residual_Sugar 1 0.18 6457.1 1303.2
## - Free_Sulfur_Dioxide 1 0.28 6457.2 1303.3
## - Chloride 1 0.52 6457.4 1303.5
## - Total_Sulfur_Dioxide 1 1.36 6458.2 1304.1
## - Density 1 2.01 6458.9 1304.6
## - Fixed_Acidity 1 2.41 6459.3 1304.9
## <none> 6456.9 1305.0
## - pH_Level 1 4.52 6461.4 1306.5
## + Sulphate 1 0.01 6456.9 1307.0
## - Volatile_Acidity 1 5.88 6462.8 1307.6
## - Citric_Acid 1 11.24 6468.1 1311.8
## - Alcohol_Level 1 30.70 6487.6 1326.8
## - Acidity_Index 1 175.21 6632.1 1437.0
## - Label_Appeal 1 1473.98 7930.9 2331.5
## - Quality_Rating 1 1815.71 8272.6 2542.6
##
## Step: AIC=1303.19
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Chloride + Free_Sulfur_Dioxide + Total_Sulfur_Dioxide + Density +
## pH_Level + Alcohol_Level + Label_Appeal + Acidity_Index +
## Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Free_Sulfur_Dioxide 1 0.28 6457.3 1301.4
## - Chloride 1 0.52 6457.6 1301.6
## - Total_Sulfur_Dioxide 1 1.35 6458.4 1302.2
## - Density 1 2.00 6459.1 1302.7
## - Fixed_Acidity 1 2.42 6459.5 1303.1
## <none> 6457.1 1303.2
## - pH_Level 1 4.50 6461.6 1304.7
## + Residual_Sugar 1 0.18 6456.9 1305.0
## + Sulphate 1 0.01 6457.1 1305.2
## - Volatile_Acidity 1 5.85 6462.9 1305.7
## - Citric_Acid 1 11.15 6468.2 1309.8
## - Alcohol_Level 1 30.70 6487.8 1324.9
## - Acidity_Index 1 175.10 6632.2 1435.0
## - Label_Appeal 1 1476.36 7933.4 2331.2
## - Quality_Rating 1 1815.58 8272.6 2540.6
##
## Step: AIC=1301.41
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Chloride + Total_Sulfur_Dioxide + Density + pH_Level + Alcohol_Level +
## Label_Appeal + Acidity_Index + Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Chloride 1 0.51 6457.8 1299.8
## - Total_Sulfur_Dioxide 1 1.35 6458.7 1300.5
## - Density 1 2.06 6459.4 1301.0
## - Fixed_Acidity 1 2.45 6459.8 1301.3
## <none> 6457.3 1301.4
## - pH_Level 1 4.55 6461.9 1302.9
## + Free_Sulfur_Dioxide 1 0.28 6457.1 1303.2
## + Residual_Sugar 1 0.18 6457.2 1303.3
## + Sulphate 1 0.01 6457.3 1303.4
## - Volatile_Acidity 1 5.81 6463.1 1303.9
## - Citric_Acid 1 11.08 6468.4 1308.0
## - Alcohol_Level 1 30.85 6488.2 1323.2
## - Acidity_Index 1 174.89 6632.2 1433.1
## - Label_Appeal 1 1476.48 7933.8 2329.4
## - Quality_Rating 1 1816.13 8273.5 2539.1
##
## Step: AIC=1299.8
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Total_Sulfur_Dioxide + Density + pH_Level + Alcohol_Level +
## Label_Appeal + Acidity_Index + Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Total_Sulfur_Dioxide 1 1.34 6459.2 1298.8
## - Density 1 2.04 6459.9 1299.4
## - Fixed_Acidity 1 2.49 6460.3 1299.7
## <none> 6457.8 1299.8
## - pH_Level 1 4.55 6462.4 1301.3
## + Chloride 1 0.51 6457.3 1301.4
## + Free_Sulfur_Dioxide 1 0.27 6457.6 1301.6
## + Residual_Sugar 1 0.18 6457.7 1301.7
## + Sulphate 1 0.01 6457.8 1301.8
## - Volatile_Acidity 1 5.82 6463.7 1302.3
## - Citric_Acid 1 11.17 6469.0 1306.5
## - Alcohol_Level 1 30.75 6488.6 1321.6
## - Acidity_Index 1 174.53 6632.4 1431.2
## - Label_Appeal 1 1476.54 7934.4 2327.8
## - Quality_Rating 1 1818.25 8276.1 2538.7
##
## Step: AIC=1298.84
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## Density + pH_Level + Alcohol_Level + Label_Appeal + Acidity_Index +
## Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Density 1 1.97 6461.1 1298.4
## - Fixed_Acidity 1 2.46 6461.6 1298.7
## <none> 6459.2 1298.8
## + Total_Sulfur_Dioxide 1 1.34 6457.8 1299.8
## - pH_Level 1 4.56 6463.7 1300.4
## + Chloride 1 0.50 6458.7 1300.5
## + Free_Sulfur_Dioxide 1 0.26 6458.9 1300.6
## + Residual_Sugar 1 0.17 6459.0 1300.7
## + Sulphate 1 0.01 6459.2 1300.8
## - Volatile_Acidity 1 5.97 6465.2 1301.5
## - Citric_Acid 1 11.27 6470.5 1305.6
## - Alcohol_Level 1 30.63 6489.8 1320.5
## - Acidity_Index 1 176.08 6635.3 1431.4
## - Label_Appeal 1 1475.21 7934.4 2325.8
## - Quality_Rating 1 1820.12 8279.3 2538.6
##
## Step: AIC=1298.36
## Cases_Purchased ~ Fixed_Acidity + Volatile_Acidity + Citric_Acid +
## pH_Level + Alcohol_Level + Label_Appeal + Acidity_Index +
## Quality_Rating
##
## Df Sum of Sq RSS AIC
## - Fixed_Acidity 1 2.43 6463.6 1298.2
## <none> 6461.1 1298.4
## + Density 1 1.97 6459.2 1298.8
## + Total_Sulfur_Dioxide 1 1.26 6459.9 1299.4
## - pH_Level 1 4.52 6465.7 1299.9
## + Chloride 1 0.48 6460.7 1300.0
## + Free_Sulfur_Dioxide 1 0.32 6460.8 1300.1
## + Residual_Sugar 1 0.16 6461.0 1300.2
## + Sulphate 1 0.01 6461.1 1300.3
## - Volatile_Acidity 1 5.89 6467.0 1300.9
## - Citric_Acid 1 11.20 6472.3 1305.0
## - Alcohol_Level 1 30.66 6491.8 1320.0
## - Acidity_Index 1 178.22 6639.4 1432.5
## - Label_Appeal 1 1477.29 7938.4 2326.3
## - Quality_Rating 1 1821.48 8282.6 2538.6
##
## Step: AIC=1298.24
## Cases_Purchased ~ Volatile_Acidity + Citric_Acid + pH_Level +
## Alcohol_Level + Label_Appeal + Acidity_Index + Quality_Rating
##
## Df Sum of Sq RSS AIC
## <none> 6463.6 1298.2
## + Fixed_Acidity 1 2.43 6461.1 1298.4
## + Density 1 1.93 6461.6 1298.7
## + Total_Sulfur_Dioxide 1 1.23 6462.3 1299.3
## - pH_Level 1 4.56 6468.1 1299.8
## + Chloride 1 0.52 6463.1 1299.8
## + Free_Sulfur_Dioxide 1 0.35 6463.2 1300.0
## + Residual_Sugar 1 0.18 6463.4 1300.1
## + Sulphate 1 0.01 6463.6 1300.2
## - Volatile_Acidity 1 6.00 6469.6 1300.9
## - Citric_Acid 1 11.10 6474.7 1304.8
## - Alcohol_Level 1 30.69 6494.3 1319.9
## - Acidity_Index 1 189.24 6652.8 1440.6
## - Label_Appeal 1 1477.19 7940.8 2325.8
## - Quality_Rating 1 1821.18 8284.8 2537.9# Summary
stargazer(linear_modelA,
type = "text",
title = "Table 3: Logistic Poisson Model A",
dep.var.labels = "Cases of Wine Purchased",
covariate.labels = c("Citric Acid", "pH Level", "Alcohol Level", "label Appeal", "Acidity Index", "Quality Rating")
)##
## Table 3: Logistic Poisson Model A
## ===============================================
## Dependent variable:
## ---------------------------
## Cases of Wine Purchased
## -----------------------------------------------
## Citric Acid -0.096**
## (0.044)
##
## pH Level 0.238***
## (0.081)
##
## Alcohol Level 0.045*
## (0.024)
##
## label Appeal 0.022***
## (0.005)
##
## Acidity Index 0.659***
## (0.020)
##
## Quality Rating -0.165***
## (0.014)
##
## Quality_Rating 0.714***
## (0.019)
##
## Constant 3.045***
## (0.154)
##
## -----------------------------------------------
## Observations 5,002
## R2 0.456
## Adjusted R2 0.455
## Residual Std. Error 1.138 (df = 4994)
## F Statistic 596.939*** (df = 7; 4994)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Comments: Using stepwise’s feature engineering method, it dropped most of the detailed records of each wine’s acidity variables. There are only 5 variables; they are “pH_Level”, “Alcohol_Level”, “Label_Appeal”, “Acidity_Index”, and “Quality_Rating”. they all are statistically significant. The AIC value of the regression made up by these five variables is 353.66 which is the lowest among the other combinations of variables used in linear regression.
In the data preparation process, I generated two versions of data because of large amounts of missing values in the variable “Quality_Rating”. One data set, called “train”, left the missing values untouched. In the other, “train_copy”, I eliminated all rows in the data set where “Quality_Rating” was missing. Potentially, the difference of its values in this variable can lead to different results in coefficients level as well as regressors chosen in a linear model. So, next I will use the same feature engineering method on the untouched data set “train”, to see if there’s any changes to our results.
Coefficient Interpretation:
Each rank increase in label appeal, consumers are more likely to purchase, roughly, 0.62 more cases of that wine.
Each unit increase in acidity index of a wine, consumers would respond to that by decrease 0.132 less case of that wine.
Multi-linear Logistic Regression Model B
# A full multi-linear regression model
linear_modelB <- lm(Cases_Purchased ~ Label_Appeal +
Quality_Rating +
Acidity_Index +
Alcohol_Level +
pH_Level,
data = wine_data)
# Summary
stargazer(linear_modelB,
type = "text",
title = "Table 3: Logistic Poisson Model B",
dep.var.labels = "Cases of Wine Purchased",
covariate.labels = c("label Appeal", "Citric Acid", "Quality Rating", "Acidity Index", "Bound Sulfur Dioxide", "Alcohol Level", "pH Level", "Intercept")
)##
## Table 3: Logistic Poisson Model B
## ================================================
## Dependent variable:
## ---------------------------
## Cases of Wine Purchased
## ------------------------------------------------
## label Appeal 0.660***
## (0.020)
##
## Citric Acid 0.715***
## (0.019)
##
## Quality Rating -0.162***
## (0.014)
##
## Acidity Index 0.022***
## (0.005)
##
## Bound Sulfur Dioxide 0.041*
## (0.024)
##
## Alcohol Level 3.072***
## (0.152)
##
## ------------------------------------------------
## Observations 5,002
## R2 0.454
## Adjusted R2 0.454
## Residual Std. Error 1.139 (df = 4996)
## F Statistic 831.087*** (df = 5; 4996)
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Comments: according to my research, at a moderate amount, sulfur helps to preserve the fruitiness of the wine as well as its freshness. However, excessive usage of sulfur will caused the taste of wine to change. Therefore, in this model, my intention is to test whether this theoretically effect holds to true or not in our data.
The coefficient value of bound sulfur dioxide is 0.0001 which suggest that each unit increase in the sulfur content will cause the consumer to purchase a little more. However, this number is not statistically significant from zero. Therefore, we did not see the theoretical effect of sulfur on consumers’ purchases.
poisson_A <- glm(data = wine_data ,
formula = Cases_Purchased ~ . ,
family = poisson(link = 'log'))
summary(poisson_A)##
## Call:
## glm(formula = Cases_Purchased ~ ., family = poisson(link = "log"),
## data = wine_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1839 -0.2705 0.0549 0.3656 1.6220
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.352e+00 2.819e-01 4.795 1.62e-06 ***
## Fixed_Acidity -1.340e-03 1.516e-03 -0.884 0.37681
## Volatile_Acidity -2.297e-02 2.035e-02 -1.129 0.25904
## Citric_Acid 6.571e-02 3.709e-02 1.772 0.07644 .
## Residual_Sugar -1.709e-05 2.974e-04 -0.057 0.95416
## Chloride 1.591e-02 3.173e-02 0.501 0.61612
## Free_Sulfur_Dioxide -1.450e-05 6.864e-05 -0.211 0.83264
## Total_Sulfur_Dioxide 3.386e-05 4.645e-05 0.729 0.46608
## Density -2.156e-01 2.766e-01 -0.780 0.43560
## pH_Level 1.090e-02 1.093e-02 0.997 0.31856
## Sulphate -2.296e-04 1.196e-02 -0.019 0.98468
## Alcohol_Level 5.591e-03 2.098e-03 2.665 0.00769 **
## Label_Appeal 1.804e-01 9.001e-03 20.045 < 2e-16 ***
## Acidity_Index -4.635e-02 6.640e-03 -6.981 2.94e-12 ***
## Quality_Rating 1.820e-01 8.418e-03 21.615 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4441.0 on 5001 degrees of freedom
## Residual deviance: 3000.8 on 4987 degrees of freedom
## AIC: 17978
##
## Number of Fisher Scoring iterations: 5countreg::rootogram(poisson_A,
xlab = "Total Number of Cases of Wine Purchased",
main = "Fit of the Poisson Logistic Model A")
Comments: The red curved line is the theoretical Poisson fit. “Hanging” from each point on the red line is a bar, the height of which represents the difference between expected and observed counts. A bar hanging below 0 indicates underfitting. A bar hanging above 0 indicates overfitting. The counts have been transformed with a square root transformation to prevent smaller counts from getting obscured and overwhelmed by larger counts
Poisson Model B
\[ \widehat{Purchases} = \beta_0 \ + \ \beta_1 \ Label \ + \ \beta_2 \ Rating \ + \ \beta_3 \ AcidityIndex \ + \ \beta_4 \ Alcohol \ + \ \epsilon \]
# Poisson logistic regression
poisson_B <- glm(Cases_Purchased ~
Label_Appeal +
Quality_Rating +
Acidity_Index +
Alcohol_Level,
family = poisson(link = 'log'),
data = wine_data)
poisson_B##
## Call: glm(formula = Cases_Purchased ~ Label_Appeal + Quality_Rating +
## Acidity_Index + Alcohol_Level, family = poisson(link = "log"),
## data = wine_data)
##
## Coefficients:
## (Intercept) Label_Appeal Quality_Rating Acidity_Index Alcohol_Level
## 1.187640 0.180776 0.182244 -0.046866 0.005472
##
## Degrees of Freedom: 5001 Total (i.e. Null); 4997 Residual
## Null Deviance: 4441
## Residual Deviance: 3008 AIC: 17970# Print results
stargazer( poisson_B,
type = "text",
title = "Table 5: Logistic Poisson Model A ",
dep.var.labels = c("Cases of Wine Purchased")
)##
## Table 5: Logistic Poisson Model A
## =============================================
## Dependent variable:
## ---------------------------
## Cases of Wine Purchased
## ---------------------------------------------
## Label_Appeal 0.181***
## (0.009)
##
## Quality_Rating 0.182***
## (0.008)
##
## Acidity_Index -0.047***
## (0.006)
##
## Alcohol_Level 0.005***
## (0.002)
##
## Constant 1.188***
## (0.059)
##
## ---------------------------------------------
## Observations 5,002
## Log Likelihood -8,977.854
## Akaike Inf. Crit. 17,965.710
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Comments: After employing the poisson logistic model, there are only three out of five regressors chosen from the stepwise method in the previous linear models shown to be statistically significant. Both alcohol level and pH level become insignificant in this model.
I think this intuitively makes sense, as acidity index provides a more direct and clear signal about the acidity of a particular wine than pH level does. Perhaps, pH level matters only to wine producers as they use scientific ways to measure, and to standardize their products. However, it requires some training for consumers to understand the pH level which some population does not have. Also, we see here that alcohol level become statistically insignificant. This might be a result of multiple factors; A) wine tend to have similar alcohol level, thus consumers aren’t considering it as a factor in their wine purchase process. B) For wine, in many cases, are viewed as a complement to some dishes. Therefore, the flavor of a wine is more of a deciding factor than alcohol level in a consumer’s wine purchase process.
Coefficient Interpretation:
Each rank increase in a wine’s label appearances, the odds of a consumer purchase one additional case of that wine would increase by 1.20%.
Each unit increase in acidity index of a wine will lead to the odds of a consumer purchase one less case of that wine by 0.973%.
# Check for over-dispersion
dispersiontest(poisson_B)##
## Overdispersion test
##
## data: poisson_B
## z = -40.986, p-value = 1
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 0.4119926Comments: We have a very slight overdispersion issue here. However, the magnitude of it so small that it can be ignored.
- Over-dispersion: In Poisson models, the variance should, in theory, equal to the mean. Over-dispersion issues arises when the variance in a Poisson model is much greater than its mean. If over-dispersion is present in a data set, the estimated standard errors and test statistics the overall goodness-of-fit will be distorted and adjustments must be made. Sometimes, variance can also be smaller than the mean and this is called a under-dispersion issue.
Negative Binomial Model A
nbm_A <- glm.nb(Cases_Purchased ~ . ,
data = wine_data)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reachedstargazer(nbm_A,
type = 'text',
title = "Table 5: Negative Binomial Model A "
)##
## Table 5: Negative Binomial Model A
## ================================================
## Dependent variable:
## ---------------------------
## Cases_Purchased
## ------------------------------------------------
## Fixed_Acidity -0.001
## (0.002)
##
## Volatile_Acidity -0.023
## (0.020)
##
## Citric_Acid 0.066*
## (0.037)
##
## Residual_Sugar -0.00002
## (0.0003)
##
## Chloride 0.016
## (0.032)
##
## Free_Sulfur_Dioxide -0.00001
## (0.0001)
##
## Total_Sulfur_Dioxide 0.00003
## (0.00005)
##
## Density -0.216
## (0.277)
##
## pH_Level 0.011
## (0.011)
##
## Sulphate -0.0002
## (0.012)
##
## Alcohol_Level 0.006***
## (0.002)
##
## Label_Appeal 0.180***
## (0.009)
##
## Acidity_Index -0.046***
## (0.007)
##
## Quality_Rating 0.182***
## (0.008)
##
## Constant 1.352***
## (0.282)
##
## ------------------------------------------------
## Observations 5,002
## Log Likelihood -8,975.046
## theta 146,680.100 (281,104.700)
## Akaike Inf. Crit. 17,980.090
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Comments: This version of the negative binomial model gives a result that is almost identical to what we’ve already gotten from the two poisson models. This is because when our sample is greater 100 (n >/= 100), both models would converge, hence, producing similar results. Also, This outcome can also be found in the rootgram below where we see the same trend of the red line as well as the same distance from 0 line across all bars. Both of these two things suggest that given the same regressors and large enough sample size, the Poisson and Negative Binary models will be exactly identical.
Additionally, it’s worth to point our that the pH level has a negative coefficient in the linear model but in both the poisson and negative binomial model, they have a positive coefficient value.
Coefficient Interpretation:
Each rank increase in a wine’s label appearances, the odds of a consumer purchase one additional case of that wine would increase by 20.8%.
Each unit increase in acidity index of a wine will lead to the odds of a consumer purchase one less case of that wine by 0.973%.
countreg::rootogram(nbm_A,
xlab = "Total Number of Cases of Wine Purchased",
main = "Fit of the Negative Binomial Model A")
Negative Binomial Model B
\[ \widehat{Purchases} = \beta_0 \ + \ \beta_1 \ Label \ + \ \beta_2 \ Rating \ + \ \beta_3 \ VAcidity \ + \ \beta_4 \ Alcohol \ + \ \beta_5 \ pH+ \ \epsilon \]
nbm_B <- glm.nb(Cases_Purchased ~
Label_Appeal +
Quality_Rating +
Volatile_Acidity +
Alcohol_Level +
pH_Level,
data = wine_data
)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reachedstargazer(nbm_B,
type = 'text',
title = "Table 5: Negative Binomial Model B ",
dep.var.labels = c("Cases of Wine Purchased")
)##
## Table 5: Negative Binomial Model B
## =============================================
## Dependent variable:
## ---------------------------
## Cases of Wine Purchased
## ---------------------------------------------
## Label_Appeal 0.176***
## (0.009)
##
## Quality_Rating 0.188***
## (0.008)
##
## Volatile_Acidity -0.031
## (0.020)
##
## Alcohol_Level 0.006***
## (0.002)
##
## pH_Level 0.016
## (0.011)
##
## Constant 0.773***
## (0.047)
##
## ---------------------------------------------
## Observations 5,002
## Log Likelihood -9,003.534
## theta 147,353.500 (284,632.000)
## Akaike Inf. Crit. 18,019.070
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01nbm_C <- glm.nb(Cases_Purchased ~
Label_Appeal +
Quality_Rating +
Acidity_Index +
Alcohol_Level,
data = wine_data
)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached5. Prediction
Test data set is already cleaned the same way as the training data set. I exclude the process in the report.
I chose Poisson Model B and Negative Binomial Model B for prediction because both has the lowest AIC value in their own category of models. Here, we excluded models from multi-linear models because the models aren’t suitable for classification predictions.
Prediction Outcome Distribution
# Poisson Model A
pdata <- predict.glm(poisson_B, newdata = test, type = 'response'
)
# Negative Binomial Model A
pdata2 <- predict.glm(nbm_B, newdata = test, type = 'response'
)
pdata3 <- predict.glm(nbm_C, newdata = test, type = 'response')par(mfrow = c(2, 2))
# Print Distribution of predictions
hist(pdata, xlab = "Number of cases purchased" , main = "Predictions of Poisson Model B")
hist(pdata2, xlab = "Number of cases purchased", main = "Predictions of Negative Binomial Model B")
hist(pdata3, xlab = "Number of cases purchased", main = "Predictions of Negative Binomial Model C")
# Round the results
rounded_poisson <- round(pdata)
rounded_negative_bi <- round(pdata2)
rounded_negative_bi2 <- round(pdata3)Confusion Matrix for Poisson Model
# Poisson Model
confusionMatrix(data = as.factor(rounded_poisson),
reference = as.factor(test$Cases_Purchased),
positive = "1"
) ## Warning in levels(reference) != levels(data): longer object length is not a
## multiple of shorter object length## Warning in confusionMatrix.default(data = as.factor(rounded_poisson), reference
## = as.factor(test$Cases_Purchased), : Levels are not in the same order for
## reference and data. Refactoring data to match.## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0
## 2 22 9 39 46 8 2 0 0 0
## 3 62 2 57 171 155 34 5 0 0
## 4 9 0 2 60 204 110 23 0 0
## 5 0 0 0 3 38 77 32 9 0
## 6 0 0 0 0 4 27 22 8 0
## 7 0 0 0 0 0 3 4 2 2
## 8 0 0 0 0 0 0 1 0 0
##
## Overall Statistics
##
## Accuracy : 0.4113
## 95% CI : (0.3839, 0.4392)
## No Information Rate : 0.3267
## P-Value [Acc > NIR] : 2.158e-10
##
## Kappa : 0.2352
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
## Sensitivity 0.00000 0.000000 0.39796 0.6107 0.4988 0.3043
## Specificity 1.00000 1.000000 0.92461 0.6759 0.7580 0.9179
## Pos Pred Value NaN NaN 0.30952 0.3519 0.5000 0.4843
## Neg Pred Value 0.92572 0.991214 0.94760 0.8577 0.7571 0.8390
## Prevalence 0.07428 0.008786 0.07827 0.2236 0.3267 0.2021
## Detection Rate 0.00000 0.000000 0.03115 0.1366 0.1629 0.0615
## Detection Prevalence 0.00000 0.000000 0.10064 0.3882 0.3259 0.1270
## Balanced Accuracy 0.50000 0.500000 0.66128 0.6433 0.6284 0.6111
## Class: 6 Class: 7 Class: 8
## Sensitivity 0.25287 0.105263 0.0000000
## Specificity 0.96652 0.992701 0.9992000
## Pos Pred Value 0.36066 0.181818 0.0000000
## Neg Pred Value 0.94542 0.986301 0.9984013
## Prevalence 0.06949 0.015176 0.0015974
## Detection Rate 0.01757 0.001597 0.0000000
## Detection Prevalence 0.04872 0.008786 0.0007987
## Balanced Accuracy 0.60970 0.548982 0.4996000Confusion Matrix for Negative Binomial Models
# Negative Binomial Model 2
confusionMatrix(data = as.factor(rounded_negative_bi),
reference = as.factor(test$Cases_Purchased),
positive = "1"
) ## Warning in levels(reference) != levels(data): longer object length is not a
## multiple of shorter object length## Warning in confusionMatrix.default(data = as.factor(rounded_negative_bi), :
## Levels are not in the same order for reference and data. Refactoring data to
## match.## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0
## 2 24 9 50 48 12 2 0 0 0
## 3 61 2 46 168 170 31 5 0 0
## 4 8 0 2 64 193 122 25 1 0
## 5 0 0 0 0 29 71 31 7 0
## 6 0 0 0 0 5 26 22 10 0
## 7 0 0 0 0 0 1 4 1 2
## 8 0 0 0 0 0 0 0 0 0
##
## Overall Statistics
##
## Accuracy : 0.4034
## 95% CI : (0.376, 0.4311)
## No Information Rate : 0.3267
## P-Value [Acc > NIR] : 7.598e-09
##
## Kappa : 0.2256
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
## Sensitivity 0.00000 0.000000 0.51020 0.6000 0.4719 0.28063
## Specificity 1.00000 1.000000 0.91768 0.6759 0.7367 0.93293
## Pos Pred Value NaN NaN 0.34483 0.3478 0.4651 0.51449
## Neg Pred Value 0.92572 0.991214 0.95664 0.8544 0.7419 0.83662
## Prevalence 0.07428 0.008786 0.07827 0.2236 0.3267 0.20208
## Detection Rate 0.00000 0.000000 0.03994 0.1342 0.1542 0.05671
## Detection Prevalence 0.00000 0.000000 0.11581 0.3858 0.3315 0.11022
## Balanced Accuracy 0.50000 0.500000 0.71394 0.6380 0.6043 0.60678
## Class: 6 Class: 7 Class: 8
## Sensitivity 0.25287 0.0526316 0.000000
## Specificity 0.96481 0.9943228 1.000000
## Pos Pred Value 0.34921 0.1250000 NaN
## Neg Pred Value 0.94533 0.9855305 0.998403
## Prevalence 0.06949 0.0151757 0.001597
## Detection Rate 0.01757 0.0007987 0.000000
## Detection Prevalence 0.05032 0.0063898 0.000000
## Balanced Accuracy 0.60884 0.5234772 0.500000# Negative Binomial Model 2
confusionMatrix(data = as.factor(rounded_negative_bi2),
reference = as.factor(test$Cases_Purchased),
positive = "1"
) ## Warning in levels(reference) != levels(data): longer object length is not a
## multiple of shorter object length## Warning in confusionMatrix.default(data = as.factor(rounded_negative_bi2), :
## Levels are not in the same order for reference and data. Refactoring data to
## match.## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0 0
## 2 22 9 39 46 8 2 0 0 0
## 3 62 2 57 171 155 34 5 0 0
## 4 9 0 2 60 204 110 23 0 0
## 5 0 0 0 3 38 77 32 9 0
## 6 0 0 0 0 4 27 22 8 0
## 7 0 0 0 0 0 3 4 2 2
## 8 0 0 0 0 0 0 1 0 0
##
## Overall Statistics
##
## Accuracy : 0.4113
## 95% CI : (0.3839, 0.4392)
## No Information Rate : 0.3267
## P-Value [Acc > NIR] : 2.158e-10
##
## Kappa : 0.2352
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
## Sensitivity 0.00000 0.000000 0.39796 0.6107 0.4988 0.3043
## Specificity 1.00000 1.000000 0.92461 0.6759 0.7580 0.9179
## Pos Pred Value NaN NaN 0.30952 0.3519 0.5000 0.4843
## Neg Pred Value 0.92572 0.991214 0.94760 0.8577 0.7571 0.8390
## Prevalence 0.07428 0.008786 0.07827 0.2236 0.3267 0.2021
## Detection Rate 0.00000 0.000000 0.03115 0.1366 0.1629 0.0615
## Detection Prevalence 0.00000 0.000000 0.10064 0.3882 0.3259 0.1270
## Balanced Accuracy 0.50000 0.500000 0.66128 0.6433 0.6284 0.6111
## Class: 6 Class: 7 Class: 8
## Sensitivity 0.25287 0.105263 0.0000000
## Specificity 0.96652 0.992701 0.9992000
## Pos Pred Value 0.36066 0.181818 0.0000000
## Neg Pred Value 0.94542 0.986301 0.9984013
## Prevalence 0.06949 0.015176 0.0015974
## Detection Rate 0.01757 0.001597 0.0000000
## Detection Prevalence 0.04872 0.008786 0.0007987
## Balanced Accuracy 0.60970 0.548982 0.49960006. Conclusion Remarks
The negative binomial model 1 has a higher sensitivity in class 1, 2, and 3. The Poisson model has a better sensitivity in class 5, 6, and 7. Overall specificity values are roughly the same throughout both models. Both models are pretty good in terms of not predicting an outcome when it’s not that number, that is why we see the specificity rate is high in both models. However, the accuracy rate is higher in the poisson model, at a level of 41% , than the negative binomial model which has a slightly lower level of 40%. Therefore, we could say that the Poisson model has a better predictive power. However, the overall performance of the two models are extremely similar. Both models also have a similar kappa value around 0.22 which suggests that they are in the acceptable range.
The negative binomial model 2 is identical to the Poisson model.