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.
optim
is your friend here $\endgroup$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, useoptim
, and get standard errors from the hessian. $\endgroup$