```julia
using CSV
using DataFrames
using PyPlot
using ScikitLearn # machine learning package
using StatsBase
using Random
using LaTeXStrings # for L"$x$" to work instead of needing to do "\$x\$"
using Printf

# (optional)change settings for all plots at once, e.g. font size
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 16

# (optional) change the style. see styles here: https://matplotlib.org/3.1.1/gallery/style_sheets/style_sheets_reference.html
PyPlot.matplotlib.style.use("seaborn-white")   
```

## classifying breast tumors as malignant or benign

source: [UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic))

> Features are computed from a digitized image of a fine needle aspirate (FNA) of a breast mass. They describe characteristics of the cell nuclei present in the image.

The mean radius and smoothness of the cell nuclei (the two features) and the outcome (M = malignant, B = benign) of the tumor are in the `breast_cancer_data.csv`.


```julia
df = CSV.read("breast_cancer_data.csv")
df[!, :class] = map(row -> row == "B" ? 0 : 1, df[:, :outcome])
first(df, 5)
```




<table class="data-frame"><thead><tr><th></th><th>mean_radius</th><th>mean_smoothness</th><th>outcome</th><th>class</th></tr><tr><th></th><th>Float64</th><th>Float64</th><th>String</th><th>Int64</th></tr></thead><tbody><p>5 rows × 4 columns</p><tr><th>1</th><td>13.85</td><td>1.495</td><td>B</td><td>0</td></tr><tr><th>2</th><td>9.668</td><td>2.275</td><td>B</td><td>0</td></tr><tr><th>3</th><td>9.295</td><td>2.388</td><td>B</td><td>0</td></tr><tr><th>4</th><td>19.69</td><td>4.585</td><td>M</td><td>1</td></tr><tr><th>5</th><td>9.755</td><td>1.243</td><td>B</td><td>0</td></tr></tbody></table>



## visualize the two classes distributed in feature space


```julia
markers = Dict("M" => "x", "B" => "o")

figure()
xlabel("mean radius")
ylabel("mean smoothness")
for df_c in groupby(df, :outcome)
    outcome = df_c[1, :outcome]
    scatter(df_c[:, :mean_radius], df_c[:, :mean_smoothness], label="$outcome", 
        marker=markers[outcome])
end
legend()
axis("equal")
```

## get data ready for classifiation in scikitlearn

scikitlearn takes as input:
* a feature matrix `X`, which must be `n_samples` by `n_features`
* a target vector `y`, which must be `n_samples` long (of course)


```julia
n_tumors = nrow(df)

X = zeros(n_tumors, 2)
y = zeros(n_tumors)
for (i, tumor) in enumerate(eachrow(df))
    X[i, 1] = tumor[:mean_radius]
    X[i, 2] = tumor[:mean_smoothness]
    y[i] = tumor[:class]
end
X # look at y too!
```




    300×2 Array{Float64,2}:
     13.85   1.495
      9.668  2.275
      9.295  2.388
     19.69   4.585
      9.755  1.243
     16.11   4.533
     14.78   2.45 
     15.78   3.598
     15.71   1.972
     14.68   3.195
     13.71   3.856
     21.09   4.414
     11.31   1.831
      ⋮           
     11.08   1.719
     18.94   5.486
     15.32   4.061
     14.25   5.373
     20.6    5.772
      8.671  1.435
     11.64   2.155
     12.06   1.171
     13.88   1.709
     14.9    3.466
     19.59   2.916
     14.81   1.677



## logistic regression

let $\mathbf{x} \in \mathbb{R}^2$ be the feature vector describing a tumor. let $T$ be the random variable that denotes whether the tumor is benign (0) or malignant (1). the logistic model is a probabilistic model for the probability that a tumor is malignant given its feature vector:

\begin{equation}
    \log \frac{Pr(T=1 | \mathbf{x})}{1-Pr(T=1 | \mathbf{x})} = \beta_0 + \boldsymbol \beta^\intercal \mathbf{x}
\end{equation}
where $\beta_0$ is the intercept and $\boldsymbol \beta \in \mathbb{R}$ are the weights for the features. 

we will use scikitlearn to learn the  $\beta_0$ and $\boldsymbol \beta$ that maximize the likelihood.


```julia
@sk_import linear_model : LogisticRegression
```

    WARNING: redefining constant LogisticRegression





    PyObject <class 'sklearn.linear_model.logistic.LogisticRegression'>




```julia
lr = LogisticRegression(penalty="none", solver="newton-cg")

lr.fit(X, y)

println("β = ", lr.coef_)
println("β₀ = ", lr.intercept_)
```

    β = [1.1686607782175527 0.9420681231447378]
    β₀ = [-19.387890643955814]


prediction of the probability that a new tumor is 0 (benign) or 1 (malignant)


```julia
lr.predict_proba([10.0 2.5])
```




    1×2 Array{Float64,2}:
     0.995256  0.00474403



## visualize the learned model $Pr(T=1|\mathbf{x})$


```julia
radius = 5:0.25:30
smoothness = 0.0:0.25:20.0

lr_prediction = zeros(length(smoothness), length(radius))
for i = 1:length(radius)
    for j = 1:length(smoothness)
        x = [radius[i] smoothness[j]]
        lr_prediction[j, i] = lr.predict_proba(x)[2] # proba for probability
    end
end
```


