# Solved – 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!

Contents

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