Showing that $\sum_{i=1}^n (y_i-\hat{y_i})(\hat{y_i} - \bar{y}) = 0$ for the generalized linear model Exercise :

Prove that for the generalized linear model, it is :
  $$\sum_{i=1}^n (y_i-\hat{y_i})(\hat{y_i} - \bar{y}) = 0$$

Question : How would one proceed with proving that for the generalized linear model ? I can prove it for the simple linear model but I seem to be stuck a bit for the generalized case. 
 A: In the linear model, $Y=X\beta +\epsilon$, $\hat \beta = (X'X)^{-1}X'Y$ and $\hat Y = X\hat \beta = X(X'X)^{-1}X'Y$
$\sum_{i=1}(y_i−\hat y_i)(\hat y_i−\bar y)=(Y−\hat Y)′(\hat Y−\bar Y) = Y'(I-H)(H-\frac 1 n J)Y = Y'(I-H)HY - Y(I-H)\frac 1 n JY $
where $H=X(X'X)^{-1}X'$ , and $J$ is matrix with elements 1.
$Y'(I-H)HY = Y'(H-HH)Y =0$
$Y'(I-H)\frac 1 n JY = Y'(I-H) \bar Y = \bar Y \sum_{i=1}(y_i−\hat y_i) = 0$ 
So $\sum_{i=1}(y_i−\hat y_i)(\hat y_i−\bar y) = 0$
A: $$\begin{array}{rcl} \sum_{i=1}^n (y_i-\hat{y_i})(\hat{y_i} - \bar{y}) =  \sum_{i=1}^n \epsilon_i(\hat{y_i} - \bar{y}) 
 = \underbrace{\sum_{i=1}^n \epsilon_i \hat{y_i}}_{\text{ equals zero if } \epsilon \perp \hat{y} } &-& \bar{y_i}\underbrace{  \sum_{i=1}^n \epsilon_i}_{  \text{equals zero} } 
\end{array}$$
I do not believe that you have $\epsilon \perp \hat{y}$ for GLM, this is only the case if you have the residual vector represent the shortest distance between $y$ and the solution space of all possible $\hat{y}$, which is the least squares solution and relates to the assumption of Gaussian distributed errors.
Neither is $ \sum_{i=1}^n \epsilon_i = 0 $ general. This means $\epsilon$ is perpendicular to the intercept term. Even in the case of a simple linear model this may not be true (simple linear regression without the intercept term).

(Computational) counter-examples:
> # some data
> x <- c(1.0, 2.0, 3.0, 4.0)
> y <- c(1.1, 2.1, 3.2, 4.0)
> 
> # OLS (minimize squared errors)
> hat_y <- predict(lm(y~x))
> sum((y-hat_y)*(hat_y-mean(y)))
[1] 6.366435e-16
> 
> # GLM (using Gamma)
> hat_y <- predict(glm(y~x, family = Gamma(link="identity")),type="response")
> sum((y-hat_y)*(hat_y-mean(y)))
[1] -0.1042444
> 
> # GLM (using Gaussian but different link)
> hat_y <- predict(glm(y~x, family = gaussian(link="log")),type="response")
> sum((y-hat_y)*(hat_y-mean(y)))
[1] 0.1504989
>
> # OLS without intercept term
> hat_y <- predict(lm(y~0+x))
> sum((y-hat_y)*(hat_y-mean(y)))
[1] -0.26

