How to derive formula for marginal probability of choosing nest in nested logit model?
As it turned out, my previous logic was wrong. Here is how it should be done.

Marginal probability of choosing nest $k$ is
$$P_{nB_k} = P\left[\max_{j\in B_k} U_{njk} \geq \max_{j\in B_s} U_{njs}, \forall s\neq k \right]\\
= P\left[W_{nk}+\epsilon_{nk}+\max_{j\in B_k}(Y_{nj}+e_{nj}) \geq W_{ns}+\epsilon_{ns}+\max_{j\in B_s}(Y_{nj}+e_{nj}), \forall s\neq k \right]$$

Then as $e_{nj}$ is iid $Gumbel(0,\lambda_t)$, $\max_{j\in B_k}(Y_{nj}+e_{nj})$ is iid Gumbel with location parameter $\lambda_kI_{nk}$ and scale parameter $\lambda_k$. The Gumbel distribution is preserved over linear transformations so 
$$ \max_{j\in B_k}(Y_{nj}+e_{nj})=\lambda_kI_{nk}+\xi_{nk}$$
where $\xi_{nk}$ is iid $Gumbel(0,\lambda_t)$. Substituting back to marginal probability, get
$$P_{nB_k} = P\left[(\epsilon_{nk}+\xi_{nk})-(\epsilon_{nk}+\xi_{nk}) \geq (W_{nk}+\lambda_kI_{nk})-(W_{ns}+\lambda_sI_{ns}), \forall s\neq k \right]\\
=\frac{e^{W_{nk}+\lambda_kI_{nk}}}{\sum_{s=1}^Ke^{W_{ns}+\lambda_sI_{ns}}} $$

the last equality follows from the fact that $\epsilon_{nk}+\xi_{nk}$ is iid $Gumbel(0,1)$ by assumption on $\epsilon_{nk}$.