# Estimating alcohol content in red wines

* Author: Martin Rožnovják
* Last edited: 2019-02-11
* Organization: Metropolia University of Applied Sciences

## What is this?

This notebook is a school assignment for a course called *Cognitive Systems - Mathematics and Methods*.

Its objective is to conduct linear regression analysis on wine properties
from the following UCI dataset https://archive.ics.uci.edu/ml/datasets/Wine+Quality.

Personally, I prefer red wine over white so I'll analyze the red wine dataset.
I'll be interested in estimating the alcohol content of the wines instead
of their quality - just out of interest and objectivity of the measurements,
not because I'm a student... :)

## Let's begin
### All the necessary imports


```python
import IPython

import scipy
import sklearn

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

from sklearn import metrics, preprocessing
from sklearn.linear_model import LinearRegression
```

Setting options, in particular those for plotting, often interferes with Jupyter's own initializations on import. Therefore, I just do it in the next cell.


```python
### plotting ###
# prettier and larger basic graphs
sns.set(rc={
    'figure.figsize':(18,8),
    'axes.titlesize':14,
})

### pandas ###
# no need to see many decimal places and makes nicer horizontal fits :-)
pd.options.display.float_format = '{:.3f}'.format
# pd.options.display.precision = 3
# make the tables more compact vertically, too
pd.options.display.max_rows = 20

### numpy ###
# same as for pandas - max. 3 decimal places
np.set_printoptions(formatter={'float_kind':'{:.3f}'.format})
# np.set_printoptions(precision=3)
```

### Fetching the dataset


