# Why does Quadratic (Normal/Laplace) Approximation fail on multilevel models?

In Statistical Rethinking, 2nd Edition, section 13.1, Richard McElreath says:

Why doesn’t simple quadratic approximation, using for example quap, work with multilevel models? When a prior is itself a function of parameters, there are two levels of uncertainty. This means that the probability of the data, conditional on the parameters, must average over each level. Ordinary quadratic approximation cannot handle the averaging in the likelihood, because in general it’s not possible to derive an analytical solution. That means there is no unified function for calculating the log-posterior. So your computer cannot directly find its minimum (the maximum of the posterior).

Corresponding rethinking code:

m12.2 <- ulam(
alist(
surv ~ dbinom(density, p),
logit(p) <- a_tank[tank],
a_tank[tank] ~ dnorm(a, sigma),
a ~ dnorm(0, 1),
sigma ~ dcauchy(0, 1)
),
data=d, iter=4000, chains=4
)


Corresponding Stan code:

data{
int<lower=1> N;
int<lower=1> N_tank;
int surv[N];
int density[N];
int tank[N];
}
parameters{
vector[N_tank] a_tank;
real a;
real<lower=0> sigma;
}
model{
vector[N] p;
sigma ~ cauchy( 0 , 1 );
a ~ normal( 0 , 1 );
a_tank ~ normal( a , sigma );
for ( i in 1:N ) {
p[i] = a_tank[tank[i]];
}
surv ~ binomial_logit( density , p );
}
generated quantities{
vector[N] p;
for ( i in 1:N ) {
p[i] = a_tank[tank[i]];
}
}


I can not figure it out how "prior is itself a function of parameters" could be a problem. It seems that they're just combined into a prior, so there is nothing differ from a single-level model?

Here's my derivation of laplace approximation, maybe I'm wrong in some steps so I can't understand Richard's statement? Please correct me.

From Bayesian Data Analysis, 3rd edition, P84.

\begin{align*} p(\theta|y) \approx N(\hat{\theta}, (I(\hat{\theta}))^{-1}) \\ I(\theta) = -\frac{d^2}{d\theta^2} \log p(\theta|y) \end{align*}

So

$$I(\theta) = -\frac{d^2}{d\theta^2} \log p(\theta|y) = -\frac{d^2}{d\theta^2} (\log p(\theta,y) - \log p(y)) = -\frac{d^2}{d\theta^2} \log p(\theta,y)$$

As the code just specified joint distribution $$\log p(\theta,y)$$, like:

function logp_theta_y(){
target = 0
target += dcauchy(sigma, 0, 1, log=TRUE);
target += dnorm(a, 0, 1, log=TRUE);
target += dnorm(a_tank, a, sigma, log=TRUE);
for ( i in 1:N ) {
target += dbinom(surv[i], density[i], logistic(a_tank[tank[i]]), log=TRUE);
}
return target;
}


In fact, I run following Julia code (Sorry, I'm not that familar with R) to get Laplace approximation and found that they're close to each other:

using CSV
using DataFrames
using Distributions
using ForwardDiff
using Optim
using LinearAlgebra

df = DataFrame(CSV.File("reedfrogs.csv")) # https://github.com/rmcelreath/rethinking/blob/master/data/reedfrogs.csv

S = df.surv
N = df.density

N_tank = length(S)

logistic(x) = 1/(1+exp(-x))

function logp_y_p(par)
a_tank = par[1:N_tank]
a = par[N_tank+1]
sigma = exp(par[N_tank+2])

target = 0
target += logpdf(Cauchy(0, 1), sigma) + par[N_tank+2]
target += logpdf(Normal(0, 1), a)
target += logpdf.(Normal.(a, sigma), a_tank) |> sum
target += logpdf.(Binomial.(N, logistic.(a_tank)), S) |> sum

target
end

f(x) = -logp_y_p(x)
x0 = randn(50)

res = optimize(f, x0, LBFGS(); autodiff = :forward)

par_est = Optim.minimizer(res)
FI = -ForwardDiff.hessian(logp_y_p, par_est)
Sigma = inv(FI)

function normal_to_lognormal(mu, sigma)
exp(mu+sigma^2/2), sqrt((exp(sigma^2)-1) * exp(2*mu+sigma^2))
end

mu_sigma, sigma_sigma = normal_to_lognormal(par_est[end], sqrt(diag(Sigma)[end]))

mu_vec = copy(par_est)
mu_vec[end] = mu_sigma

sigma_vec = sqrt.(diag(Sigma))
sigma_vec[end] = sigma_sigma

[mu_vec sigma_vec]