3
$\begingroup$

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,
              start=c(log(mean(quine$Days)),0,0,0,0,0,0))

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ 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.) $\endgroup$
    – whuber ♦
    Dec 12, 2012 at 18:27
  • 1
    $\begingroup$ It is the variance. $\endgroup$
    – scalar
    Dec 12, 2012 at 18:29
  • $\begingroup$ optim is your friend here $\endgroup$
    – conjugateprior
    Dec 13, 2012 at 21:26
  • $\begingroup$ 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. $\endgroup$
    – conjugateprior
    Dec 13, 2012 at 21:45

2 Answers 2

Reset to default
2
$\begingroup$

Lok,

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: http://data.princeton.edu/R/glms.html.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – scalar
    Dec 13, 2012 at 16:18
  • $\begingroup$ 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. $\endgroup$
    – RobertF
    Dec 13, 2012 at 22:05
  • $\begingroup$ 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. $\endgroup$
    – RobertF
    Dec 13, 2012 at 22:13
  • $\begingroup$ 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? $\endgroup$
    – scalar
    Dec 13, 2012 at 22:23
  • $\begingroup$ Just use family=gamma(link="log") - I don't think you need to specify variance. Yes, gamma.shape returns the theta (=shape) estimate. $\endgroup$
    – RobertF
    Dec 14, 2012 at 16:00
3
$\begingroup$

  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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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