```python
wines_df = pd.read_csv(
    r'https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv',
    sep=';',
)

# a peek whether it went alright
wines_df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>fixed acidity</th>
      <th>volatile acidity</th>
      <th>citric acid</th>
      <th>residual sugar</th>
      <th>chlorides</th>
      <th>free sulfur dioxide</th>
      <th>total sulfur dioxide</th>
      <th>density</th>
      <th>pH</th>
      <th>sulphates</th>
      <th>alcohol</th>
      <th>quality</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>7.400</td>
      <td>0.700</td>
      <td>0.000</td>
      <td>1.900</td>
      <td>0.076</td>
      <td>11.000</td>
      <td>34.000</td>
      <td>0.998</td>
      <td>3.510</td>
      <td>0.560</td>
      <td>9.400</td>
      <td>5</td>
    </tr>
    <tr>
      <th>1</th>
      <td>7.800</td>
      <td>0.880</td>
      <td>0.000</td>
      <td>2.600</td>
      <td>0.098</td>
      <td>25.000</td>
      <td>67.000</td>
      <td>0.997</td>
      <td>3.200</td>
      <td>0.680</td>
      <td>9.800</td>
      <td>5</td>
    </tr>
    <tr>
      <th>2</th>
      <td>7.800</td>
      <td>0.760</td>
      <td>0.040</td>
      <td>2.300</td>
      <td>0.092</td>
      <td>15.000</td>
      <td>54.000</td>
      <td>0.997</td>
      <td>3.260</td>
      <td>0.650</td>
      <td>9.800</td>
      <td>5</td>
    </tr>
    <tr>
      <th>3</th>
      <td>11.200</td>
      <td>0.280</td>
      <td>0.560</td>
      <td>1.900</td>
      <td>0.075</td>
      <td>17.000</td>
      <td>60.000</td>
      <td>0.998</td>
      <td>3.160</td>
      <td>0.580</td>
      <td>9.800</td>
      <td>6</td>
    </tr>
    <tr>
      <th>4</th>
      <td>7.400</td>
      <td>0.700</td>
      <td>0.000</td>
      <td>1.900</td>
      <td>0.076</td>
      <td>11.000</td>
      <td>34.000</td>
      <td>0.998</td>
      <td>3.510</td>
      <td>0.560</td>
      <td>9.400</td>
      <td>5</td>
    </tr>
  </tbody>
</table>
</div>



No need for the "quality" column, let's forget it...


```python
wines_df.drop(columns=["quality"], inplace=True)
```

### Basic stats and info


```python
print('Summary about the DataFrame:')
wines_df.info();
```

    Summary about the DataFrame:
    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 1599 entries, 0 to 1598
    Data columns (total 11 columns):
    fixed acidity           1599 non-null float64
    volatile acidity        1599 non-null float64
    citric acid             1599 non-null float64
    residual sugar          1599 non-null float64
    chlorides               1599 non-null float64
    free sulfur dioxide     1599 non-null float64
    total sulfur dioxide    1599 non-null float64
    density                 1599 non-null float64
    pH                      1599 non-null float64
    sulphates               1599 non-null float64
    alcohol                 1599 non-null float64
    dtypes: float64(11)
    memory usage: 137.5 KB


All looks good so far, nothing missing (just as promised in the dataset description), correct data types...


```python
print('Basic descriptive statistics:')
wines_df.describe()
```

    Basic descriptive statistics:





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>fixed acidity</th>
      <th>volatile acidity</th>
      <th>citric acid</th>
      <th>residual sugar</th>
      <th>chlorides</th>
      <th>free sulfur dioxide</th>
      <th>total sulfur dioxide</th>
      <th>density</th>
      <th>pH</th>
      <th>sulphates</th>
      <th>alcohol</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
      <td>1599.000</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>8.320</td>
      <td>0.528</td>
      <td>0.271</td>
      <td>2.539</td>
      <td>0.087</td>
      <td>15.875</td>
      <td>46.468</td>
      <td>0.997</td>
      <td>3.311</td>
      <td>0.658</td>
      <td>10.423</td>
    </tr>
    <tr>
      <th>std</th>
      <td>1.741</td>
      <td>0.179</td>
      <td>0.195</td>
      <td>1.410</td>
      <td>0.047</td>
      <td>10.460</td>
      <td>32.895</td>
      <td>0.002</td>
      <td>0.154</td>
      <td>0.170</td>
      <td>1.066</td>
    </tr>
    <tr>
      <th>min</th>
      <td>4.600</td>
      <td>0.120</td>
      <td>0.000</td>
      <td>0.900</td>
      <td>0.012</td>
      <td>1.000</td>
      <td>6.000</td>
      <td>0.990</td>
      <td>2.740</td>
      <td>0.330</td>
      <td>8.400</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>7.100</td>
      <td>0.390</td>
      <td>0.090</td>
      <td>1.900</td>
      <td>0.070</td>
      <td>7.000</td>
      <td>22.000</td>
      <td>0.996</td>
      <td>3.210</td>
      <td>0.550</td>
      <td>9.500</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>7.900</td>
      <td>0.520</td>
      <td>0.260</td>
      <td>2.200</td>
      <td>0.079</td>
      <td>14.000</td>
      <td>38.000</td>
      <td>0.997</td>
      <td>3.310</td>
      <td>0.620</td>
      <td>10.200</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>9.200</td>
      <td>0.640</td>
      <td>0.420</td>
      <td>2.600</td>
      <td>0.090</td>
      <td>21.000</td>
      <td>62.000</td>
      <td>0.998</td>
      <td>3.400</td>
      <td>0.730</td>
      <td>11.100</td>
    </tr>
    <tr>
      <th>max</th>
      <td>15.900</td>
      <td>1.580</td>
      <td>1.000</td>
      <td>15.500</td>
      <td>0.611</td>
      <td>72.000</td>
      <td>289.000</td>
      <td>1.004</td>
      <td>4.010</td>
      <td>2.000</td>
      <td>14.900</td>
    </tr>
  </tbody>
</table>
</div>




```python
print('Almost the same thing again, this time graphically:')
wines_df.plot.box(subplots=True, layout=(2, 6), figsize=(20, 10));
```

There seem to be several outliers, I will get back to them soon.


```python
print('Yet different perspective, feature histograms:')
wines_df.hist(bins=20, figsize=(18, 14))
plt.tight_layout()
```


```python
ax = sns.heatmap(
    wines_df.corr(),
    annot=True,
    vmin=-1,
    vmax=1,
    cmap="coolwarm",
    fmt='0.2f',
    linewidths=1,
)
ax.set_title('Correlation Matrix');
```

### Training-testing split

Error estimates are more credible when done on "unseen" data, therefore, I will set aside randomly selected 30% of the data for testing. 

I will also make a separate training set with removed outliers for inspecting their role in the regression - how they skew the regression and influence generalization.
I will remove rows that contain a value more than 3 standard deviations in any of the columns.


```python
train_df, test_df = sklearn.model_selection.train_test_split(
    wines_df,
    test_size=0.3,
    shuffle=True,
)

