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

Discussion:

Danson, Bryan

2012-03-19 20:06:57 UTC

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120319/b28f4b28/attachment-0001.pl>

Ben Bolker

2012-03-19 21:18:11 UTC

Permalink

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

Bill Pikounis

2012-03-19 21:56:16 UTC

Permalink

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

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

Zhang,Yanwei

2012-03-19 22:16:35 UTC

Permalink

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

Is there a way to run a GLMM with a tweedie distribution?

_______________________________________________
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.

Rubén Roa

2012-03-20 07:34:59 UTC

Permalink

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

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

Ben Bolker

2012-03-20 15:01:44 UTC

Permalink

Post by 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.

For what it's worth, the upcoming/development version of
lme4 (now on R-forge) should work with custom family arguments,
so this approach *should* be possible with glmer as well as
glmmPQL ... (but I would also definitely give a thumbs-up
to cplm, which looks quite powerful).

Anyone who wants to do some http://glmm.wikidot.com/faq - editing
is welcome ...

Zhang,Yanwei

2012-03-22 03:34:50 UTC

Permalink

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

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

Ben Bolker

2012-03-20 15:01:44 UTC

Permalink

Zhang,Yanwei

2012-03-22 03:34:50 UTC

Permalink

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

Rubén Roa

2012-03-20 07:34:59 UTC

Permalink

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

Bill Pikounis

2012-03-19 22:18:16 UTC

Permalink

And I overlooked that Ben Bolker was part of that sig-eco thread and
was the one who reported the success... Sorry for my lack of
attribution, Ben....

Bill

Post by Bill Pikounis
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
https://stat.ethz.ch/pipermail/r-sig-ecology/2011-May/002158.html
Hope that helps.
Bill

Is there a way to run a GLMM with a tweedie distribution?

Zhang,Yanwei

2012-03-19 22:16:35 UTC

Permalink

Is there a way to run a GLMM with a tweedie distribution?

Bill Pikounis

2012-03-19 22:18:16 UTC

Permalink

And I overlooked that Ben Bolker was part of that sig-eco thread and
was the one who reported the success... Sorry for my lack of
attribution, Ben....

Bill

Is there a way to run a GLMM with a tweedie distribution?

Bill Pikounis

2012-03-19 21:56:16 UTC

Permalink

Is there a way to run a GLMM with a tweedie distribution?

Danson, Bryan

2012-03-20 16:00:08 UTC

Permalink

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120320/f0203932/attachment-0001.pl>

Ben Bolker

2012-03-21 01:59:16 UTC

Permalink

Danson, Bryan <Bryan.Danson at ...> writes:
[snip]

glm.fit: algorithm did not converge

This is probably harmless -- it means that an intermediate
GLM step didn't quite work, probably because you have strongly
separated data (i.e. some places/factor combinations
etc. with all-zero or all-one data)

Does anyone have any suggestions on how to address this? My data is
fairly simple, a distance measurement for each of several trap types
on different dates.
Also, I was trying to produce a summary of the model, however I cannot get the

commands to work. I tried

Error in UseMethod("fixef") : no applicable method for 'fixef'
applied to an object of class "c('cpglmm', 'mer', 'cplm')"
summary signature(object = "cpglm")

Can we see the results of sessionInfo()? I suspect you have a
problem with methods from some packages masking others. If you have
installed lme4 from r-forge, I suspect you should re-install it from
CRAN ...

Ben

Ben Bolker

2012-03-21 01:59:16 UTC

Permalink

Danson, Bryan <Bryan.Danson at ...> writes:
[snip]

glm.fit: algorithm did not converge

commands to work. I tried

Error in UseMethod("fixef") : no applicable method for 'fixef'
applied to an object of class "c('cpglmm', 'mer', 'cplm')"
summary signature(object = "cpglm")

Danson, Bryan

2012-03-21 16:08:51 UTC

Permalink

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120321/b5db4a41/attachment-0001.pl>

Zhang,Yanwei

2012-03-22 02:58:22 UTC

Permalink

Hi Bryan,

The "cpglmm" uses the "glm.fit" function to generate initial values. So the "glm.fit does not converge" message means that when generating the initial values using GLM, the model does not converge. But this should not be a problem as long as you get a converged "cpglmm" estimate - the final estimates are independent of the initial values. I suspect if you supply initial values, this message will go away. But thank you for reporting this - I will suppress this kind of message in the new release to make it less confusing.

I believe the "summary" function does not work because you have the "nlme" package in front of the "cplm" package in the search path. If you just detach the "nlme" package, it should work.

You might want to use the Bayesian tweedie mixed models to assess the "p-values". The function "bcpglmm" does that, but the released version is not quite stable. I was using a block Metropolis update for the random effects, but this oftentimes leads to poor convergence because of the difficulty in tuning the proposal covariance matrix. In the development version, I have added an option to perform the na?ve Gibbs sampler, which proves to be much easier to tune, although at the cost of slower speed. The new version should be released within a month.

