Lagrange's equation is form invariant under EVERY coordinate transformation. Hamilton's equations are not under EVERY phase space transformation. Why? When we make an arbitrary invertible, differentiable coordinate transformation $$s_i=s_i(q_1,q_2,...q_n,t),\forall i,$$ the Lagrange's equation in terms of old coordinates $$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{q}_i}\right)-\frac{\partial L}{\partial q_i}=0,\forall i$$ changes to $$\frac{d}{dt}\left(\frac{\partial \hat{L}}{\partial\dot{s}_i}\right)-\frac{\partial \hat{L}}{\partial s_i}=0, \forall i$$ where $\hat{L}(s,\dot{s},t)$ is obtained from $L(q,\dot{q},t)$ by $$\hat{L}(s,\dot{s},t)=L(q(s,t),\dot{q}(s,\dot{s},t),t).$$
When we go to the Hamiltonian framework, things do not remain so simple. Under an arbitrary transformation in phase space, $$~~~Q_i\to Q_i(q_1,q_2...,p_1,p_2,..,t),\forall i,\\
P_i\to P_i(q_1,q_2...,p_1,p_2,..,t), \forall i$$ the Hamilton's equations do not remain form invariant. This only happens for a restricted class of transformations, called the canonical transformations.  Also the new Hamiltonian is not obtained from the old Hamiltonin by $$\hat{H}(Q,P,t)= H(q_i(Q,P,t),p_i(Q,P,t),t)$$ even when the transformation is canonical, unless that canonical transformation is also a symmetry.
Why is the reason for this?
 A: The other answers have explained the issue in a mostly mathematically formal and rigorous way. I want to add to the discussion by instead trying to explain why the Lagrange equations of motion stay the same in an informal yet hopefully intuitive way and then discuss what is different for Hamiltons equations of motion.