# zcore = (x - x.avg)/x.std
good_zscore = np.abs(scipy.stats.zscore(train_df)) <= 3
trim_index = good_zscore.all(axis=1) # all values in a row have good zscore
train_df_trimmed  = train_df[trim_index]
```

Splitting into inputs (explanatory features) and results (responses).


```python
def xy_split(df):
    x_labels = wines_df.columns[:-1]
    y_label = wines_df.columns[-1]
    return df[x_labels], df[y_label]

x_train,      y_train      = xy_split(train_df)
x_train_trim, y_train_trim = xy_split(train_df_trimmed)
x_test,       y_test       = xy_split(test_df)
```


```python
# a simple utility
def compare_distributions(
    title='Kernel densities (distribution) of response values',
    ax=None,
    kde=True,
    hist=False,
    **label_value_pairs
):
    for label, values in label_value_pairs.items():
        ax = sns.distplot(
            values,
            kde=kde,
            hist=hist,
            label=label,
            ax=ax,
            axlabel=False
        )
        
    if title: ax.set_title(title)
    # if hist it won't set the legend automatically
    if hist: ax.legend(label_value_pairs.keys())
    
    return ax
```


```python
# ordinary histogram would be quite confusing...
compare_distributions(
    Training=y_train,
    Trimmed_training=y_train_trim,
    Validation=y_test
);
```


```python
print('Training data:')
train_df.plot.box(subplots=True, figsize=(18,5))
plt.tight_layout()
```


```python
print('Testing data:')
test_df.plot.box(subplots=True, figsize=(18,5))
plt.tight_layout()
```

### Model

The model and regression method are perhaps the simplest possible - i.e. a multilinear model
$\theta_0 + \theta_1 a_1 + \theta_2 a_2 + \dots + \theta_n a_n$
where $a_i$ is the $i$-th attribute (explanatory feature) and $\theta_i, i \in \{0, 1, \dots, n\}$
are the regression coefficients obtained using the least squares method.

If the data has unit variance and zero mean, the absolute value of the resulting coefficients can be
interpreted as the "importance/usefulness" of the corresponding attribute for the problem.
Therefore, I will shift and scale the training data accordingly.

<em>Note: it does not really matter whether the response (results) column is normalized.
If it is not, then the regression coefficients will be multiplied by the standard deviation
of the response column and intercept (the constant term) will be its mean.
I find working with the original values easier in this case (less to worry about).</em>


```python
# some utils

def print_regr_coefs(lin_reg):
    print(f'Intercept: {lin_reg.intercept_:.3f}')
    print('Coefs:', lin_reg.coef_)
    
def print_metrics(y_true, y_pred, comment=None, show_dists=True):
    mae = metrics.mean_absolute_error(y_true, y_pred)
    mse = metrics.mean_squared_error(y_true, y_pred)
    r2 = metrics.r2_score(y_true, y_pred)
    
    desc = f'MAE: {mae:.3f}  MSE: {mse:.3f}  R2: {r2:.3f}'
    if comment:
        desc += f'  [{comment}]'
        
    if show_dists:
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 3))
        fig.suptitle(desc)
        
        residuals = (y_true - y_pred).values
        abs_residuals = np.abs(residuals)
        
        # distribution of residuals
        sns.distplot(residuals, ax=ax1)
        sns.distplot(abs_residuals, ax=ax1)
        start = np.floor(min(residuals))
        end = np.ceil(max(abs_residuals))
        ax1.xaxis.set_ticks(np.round(np.linspace(start, end, 11),2))
        ax1.legend(['Residuals','Absolute residuals'])
        
        # dist. of response values
        compare_distributions(
            title=None,
            ax=ax2,
            hist=True,
            Measured=y_true,
            Predicted=y_pred,
        )
        plt.show()
    else:
        print(desc)
    
    
