First of all, it is a bit unclear what do you want exactly. As in:
Do you want an approach that doesn't use coordinates at all, but can, for example, use frames?
Do you want a completely "invariant" approach that absolutely doesn't use any local trivializations of any fiber bundles whatsoever? If so, are you willing to work on the total spaces of said bundles?
Do you want an approach that is global, regardless of any other considerations?
Do you want an approach that doesn't use indices, but everything else is game?
This is not clear to me. Especially, that despite the fact that your issue with Arnold Neumaier's answer was that it uses frames, but the electrodynamics example you have given essentially also uses frames.
Why? Because if $P\rightarrow M$ is a principal $\text{U}(1)$ bundle with principal Ehresmann connection $\omega$, then what you call $A$ is essentially defined as follows. If $(U,\varphi)$ is a local section of $P$ ($U\subseteq M$ is the domain and $\varphi:U\rightarrow P$ is the map), then we define $A^{(U)}=-i\varphi^\ast\omega$, valid in $U$. Since $\omega$ is a pseudotensorial $1$-form of type $\text{Ad}$ on $P$, these pullbacks in overlapping neighborhoods will not fit into one well-defined global object on $M$. Essentially, using $A$ is the same as using $\omega^a_{\ b}$ (the connection forms) in GR. Also note that thus $A$ is well-defined globally if and only if $P$ admits global sections.
With that said, there are several approaches one can take to vary the Einstein-Hilbert action in an "invariant" manner.
The global geometrical approach:
The biggest difficulty is to derive the variation of the volume element $đ\mu=\sqrt{-g}d^4x$. The easiest is of course to expand in local coordinates, but we want to avoid that for now. Aside from working directly on some fiber bundle, I think one needs to use at least frames here.
One reason why you need to use frames here is that a differential form vanishes when you feed linearly dependent vectors into it. So if $đ\mu$ is a volume form, and $X_1,...,X_n$ are vector fields, then $$ đ\mu(X_1,...,X_n)\neq 0\Leftrightarrow X_1,...,X_n\ \text{is a linearly independent set.} $$ But if these $n$ vector fields are linearly independent, then they form a frame, essentially.
The best thing we can do is to use frames, but only use them "passively". Eg. the volume form is not defined in terms of the frame, but the frame is used to derive relations.
We know that if $\Omega$ is any $n$-form and $đ\mu$ is a volume form, then there is a function $\alpha$ such that $$ \Omega=\alpha\ đ\mu, $$ in particular, since the variation of $đ\mu$ is an $n$-form, there is an $\alpha$ such that $$ \delta đ\mu=\alpha đ\mu. $$
If $(e_1,...,e_n)$ is a positively oriented orthonormal frame, then we have $đ\mu(e_1,...,e_n)=1$, so we have $$ \alpha=(\delta đ\mu)(e_1,...,e_n). $$
Although it is not the most rigorous of approaches, it is not immoral or incorrect to obtain this by considering a first order Taylor expansion. Let $g_\epsilon$ be a one-parameter family of metrics, $đ\mu_\epsilon$ the corresponding 1-parameter family of volume forms, and $e_a^\epsilon$ the corresponding 1-parameter family of orthonormal frames.
We have $g_\epsilon(e^\epsilon_a,e^\epsilon_b)=\eta_{ab}$ for all $\epsilon$, so the $\epsilon$-derivative at $0$ is $$0= \delta g(e_a,e_b)+2g(\delta e_a,e_b). $$ Since the $e_b$ are basis vectors, we have $$ g(\delta e_a,\bullet)=-\frac{1}{2}\delta g(e_a,\bullet), $$ so $$ \delta e_a=-\frac{1}{2}\delta g(e_a,\bullet)^\sharp, $$ where $\sharp$ is "raising an index". In index notation this result is $\delta e^\mu_a=\delta g_{\nu}^{\ \mu}e^\nu_a$. Clearly this expression is linear in $e_a$, so let us define $A(e_a)=\delta g(e_a,\bullet)^\sharp$ to be this linear map.
Now, we know that $$ 1=đ\mu_\epsilon(e_1^\epsilon,...,e^\epsilon_n) $$ for all $\epsilon$s. If we differentiate at $\epsilon=0$, we obtain $$ 0=(\deltađ\mu)(e_1...e_n)+đ\mu(\delta e_1,...,e_n)+...+đ\mu(e_1,...,\delta e_n). $$ I personally find the evaluation of the terms other than the first one (which we wish to express) difficult, so instead of directly differentiating at $\epsilon=0$, we expand in a truncated Taylor series. To first order in $\epsilon$, we have $$ e_a^\epsilon=e_a+\epsilon A(e_a)=(1+\epsilon A)e_a, $$ so to first order in $\epsilon$ we also have $$ 1=đ\mu_\epsilon(e^\epsilon_1,...,e^\epsilon_n)=đ\mu_\epsilon(e_1,...,e_n)\det(1+\epsilon A)=đ\mu_\epsilon(e_1...e_n)(1+\epsilon\text{Tr}A)\\ =(1+\epsilon\alpha)đ\mu(e_1...e_n)(1+\epsilon\text{Tr}A)=(1+\epsilon \alpha)(1+\epsilon\text{Tr}A)=1+\epsilon(\alpha+\text{Tr}A)+O(\epsilon^2), $$ so we obtain $$\alpha=-\text{Tr}A=\frac{1}{2}\text{Tr}_g\delta g=-\frac{1}{2}\text{Tr}_g\delta g^\ast, $$ where $\text{Tr}_g$ is the metric trace (eg. $h_{\mu\nu}\mapsto h_{\mu\nu}g^{\mu\nu}$), and $g^\ast$ is the inverse/dual metric. We have used the well known relation that the variations of the direct and the inverse metrics are negatives of one another.
We have thus obtained in a reasonably invariant manner that $$\delta đ\mu=-\frac{1}{2}\text{Tr}_g\delta g^\ast đ\mu=-\frac{1}{2}g_{\mu\nu}\delta g^{\mu\nu}_\ast đ\mu.$$
For the rest of the derivation, I will use abstract index notation, which is global and coordinate-free.
We know that (see Wald - General Relativity for a derivation) if $\nabla^\prime$ and $\nabla$ are two linear connections, with difference tensor $C:\ \nabla^\prime=\nabla+C$, then the corresponding curvature tensors are related as $$ R^{\prime\rho}_{\ \sigma\mu\nu}=R^\rho_{\ \sigma\mu\nu}+2\nabla_{[\mu} C^\rho_{\nu]\sigma}+2C_{[\mu|\lambda|}^{\rho}C^\lambda_{\nu]\sigma}. $$
In particular, let $\nabla^\epsilon$ be the 1-parameter family of Levi-Civita connections induced by the 1-parameter family of metrics, and $C^\epsilon$ be the 1-parameter family of difference tensors between $\nabla^\epsilon$ and the unperturbed LC-connection. Then $$ \delta\nabla=\frac{d}{d\epsilon}\nabla^\epsilon\ |_{\epsilon=0}=\frac{d}{d\epsilon}\left(\nabla+C^\epsilon\right)_{\epsilon=0}=\frac{d}{d\epsilon}C^\epsilon |_{\epsilon=0}\equiv\delta\Gamma, $$ and the $\epsilon$-family of curvature tensors is $$ (R^{\epsilon}){}^\rho_{\ \sigma\mu\nu}=R^\rho_{\ \sigma\mu\nu}+2\nabla_{[\mu}(C^\epsilon)^\rho_{\nu]\sigma}+2(C^\epsilon)^\rho_{[\mu|\lambda|}(C^\epsilon)^\lambda_{\nu]\sigma}, $$
so $$ \delta R^\rho_{\ \sigma\mu\nu}=\nabla_\mu\delta\Gamma^\rho_{\nu\sigma}-\nabla_\nu\delta\Gamma^\rho_{\mu\sigma}. $$
We have thus resolved the calculation of the variations of the connection, curvature tensor and volume form without using local coordinate or frame expansions. You can fill in the rest of the details yourself.
Further reading:
"Besse: Einstein Manifolds".
This is especially recommended if you really really hate indices, as Besse doesn't use abstract index notation, he uses standard mathematician notation. He does not derive the variation of the volume form however, just leaves it to the reader as an excercise. Relevant sections are Chapter 1 - K: First variations of Curvature Tensor Fields and Chapter 4 - C: Total Scalar Curvature: First Order Properties.
Orthonormal frame approach:
This has been covered by Arnold Neumaier's answer, and this is the approach that it closest to your electrodynamics example in spirit. I will only note one thing, that if we wish for the initial value problem in GR to be well-defined, we want spacetime to be globally hyperbolic, which implies that topologically $M=\mathbb R\times\Sigma$. Since $\mathbb R$ is parallelizable, the parallelizability of $M$ hinges upon whether the 3-manifold $\Sigma$ is parallelizable or not.
On physical grounds, we wish space to be orientable, so $\Sigma$ should be orientable. However it is known that every orientable 3-manifold is parallelizable, so if these physically reasonable requirements are met, then 4-dimensional spacetimes are parallelizable.
Therefore, these orthonormal-frame based approaches to GR are actually global!
Further reading:
"Thirring: Course on Mathematical Physics Vol 2 "
"Straumann: General Relativity"
Books by the Loop Quantum Gravity crowd, like Thiemann, Rovelli, Gambini etc. also tend to treat the action principle for GR via orthonormal frames.
Principal bundle approach:
You can regard every tensor field (on $M$) as a certain array of functions on the frame bundle of $M$ that satisfy certain equivariance properties.
For example, in the usual, frame-based approaches, a vector field is something like $V^a(x)$. But these components do not only depend on the manifold point $x$ but also on a frame chosen at $x$. So these $V^a$ are actually functions on the frame bundle: $$ V^a(x,e), $$ where points of the frame bundle $F(M)\rightarrow M$ are denoted as pairs $(x,e)$, where $e$ is a frame at $x$.
These functions define a vector field if and only if for any $\Lambda\in\text{GL}(n,\mathbb R)$, we have the $$ V^a(x,e\Lambda)=(\Lambda^{-1})^a_{\ b}V^b(x,e) $$ equivariance property. Similarly for other tensor fields.
It is important to note that despite the indices and components, these functions are global and completely invariant.
One can develop a tensor calculus on the frame bundle that essentially mirrors the usual index-based local tensor calculus on the base manifold, but is global and totally frame and coordinate independent.
Further reading:
"David Bleecker: Gauge Theory and Variational Principles"
This book treats gauge theories, gravity and gauge theories + gravity in a completely invariant and mathematically precise manner using principal fiber bundles, action principles included.
