How to explain covariance in logistic regression + analogy to linear regression Introduction
Linear model
In linear regression we predict continuous variable $Y \in R^n$ with use of $n \times p$ deterministic plan matrix $X$ and theoretical model (let's ignore intercept coefficient $\beta_0$ for simplicity)
$$Y_i = \beta_1 X_1 + ... \beta_p X_p + \epsilon_i,$$
where random error $\epsilon_i \sim N(0,\sigma^2).$ The predictor of $Y_i$ is given by
$$\hat Y_i = E(Y_i) = \hat\beta_1 X_1 + ... \hat\beta_p X_p.$$
Because $Y$ is a linear transformation of $\epsilon$ it has distribution $N(\beta X, \sigma^2)$. We know that OLS estimator $\hat \beta$ is given by
$$\hat \beta = (X^T X)^{-1} X^T Y$$
and it can be shown (and I understand the derivation) that $\hat \beta \sim N(\beta, \sigma^2 (X^TX)^{-1})$.
Logistic model
In logistic regression we predict the probability
$\mu_i := P(Y_i = 1)$
of binary variable $Y \in \{0,1\}^n$ with the same matrix $X$ and model
$$\log\left(\frac{\mu_i}{1-\mu_i}\right) = \beta_1 X_1 + ... \beta_p X_p$$
(in my lecture slides there is no $\epsilon_i$ added for some reason, so let's assume that there's no $\epsilon_i$.) It is known that there is no closed form of $\hat \beta$ in this case, instead we need to use some optimization algorithms like Newton's method to find it.
The problem
From my lecture slides $\hat \beta$ has asymptiotic distribution $$N(\beta, (X^T S X)^{-1}),$$
where $S_{i,i} = \mu_i(1 - \mu_i)$ for certain $\hat \beta$ and as far as I understand, it is simply a covariance matrix of $Y_i$ (is it true?). Similarly in linear model $I_{n \times n} \sigma^2$ was covariance matrix of $Y_i$ (was it?). We are not given the explanation though. The slides just point out that there is some kind of analogy between linear and logistic regression estimator distributions.

*

*I want to fully understand the analogy and get  to know where does the logistic distribution come from (I know it is based on some theorem about asymptotic estimators, but it's all I know).

*What bothers me is that in linear model the covariance matrix of $Y$ is not inverted, while in logistic it is ($S^{-1}$). What is wrong with that?

 A: $\newcommand{\E}{\operatorname{E}}$
$\newcommand{\var}{\operatorname{Var}}$
$\newcommand{\D}{\operatorname{D}\!}$
Unlike in ordinary least squares, the estimator $\widehat{\beta}$ for logistic regression is found via maximum likelihood estimation. Under certain regularity conditions, it can be shown that the MLE $\widehat{\beta}$ is asymptotically normal, with
\begin{equation*}
\widehat{\beta} \overset{d}{\longrightarrow} N(\beta, \mathcal{I}_T^{-1}).
\end{equation*}
Above, $\mathcal{I}_T$ is the total Fisher information matrix
\begin{equation*}
\mathcal{I}_T(\beta) = \E\left[-\frac{\partial^2}{\partial\beta\partial\beta^\top}\ell(\beta \,|\, Y)\right]
\end{equation*}
where $\ell$ is the log-likelihood of the data.
Let $\mu := \E[Y]$. Then it can be shown that for the generalized linear model $g(\mu) = X\beta$, the total Fisher information matrix is given by
\begin{equation*}
\mathcal{I}_T(\beta) = (\D\mu)^\top\var[Y]^{-1}\D\mu
\end{equation*}
where $\D\mu$ is the Jacobian matrix of $\mu$ with respect to $\beta$.
To illustrate this, we may consider the ordinary least squares case $Y \sim N(0, \sigma^2I)$. Here, $\mu = X\beta$, so $\D\mu = X$. Furthermore, we have $\var[Y] = \sigma^2I$. Hence,
\begin{equation*}
\mathcal{I}_T^{-1}(\beta) = \left[X^\top(\sigma^2I)^{-1}X\right]^{-1} = \sigma^2(X^\top X)^{-1}
\end{equation*}
That is to say that the asymptotic variance of $\widehat{\beta}$ in OLS is $\sigma^2(X^\top X)^{-1}$. Indeed, we know that this variance actually happens to be exact.
In logistic regression, we proceed similarly. We have that each $Y_i \sim \operatorname{Bernoulli}(\mu_i)$. Hence, $\var[Y_i] = \mu_i(1-\mu_i)$. Our model is $\operatorname{logit}(\mu) = X\beta$. Calculating the derivative, we find that
\begin{equation*}
\D\mu = \begin{bmatrix}
\mu_1(1-\mu_1) & 0 & \cdots & 0 \\
0 & \mu_2(1-\mu_2) & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \mu_n(1-\mu_n)
\end{bmatrix}X = \var[Y]X
\end{equation*}
It follows that
\begin{equation*}
\mathcal{I}_T^{-1}(\beta) = \left(X^\top\var[Y] \cdot \var[Y]^{-1} \cdot \var[Y]X\right)^{-1} = (X^\top \var[Y]X)^{-1}
\end{equation*}
as was stated in your class notes.
