%% Chapter 12: Logistic Regression
% So far we used linear models to predict a continuous reponse variable y using a set of continuous or discrete predictors x_1, \ldots,x_n. Now we turn our attention to predicting a discrete reponse variable using both continuous and discrete predictors. There are many applications where the response variable is discrete, such as predicting 1.) if a team will win a sporting event, 2.) if a patient will have a heart attack, or 3.) if a medical device will fail.
% Let's assume the response variable y has two possible outcomes, 0 and 1. (This is not really a limitation, as response variable with multiple outcomes can be modeled as a series of binary outcomes.) Linear models are unbounded, real-valued functions. The responses predicted by linear models can therefore range from negative to positive infinity. There is no way to use only field algebra or linear algebra to fit a function that predicts only two events. Instead, we use a link function to translate the continuous output of a linear model into a binary prediction.
% Rather than predict the binary response, we could try predicting the probability that the response is equal to 1. Probabilities are continous, but they are bounded to the interval [0,1]. We don't have a method to force a linear models to only make predictions in this range. We could instead predict the odds of the response variable, which is the probability that the response equals 1 divided by the probability that the response equals 0.
% \mathrm{odds}(y) = \frac{P(y=1)}{P(y=0)}
% If the \mathrm{odds}(y)=2, then the probability that y equals 1 is twice as large as the probability that y equals 0. (Odds are usually expressed as a proportion, so an odds of 2 is written as 2:1, or "two to one".) We can convert between probabilities and odds by remembering that probabilities sum to 1, or P(y=0)+P(y=1)=1. Then
% \mathrm{odds}(y) = \frac{P(y=1)}{P(y=0)}=\frac{P(y=1)}{1-P(y=1)}\Rightarrow P(y=1)=\frac{\mathrm{odds}(y)}{1+\mathrm{odds}(y)}
% The odds of a discrete variable are always nonnegative and tend toward infinity as P(y=0) approaches 0. The odds function is therefore unbounded, but only in the positive direction.

axargs = {'Interpreter','latex','FontSize',24};
pargs = {'LineWidth',4};
P = linspace(0,1);
plot(P,P./(1-P),pargs{:})
xlabel('$$P(y=1)$$',axargs{:})
ylabel('$$\mathrm{odds}(y)$$',axargs{:})

% The natural logarithm of the odds, however, is a continuous variable. As shown below, the "log odds" is a sigmoid function with horizontal asymptotes at 0 and 1.

P = linspace(1e-4,1-1e-4,100);
plot(log(P./(1-P)),P,pargs{:})
ylabel('$$P(y=1)$$',axargs{:})
xlabel('$$\log\left(\mathrm{odds}(y)\right)$$',axargs{:})