def test_report(y_true, y_pred):
    err_df = pd.DataFrame()
    err_df['measured'] = y_true
    err_df['predicted'] = y_pred
    err_df['residual'] = err_df.measured - err_df.predicted
    err_df['abs_residual'] = np.abs(err_df.residual)

    # show statistics
    print('Response statistics:')
    IPython.display.display(err_df.describe())
    
    ## plot results and residuals
    # true vs. predicted
    ax1 = plt.subplot(1,2,1)
    ax1.set_title('Measured values vs. predicted')
    ax1.plot(y_pred, y_pred)
    err_df.plot.scatter('measured', 'predicted', c='orange', ax=ax1)
    
    # residuals 
    ax2 = plt.subplot(1,2,2)
    ax2.set_title('Residuals')
    # y=0 line
    x_infimum = np.floor(err_df.measured.min())
    x_supremum = np.ceil(err_df.measured.max())
    ax2.plot([x_infimum, x_supremum], [0, 0])
    err_df.plot.scatter('measured', 'residual', c='orange', ax=ax2)

    plt.tight_layout()
    plt.show()
```


```python
print('Regressing on all features, whole training set\n')

# scaling
x_scaler = sklearn.preprocessing.StandardScaler()
x_scaler.fit(x_train)
x_train_scaled = x_scaler.transform(x_train)
x_train_trim_scaled = x_scaler.transform(x_train_trim)
x_test_scaled = x_scaler.transform(x_test)

# regression
reg = LinearRegression().fit(x_train_scaled, y_train)

# reports
print_regr_coefs(reg)

y_pred_train = reg.predict(x_train_scaled)
print_metrics(y_train, y_pred_train, 'training')

y_pred_test = reg.predict(x_test_scaled)
print_metrics(y_test, y_pred_test, 'validation')
print()
test_report(y_test, y_pred_test)
```

The mean average error (same as mean absolute residual) tells us
that the model is on average off by 0.48% off from the actual percentual content of the alcohol,
considering that the range of the alcohol volume is 8.4% - 14.9%,
this can be translated to about 7.4% relative error.
From the "Response statistics" table, based on the quantiles of absolute residual,
we can see that a quarter of the results are predicted within 0.2% (absolute error),
half of them within 0.4%, three quarters only within 0.7%
and the worst error is 2.4% off.


```python
print('Features by "importance" for the last regression:')
coef_df = pd.DataFrame()
coef_df['feature'] = x_train.columns
coef_df['coef'] = reg.coef_
coef_df['abs_coef'] = abs(coef_df['coef'])
coef_df['natural_coef'] = x_scaler.inverse_transform(reg.coef_)
coef_df.sort_values('abs_coef', ascending=False, inplace=True)
coef_df.reset_index(drop=True, inplace=True)
coef_df
```

    Features by "importance" for the last regression:





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>feature</th>
      <th>coef</th>
      <th>abs_coef</th>
      <th>natural_coef</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>density</td>
      <td>-1.118</td>
      <td>1.118</td>
      <td>0.995</td>
    </tr>
    <tr>
      <th>1</th>
      <td>fixed acidity</td>
      <td>0.904</td>
      <td>0.904</td>
      <td>9.912</td>
    </tr>
    <tr>
      <th>2</th>
      <td>pH</td>
      <td>0.571</td>
      <td>0.571</td>
      <td>3.395</td>
    </tr>
    <tr>
      <th>3</th>
      <td>residual sugar</td>
      <td>0.397</td>
      <td>0.397</td>
      <td>3.093</td>
    </tr>
    <tr>
      <th>4</th>
      <td>sulphates</td>
      <td>0.200</td>
      <td>0.200</td>
      <td>0.690</td>
    </tr>
    <tr>
      <th>5</th>
      <td>citric acid</td>
      <td>0.170</td>
      <td>0.170</td>
      <td>0.305</td>
    </tr>
    <tr>
      <th>6</th>
      <td>chlorides</td>
      <td>-0.086</td>
      <td>0.086</td>
      <td>0.083</td>
    </tr>
    <tr>
      <th>7</th>
      <td>total sulfur dioxide</td>
      <td>-0.076</td>
      <td>0.076</td>
      <td>43.908</td>
    </tr>
    <tr>
      <th>8</th>
      <td>volatile acidity</td>
      <td>0.076</td>
      <td>0.076</td>
      <td>0.543</td>
    </tr>
    <tr>
      <th>9</th>
      <td>free sulfur dioxide</td>
      <td>-0.025</td>
      <td>0.025</td>
      <td>15.461</td>
    </tr>
  </tbody>
</table>
</div>




```python
# I will leave this here, perhaps I had too long weekend...
print('Formula for estimating alcohol content in red wine:')
print('(features by contribution, no need to scale)')
IPython.display.Math(
    r'\begin{align}' +
    r'{:.3f}\ &+\ '.format(reg.intercept_) + 
    ''.join([
        r'\\& {:+.3f} \cdot [\verb|{}|]'.format(coef, ftr)
        for coef, ftr in coef_df[['natural_coef', 'feature']].values
    ]) +
    r'\end{align}'
)
```

    Formula for estimating alcohol content in red wine:
    (features by contribution, no need to scale)





$\displaystyle \begin{align}10.396\ &+\ \\& +0.995 \cdot [\verb|density|]\\& +9.912 \cdot [\verb|fixed acidity|]\\& +3.395 \cdot [\verb|pH|]\\& +3.093 \cdot [\verb|residual sugar|]\\& +0.690 \cdot [\verb|sulphates|]\\& +0.305 \cdot [\verb|citric acid|]\\& +0.083 \cdot [\verb|chlorides|]\\& +43.908 \cdot [\verb|total sulfur dioxide|]\\& +0.543 \cdot [\verb|volatile acidity|]\\& +15.461 \cdot [\verb|free sulfur dioxide|]\end{align}$



## Conclusion

I experimented with leaving some columns and with the threshold for trimming outliers,
the results were, however, insignificant and inconsistent, so I won't clutter the notebook
with examples of those, I will just leave one quite general cell for quickly playing
with  these adjustments (just comment out columns, change training set, etc.).

This approximation method is trying to minimize the MSE metric.
In my opinion, it would be better to minimize error so that it favors
the majority of the data and not being influenced much by edge cases.


```python
### tinkering cell ###

