# Parameter estimation for a model defined by $Y_i \sim N(\mu_i,\theta\mu^{2}_{i})$

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. • 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, 2012 at 18:27 • It is the variance. Dec 12, 2012 at 18:29 • optim is your friend here Dec 13, 2012 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. Dec 13, 2012 at 21:45 ## 2 Answers 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. • 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. Dec 13, 2012 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. Dec 13, 2012 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. Dec 13, 2012 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? Dec 13, 2012 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. Dec 14, 2012 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