The application is that I want to know how $X$ maps not only to $y$, but to the variance of $y$. I think I've worked out a reasonable solution for doing so using the gamma distribution and the inverse square link, but I've got a couple of questions. The setup is as follows: $$ y = f(X) + \epsilon\\ \epsilon = \phi h(X)\sigma $$ $\phi$ is an unknown scale parameter, $h(X)$ is a deterministic function, and $\sigma$ is a unit-variance normal random variable.
I fit the first part using OLS. Coefficients are consistent in the presence of heteroskedasticity, and I'm not worried about intervals yet.
Back to the variance. I square everything and take the expectation, subbing in the estimated residuals:
$$ E[\hat\epsilon^2] = \phi^2h(X)^2 $$ The sigma drops out because the expectation of a squared normal is the variance of the normal, which is 1.
Here is where I am less sure of the way forward. If I take $h(X)^2$ to be $1/\mu$, and use an inverse square link function, then I have $\mu= h(X)$. This I can estimate with a GLM (and when I do so the diagnostic plots look pretty good).
If I estimate the GLM using ML (as opposed to the GCV or REML options in GAM), I can extract the ML estimates of the scale parameter, which is $1/\phi$.
Going back to the original model, I can get the expected value of the residual at a given $X$ by taking the square root of $\phi^2\hat h(X)^2$. And if I want to improve the efficiency of my OLS estimates, I can use the estimates as weights and go GLS, and then iterate further until approximate convergence.
Anything wrong with this? Gamma is confusing and there are two different parameterizations that arbitrarily involve inverting one thing and not another.
Is the use of the inverse square link valid here, in the sense of, "am I understanding the role of the link function correctly"? Seems a lot less intuitive than with e.g. logit for binomial data.
Also, practically, when I fit the gamma GLM (or GAM) in
mgcv and then ask for
mod$scale, am I getting $\phi$ or $1/\phi$?
Lastly, is there any way to know whether the iterative GLS/gamma GLM procedure that I propose converges to truth after some fashion?