# Lesson 10 | 簡易數據分析
## 第一節:描述性統計(1)
* 這節課我們終於要開始進行數據分析了,我們先簡單介紹幾個描述性統計的函數
* 請至[這裡](https://linchin.ndmctsgh.edu.tw/data/Example_data.csv)下載本週的範例資料,另外,請同學們先幫我安裝```pandas```以及```numpy```。
```python=
# 載入套件
import pandas as pd
import numpy as np
```
```python=
data = pd.read_csv("Example_data.csv")
data
```
* 以下幾個```numpy```的函式為常見的描述性統計,分別是:
* 函式「.mean()」可以幫助我們計算算數平均值。
```python=
np.mean(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「.std()」可以幫助我們計算標準差。
```python=
np.std(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「.var()」可以幫助我們計算變異數。
```python=
np.var(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「.median()」可以幫助我們計算中位數。
```python=
np.median(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「quantile()」可以幫助我們計算百分位數。
```python=
np.quantile(data['eGFR'][pd.isna(data['eGFR']) == False], [0, 0.25, 0.50, 0.75, 1])
```
```python=
# 50百分位數
np.quantile(data['eGFR'][pd.isna(data['eGFR']) == False], 0.50)
```
```python=
# 95百分位數
np.quantile(data['eGFR'][pd.isna(data['eGFR']) == False], 0.95)
```
* 函數「min()」可以幫助我們找出最小值。
```python=
np.min(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「max()」可以幫助我們找出最大值。
```python=
np.max(data['eGFR'][pd.isna(data['eGFR']) == False])
```
* 函數「table()」可以幫助我們計數各類別的頻率。
```python=
np.unique(data['Disease'], return_counts = True)
```
* 稍微加工一下就可以產生各類別的百分比。
```python=
np.unique(data['Disease'], return_counts = True)[1]/len(data['Disease'])
```
## 第一節:描述性統計(2)
* 利用索引,我們將能得到某種條件下的描述性統計,舉例來說我們想獲得Case組的eGFR平均數,可以使用下列語法。
```python=
clean_data = data[(pd.isna(data['eGFR']) == False) & (pd.isna(data['Disease']) == False)]
np.mean(clean_data['eGFR'][clean_data['Disease'] == 1])
```
* Control組的亦同
```python=
clean_data = data[(pd.isna(data['eGFR']) == False) & (pd.isna(data['Disease']) == False)]
np.mean(clean_data['eGFR'][clean_data['Disease'] == 0])
```
* 如果我們想要把平均數±標準差表示出來,在Python裡的String可以直接使用```+```號進行處理:
```python=
clean_data = data[(pd.isna(data['eGFR']) == False) & (pd.isna(data['Disease']) == False)]
str(np.mean(clean_data['eGFR'][clean_data['Disease'] == 1])) + ' ± ' + str(np.std(clean_data['eGFR'][clean_data['Disease'] == 1]))
```
* 我們發現這樣的呈現相當醜,我們可以使用函數「round()」來指定我們想要的小數點位數:
```python=
clean_data = data[(pd.isna(data['eGFR']) == False) & (pd.isna(data['Disease']) == False)]
m = np.mean(clean_data['eGFR'][clean_data['Disease'] == 1])
s = np.std(clean_data['eGFR'][clean_data['Disease'] == 1])
str(round(m, 2)) + ' ± ' + str(round(s, 2))
```
## 練習1:
* 現在請利用自訂函數功能,寫出兩個函數,規格如下:
* 函數之input必須為一類別變項及一連續變項。
* 函數中必須設定一個變項,讓使用者能決定顯示小數點的位數。
* 問題1:第一個函數必須輸出在該類別變項為不同類別時,該連續變項之『平均數±標準差』,例如:
```python=
func_1(x = data['Disease'], y = data['eGFR'])
```
```python=
func_1(x = data['Education'], y = data['SBP'], dig = 1)
```
<details>
<summary>解答</summary>
```python=
def func_1(x, y, dig = 3):
x = x[(pd.isna(x) == False) & (pd.isna(y) == False)]
y = y[(pd.isna(x) == False) & (pd.isna(y) == False)]
levels = list(set(x))
out_list = list()
for i in levels:
m = np.mean(y[x == i])
s = np.std(y[x == i])
m_s_str = str(round(m, dig)) + ' ± ' + str(round(s, dig))
out_list.append(m_s_str)
return out_list
```
</details>
* 問題2:第二個函數必須輸出在該類別變項為不同類別時,該連續變項之『中位數(25百分位-75百分位)』,例如:
```python=
func_2(x = data['Disease'], y = data['eGFR'])
```
```python=
func_2(x = data['Education'], y = data['SBP'], dig = 1)
```
<details>
<summary>解答</summary>
```python=
def func_2(x, y, dig = 3):
x = x[(pd.isna(x) == False) & (pd.isna(y) == False)]
y = y[(pd.isna(x) == False) & (pd.isna(y) == False)]
levels = list(set(x))
out_list = list()
for i in levels:
q_list = np.quantile(y[x == i], [0.25, 0.50, 0.75])
q_str = str(round(q_list[1], dig)) + ' (' + str(round(q_list[0], dig)) + '-' + str(round(q_list[2], dig)) + ')'
out_list.append(q_str)
return out_list
```
</details>
## 第二節:線性迴歸(1)
* 線性迴歸(linear regression)是高中已經教過的基本分析技術,目標是從一群數據點中找出$x_i$與$y_i$之間的線性關係,其方程式如下:
$$\hat{y_i} = b_{0} + b_{1}x_i$$
* 這裡在高中上課的時候有提到我們要如何解這個$b_0$與$b_1$的係數,我們是透過先定義誤差為$ϵ$:
$$y_i = b_{0} + b_{1}x_i + \epsilon$$
$$\epsilon = y_i - b_{0} - b_{1}x_i$$
* 接著再定義找出一組$b0$與$b1$使$\sum{\epsilon}^2$最小化:
$$\mbox{min}(\sum{\epsilon}^2) = \mbox{min}(\sum\limits_{i=1}^{n}(y_i - b_{0} - b_{1}x_i)^2)$$
## 第二節:線性迴歸(2)
* 若我們要求解線性回歸,我們可以分別對$b_0$與$b_1$進行偏微分得到公式,求解其係數,我們也可以透過第三方套件來完成這件事情。
* 我們接著將使用```scikit-learn```套件,在Python裡[scikit-learn](https://scikit-learn.org/stable/)是最常用來做數據分析以及機器學習的功能強大套件請先安裝:

* 我們可以在[這裡](https://scikit-learn.org/stable/modules/linear_model.html)看到關於線型回歸的相關說明。
```python=
from sklearn import linear_model
reg = linear_model.LinearRegression()
reg.fit([[0], [1], [2]], [0, 1, 2])
```
```python=
print(reg.intercept_)
print(reg.coef_)
```
## 第三節:線性迴歸(3)
* 假設我們使用的是剛剛的資料,讓我們看看使用「linear_model」的方法:
* 記得要把```numpy array```的shape改成2維
```python=
SBP = data['SBP'][(pd.isna(data['SBP']) == False) & (pd.isna(data['DBP']) == False)].to_numpy()
DBP = data['DBP'][(pd.isna(data['DBP']) == False) & (pd.isna(data['SBP']) == False)].to_numpy()
DBP = DBP.reshape((DBP.shape[0], 1))
reg = linear_model.LinearRegression()
reg.fit(DBP, SBP)
# 讓我們看一下預測出來線性回歸的參數。
print(reg.intercept_, reg.coef_)
```
* 這樣的線性回歸方程式可以放入多個變項。
* 如果我們要用SBP以及DBP同時預測eGFR,可以用下列程式碼:
```python=
na_idx = (pd.isna(data['SBP']) == False) & \
(pd.isna(data['DBP']) == False) & \
(pd.isna(data['eGFR']) == False)
clean_data = data[na_idx]
x = clean_data[['SBP', 'DBP']].to_numpy()
y = clean_data['eGFR'].to_numpy()
reg = linear_model.LinearRegression()
reg.fit(x, y)
# 也是一樣,讓我們看一下預測出來線性回歸的參數。
print(reg.intercept_, reg.coef_)
```
## 練習2:
* 現在我們想要使用某病人的SBP以及DBP來預測eGFR,所以希望你能寫出一個自訂函數,讓使用者輸入單筆SBP與DBP的值,並為他計算出eGFR的值。
* 但是這時候我們遇到一個問題,那就是預測公式鐵定是以這樣的型式表現,因此假設有一個新的個案缺少SBP或DBP其中一項,那他將無法使用這個式子:
$$\hat{eGFR_i} = b_{0} + b_{1}SBP_i + b_{2}DBP_i$$
* 因此這個函數需要有一個功能,就是當使用者沒有輸入SBP而只輸入DBP時,你必須用這個式子協助他進行預測(依此類推):
$$\hat{eGFR_i} = b_{0} + b_{1}DBP_i$$
```python=
print(pred_func(SBP = 100))
print(pred_func(DBP = 100))
print(pred_func(SBP = 100, DBP = 100))
```
<details>
<summary>解答</summary>
```python=
def pred_func(SBP = None, DBP = None):
na_idx = (pd.isna(data['SBP']) == False) & \
(pd.isna(data['DBP']) == False) & \
(pd.isna(data['eGFR']) == False)
clean_data = data[na_idx]
if SBP is None and DBP is None:
return("Please enter at least one variable.")
elif SBP is None:
x = clean_data['DBP'].to_numpy()
y = clean_data['eGFR'].to_numpy()
x = x.reshape((x.shape[0], 1))
reg = linear_model.LinearRegression()
reg.fit(x, y)
return reg.intercept_ + reg.coef_[0] * DBP
elif DBP is None:
x = clean_data['SBP'].to_numpy()
y = clean_data['eGFR'].to_numpy()
x = x.reshape((x.shape[0], 1))
reg = linear_model.LinearRegression()
reg.fit(x, y)
return reg.intercept_ + reg.coef_[0] * SBP
else:
x = clean_data[['SBP', 'DBP']].to_numpy()
y = clean_data['eGFR'].to_numpy()
reg = linear_model.LinearRegression()
reg.fit(x, y)
return reg.intercept_ + reg.coef_[0] * SBP + reg.coef_[1] * DBP
```
</details>
* 同學們是否可以想到,我們其實應該以物件導向的設計來撰寫這個功能?
## 第三節:在線性迴歸中使用類別變項(1)
* 你可能有注意到剛剛的資料集中,其中有些變項是被編碼成0、1以及0、1、2,這樣的變項可以放入線性迴歸中嗎?
* 答案是可以的,對於已經編碼成0、1的變項我們可以直接把他放入線性迴歸,但是0、1、2的變項我們就不能直接放入線性迴歸之中。
* 那我們該怎樣處理這種有3個類別以上的類別變項呢?方法很簡單就是要把他改成這樣的編碼格式(他們稱之為one-hot encoding):
```python=
x = [0, 1, 2, 1, None, 2]
np.array([0, 1, 0, 1, None, 0,
0, 0, 1, 0, None, 1]).reshape((6, 2), order='F')
```
## 第三節:在線性迴歸中使用類別變項(2)
* 至於要怎樣把他轉換成這個格式呢?我們可以用這種方式:
```python=
cato = list(set(data['Education'][pd.isna(data['Education']) == False]))
print(cato)
data['education_1'] = 0
data['education_2'] = 0
for i in range(len(data['Education'])):
if data.loc[i, 'Education'] == cato[1]:
data.loc[i, 'education_1'] = 1
elif data.loc[i, 'Education'] == cato[2]:
data.loc[i, 'education_2'] = 1
else:
data.loc[i, 'education_1'] = None
data.loc[i, 'education_2'] = None
print(data)
```
* 當然,你也可以在第三方套件裡找到這個功能的函式:
```python=
pd.get_dummies(data['Education'], dummy_na = True)
```
* 你應該會注意到「get_dummies」與我們想要的不太一樣(NA變新的一欄),所以我們還要稍微修正一下他:
```python=
def one_hot_coding(cato_df):
dummies = pd.get_dummies(cato_df, dummy_na = True)
nan_col = dummies[np.nan]
dummies[nan_col == 1] = None
col_name = dummies.columns[pd.isna(dummies.columns) == False]
col_name = col_name[1:]
final_dummies = dummies[col_name]
return final_dummies
```
* 稍微看一下結果。
```python=
one_hot_coding(data['Education'])
```
* 有了這項技術以後,那就可以把剛剛的Income以及Education都納入回歸之內了:
* 先整理資料:
```python=
data_edu = one_hot_coding(data['Education'])
data_income = one_hot_coding(data['Income'])
data_edu = data_edu.rename(columns = {1.0: 'educa_1', 2.0 : 'educa_2'})
data_income = data_income.rename(columns = {1.0: 'income_1', 2.0 : 'income_2'})
data1 = pd.concat([data, data_edu, data_income], axis = 1)
```
```python=
from sklearn import linear_model
# 進行線性回歸模型的建立。
na_idx = (pd.isna(data1['SBP']) == False) & \
(pd.isna(data1['DBP']) == False) & \
(pd.isna(data1['Income']) == False) & \
(pd.isna(data1['Education']) == False) & \
(pd.isna(data1['eGFR']) == False)
clean_data1 = data1[na_idx]
x_var = ['SBP', 'DBP', 'income_1', 'income_2', 'educa_1', 'educa_2']
x = clean_data1[x_var].to_numpy()
y = clean_data1['eGFR'].to_numpy()
reg = linear_model.LinearRegression()
reg.fit(x, y)
# 也是一樣,讓我們看一下預測出來線性回歸的參數。
print(reg.intercept_, reg.coef_)
```
## 總結
* 有了關於```pandas```以及```numpy```的熟悉使用,再加上```scikit-learn```我們可以輕易地使用python進行資料分析。
* 由於課程設計的限制,我們並沒有完整的介紹Python裡的統計分析功能,在後面的課程將會著重在機器學習以及深度學習等演算法上。
* 現在你應該有能力獨自完成從蒐集資料至資料分析的整個過程,程式語言將能協助你建立pipeline,讓你能隨時看到最新結果!