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

[Rubén Roa]

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