Thanks.

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 Danson, Bryan
Sent: Wednesday, March 21, 2012 11:09 AM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] glmm with a tweedie distribution

This is probably harmless -- it means that an intermediate GLM step didn't quite work,
probably because you have strongly separated data (i.e. some places/factor combinations etc.
with all-zero or all-one data)

I do have strongly separated data, as some trap types did not move at all (therefore all zeros), and others moved a lot.

Can we see the results of sessionInfo()? I suspect you have a problem with methods from some
packages masking others. If you have installed lme4 from r-forge, I suspect you should re-
install it from CRAN ...

The results from sessionInfo(model) are:

Error in as.list.default(X) :
no method for coercing this S4 class to a vector

The packages in my workspace were (installed in order):

car
ggplot2
sos
glmmADMB
lme4
nlme
plotrix
cplm

I tried removing all that had to do with GLMMs (glmmADMB, lme4, nlme, cplm) and reinstalling only 'cplm' and have been able to get the summary() function to work. This results in a list of t-values for each trap type. From my background reading, I understand that p-values are difficult to determine in GLMMs, however, I am not sure how to interpret the t-values to estimate which trap type is different from the wood traps. Here are the resulting t-values:

Estimate Std. Error t value
(Intercept) 0.05134 0.68273 0.075
trap_type5-slat -1.59469 0.44568 -3.578
trap_typevw-6 -3.79851 0.71923 -5.281
trap_typevw-sponge -2.41591 0.52667 -4.587
trap_typewire basket -27.93281 386.0848 -0.072
trap_typewire on frame -4.77199 0.90967 -5.246
trap_typewood-6 0.29933 0.32652 0.917
trap_typewood-thick 0.27437 0.34785 0.789
trap_typeyb-6 -4.27295 0.80548 -5.305
trap_typeyf-6 -2.88076 0.58270 -4.944

According to the exploratory plot, it is likely that the wood-6 and wood-thick traps are not different from the standard control. The others are probably different.

Is there a way to know for sure?

Thank you again,

Bryan

_____________________
Bryan Danson
Biological Scientist I
Fish and Wildlife Research Institute
Florida Fish and Wildlife Conservation Commission

"The significant problems we have cannot be solved at the same level of thinking with which we created them."
~Albert Einstein

[[alternative HTML version deleted]]

_______________________________________________
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.

Danson, Bryan

2012-03-22 15:40:25 UTC

Permalink

Hi Wayne,

Thank you for the help. The summary function does work once the "nlme" package is removed. I was able to get the bcpglmm() function to work as well. However, the results give a "Mean", "SD", "Na?ve SE", and "Time-series SE" for each trap type. And I am confused to how to interpret this.

I supposed I should have prefaced all of these emails with the fact that I have very little statistical background. I have taken courses in undergrad and graduate school, however they are introductory courses and fall far short of this level of mixed-modeling. I am therefore trying to teach myself this modeling. I have ordered a few books that have been mentioned on this list to help, I am just waiting on them to come in.

But with that in mind, thank you all so much for your help. I have learned a great deal so far.

Bryan

-----Original Message-----
From: Zhang,Yanwei [mailto:Yanwei.Zhang at cna.com]
Sent: Wednesday, March 21, 2012 10:58 PM
To: Danson, Bryan; r-sig-mixed-models at r-project.org
Subject: RE: Re: [R-sig-ME] glmm with a tweedie distribution

Hi Bryan,

The "cpglmm" uses the "glm.fit" function to generate initial values. So the "glm.fit does not converge" message means that when generating the initial values using GLM, the model does not converge. But this should not be a problem as long as you get a converged "cpglmm" estimate - the final estimates are independent of the initial values. I suspect if you supply initial values, this message will go away. But thank you for reporting this - I will suppress this kind of message in the new release to make it less confusing.

I believe the "summary" function does not work because you have the "nlme" package in front of the "cplm" package in the search path. If you just detach the "nlme" package, it should work.

You might want to use the Bayesian tweedie mixed models to assess the "p-values". The function "bcpglmm" does that, but the released version is not quite stable. I was using a block Metropolis update for the random effects, but this oftentimes leads to poor convergence because of the difficulty in tuning the proposal covariance matrix. In the development version, I have added an option to perform the na?ve Gibbs sampler, which proves to be much easier to tune, although at the cost of slower speed. The new version should be released within a month.

Thanks.

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 Danson, Bryan
Sent: Wednesday, March 21, 2012 11:09 AM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] glmm with a tweedie distribution

This is probably harmless -- it means that an intermediate GLM step didn't quite work,
probably because you have strongly separated data (i.e. some places/factor combinations etc.
with all-zero or all-one data)

I do have strongly separated data, as some trap types did not move at all (therefore all zeros), and others moved a lot.

Can we see the results of sessionInfo()? I suspect you have a problem
with methods from some
packages masking others. If you have installed lme4 from r-forge, I
suspect you should re-
install it from CRAN ...

The results from sessionInfo(model) are:

