Why do we use GLM? I'm trying to justify using of GLM model in my project instead of a simple linear regression. A lot of sources that I've seen contain the statement that "GLM allow us to build regression models for response variables that are not normally distributed". I can't understand what that is supposed to mean since we don't make any assumptions about variable distribution. Could you, please, explain  why it's often more appropriate (theoretically)  to use GLM if your response variable is,  say, count. 
 A: You do make assumptions about response variable distribution in normal linear models. Normal linear models can be written as
$$y_i = x_i'\beta + \varepsilon_i,\qquad \varepsilon_i\sim\mathcal{N}(0,\sigma^2)$$
but that means that $y_i$ is a linear transformation of a normally distributed random variable, thus its distribution is
$$y_i\sim\mathcal{N}(x_i'\beta,\sigma^2)$$
BTW, it is this normal distribution that makes $t$- and $F$-tests possible: $\hat\beta$ is normally distributed because it is a linear transformation of $y$, $\hat\beta=(X'X)^{-1}X'y$.
Further, that means that $y_i$ is continuous and can be positive, zero, or negative.
When you have a binary response variable, it can't be normally distributed. You are actually interested in the probability $P(y=1)$ vs $P(y=0)$. You could try a linear probability model, i.e.
$$P(y=1\mid x)=p(x)=x'\beta$$
where
$$\begin{align}E(y\mid x)&=1\cdot p(x)+0\cdot(1-p(x))=p(x)=x'\beta\\
V(y\mid x) &= p(x)(1-p(x))=x'\beta(1-x'\beta)\end{align}$$
but


*

*the error is not normally distributed, since it can just take two vales:
$$\varepsilon=y-x'\beta\quad\Rightarrow\quad
\varepsilon=\begin{cases} 1-x'\beta & \text{if } y=1\\-x'\beta & \text{if } y=0\end{cases}$$

*the error is not homoscedastic:
$$V(\varepsilon)=E(\varepsilon^2)=x'\beta(1-x'\beta)$$
thus the error variance depends on $x$

*$E(y\mid x)=p(x)$ is a probability, thus it must be greater than 0 and lesser than 1, but the linear predictor $E(y\mid x)=x'\beta$ can take any real value; you could try:
$$p(x)=\begin{cases} 0 & \text{if } x'\beta<0 \\
x'\beta & \text{if } 0\le x'\beta \le 1 \\
1 & \text{if } x'\beta > 1 \end{cases}$$
but the graph of $p(x)$ would be a broken line. How could you interpret it? What could its vertices mean? (BTW: interpretation is even more important than p-values: high p-values which make sense are better than low p-values which do not make sense.)


This is why one looks for a smooth function which maps $[0,1]$ to $\mathbb{R}$:
$$x'\beta=g(E(y\mid x))$$
and you get:
$$P(y=1\mid x)=G(x'\beta),\qquad G(x'\beta)=g^{-1}(x'\beta)$$
In logit models, 
$$g(x)=\ln\left(\frac{x}{1-x}\right),\qquad g^{-1}(x)=\frac{1}{1+e^{-x}}$$
where $g^{-1}$ is the cdf of a logistic random variable.
When you have to model count data, a natural (not always perfect) solution is:
$$P(y=k)=\frac{e^{-\theta}\theta^k}{k!}\quad\text{i.e.}\quad y\sim\text{Poisson}(\theta)$$
so that
$$E(y\mid x)=\theta=\exp(x'\beta)$$
thus $E(y\mid x) > 0$ even if $x'\beta < 0$. Other solutions are zero-inflated or negative-binomial models.
