manyglm (package mvabund) giving error when using family = Gamma(link = log)

Question

I am trying to test if my three experimental factors (see example_data1.csv) affect a multivariate dataset (5 variables measured on the same experimental unit, see example_data2.csv), using the manyglm function form the mvabund package. All the data are positive, continuous numeric variables. The explanatory variables are factors. Here are the example datasets: (https://github.com/dwarton/mvabund/issues/117)

the code I am running is:

data1<-read.csv("example_data1.csv")%>%
column_to_rownames(var = "obs")%>%
as.matrix()
data2<-read.csv("example_data2.csv")

manyglm

set.seed(210)
mod1 <- manyglm(data1 ~ factor1 * factor2 * factor3,
family=Gamma(link = "log"),
data=data2
)

with family = Gamma, I get the following error:

Error in if (any(z$var.est == 0)) { :
missing value where TRUE/FALSE needed

All univariate relationships work fine:

model<-glm(variable1~ factor1 x factor2 x factor3, family=Gamma(link = "log"),data = data).

I think this could be related to the estimation as detailed here: https://github.com/dwarton/mvabund/issues/91

And I was wondering if anyone had encountered a similar issue and how they solved it.

Other distributions (negative binomial, possion) work/run but they are not the correct distribution for my data. I have also tried family = "tweedie" (which could be appropriate for my data) as this should work according to the mvabund documentation here: https://environmentalcomputing.net/statistics/glms/glm-2/, but it is not supported according to the CRAN documentation.

I have also tried the manyglm() model with each factor individually, and it seems that factor 2 is causing the issue.

Answer 1

I am not familiar with mvabund or manyglm, so you will have to check if this is doing what you want:

# Load packages
library(tidyr)
library(mvabund)


# Import data
data1<-read.csv("example_data1.csv") %>% 
  select(-obs)
data2<-read.csv("example_data2.csv")


# Process data
data1 <- mvabund(data1)
is.mvabund(data1)
data2 <- as.data.frame(data2)

#Merge 
data <- merge(data1, data2, by.x = "row.names", by.y = "row.names")


# Fit a multivariate model using mvabund with multiple response variables
set.seed(210)
mod1 <- manyglm(data1 ~ factor1 * factor2 * factor3, 
                family = Gamma(link = "log"),
                data = data
)


# See Model
summary(mod1)
plot(mod1)