Okay, I checked and you got it mostly right, except the factor $\mu(u)$ in front of the power of $u$ at the very end, in the main equation.

Let's simplify slightly your notations and write $R(\Omega, u)$ and $R(\Omega, u + \varphi)$ for the Rayleigh quotients. I do not compute using a real parameter ; rather, I perform some generalized Taylor exapnsions to get the full derivative of the quotient, which means that $\varphi$ is to be thought of as small itself. Denote by $\|u\|_q^2$ the denominator of $R$.

The numerator of $R(\Omega, u + \varphi)$ expands as  
\begin{equation}
\int_{\Omega} |\nabla u |^2 + \alpha \int_{\partial \Omega} |u|^2 + 2 \int_{\Omega} \nabla u \cdot \nabla \varphi + 2 \alpha \int_{\partial \Omega} u \varphi + \mathcal{O}(\varphi^2).
\end{equation}

Its denominator writes 

\begin{align}
& \|u\|_q^2 \left( 1 + \frac{q}{\|u\|_q^q} \int_{\Omega} |u|^{q-2} u \varphi + \mathcal{O}(\varphi^2) \right)^{\frac{2}{q}} \\
& = \|u\|_q^2 \left( 1 + \frac{2}{\|u\|_q^q} \int_{\Omega} |u|^{q-2} u \varphi + \mathcal{O}(\varphi^2) \right) .
\end{align}

Thanks to the usual $(1+x)^{-1} = 1-x + \mathcal{O}(x^2)$, we get :

\begin{equation}
R(\Omega, u + \varphi) = R(\Omega, u) - \frac{2}{\|u\|_q^q} R(\Omega,u) \int_{\Omega} |u|^{q-2} u \varphi + \frac{2}{\|u\|_q^2} \left( \int_{\Omega} \nabla u \cdot \nabla \varphi + \alpha \int_{\Omega} u \varphi \right) + \mathcal{O}(\varphi^2).
\end{equation}

Using integration by parts, this rewrites 

\begin{equation}
R(\Omega, u + \varphi) = \lambda - \frac{2}{\|u\|_q^q} \lambda \int_{\Omega} |u|^{q-2}u \varphi  + \frac{2}{\|u\|_q^2} \left( - \int_{\Omega} \Delta u \varphi + \int_{\partial \Omega} (\partial_{\nu} u  + \alpha u) \varphi\right).
\end{equation}

Finally we are done. The boundary condition that you got was the right one, whereas the main equation is now

\begin{equation}
- \frac{\Delta u}{\|u\|_q^2} = \lambda \frac{|u|^{q-2}u}{\|u\|_q^q},
\end{equation}

i.e.

\begin{equation}
- {\Delta u} = \lambda \frac{|u|^{q-2}u}{\|u\|_q^{q-2}}.
\end{equation}
Notice that if we set $q=2$ and $\alpha = 0$, we indeed get back to the usual eigenvalue problem for the Neumann laplacian, which is your first case.