[R-sig-ME] glmm with a tweedie distribution

[Zhang,Yanwei]

Hi Ruben, 

Thank you for the input and the reference. The profile likelihood approach is of course a standard way to estimate the variance function. Because it's so fast, I'm planning to use this to generate starting values in the new version of the cplm package. Indeed, there is another much easier way - since the p parameter is only related to the shape parameter in the underlying gamma distribution and it changes little as the mean varies, you can fit a Gamma regression to the severity data (the positive catch rates only) to get a rough estimate of the shape parameter and then derive the value of p. 

That being said, I still prefer the likelihood-based approach for the Tweedie mixed models.  There is some  philosophical gap when using the PQL approach with the profile likelihood. 1) it's an approximation to the model rather than the underlying likelihood, but the profile likelihood approach depends on the approximated likelihood. 2) it's a quasi-likelihood - glmmPQL reports NA on the likelihood - so a more appropriate way is to use the extended quasi-likelihood. But this will need some ad hoc adjustment because the extended quasi-likelihood does not allow zeros. See the argument in Dunn and Smyth (2005): Series evaluation of Tweedie exponential dispersion models densities. Statistics and Computing, 15, 267-280.    
Doing this via the new glmer function may be a better way. Although the results may not be quite different, the latter has more conceptual clarity. 

Wayne
 

-----Original Message-----
From: Rubén Roa [mailto:rroa at azti.es] 
Sent: Tuesday, March 20, 2012 2:35 AM
To: Zhang,Yanwei; Bill Pikounis; r-sig-mixed-models at r-project.org
Subject: RE: [R-sig-ME] glmm with a tweedie distribution

I wouldn't call it ad hoc. The power parameter p in the variance function that defines the Tweedie family of exponential distributions, v(mu)=phi*mu^p, can be estimated via profile likelihood, and then the maximum profile likelihood estimate of the p parameter can be inserted in the glmm, essentially estimating the glmm by an estimated likelihood. So there are two stages of approximation but the approximation methods are not ad hoc, they are pretty much mainstream approximation methods to complex likelihoods. Here is a pseudo code using the tweedie package and glmmPQL from MASS (plus msm). For a published application you can see Tascheri, Saavedra-Nievas, Roa-Ureta. 2010. Statistical models to standardize catch rates in the multi-species trawl fishery for Patagonian grenadier (Macruronus magellanicus) off Southern Chile. Fisheries Research 105:200-214.

HTH

Ruben

--
Dr. Ruben H. Roa-Ureta
Senior Researcher, AZTI Tecnalia,
Marine Research Division,
Txatxarramendi Ugartea z/g, 48395, Sukarrieta,
Bizkaia, Spain

##################################  Tweedie ####################################
#estimating variance power parameter

#libraries

library(tweedie)

library(MASS)

library(msm)

MyCaseStudy.Tweedie.prof <- tweedie.profile(MyResponseVar ~ log(Y) + log(Z) + ... + Factor1 + Factor2 + ... ,
p.vec=seq(1.0,2.0,by=0.1), data=MyData, link.power=0, fit.glm=FALSE, do.smooth=TRUE, do.plot=TRUE,
method="inversion",conf.level=0.95, phi.method= "mle",verbose=TRUE)

#distributional plot

y <- rtweedie(1000, p= MyCaseStudy.Tweedie.prof $p.max, mu=1, phi= MyCaseStudy.Tweedie.prof $phi.max)
tweedie.plot(seq(0, max(y), length=100), mu=mean(y), p= MyCaseStudy.Tweedie.prof $p.max, phi= MyCaseStudy.Tweedie.prof $phi.max)

#fitting the glmm
MyCaseStudy..Tweedie.mix <- glmmPQL(fixed = MyResponseVar ~ log(Y) + log(Z) + ... + Factor1 + Factor2 + ... ,
random = list(~1|RE), family=tweedie(var.power = MyCaseStudy.Tweedie.prof$p.max, link.power=0), data= MyData)

MyCaseStudy.Tweedie.mix.sum  <- summary(MyCaseStudy.Tweedie.mix)

#estimated covariance of estimates of a subset of coefficients, [3:11]

MyCaseStudy.Tweedie.mix.year.cov <- round(deltamethod(g=list(~exp(x1),~exp(x2),~exp(x3),~exp(x4),~exp(x5),
~exp(x6),~exp(x7),~exp(x8),~exp(x9)), mean= MyCaseStudy.Tweedie.mix.sum$coef$fixed[3:11], cov=vcov(MyCaseStudy.Tweedie.mix)[3:11,3:11],
ses=FALSE),5)

MyCaseStudy.Tweedie.mix.year.se  <- sqrt(diag(MyCaseStudy.Tweedie.mix.year.cov))

MyCaseStudy.Tweedie.mix.year.cor <- cov2cor(MyCaseStudy.Tweedie.mix.year.cov)

-----Mensaje original-----
De: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] En nombre de Zhang,Yanwei
Enviado el: lunes, 19 de marzo de 2012 23:17
Para: Bill Pikounis; r-sig-mixed-models at r-project.org
Asunto: Re: [R-sig-ME] glmm with a tweedie distribution

One problem with the "glmmPQL" is that the variance function can not be estimated - you need to pre-specify it in an ad hoc way. The "cpglmm" function in the "cplm" package estimates it directly from the data along with other parameters using MLE. But of course, you can use glmmPQL to generate starting values that are fed to cpglmm.   

Regards,
Wayne 

-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Bill Pikounis
Sent: Monday, March 19, 2012 4:56 PM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] glmm with a tweedie distribution

Bryan,
You might also wish to try the glmmPQL function in Venables and Ripley's MASS package. Someone reported success with it on the SIG-ECO R list nearly a year ago:

https://stat.ethz.ch/pipermail/r-sig-ecology/2011-May/002158.html

Hope that helps.

Bill

On Mon, Mar 19, 2012 at 17:18, Ben Bolker <bbolker at gmail.com> wrote:
> Danson, Bryan <Bryan.Danson at ...> writes:
>
>>
>> Is there a way to run a GLMM with a tweedie distribution?
>>
>
>  Yes, using the 'cplm' package.  I expanded the section on ZIGLMMs on 
> http://glmm.wikidot.com/faq to include a bit more stuff on ZIGLMMs, 
> including a bit on ZIGLMMs with a continuous response for the non-zero 
> data.
>
>  Ben
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

NOTICE:  This e-mail message, including any attachments and appended messages, is for the sole use of the intended recipients and may contain confidential and legally privileged information.
If you are not the intended recipient, any review, dissemination, distribution, copying, storage or other use of all or any portion of this message is strictly prohibited.
If you received this message in error, please immediately notify the sender by reply e-mail and delete this message in its entirety.

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models