# scaling
columns = [
    'density',
    'fixed acidity',
    'pH',
    'residual sugar',
    'sulphates',
    'citric acid',
    'volatile acidity',
    'chlorides',
    'total sulfur dioxide',
    'free sulfur dioxide'
]
x_scaler = sklearn.preprocessing.StandardScaler()
x_scaler.fit(x_train_trim[columns])
x_train_scaled = x_scaler.transform(x_train[columns])
x_train_trim_scaled = x_scaler.transform(x_train_trim[columns])
x_test_scaled = x_scaler.transform(x_test[columns])

# regression
reg = LinearRegression().fit(x_train_trim_scaled, y_train_trim)
# reg.intercept_ = 10.2

# reports
print_regr_coefs(reg)

y_pred_train_trim = reg.predict(x_train_trim_scaled)
print_metrics(y_train_trim, y_pred_train_trim, 'training trimmed')

y_pred_train = reg.predict(x_train_scaled)
print_metrics(y_train, y_pred_train, 'whole training data')

y_pred_test = reg.predict(x_test_scaled)
print_metrics(y_test, y_pred_test, 'validation')
print()
test_report(y_test, y_pred_test)
```

9c2bd69d9657c7aec08e2234d834dc1893a0dd46
ipynb
Jupyter Notebook
Cognitive_Systems-Mathematics_and_Methods/week04/Roznovjak_Assignment_4-Linear_regression_on_red_wines.ipynb
rozni/uni-ml
0667c7504927ea3bd1850d118708ea72b4b43430
MIT
Cognitive_Systems-Mathematics_and_Methods/week04/Roznovjak_Assignment_4-Linear_regression_on_red_wines.ipynb
rozni/uni-ml
0667c7504927ea3bd1850d118708ea72b4b43430
MIT
Cognitive_Systems-Mathematics_and_Methods/week04/Roznovjak_Assignment_4-Linear_regression_on_red_wines.ipynb
rozni/uni-ml
0667c7504927ea3bd1850d118708ea72b4b43430
MIT
Qwen/Qwen-72B
1. YES
2. YES