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**hide

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

#### Best Answer

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.

### Similar Posts:

- Solved – Gamma vs tweedie distribution for large productivity dataset
- Solved – Does the dependent variable in a GLM have to be transformed before running the model or does the model do it
- Solved – How to calculate the Tweedie prediction based on model coefficients
- Solved – GLM Tweedie dispersion parameter
- Solved – Fitting a glm to a zero inflated positive continuous response