Accounting for metric tensor derivatives in Einstein-Hilbert action I'm puzzling over the canonical derivation of GR from the Einstein-Hilbert action; getting the derivation to gel with an explicit treatment of the functional derivative isn't working out. So the derivation (drawn here from Wikipedia, though other literature is similar) begins,
$$
I = \int{\sqrt{-g}d^4x} \left[\frac{1}{2 \kappa}R + \mathcal{L}\right]
$$
and immediately proceeds to 
$$
\delta I = 0 =\int d^4x \delta g_{\mu \nu}\left[\frac{1}{2\kappa}\frac{\delta\left(\sqrt{-g}R\right)}{\delta g_{\mu \nu}}+\frac{\delta\left(\sqrt{-g}\mathcal{L}\right)}{\delta g_{\mu \nu}}\right].
$$
But the Ricci scalar depends on the first and second derivatives of the metric tensor, so why do we not have factors
$$
\delta g_{\mu \nu, \alpha}~, \qquad \delta g_{\mu \nu,\alpha \beta}~,
$$
against which we vary as well? Maybe there is some identity that in this case causes these terms to vanish, but I don't see it.
 A: You absolutely do have these terms.  Most people just always implicitly integrate by parts, and realistically, they hide these terms in other terms, because the algebra blows up into tons of terms very quickly.  A very laborious version of this is worked out in the Classical Field Theory book in the Landau and Lifschitz series.  
Alternately, you can use the Palantini form of the variation, and vary the Christoffel symbols, rather than the metric.
A: I) Even before one varies the action, recall that the Einstein-Hilbert (EH) Lagrangian density is 
$$\tag{1} {\cal L}_{EH}~\sim~\sqrt{-\det(g_{\cdot\cdot})} \left\{g^{\mu\nu}~ R_{\mu\nu}(\Gamma_{LC},\partial\Gamma_{LC})-2\Lambda\right\}  $$ 
where $\Gamma_{LC}$ refer to the Levi-Civita (LC) Christoffel symbols, which in turn depend on up to first-order derivatives $\partial g_{\cdot\cdot}$ of the metric $g_{\mu\nu}$. 
By inspection of eq. (1), we see that the EH Lagrangian density (1) is linear in second-order derivatives
$$\tag{2}  {\cal L}_{EH} ~=~ {\cal F}^{\mu\nu\lambda\sigma}(g_{\cdot\cdot})~\partial_{\mu} \partial_{\nu}g_{\lambda\sigma}   +  {\cal F}(g_{\cdot\cdot},\partial g_{\cdot\cdot}) , $$
which we can rewrite as a function $F(g_{\cdot\cdot},\partial g_{\cdot\cdot})$ which depends on up to first-order derivatives, plus a total divergence term: 
$$\tag{3}  {\cal L}_{EH} ~=~F(g_{\cdot\cdot},\partial g_{\cdot\cdot}) + d_{\mu} \left[ F^{\mu}(g_{\cdot\cdot},\partial g_{\cdot\cdot}) \right]. $$
Total divergence terms in the Lagrangian density give rise to boundary terms in the action. If we just want to derive the EFE in the interior bulk, away from boundaries, then total divergence terms can be neglected. See also this related Phys.SE post.
II) However, the full story is more complicated. Bear in mind that to have a consistent stationary action principle with well-defined functional/variational derivatives, appropriate boundary conditions should be assigned. This is not possible without adding the Gibbons–Hawking–York (GHY) boundary term to the EH action.