```julia
figure()
pcolor(radius, smoothness, lr_prediction, cmap="viridis", 
    vmin=0.0, vmax=1.0) # vmin,vmax enforce colorbar scale
colorbar()

for df_c in groupby(df, :outcome)
    outcome = df_c[1, :outcome]
    scatter(df_c[:, :mean_radius], df_c[:, :mean_smoothness], label="$outcome", 
        marker=markers[outcome])
end
legend()
```

## making decisions: the ROC curve

this depends on the cost of a false positive versus false negative. (here, "positive" is defined as testing positive for "malignant")

> "I equally value minimizing (1) false positives and (2) false negatives."

$\implies$ choose $Pr(T=1|\mathbf{x})=0.5$ as the decision boundary.

> "I'd rather predict that a benign tumor is malignant (false positive) than predict that a malignant tumor is benign (false negative)."

$\implies$ choose $Pr(T=1|\mathbf{x})=0.2$ as the decision boundary. Even if there is a relatively small chance that the tumor is malignant, we still take action and classify it as malignant...

the receiver operator characteristic (ROC) curve is a way we can evaluate a classification algorithm without imposing our values and specifying where the decision boundary should be.


```julia
@sk_import metrics : roc_curve
@sk_import metrics : auc
```

    WARNING: redefining constant roc_curve





    PyObject <function auc at 0x7fe80b1fb7b8>




```julia
scores = lr.predict_proba(X)[:, 2]
```




    300-element Array{Float64,1}:
     0.14263839860751615  
     0.0026092670076477047
     0.0018782588829001088
     0.9996447815854005   
     0.0010942250321098578
     0.9760986635528334   
     0.5480964640686227   
     0.9200581620768566   
     0.6962552287733829   
     0.6852396841339906   
     0.5663718098629624   
     0.9999187138695917   
     0.011596238415532208 
     ⋮                    
     0.008004504524685206 
     0.9996348109797668   
     0.9122747253539892   
     0.9111094723703275   
     0.9999599017831867   
     0.0003696570017181287
     0.022876063860596256 
     0.014910308332017176 
     0.17409413839661506  
     0.7842086355766936   
     0.9980794947953938   
     0.37749924416866815  




```julia
fpr, tpr, thresholds = roc_curve(y, scores)
```




    ([0.0, 0.0, 0.0, 0.005555555555555556, 0.005555555555555556, 0.011111111111111112, 0.011111111111111112, 0.027777777777777776, 0.027777777777777776, 0.03333333333333333  …  0.37222222222222223, 0.43333333333333335, 0.43333333333333335, 0.4722222222222222, 0.4722222222222222, 0.5055555555555555, 0.5055555555555555, 0.7666666666666667, 0.7666666666666667, 1.0], [0.0, 0.008333333333333333, 0.7416666666666667, 0.7416666666666667, 0.7583333333333333, 0.7583333333333333, 0.7666666666666667, 0.7666666666666667, 0.8083333333333333, 0.8083333333333333  …  0.9583333333333334, 0.9583333333333334, 0.975, 0.975, 0.9833333333333333, 0.9833333333333333, 0.9916666666666667, 0.9916666666666667, 1.0, 1.0], [1.9999999999999254, 0.9999999999999254, 0.7959831048160733, 0.7882403199356263, 0.7732403837747389, 0.7622568750185913, 0.7614526612308641, 0.6962552287733829, 0.6282401802701606, 0.5888691371793174  …  0.09889839902698462, 0.065245657335008, 0.0601149551828247, 0.05261077303643403, 0.0518823253258252, 0.043527327155129115, 0.04345926758794682, 0.008416240229296301, 0.008004504524685206, 0.00011874278075119236])




```julia
figure()
plot([0, 1], [0, 1], color="k")
plot(fpr, tpr, color="darkorange",
         label=@sprintf("ROC curve (area = %0.2f)", auc(fpr, tpr))
    )
legend()
xlabel("false positive rate")
ylabel("true positive rate")
```

tradeoff:
* threshold too small: classify all of the tumors as malignant, false positive rate very high
* threshold too large: classify all of the tumors as benign, false negative rate very high

somewhere in the middle (but still depending on the cost of a false positive versus false negative) is where we should operate.

the `auc`, area under the curve, has a probabilistic interpretation:
>  the area under the curve (often referred to as simply the AUC) is equal to the probability that a classifier will rank a randomly chosen positive instance higher than a randomly chosen negative one (assuming 'positive' ranks higher than 'negative') -[Wikipedia](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)

**warning**: always split your data into test or train or do cross-validation to assess model performance. we trained on all data here to see the mechanics of fitting a logistic regression model to data, visualizing the model, and creating an ROC curve.

ec08332d2b3304ba2eae5df9a5ca8baa28969650
ipynb
Jupyter Notebook
CHE599-IntroDataScience/lectures/logistic_regression/logistic regression.ipynb
leanth/OSUCoursework
ccfbf5f9daa8f6d3818bb5e4cc8df7c5135a5f34
MIT
CHE599-IntroDataScience/lectures/logistic_regression/logistic regression.ipynb
leanth/OSUCoursework
ccfbf5f9daa8f6d3818bb5e4cc8df7c5135a5f34
MIT
CHE599-IntroDataScience/lectures/logistic_regression/logistic regression.ipynb
leanth/OSUCoursework
ccfbf5f9daa8f6d3818bb5e4cc8df7c5135a5f34
MIT
Qwen/Qwen-72B
1. YES
2. YES