# Tweedie distribution GLM for manyany() in {mvabund} package

My data follow a Tweedie distribution, and I'm working with multivariate abundances. So I'm trying to use the manyany() function in the mvabund package. manyany() fits many univariate models to multivariate abundance data, in this case I'm trying to fit Generalized Linear Models (GLMs).

Using the example data set in mvabund, I've given my code below. Although the example data set is Poisson distributed, the same error message appears for my own data set as it does for the example below.

data(spider)
abund <- spider$abund X <- as.matrix(spider$x)

## To fit a log-linear model assuming counts are poisson:
spidPois <- manyany("glm",abund,data=X,abund~X,family=poisson())

## Fitting tweedie GLMs, error message given.
spidTwee <- manyany("glm",abund,data=X,abund~X,family=tweedie(var.power=1.2))

## Error in if (xi != power) { : missing value where TRUE/FALSE needed


I set the var.power=1.2 because a value of 1 would be a Poisson. When looking at the command rtweedie() I noticed that xi is said to be a synonym for power. I'm stuck now though and unsure where to go from here.

One option would be to just use a negative binomial GLM, which I've tried and it works fine. This answer indicates that a Tweedie would be better, but I guess an extension would be how to determine how poorly a negative binomial GLM would be fitting your data.

Any ideas or thoughts would be greatly appreciated!

## Update:

I tried using gam() with manyany(), specifically with the Tweedie using tw(), as recommended.

library(mgcv)
gam <- mgcv::gam
tw <- mgcv::tw
spidTweeGAM <- manyany("gam", abund, data = X, abund~X, family=tw())

#### Error in manyany("gam", abund, data = X, abund ~ X, family = tw()) :
family argument has length more than one but not equal to the number of columns in yMat (!?)


I don't really know how to interpret this error message though! Again, any help is appreciated.

You could do a manyany() using gam() from package mgcv. From version 1.8-0 this has had a tw() family function (so you use it as family = tw etc). This will also do a search for the Tweedie parameter as part of the model fitting (I suppose this is your var.power) so you don't even need to state this value or assume it is the same for all species.
The same package has an nb family function for the negative binomial distribution with search for $\theta$, so you don't need to specify this either. You could assess the relative fitness of the two approaches using prediction error (MSE or RMSE) for example either on hold-out data or in a cross-validation process.
• A GAM as per mgcv is a penalised GLM; try it with gam() but just don't include any smooth terms (for a single species) and see if it works. I think gam() is OK without having any smooth terms at all in the model, in which case it is just a GLM. – Gavin Simpson Feb 17 '15 at 20:32
Just found out the answer. When you use manyany with a Tweedie distribution, there is an argument within the tweedie() function, but you also have to specify the var.power parameter in the manyany function.