Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a set of data I need to model using the following:

$$Y_i \sim N(\mu_i,\theta\mu^{2}_{i}) \quad \text{and}\quad \log \mu_i = \beta^{T}X_{i}$$

($\theta\mu_i^2$ is the variance of $Y_i$). I need to estimate the model parameters $\beta$ and $\theta$. I am not sure how to go about that. It was suggested to try to integrate out $\theta$ using joint MLE but once again I am not sure how to go about that.

I am also trying to implement this model in R, but as I am new to R and linear models I wasn't quite sure how to go about this. I tried:

glm(Days ~ ., family=gaussian(link="log"), data=quine,

But I am not sure if this is the right model to use. Should I be looking into nls or gls?

Any suggestions would be appreciated.

Edit: So thinking about this some more, would this be equivalent to a weighted poisson with weight $\frac{1}{\theta\mu}$. So in R I would implement maybe as:

glm(Days ~ ., family=poisson, data=quine,weights = 1/(theta*predict(model1,type="pear")))

Where model1 is an un-weighted poisson.

I still don't know how to estimate theta though.

share|improve this question
In your notation, is $\theta \mu_i^2$ the standard deviation or the variance? (The solution appears to be simpler if it's the latter.) – whuber Dec 12 '12 at 18:27
It is the variance. – scalar Dec 12 '12 at 18:29
optim is your friend here – conjugateprior Dec 13 '12 at 21:26
If you're going down the 'define your own GLM' route then you should take a look at the family help page. There is a parameterised link function example there that you could adjust to get your theta in. Personally I'd still just code the log likelihood, use optim, and get standard errors from the hessian. – conjugateprior Dec 13 '12 at 21:45
up vote 2 down vote accepted


I'm not sure a Poisson regression would be appropriate here, since by definition E(X)=Var(X) for a Poisson distribution, and you're looking to fit the data to a $N(\mu, \theta\mu^{2})$ distribution (correct me if I'm wrong).

The glm() function in R does include an option for user-defined variance and link functions using maximum quasi-likelihood, see the following site:

share|improve this answer
Using your suggestion, I tried the following in R: glm(Days ~., family=quasi(variance="mu^2", link="log"),data=quine) I would still need to figure out a way to estimate theta. – scalar Dec 13 '12 at 16:18
Hmm, I wonder if the gamma family of distributions is what you're looking for, where $Y \sim Gamma(k\theta, k\theta^{2})$. Y is restricted to positive values with the gamma distribution, which may or may not match your data. For large values of the scaling parameter k, the gamma dbn approaches the normal dbn, which is in line with your initial model. – RobertF Dec 13 '12 at 22:05
Yea, I think gamma should work - $\mu= k\theta$, so then $Var(Y) = (1/k)\mu^{2}$, where $1/k$ is the additional variance parameter. – RobertF Dec 13 '12 at 22:13
Thank you for the help. So looking at the ?glm in R, a gamma would be fit using the following: model1 <- glm(Days ~.,quasi(link=log, variance="mu^2"), data=quine Would I use gamma.shape(gm, verbose = TRUE) to extract my estimate of theta? – scalar Dec 13 '12 at 22:23
Just use family=gamma(link="log") - I don't think you need to specify variance. Yes, gamma.shape returns the theta (=shape) estimate. – RobertF Dec 14 '12 at 16:00
  1. Set up the Log Likelihood function using the Normal density function provided. This is a function of $θ$ and $\mu_{i}$ .

  2. Next substitute $\mu_{i}$ in terms of $β$ and $X_i$.

  3. You now have the Log Likelihood function written in terms of the parameters $β$ and $θ$.

Write R function to calculate MLE

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.