I was following along with the MIT Opencourseware "Statistics for Applications" and their analysis of GLMs only covers discussion of GLMs whose dependent variable $y$ is distributed according to a distribution in the canonical exponential family (i.e. of the form below): $$ P(y\ |\ \theta) = h(y)\exp(\theta^Ty - A(\theta)) $$ If this is the case, then one can easily determine that $\mu = \mathbb{E}[Y] = A'(\theta)$, and since $\mu$ is related to the linear predictor via the link function $g$, we have $$ X^T\beta = g(A'(\theta)),\qquad \theta = (g\circ A')^{-1}(X^T\beta) $$
However, I'm curious about the more general case where the distribution belongs to the overdispersed exponential family: $$ P(\mathbf y\ |\ \theta,\tau) = h(y,\tau)\exp\left(\frac{\mathbf b(\mathbf \theta)^T\mathbf{T}(\mathbf y) - A(\theta)}{d(\tau)}\right) $$
Using the identity $\mathbb{E}[\nabla_\theta\ell] = 0$ where $\ell$ is the log-likelihood function I was able to arrive at $$ [D_\theta \mathbf b]^T\mathbb{E}[\mathbf T(\mathbf y)] = \nabla_\theta A $$ but there are two problems with this. For one, the dimension of $\mathbf b$ might not equal the dimension of $\theta$, and even if it did this would not guarantee the matrix $D_\theta\mathbf b$ (the Jacobian of $\mathbf b$) is invertible. Second, even if it were invertible, we would only arrive at $$ \mathbb{E}[\mathbf T(\mathbf y)] = ([D_\theta \mathbf b]^T)^{-1}\nabla_\theta A $$ which doesn't give us any direct correspondence between the parameter $\theta$ and the linear predictor $\mathbf X^T\beta$, since the link function links $\mathbb{E}[\mathbf y]$ to $\mathbf X^T\beta$, not $\mathbb{E}[\mathbf T(\mathbf y)]$.
So, how are GLMs of this form implemented? Where can I go to read more about how this works?