First let us deal with Lagrange. Let us remember what the actual problem is. We have an action of the form $S=\int dt L$ and we want to find the trajectory for which this action is extremal. Here the Lagrangian $L$ is (for a given point in time) a function of properties of the trajectory like the actual physical location and velocity. To actually do a computation we choose coordinates to describe the physical space, call them $q(t)$, where the dependence of the Lagrangian on the physical properties of the trajectory (and time but let us assume from now on that we do not have an explicit time dependence) can now be written down as $L(q(t),\dot{q}(t),t)$ (I only look at a 1D problem, the generalization to more dimensions is conceptually easy and this way I have less to type). Of yourse our choice of coordinates does not influence the solution. If we instead take different coordinates $Q(t)$, we will still get the same trajectory, just in these different coordinates, because the problem we are solving is the same. We just have to know the changes of coordinates $Q(q)$ and $q(Q)$ and the Lagrangian must now be written as $L(q(Q(t)),\dot{q}(Q(t),\dot{Q}(t)))$ where the dependence $\dot{q}(Q(t),\dot{Q}(t))$ can be found by taking the time derivative of $q(Q(t))$.
Actually we even know how to solve the problem: The Lagrange equation of motion which as a mathematical fact is also known as the Euler-Lagrane equation. This equation states that the integral $\int dt G(x(t),\dot{x}(t))$ becomes extremal if $$\frac{d}{dt}\frac{\partial G}{\partial \dot{x}}-\frac{\partial G}{\partial x}=0$$
It has been derived by only using the information that the integral (which in our case is the action) is extremal. So if two people solve the problem, person 1 with the $q$-coordinates, person 2 with $Q$, both know that they can use the Euler Lagrange equation in their respective coordinates and both will get the same result, just each in their coordinates, because both solved the same problem.
But why isn't that the case for Hamiltons equations of motion? To see that let us have a look on how we go from Lagrange to Hamilton. The very basic idea is that we introduce a "new" variable $p$. This is of course not really new, it is in fact definied to be $p=\frac{\partial L}{\partial \dot{q}}$, so at each point in time it is a function of $q(t)$ and $\dot{q}(t)$, i.e. $p(q,\dot{q})$. We assume that this dependence can be inverted to find $\dot{q}=f(q,p)$ for some function $f$. Now we can just write our action as $S=\int dt L(q(t),f(q(t),p(t)))$ and again want to find functions $q(t)$ and $p(t)$ that make it extremal, right? No, that would be wrong. At least if we now vary $p$ and $q$ independently, we will not get the correct result. Instead we have to respect the condition that $\dot{q}=f(q,p)$. So now we don't just have an optimization problem, but an optimization problem with constraints, namely that $\dot{q}=f(q,p)$. To implement this constraint one could use the method of Lagrange multipliers, i.e. one adds to $L$ a term like $\lambda(\dot{q}-f(q,p))$ where $\lambda$ is this Lagrange multiplier. But that is not the only possibility. Instead we add a term $p(\dot{q}-f(q,p))$. One can verify that this results in the correct equation of motion for $q$ when varying $p$ and $q$ independently, i.e. using the Euler-Lagrange equation for $$\int dt (L(q(t),f(q(t),p(t)))+p(t)(\dot{q}(t)-f(q(t),p(t))))$$
Reordering the terms and defining $H(q,p)=pf(q,p)-L(q,f(q,p))$ we have rewritten our initial problem now as the optimization of $$\int dt (p(t)\dot{q}(t)-H(q(t),p(t))$$
Here $q$ and $p$ get varied independently. If we do that we get Hamiltons equation of motion.
But where does our argument from above now fail when we consider a "change of coordinates" $Q(q,p)$, $P(q,p)$? We still have an optimization problem as before. But now the difference is that, where before we had a completely arbitrary function in the integral and just plugged that into the Euler-Lagrange equation to get Lagranges equation of motion, we now rely on the integrand having a specific form. This specific form will not generally be preserved when changing coordinates. So of course we can still use the Euler-Lagrange equations to extremize the integral, but we will not necessarily get Hamiltons equations in the new coordinates. Let's see what happens:
As before the change of coordinates can be implemented by replacing the old variables by the corresponding expressions depending on the new variables, i.e. $q\to q(Q,P)$ and $p\to p(Q,P)$, leading to $$ \int dt (p(Q(t),P(t))\dot{q}(Q(t),P(t))-H(q((Q(t),P(t))),p((Q(t),P(t))))$$ having to be optimized. Observe that it does not have the form $$\int dt (P(t)\dot{Q}(t)-H(q((Q(t),P(t))),p((Q(t),P(t))))$$ which would be the expression leading to Hamiltons equations in the new coordinates. In fact the difference between the integrand we have and the integrand giving us Hamiltons equations is $p(Q(t),P(t))\dot{q}(Q(t),P(t))-P(t)\dot{Q}(t)$.
There are of course cases where this difference does not impact the actual equations. A sufficient condition for that is that it is a total time derivative, i.e. something like $\frac{d}{dt} F(Q(t),P(t))$ for some function $F$. This leads to the same equation of motion because when integrating a time derivative one just gets the function $F$ evaluated at the boundaries of the intervall of integration (which I omitted everywhere) and the variations of $q,p$ are zero at the boundary (this is used in the derivation of Lagranges equation of motion, if you don't remember it, review that part!), so the integral of the time derivative does not change under the variation.
By the way: This condition is sufficient, but not necessary. For example you could imagine adding a term to the integrand which is proportional to the integrand itself. This would of course also leave the equations of motion invariant. But in the case of canonical transformations it is generally assumed that the difference is a total time derivative.
So now we know that if $p(Q(t),P(t))\dot{q}(Q(t),P(t))-P(t)\dot{Q}(t)$ ist a total time derivative, i.e. equal to $\frac{\partial F}{\partial Q}\dot{Q}(t)+\frac{\partial F}{\partial P}\dot{P}(t)$, then we definitely get Hamiltons equations in both coordinate systems. One can also see that by plugging $\frac{\partial F}{\partial Q}\dot{Q}(t)+\frac{\partial F}{\partial P}\dot{P}(t)$ into the Euler-Lagrange equation for $P,Q$ which gives trivially zero, without restricting $P,Q$ in any way.
Because the mixed second derivatives have to be equal, i.e. $\frac{\partial^2 F}{\partial Q\partial P}=\frac{\partial^2 F}{\partial P\partial Q}$ one obtains from $p(Q(t),P(t))\dot{q}(Q(t),P(t))-P(t)\dot{Q}(t)=\frac{\partial F}{\partial Q}\dot{Q}(t)+\frac{\partial F}{\partial P}\dot{P}(t)$ the condition $\frac{\partial q}{\partial Q}\frac{\partial p}{\partial P}-1=\frac{\partial p}{\partial Q}\frac{\partial q}{\partial P}$ which can be rewritten in the more familiar form $\{q,p\}_{Q,P}=1$. This shows that this Poisson bracket being equal to 1 (i.e. the change of coordinates being a canonical transformation) is a sufficient condition for Hamiltons equations to be preserved.
A: You should really think about the variables we use as being like coordinates on some manifold, the configuration space (roughly the same as the phase space, I won't be careful about the distinction). In this language, changing variables is equivalent to changing coordinates on this manifold. The action is some scalar function on this space, and we can take coordinate derivatives, $\frac{\delta S}{\delta q_a}$ with respect to whatever coordinates $q_a$ we are using on the space. As in multivariable calculus, we can form the directional derivatives $D_v S = v_a\frac{\delta S}{\delta q_a}$ for any vector with components $v_a$ we like. If you want something more formal and geometric, the directional derivative is a Lie derivative on the configuration manifold.
Now, remember that when we vary the action, we demand that the variation be stationary. That is, we demand that all directional derivatives vanish, meaning $D_v S=0$ for all vectors $v$. This statement about the vanishing of directional derivatives, you'll note, it entirely agnostic to the coordinates we use, but nonetheless implies that if we use coordinates $q_a$, that $\frac{\delta S}{\delta q_a}=0$ for each $a$. But any coordinate system would result in the same condition that all the coordinate derivatives vanish. This should also be in some sense familiar from multivariable calculus.
So this is casting the invariance of the Euler-Lagrange equations in a geometric language. Aside from being nice, this is also going to be the right language to understand what's going on in the Hamiltonian picture.
The phase space is normally coordinatized by pairs of coordinates $(p^a,q_a)$, but this is really not necessary. At the end of the day, the phase space is again a manifold and the $(p,q)$ are simply a special coordinate system on that manifold (Darboux's theorem implies that we can always, at least locally, find such a coordinate system). The thing that defines these special coordinates, really, is that the symplectic form takes a very special form.
In case you are less than familiar with symplectic forms, let me do the following to motivate the idea. Instead of using the coordinates $(p^a,q_a)$, instead use a collective coordinate $\xi^a = \langle p^a, q_a\rangle$, so all I've really done is put the $p$'s and $q$'s into one big vector. Just to be clear, if $q_a$ and $p^a$ were $n$-dimensional vectors, then $\xi^a$ is a $2n$-dimensional vector formed by concatinating the components (well, any way of putting the components together will do...this will just change the precise form of the $\Omega$ I introduce in a moment by permuting it's rows and columns appropriately). In terms of this, the Hamilton's equations may now be written
$$
\frac{d\xi^a}{d t}=\Omega^{ab}\frac{\partial H}{\partial \xi^b}
$$
where $\Omega$ is some matrix. Usually, this looks something like
$$
\Omega=\left(\begin{array}{cc}
0 & -1\\
1 & 0\end{array}\right).
$$
This $\Omega$ is typically known as the inverse of the symplectic form. Though sometimes you'll just hear it called the symplectic form (an abuse of language, but not uncommon) or, more accurately, a Poisson bivector. These names are not so important to what I want to say, but I figure I may as well mention the correct terminology for anyone who wants to try searching online themselves.
Now, the symplectic form does, in fact, transform under coordinate changes just like you might expect any tensorial object over a manifold to do. If we take on faith that the symplectic form should transform as a tensor under coordinate changes, then we already know how the right-hand side of the rewritten Hamilton's equation transforms if we were to transform to some other coordinate system. But let's not worry about this assumption and start by investigating the left-hand side.
Suppose we perform some transformation $\xi^a=\xi^a(\zeta)$ to a new coordinate system $\zeta$. Then by chain rule,
$$
\frac{d \xi^a}{d t}=\frac{\partial \xi^a}{\partial \zeta ^\mu}\frac{d\zeta^\mu}{d t},
$$
so we see, perhaps unsurprisingly the time derivative also transforms as a tensor under a coordinate change (I used $\mu$ for the indices of the new coordinates just to keep things visually distinct).
So in the end, we find that Hamilton's equations transform as
$$
\frac{d \xi^a}{d t}=\Omega^{ab}\frac{\partial H}{\partial \xi^b}\implies \frac{\partial \xi^a}{\partial \zeta ^\mu}\frac{d\zeta^\mu}{d t} = \Omega^{ab}\frac{\partial \zeta^\nu}{\partial \xi^b}\frac{\partial H}{\partial \zeta^\nu},
$$
which if we move the Jacobian on the left to the right, we find
$$
\frac{d\zeta^\mu}{dt}=\left(\frac{\partial\zeta^\mu}{\partial \xi^a}\Omega^{ab}\frac{\partial\zeta^\nu}{\partial\xi^b}\right)\frac{\partial H}{\partial \zeta^\nu}
$$
and hence we see that if we define a new symplectic form $\Omega^\prime$ by
$$
\Omega^{\prime\, \mu\nu}=\frac{\partial\zeta^\mu}{\partial \xi^a}\Omega^{ab}\frac{\partial\zeta^\nu}{\partial\xi^b},
$$
(which is equivalent to my assertion that $\Omega$ should be tensorial under coordinate changes) Hamilton's equation actually are invariant in the sense that we still have equations of the form
$$
\frac{d \zeta^\mu}{dt}=\Omega^{\prime\, \mu\nu}\frac{\partial H}{\partial\zeta^\nu}.
$$
The only difference is that the components of $\Omega$ and $\Omega^\prime$ might not be the same. But really this shouldn't be such a shock after all this setup since we wouldn't expect the components of a tensor to remain the same after a coordinate change.
Consider as an example the Minkowski metric. We know what this looks like in Cartesian coordinates. If we changed to polar coordinates, for example, of course the component entries in the metric change, but it's still, in a very real sense, the same metric. We just have a new representation thereof.
So where to canonical transformations fit into all this? They are simply the very special coordinate transformations which actually leave the components of the symplectic form invariant. Formally, these are coordinate transformations generated by vector fields over phase space whose Lie derivative of the symplectic form vanishes. This is very similar in many respects to a vector field being a Killing vector of some metric.
Finally, I should point out that the way I have framed the entire discussion above, it may seem strange why we should consider canonical transformations at all. After all, we can use any transformation at the cost of a nice form for the symplectic form. Perhaps non-canonical transforms can put the equations in a nice form.
In principle this is of course true. However canonical transformations play a very essential role which is intimately tied to Noether's theorem and symmetry. Essentially one can guarantee that every symmetry of the action corresponds to precisely a canonical transformation. Furthermore, only vector fields which correspond to canonical transformations are guaranteed to have a charge associated to them (like the Hamiltonian is the charge associated to time evolution (which is itself a canonical transformation)).
A: *

*A more geometric approach to the Hamiltonian formulation is to consider the $(2n+1)$-dimensional contact manifold ${\cal M}$ with coordinates $(q^i,p_j,t)$. The Hamiltonian action functional is
$$S_H[\gamma]~=~\int_I \gamma^{\ast} \Theta, \qquad \Theta~=~p_j \mathrm{d}q^j -H \mathrm{d}t,  \tag{1}$$
where $\gamma:I\to {\cal M}$ is a curve.
From this it becomes clear that the Hamiltonian $H$ is not a scalar object but rather the $t$-component of a 1-form/co-vector, and therefore transforms non-trivially under coordinate transformations. This answers OP's main question. See also e.g. my related Phys.SE answer here.


*Furthermore, the Hamiltonian formulation can be generalized to non-canonical coordinates, cf. e.g. my Phys.SE answer here.
