Try this explanation on for size. I emphasise that it is my own way of understanding the $U(1)$ Gauge Invariance of electrodynamics and I haven't seen it elsewhere in exactly the same words.
The derivation you cite is wontedly given in the context of the semi-classical (i.e. first quantized, or before the quantum field is introduced) Dirac or Schrödinger equation for the electron: the same reasoning applies to both. These equations describe a fermion, so one cannot re-interpret their particle fields $\psi$ as a macroscopic, classically-measureable field: the Pauli exclusion principle means that you cannot copy your fermion and have $N$ (where $N$ is very big) particles in the same quantum state. Otherwise, you could in principle measure the full complex value of $\psi(x,y,z,t)$ to arbitrary accuracy by copying $\psi$ in this way and then doing a classical measurement - something that you can do with bosons (see below).
What does this mean at the one particle level? It means that only $|\psi|^2$ is experimentally meaningful for one particle: you can measure the probability of finding the particle at a position in space, but the phase of $\psi$ has no such meaning. So we should be able to multiply $\psi$ by an arbitrary phasing function $e^{i\,\alpha(x,y,z,t)}$ and get something that means physically the same thing. But, in the light of the structure of the Dirac or Schrödinger equations this seems absurd: most assuredly:
$\partial_j \left[e^{i\,\alpha(x,y,z,t)}\,\psi(x,y,z,t)\right]\neq e^{i\,\alpha(x,y,z,t)}\, \partial_j\,\psi(x,y,z,t)$
but rather
$\partial_j\left[ e^{i\,\alpha(x,y,z,t)}\,\psi(x,y,z,t)\right]= e^{i\,\alpha(x,y,z,t)}\, \left[\partial_j\,\psi(x,y,z,t) + i\, \psi(x,y,z,t)\,\partial_j\,\alpha(x,y,z,t)\right]$
and so such an assumption of invariance with respect to arbitrary phasing is invalidated because phasing patently "ruins" the structure of the Dirac or Schrödinger equations unless the phase factor $\alpha$ is globally constant. This is reasonable: the phase of $\psi$ is most definitely part of the solutions of the equations and plays a definite role in diffraction and other wave effects that have a bearing on the intensity field $|\psi|^2$.
But one can "retrieve" the situation by postulating that the single particle is coupled to some outside field, so that the Dirac or Schrödinger equations now have the terms involving this coupled-in field: if so, we can add an arbitrary phase $\alpha(x,y,z,t)$ and keep the global equation the same by saying that, whenever we do this, we must take away a balancing $i \partial_j\,\alpha(x,y,z,t)$ from the coupled in field. Thus, we conclude, if the coupled-in field $A_j$ is at all physical, it has to give the same measurements as the field $A_j + \partial_j\,\alpha(x,y,z,t)$: the original plus any arbitrary field of the form $\partial_j\,\alpha(x,y,z,t)$, where $\alpha(x,y,z,t)$ is some suitably well-defined scalar field. But there is a field that behaves exactly like this: the electromagnetic four-potential. We can add a spatial gradient $\nabla \alpha$ to the vector part and at the same time add the scalar $\partial_t\alpha$ to the scalar electric potential without affecting the electric $\mathbf{E}$ and magnetic $\mathbf{B}$ fields.
So there we have it: the assumed "gauge invariance" suggests the electromagnetic field, because the balancing field behaves the same way as gauge-transformed vector magnetic and electric potentials. So the essential reasoning here is really a hunch: if it looks like a duck and quacks, maybe it is a duck. So too with the outside field coupling into the Dirac equation. It behaves like the potential fields of Maxwell's electromagnetism, so "we" (or rather the physicist who first thought of this - my ignorance sadly hinders my telling you who) go with our hunch and see what happens when we assume the field is the electromagnetic field.
Note that the same ideas do not apply to bosons. We can think of Maxwell's equations as the first-quantized equations for the photon. No one speaks of making Maxwell's equations invariant with respect to multiplication by an arbitrary phase function $e^{i\,\alpha(x,y,z,t)}$ in the same way as is done for the Dirac equation. This is even though the Maxwell equations can be cast in a quaternionic form that is identical to a zero-mass Dirac equation when the latter is thought of as two quaternionic equations coupled by a mass term - so the same mathematical trick would be just as valid with the Maxwell equations as it would with the Dirac equation. My interpretation is this: the photon's phase is classically meaningful: photons are bosons, so in principle we can get as many as we like in the same state: so many indeed that we can measure the phase of their common "wave function" (which is now the electromagnetic field - see caveats in the afternotes below) to arbitrary accuracy with classical measurement devices such as interferometers. Every one-particle photon state corresponds EXACTLY to a macroscopic, classical electromagnetic field that we can set up and measure in detail to arbitrary accuracy in the laboratory. So, even though the phase of one photon is "hidden" just as is the phase of one electron above - there are no "phase" eigenstates and no "phase" observable - it must still be "absolute" in the "prolifically-copy-and-classically-measure" sense just described.
This ends my answer, but I add some interesting related things below.
Maxwell's Equations from $U(1)$ Gauge Invariance
Incidentially, one can use this gauge invariance thinking to motivate or "derive" the Maxwell equations. Given suitable differentiability assumptions on the field, the simplest way to derive fields that are exactly unaffected by the gauge transformations is to form the tensor curl:
$F_{\mu\,\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$
(try this if you've not before: it leaves nonzero $F_{\mu\,\nu}$ that are unaffected by the gauge transformation). Now we postulate a Lorentz-invariant wave equation for a massless field: $\Box \mathbf{A} = \mathbf{0}$ is the obvious one (let's just go all the way back to D'Alembert!). You straight away have get the freespace Maxwell equations as equations fulfilled by the "physical", unaffected by gauge $F_{\mu,\nu}$.
How ironic that the Pauli Exclusion Principle, which was first postulated to "stop" electrons in the Bohr atom from radiating and thus told against behavior seemingly foretold by the Maxwell equations can be used to motivate a gauge transformation that pretty much leads from the Dirac equation back to the Maxwell equations!
An intriguing and uncluttered "first quantized" or "baby" formulation of QED (appealing especially to non quantum field theorists like me) can be gotten by postulating the Dirac-Maxwell equation:
$\gamma^\mu\left(i \partial_\mu - q A_\mu\right) \psi + V \psi - \psi = 0$
$\partial_\nu F^{\nu\,\mu} = q\,\bar{\psi} \gamma^\mu \psi$
with the Lorenz Gauge $\partial_\mu A^\mu = 0$
i.e. intuitively, the source for the Maxwell equations is $q$ times the probability current density (here $\bar{\psi}$ stands for charge conjugated $\psi$). A first quantized electron field is now nonlinearly coupled to a first quantized photon field. This nonlinear system can be solved exactly for the hydrogen atom potential $V$ and in other situations
A. O. Barut and J. Kraus, "Nonperturbative Quantum Electrodynamics: The Lamb Shift", Foundations of Physics, Vol. 13, No. 2, 1983
and this solution can indeed model the Lamb shift and spontaneous emission. The series solution comes out to something very like the standard QED perturbation terms and indeed one must step in and "renormalize" this solution too.
See also the works of Hilary Booth in the late 1990s.
Maxwell's Equations as the Propagation Equations Photon Wave Function
For various reasons, a first quantized "photon wave function" has some different behaviours and interpretations from the fermion field $\psi$ in the Dirac equation. It is perfectly valid to treat the solutions of Maxwell's equations as one-photon quantum states, but one must understand that there are difficulties in handling localized position eigenstates for the photon precisely analogous to those for fermions. See the works of Iwo Bialynicki-Birula, Margaret Hawton in the late 1990s and early 2000s, for example:
Iwo Bialynicki-Birula, "On the Wave Function of the Photon" Acta Physica Polonica 86, 97-116 (1994)
Margaret Hawton and William E. Baylis, "Angular momentum and the geometrical gauge of localized photon states", Phys. Rev. A 71, 033816 (2005)
