dwarton / mvabund

mvabund updates
10 stars 14 forks source link

binstruc + binomial manyglm matrix size error #66

Open Stu-dem opened 6 years ago

Stu-dem commented 6 years ago

Hello!

I am wondering if this is an issue under development or an error. I have used the binstruc function to build a dataframe for a binomial manyglm to calculate the proportional rather than discrete composition. The new binomial structured matrix builds successfully and manyglm runs successfully with family = "negative.binomial" but with family = "binomial" I get this fatal error:

gsl: ../../gsl-2.1/matrix/copy_source.c:31: ERROR: matrix sizes are different
Default GSL error handler invoked.

I have produced a (hopefully) reproducible example using the spider dataset:

library(mvabund)

data("spider")

abund <- spider$abund

abund$Total <- rowSums(abund) #calculate row totals for total number of individuals
abund_fails <- abund[,13] - abund[1:12] #calculate failures as difference between row totals and individual species abundance 
abund <- abund[1:12] #cut out total column
X <- as.data.frame(spider$x) #store environmental vars

##build binstruc function (copied from https://github.com/cran/mvabund/blob/master/R/binstruc.R)
binstruc <- function(y1, y2, order =1){

  if(any(y1<0) | any(y2<0)) stop("binomial data must consist of integers only")

  y1 <- as.matrix(y1)

  if(nrow(y1)!=NROW(y2)) stop("'y1' must have the same number of rows as 'y2'")

  if(order==1){

    nvars <- ncol(y1)
    y2 <- as.matrix(y2)
    if(is.null(colnames(y2))) { colnames(y2) <- colnames(y1)
    } else if(is.null(colnames(y1))) colnames(y1) <- colnames(y2)
    if(nvars!=ncol(y2)) stop("'y1' must have the same number of columns as 'y2'")
    y <- cbind(y1,y2)
    colnam <- colnames(y)

    if(is.null(colnam)) {
      colnames(y) <- c(paste("succ",1:nvars, sep=""), paste("fail",1:nvars, sep=""))
      struc <-  c( rep("succ", nvars), rep("fail", nvars) )
    } else {
      struc <- c( rep("succ",times = nvars), rep("fail", times=nvars) )
      colnames(y) <- paste(struc , ".", colnam, sep = "")
    }

    n <- y[ ,  struc=="succ", drop=FALSE ] + y[, struc=="fail", drop=FALSE]
    if(any(n - n[,1] != 0))
      stop("the number of trials must be the same for each variable")

  } else  if(order==2){

    nvars <- ncol(y1)
    if(NCOL(y2)==1){
      y2 <- c(y2)
    } else if(nvars!=ncol(y2)) {
      stop("'y2' must either have the same number of columns as 'y1' or be a vector")
    } else if(any(y2- y2[,1]!=0))
      stop("the number of trials must be the same for each variable")

    if(any(y2-y1<0)) stop("the number of trials must be larger than the number of successses")
    y <- cbind(y1,y2-y1)
    colnam <- colnames(y)

    if(is.null(colnam)) {
      colnames(y) <- c(paste("succ",1:nvars, sep=""), paste("fail",1:nvars, sep=""))
    } else {
      struc <- c( rep("succ",times = nvars), rep("fail", times=nvars) )
      colnames(y) <- paste(struc , ".", colnam, sep = "")
    }

  } else  stop("'order' must be in either '1' or '2'")

  return(y)

}

abund_binstruc <- binstruc(abund, abund_fails, 1) #build successes-failures matrix
abund_mvabund <- mvabund(abund_binstruc) #convert to mvabund object

abund_manyglm <- manyglm(abund_mvabund ~ soil.dry, data = X, family = "binomial") #run manyglm

This functionality looks incredibly promising! But if this is still in development could you please suggest another way of computing proportional compositions?