Error in as.list.default(X) :
no method for coercing this S4 class to a vector

The packages in my workspace were (installed in order):

car
ggplot2
sos
glmmADMB
lme4
nlme
plotrix
cplm

I tried removing all that had to do with GLMMs (glmmADMB, lme4, nlme, cplm) and reinstalling only 'cplm' and have been able to get the summary() function to work. This results in a list of t-values for each trap type. From my background reading, I understand that p-values are difficult to determine in GLMMs, however, I am not sure how to interpret the t-values to estimate which trap type is different from the wood traps. Here are the resulting t-values:

Estimate Std. Error t value
(Intercept) 0.05134 0.68273 0.075
trap_type5-slat -1.59469 0.44568 -3.578
trap_typevw-6 -3.79851 0.71923 -5.281
trap_typevw-sponge -2.41591 0.52667 -4.587
trap_typewire basket -27.93281 386.0848 -0.072
trap_typewire on frame -4.77199 0.90967 -5.246
trap_typewood-6 0.29933 0.32652 0.917
trap_typewood-thick 0.27437 0.34785 0.789
trap_typeyb-6 -4.27295 0.80548 -5.305
trap_typeyf-6 -2.88076 0.58270 -4.944

According to the exploratory plot, it is likely that the wood-6 and wood-thick traps are not different from the standard control. The others are probably different.

Is there a way to know for sure?

Thank you again,

Bryan

_____________________
Bryan Danson
Biological Scientist I
Fish and Wildlife Research Institute
Florida Fish and Wildlife Conservation Commission

"The significant problems we have cannot be solved at the same level of thinking with which we created them."
~Albert Einstein

[[alternative HTML version deleted]]

_______________________________________________
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.

Danson, Bryan

2012-03-22 15:40:25 UTC

Permalink

This is probably harmless -- it means that an intermediate GLM step didn't quite work,
probably because you have strongly separated data (i.e. some places/factor combinations etc.
with all-zero or all-one data)

I do have strongly separated data, as some trap types did not move at all (therefore all zeros), and others moved a lot.

Can we see the results of sessionInfo()? I suspect you have a problem
with methods from some
packages masking others. If you have installed lme4 from r-forge, I
suspect you should re-
install it from CRAN ...

The results from sessionInfo(model) are:

Error in as.list.default(X) :
no method for coercing this S4 class to a vector

The packages in my workspace were (installed in order):

car
ggplot2
sos
glmmADMB
lme4
nlme
plotrix
cplm

I tried removing all that had to do with GLMMs (glmmADMB, lme4, nlme, cplm) and reinstalling only 'cplm' and have been able to get the summary() function to work. This results in a list of t-values for each trap type. From my background reading, I understand that p-values are difficult to determine in GLMMs, however, I am not sure how to interpret the t-values to estimate which trap type is different from the wood traps. Here are the resulting t-values:

Estimate Std. Error t value
(Intercept) 0.05134 0.68273 0.075
trap_type5-slat -1.59469 0.44568 -3.578
trap_typevw-6 -3.79851 0.71923 -5.281
trap_typevw-sponge -2.41591 0.52667 -4.587
trap_typewire basket -27.93281 386.0848 -0.072
trap_typewire on frame -4.77199 0.90967 -5.246
trap_typewood-6 0.29933 0.32652 0.917
trap_typewood-thick 0.27437 0.34785 0.789
trap_typeyb-6 -4.27295 0.80548 -5.305
trap_typeyf-6 -2.88076 0.58270 -4.944

According to the exploratory plot, it is likely that the wood-6 and wood-thick traps are not different from the standard control. The others are probably different.

Is there a way to know for sure?

Thank you again,

Bryan

_____________________
Bryan Danson
Biological Scientist I
Fish and Wildlife Research Institute
Florida Fish and Wildlife Conservation Commission

"The significant problems we have cannot be solved at the same level of thinking with which we created them."
~Albert Einstein

[[alternative HTML version deleted]]

_______________________________________________
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.

Zhang,Yanwei

2012-03-22 02:58:22 UTC

Permalink

This is probably harmless -- it means that an intermediate GLM step didn't quite work,
probably because you have strongly separated data (i.e. some places/factor combinations etc.
with all-zero or all-one data)

I do have strongly separated data, as some trap types did not move at all (therefore all zeros), and others moved a lot.

Can we see the results of sessionInfo()? I suspect you have a problem with methods from some
packages masking others. If you have installed lme4 from r-forge, I suspect you should re-
install it from CRAN ...

Danson, Bryan

2012-03-19 20:06:57 UTC

Permalink

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120319/b28f4b28/attachment-0002.pl>

Ben Bolker

2012-03-19 21:18:11 UTC

Permalink

Is there a way to run a GLMM with a tweedie distribution?

Danson, Bryan

2012-03-20 16:00:08 UTC

Permalink

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120320/f0203932/attachment-0002.pl>

Danson, Bryan

2012-03-21 16:08:51 UTC

Permalink

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20120321/b5db4a41/attachment-0002.pl>