How does one know in general that the toric fibers of a Kahler toric manifold are Lagrangian submanifolds? I can roughly fathom this for e.g. a circle in $\mathbb{C}P^1$, but how about more general toric manifolds of higher dimension?

By definition, a toric variety of complex dimension $n$ comes with an action of a complex torus $(\mathbb{C}^{*})^n$. I guess that by Kähler toric manifold, you mean a toric variety $X$ equipped with a Kähler form $\omega$ which is preserved by a compact subgroup $U(1)^n$ of $(\mathbb{C}^{*})^n$. The toric fibration $X \rightarrow \mathbb{R}^n$ is then the moment map of this action of $U(1)^n$ on $(X,\omega)$. More precisely, choosing Hamiltonian functions $H_1$,...,$H_n$ generating the action of $U(1)^n$ (i.e. such that the corresponding Hamiltonian vector fields generate the action of $U(1)^n$), then the toric fibration $X \rightarrow \mathbb{R}^n$ is given by $x \mapsto (H_1(x),...,H_n(x))$. Because $U(1)^n$ is abelian, the vector fields generating the action commute and so the Hamiltonian functions Poisson commute. It is a general fact that if one has a symplectic manifold of real dimension $2n$, with $n$ functions $H_1,...,H_n$ Poisson commuting (and with linearly independent differentials), then the common vanishing locus of these functions is a Lagrangian submanifold (for example, by Darboux lemma, you can locally complete these $n$ functions by $n$ other functions $\theta_1$,...,$\theta_n$ such that $H_i$, $\theta_j$ are local coorrdinates such that $\omega=\sum_{i=1}^n dH_i \wedge d\theta_i$ and the result is then obvious).