It's been some time since this question has been asked, but allow me to post a new answer. Its piece of my thesis, where I've put the derivation (since its so hard to find). I hope it'll help anyone, who encounters this problem. It's mostly based on the introduction of the paper "The decay of quantum systems with a small number of open channels" by F. M. Dittes (http://www.sciencedirect.com/science/article/pii/S037015730000065X).
Consider a system described by a Hamiltonian $H=H_0+V$, where $V$ describes the perturbation upon which the scattering one wishes to study -- in our case it will be the coupling between the leads and the device. The system set in the state $|i\rangle$ at the initial time $t=-\infty$ is subject to evolution governed by the Schrödinger equation in the interaction picture. At the final time $t=+\infty$ its state is given by $\mathcal{T}\mathrm{e}^{-\mathrm{i} \int_{-\infty}^{\infty} \mathrm{d} t\,V_{\rm I}(t)}|i\rangle = U|i\rangle$, where $V_{\rm I}(t) = \mathrm{e}^{\mathrm{i} H_0 t}V\mathrm{e}^{-\mathrm{i} H_0 t}$ is the perturbation in the interaction picture, and $\mathcal{T}$ denotes time-ordering. The scattering amplitude to some state $|f\rangle$ is given by the projection of the final state:
\begin{multline}
\langle f|U|i\rangle = \langle f | i \rangle+\\+ \sum_{j_1,j_2,\dots,j_{n-1}} \sum_{n=1}^\infty (-\mathrm{i})^n
\int_{-\infty}^\infty\mathrm{d} {t_1}
\int_{-\infty}^{t_1}\mathrm{d} t_2 \dots
\int_{-\infty}^{t_{n-1}}\mathrm{d} t_n
\langle f|V_{\rm I}(t_1)|j_1\rangle\langle j_1|
V_{\rm I}(t_2)|j_2\rangle\dots \langle j_{n-1}|V_{\rm I}(t_n)|i\rangle,
\end{multline}
where $j_k$-s are summed over the eigenbasis of $H_0$. Performing the last integral:
\begin{align}
\int_{-\infty}^{t_{n-1}}\mathrm{d} t_n \langle j_{n-1}|V_{\rm I}(t_n)|i\rangle &=
\int_{-\infty}^{t_{n-1}}\mathrm{d} t_n
\mathrm{e}^{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+) t_n}
\langle j_{n-1}| V|i\rangle =
\frac{\mathrm{e}^{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+) t_{n-1}}}
{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+)}
\langle j_{n-1}| V|i\rangle.
\end{align}
Performing the second last integral:
\begin{multline}
\int_{-\infty}^{t_{n-2}}\mathrm{d} t_{n-1} \langle j_{n-2}| V_{\rm I}(t_{n-1})|j_{n-1}\rangle
\frac{\mathrm{e}^{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+) t_{n-1}}}
{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+)}
\langle j_{n-1}| V|i\rangle =\hfill\\
\begin{aligned}
&=\int_{-\infty}^{t_{n-2}}\mathrm{d} t_{n-1}
\mathrm{e}^{-\mathrm{i} (E_i-E_{j_{n-2}} + \mathrm{i} 0^+) t_{n-1}}
\langle j_{n-2}| V|j_{n-1}\rangle
\frac{1}{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+)}
\langle j_{n-1}| V|i\rangle\\
&=\frac{\mathrm{e}^{-\mathrm{i} (E_i-E_{j_{n-2}} + \mathrm{i} 0^+) t_{n-2}}}
{-\mathrm{i} (E_i-E_{j_{n-2}} + \mathrm{i} 0^+)}
\langle j_{n-2}| V|j_{n-1}\rangle
\frac{1}{-\mathrm{i} (E_i-E_{j_{n-1}} + \mathrm{i} 0^+)}
\langle j_{n-1}| V|i\rangle.
\end{aligned}
\end{multline}
It is straightforward to see that the remaining integrals (except the first one) will follow identical pattern. Summing over $j_k$-s one can contract the notation: $\sum_{j_k}|j_k\rangle\langle j_k|(E_i-E_{j_k}+\mathrm{i} 0^+)^-1 = (E_i-H_0+\mathrm{i} 0^+)^{-1} \equiv G(E_i+\mathrm{i}0^+)$, which yields:
\begin{align}
\langle f|U|i\rangle &= \langle f | i \rangle + \sum_{n=1}^\infty \frac{(-\mathrm{i})^n}{(-\mathrm{i})^{n-1}}
\int_{-\infty}^\infty\mathrm{d} t_1 \mathrm{e}^{-\mathrm{i}(E_i-E_f)t_1}
\langle f|V(GV)^{n-1}|i\rangle \\
&= \langle f | i \rangle -2\pi\mathrm{i}\delta(E_i-E_f)
\langle f|V\sum_{n=0}^\infty(GV)^n|i\rangle.
\label{eq:mwe1}
\end{align}
Note that in the case of our interest $G$ acts separately in the Hilbert spaces of the leads ($\mathscr{H}_{\rm L}$) and the device ($\mathscr{H}_{\rm D}$), whereas $V$ mediates between the two. Since both $|i\rangle$ and $|f\rangle$ belong to the $\mathscr{H}_{\rm L}$, only terms with even number of $V$ operators (including the one in front of the sum) in the equation above will contribute to the scattering amplitude. Writing explicitly $G=G_{\rm L} + G_{\rm D}$ and $V = V_{\rm LD} + V_{\rm DL}$, where $G_{\rm L}$ ($G_{\rm D}$) acts within $\mathscr{H}_{\rm L}$ ($\mathscr{H}_{\rm D}$) space and $V_{DL}$ ($V_{LD}$) acts from $\mathscr{H}_{\rm L}$ to $\mathscr{H}_{\rm D}$ (from $\mathscr{H}_{\rm D}$ to $\mathscr{H}_{\rm L}$) one gets:
\begin{align}
\langle f|U|i\rangle &= \langle f | i \rangle -2\pi\mathrm{i}\delta(E_i-E_f)
\langle f|V_{\rm LD}\sum_{n=0}^\infty(G_{\rm D}V_{\rm DL}G_{\rm L}V_{\rm LD})^n G_{\rm D}V_{\rm DL}|i\rangle \\
&= \langle f | i \rangle -2\pi\mathrm{i}\delta(E_i-E_f)
\langle f|V_{\rm LD}\frac{1}{1-G_{\rm D}V_{\rm DL}G_{\rm L}V_{\rm LD}} G_{\rm D}V_{\rm DL}|i\rangle\\
&= \langle f | i \rangle -2\pi\mathrm{i}\delta(E_i-E_f)
\langle f|V_{\rm LD}\frac{1}{G_{\rm D}^{-1}-V_{\rm DL}G_{\rm L}V_{\rm LD}} V_{\rm DL}|i\rangle.
\end{align}
The $S$-matrix by summing over $|f\rangle$ and $|i\rangle$. Note that $\sum_{f} \delta(E_i-E_f)|f\rangle\langle f| = \sum_{f:E_f=E_i} \rho(E_i)|f\rangle\langle f| = \rho(E_i)\Pi(E_i)$, where $\rho(E_i)$ is the density of states at energy $E_i$ and $\Pi(E_i)$ is the projector onto the eigensubspace of $\mathscr{H}_{\rm L}$ of energy $E_i$. Then:
\begin{align}
S(E) &= \sum_{\substack{f,i:\\E_i=E}} |f\rangle\langle f|U|i\rangle\langle i| = \Pi(E) -2\pi\mathrm{i}\rho(E)\Pi(E)V_{\rm LD}\frac{1}{E-H_{\rm D}-V_{\rm DL}G_{\rm L}V_{\rm LD}} V_{\rm DL}\Pi(E).
\end{align}
Furthermore, the self energy $V_{\rm DL}G_{\rm L}V_{\rm LD}$ can be expressed in terms of its Hermitian and antihermitian parts:
\begin{align}
V_{\rm DL}G_{\rm L}V_{\rm LD} &= \int\mathrm{d} E'\,
V_{\rm DL}\frac{\rho(E')\Pi(E')}{E-E'+\mathrm{i} 0^+}V_{\rm LD} \\
&= \int\mathrm{d} E'\,
V_{\rm DL}\rho(E')\Pi(E')V_{\rm LD}
\left(\mathcal{P}\frac{1}{E-E'} - \mathrm{i}\pi\delta(E-E')\right) \\
&= \mathcal{P}\int\mathrm{d} E'\,
\frac{W(E')W^\dagger(E')}{E-E'}
-\mathrm{i}\pi W(E)W^\dagger(E),
\end{align}
where $W(E)\equiv \sqrt{\rho(E)} V_{\rm DL}\Pi(E)$ and $\mathcal{P}$ denotes the principal value of the integral. The $S$ matrix finally reads:
\begin{equation}
S(E) = 1 - 2\pi\mathrm{i} W^\dagger(E)\frac{1}{E-H_{\rm D}-\mathcal{P}\int\mathrm{d} E'\,
\frac{W(E')W^\dagger(E')}{E-E'}
+\mathrm{i}\pi W(E)W^\dagger(E)} W(E).
\end{equation}
where we replaced the first projector with identity, sice the $S$-matrix acts in the relevant eigenspace anyway. The effective Hamiltonian is given by $H_{\rm D}+\mathcal{P}\int\mathrm{d} E'\,\frac{W(E')W^\dagger(E')}{E-E'}$, but the second term is often omited if the coupling depends weakly on the energy.