% This is exactly the type of function we're looking for. We can interpret the output of a linear model as the log of the odds of the response variable. Our linear model becomes
% \log\left(\mathrm{odds}(y)\right) = \beta_0 + \beta_1x_1 + \cdots+\beta_nx_n
% Looking at the log odds function above, we see that it transitions from 0 to 1 when \log\left(\mathrm{odds}(y)\right)=0. Thus if our continuous linear model predicts a negative value, we say the binary response variable is 0. If the linear model predicts a positive value, we say the binary response variable is 1. (If the linear model predicts exactly zero, we are equally sure that the response is 0 and 1.) Using a linear model to predict the log-odds of a discrete variable is called logistic regression. The function L(y) = \log\left(\mathrm{odds}(y)\right) is called the logit function. Because it links the response variable to the linear models, we refer to the logit (and other similar functions) as link functions.
%% Example: Predicting risk of Huntington's Disease
% Huntington's Disease is an inherited genetic condition caused by repeated CAG sequences in the Huntingtin (HTT) gene. Too many CAG repeats create a "glutamine knot" in the protein, causing toxic protein aggregates in neurons. Symptoms of Huntington's appear later life, and an individual's risk for developing the disease correlates with the number of CAG repeats.
%
% Let's build a model to predict the probability of developing Huntington's based on the number of CAG repeats. The response variable is binary (Huntington's disease or not) and the predictor variable is continuous (the number of CAG repeats in the HTT gene). To train the model, we counted the number of CAG repeats in 50 individuals with and without the disease.

load huntington.mat
scatter(hunt.CAGs,categorical(hunt.disease))
xlabel('CAG repeats',axargs{:})
ylabel("Huntington's disease",axargs{:})

% We see from these data that predicting disease status with low (<25) or high (>35) CAG repeats is straightforward. However, there is a region between 25 and 35 CAG repeats where disease status can be ambiguous. Let's build a logistic regression model to predict Huntington's status. We use the Matlab function fitglm, for "fit generalized linear model". The fitglm function is similar to fitlm; the first argument is a table of data, and the second argument is a formula describing the model. However, fitglm can use a wide range of link functions and datatypes when fitting linear models. For logistic regression using binary responses we need to specify the logit link function and a binomial distribution.

model = fitglm(hunt,'disease ~ CAGs','link','logit','Distribution','binomial')

% Remember that the model we're fitting is
% \log\left(\mathrm{odds}(y)\right) = \beta_0 + \beta_1x_1 + \cdots+\beta_nx_n
% so the parameters are not directly interpretable as probabilities. We can rearrange the model to find a formula for the odds of the response:
% \mathrm{odds}(y) = e^{\beta_0 + \beta_1x_1 + \cdots+\beta_nx_n}
% Since the odds of y is equal to P(y=1)/\left(1-P(y=1)\right), we can also solve for the probability that y=1.
% P(y=1) = \frac{1}{1+e^{-t}}, \quad t=\beta_0 + \beta_1x_1 + \cdots+\beta_nx_n
% For our Huntington's example
% P(\mathrm{disease}) = \frac{1}{1+e^{-(\beta_0+\beta_1[\mathrm{CAGs}])}} = \frac{1}{1+e^{14.03-0.51[\mathrm{CAGs}]}}

scatter(hunt.CAGs,hunt.disease)
hold on
cag_range = linspace(5,50,100);
beta = model.Coefficients.Estimate;
plot(cag_range, 1./(1+exp(-(beta(1)+beta(2)*cag_range))))
hold off
xlabel('CAG repeats',axargs{:});
ylabel('$$P(\mathrm{disease})$$',axargs{:});

% We are often interested in when P(\mathrm{disease}) = 1/2, as this is the threshold number of CAG repeats where a person is equally likely to have or not have Huntington's. We could solve the above equation, or we could recall that the logistic function reaches its midpoint when the linear model moves from negative to positive. Thus we can simply solve for when \beta_0+\beta_1[\mathrm{CAGs}]=0.
% -14.03+0.51[\mathrm{CAGs}]=0 \Rightarrow [\mathrm{CAGs}] = 14.03/0.51 \approx 28 CAG repeats
%% Interpreting Individual Coefficients
% The entire linear model (\beta_0 + \beta_1x_1 + \cdots+\beta_nx_n) is useful for predicting the probability of the reponse variable. To interpret individual coefficients we can calculate the odds ratio, or the change in odds per unit change in a single predictor variable. For our Huntington's example, let's calculate the odds ratio for adding one more CAG repeat.
% odds ratio([CAGs]) = \frac{\mathrm{odds}([\mathrm{CAGs}]+1)}{\mathrm{odds}([\mathrm{CAGs}])} = \frac{e^{\beta_0+\beta_1([\mathrm{CAGs}]+1)}}{e^{\beta_0+\beta_1[\mathrm{CAGs}]}} = \frac{e^{\beta_0}e^{\beta_1[\mathrm{CAGs}]}e^{\beta_1}}{e^{\beta_0}e^{\beta_1[\mathrm{CAGs}]}} = e^{\beta_1}
% Since \beta_1=0.51, having one more CAG repeat increases the odds of developing Huntington's by e^{0.51}=1.67 fold. For any logistic regression model, the odds ratio for the ith predictor variable is the exponential of the ith coefficient:
% \mathrm{odds}\,\,\mathrm{ratio}(x_i)=\frac{\mathrm{odds}(x_i+1)}{\mathrm{odds}(x_i)}=e^{\beta_i}