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?

  • 2
    $\begingroup$ "I can use the estimates as weights" -- well, you would use the inverses of the variances as weights, not the variances themselves. I presume you already knew that though, but its important to make the distinction clear. $\endgroup$ – Glen_b Dec 10 '13 at 22:28
  • 1
    $\begingroup$ You are not modelling conditional heteroskedasticity - your error term itself is dependent on the regressors (your second equation in the original setup), and not just its conditional variance . Conditional heteroskedasticity is directly defined at the level of conditional moments of the variable (the second moment specifically), not the variable itself : it is from the start $Ε(u^2\mid X) = g(X)$ - we do not build this relation starting from $u= k(X)$ and then squaring and taking the expected value. This way that you do it creates "a lot more of dependence" between regressors and error term. $\endgroup$ – Alecos Papadopoulos Dec 10 '13 at 23:08
  • 1
    $\begingroup$ @AlecosPapadopoulos: OK, it is possible that there is some better method; I invented the above (with a little help from google). There will definitely be dependence between the regressors and the estimate of h(x). I did it the way i did it because it seems like a way of capturing the degree to which the amplitude of the noise changes at different levels of the regressors. Do you have references for better methods? Again, all I want to do is map both the mean and the variance of $y$ to $X$. Also I don't follow your math; is not $u = f(X)$, which is unknown? $\endgroup$ – generic_user Dec 10 '13 at 23:38
  • $\begingroup$ @Glen_b: yes indeed. $\endgroup$ – generic_user Dec 10 '13 at 23:39
  • 1
    $\begingroup$ No, we do not define the error term (the noise) as a function of the regressors - this defies the very purpose of estimation (if they were directly related the "noise" would not really be noise). We define its (unknown) variance conditional on the regressors, as a function of the regressors. This is standard heteroskedasticity treatment... I have a feeling that you are mixing this "traditional" form of heteroskedasticity with the "autoregressive conditional heteroskedasticity" which is a different class of heteroskedasticity models. $\endgroup$ – Alecos Papadopoulos Dec 10 '13 at 23:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.