Function manyglm() on species abundance and trait data

Question

I am attempting to run a glm analysis using the following code:

ftP <- manyglm(AbundVec ~ as.matrix(traitVec)*trV, family = "poisson")

I currently receive the following error message:

Error in [[<-.data.frame(*tmp*, i, value = c(54L, 54L, 54L, 54L, 54L, : replacement has 60326 rows, data has 8618

As a first step, I create a species x plot matrix for the abundance data based on:

AbundMat <- acast(Abund, plot~species, value.var="occurrence", fun.aggregate=length)

Output: class(AbundMat)=matrix, typeof(AbundMat)=integer, storage.mode(AbundMat)=integer, length(AbundMat)=8618, attributes(AbundMat)= 139(plots) 62(species)

I then import files:

1) Traits: data.frame with 62obs. and 7variables (includes Factors and num); includes ; typeof(Traits)=list and class(Traits)=df

2) Treatment: data.frame with 139obs. and 2variables ($plot and $site as Factors); typeof(Treatment)=list and class(Treatment)=df

I then vectorize the abundance data following an example described in the appendix of the paper: CATS regression - a model-based approach to studying trait-based community assembly (http://onlinelibrary.wiley.com/store/10.1111/2041-210X.12280/asset/supinfo/mee312280-sup-0003-AppendixS3.pdf?v=1&s=0037b03799e63903d896ba208fb2a6cc1b5605d0):

AbundVec <- as.vector(AbundMat)
    n.rows <- NROW(AbundMat)
    n.vars <- NCOL(AbundMat)
    sppVec <- rep(row.names(Traits), each=n.rows)
    traitVec <- Traits[rep(1:n.vars, each = n.rows),]
    siteVec <- rep(row.names(AbundMat), n.vars)
    Treatment <- as.factor(Treatment)
    levels(Treatment) <- c("A", "B", "C")
    trV <- rep(Treatment, n.vars)

Finally, I execute a glm (requires package mvabund):

ftP <- manyglm(AbundVec ~ as.matrix(traitVec)*trV, family = "poisson")

A this stage I get the error message mentioned above: Error in [[<-.data.frame(*tmp*, i, value = c(54L, 54L, 54L, 54L, 54L, : replacement has 60326 rows, data has 8618

Please note that 8618*7=60326 (I have 7 trait variables)

When I compare my data files to the example data file for abundance described in the CATS regression paper, I notice that the example data file has the following output for str():

str(GrazedPlants$abundMat) num [1:32, 1:68] 0 210 195 1 617 0 1 75 1 0 ... - attr(, "dimnames")=List of 2 ..$ : chr [1:32] "FC01" "FC02" "FC03" "FC04" ... ..$ : atomic [1:68] CAPSBURS CREPVESI.HAE GALICORR MUSCNEGL ... .. ..- attr(, "levels")= chr [1:83] "ALYSALYS" "ANTHVULN" "APHYMONS" "ARENSERP" ...

→ typeof = double, class = matrix

When I run str() on my data I instead get:

str(AbundMat) int [1:139, 1:62] 0 0 0 0 0 0 0 0 0 1 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:139] "A01" "A02" "A03" "A04" ... ..$ : chr [1:62] "sp01" "sp02" "sp03" "sp04" ...

→ typeof = integer, class = matrix

The format of my data for Treatment and Traits appears to match that of the example dataset except that most of my traits are factors (while in the example they are all numerical).

So, I wonder whether the error message is due to 'AbundMat' not being of the required form class=double (but integer instead)?

Also, the attributes(levels) seem to be missing in 'AbundMat' according to the str() output.

Would it be possible to produce the matrix to be of the type double (to resemble the example data matrix)?

============ Additional note: ==========================

As an additional note, this code can be used to coerce data to use with function manyglm(). It uses just one row of data (a single plot). I wonder how this code can be adjusted for multiple sites, i.e. the 139 plots in my dataset?

y2 = c(t(grazedPlants$abundMat[1,])) # coercing the first row into a vector
y2
singleSite2 = data.frame(y2, grazedPlants$traits, check.names = F) # make a data frame from site abundances and traits
singleSite2
ft5 = manyglm(y2 ~ ., data = singleSite2, family = "negative.binomial")
ft5
summary(ft5, nBoot=1) # find which traits are important
plot(ft5)

====================================================

Many thanks in advance for any advice.

These are my datasets:

1) species x plot abundance (count) data (file name: Abund)

structure(list(plot = structure(c(47L, 47L, 48L, 51L, 53L, 54L), .Label = c("A01", "A02", "A03", "A04", "A05", "A06", "A07", "A08", "A09", "A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18", "A19", "A20", "A21", "A22", "A23", "A24", "A25", "A26", "A27", "A28", "A29", "A30", "A31", "A32", "A33", "A34", "A35", "A36", "A37", "A38", "A39", "A40", "A41", "A42", "A43", "A44", "A45", "A46", "A47", "A48", "A49", "B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B09", "B10", "B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20", "B21", "B22", "B23", "B24", "B25", "B26", "B27", "B28", "B29", "B30", "B31", "B32", "B33", "B34", "B35", "B36", "B37", "B38", "B39", "B40", "B41", "B42", "B43", "B44", "B45", "C01", "C02", "C03", "C04", "C05", "C06", "C07", "C08", "C09", "C10", "C11", "C12", "C13", "C14", "C15","C16", "C17", "C18", "C19", "C20", "C21", "C22", "C23", "C24", "C25", "C26", "C27", "C28", "C29", "C30", "C31", "C32", "C33", "C34", "C35", "C36", "C37", "C38", "C39", "C40", "C41", "C42", "C43", "C44", "C45"), class = "factor"), 
species = structure(c(46L, 47L, 47L, 47L, 47L, 47L), .Label = c("sp1", 
"sp10", "sp11", "sp12", "sp13", "sp14", "sp15", "sp16", "sp17", 
"sp18", "sp19", "sp2", "sp20", "sp21", "sp22", "sp23", "sp24", 
"sp25", "sp26", "sp27", "sp28", "sp29", "sp3", "sp30", "sp31", 
"sp32", "sp33", "sp34", "sp35", "sp36", "sp37", "sp38", "sp39", 
"sp4", "sp40", "sp41", "sp42", "sp43", "sp44", "sp45", "sp46", 
"sp47", "sp48", "sp49", "sp5", "sp50", "sp51", "sp52", "sp53", 
"sp54", "sp55", "sp56", "sp57", "sp58", "sp59", "sp6", "sp60", 
"sp61", "sp62", "sp7", "sp8", "sp9"), class = "factor"), 
occurrence = c(1L, 1L, 1L, 1L, 1L, 1L)), .Names = c("plot", "species", "occurrence"), row.names = c(NA, 6L), class = "data.frame")

2) species x trait data (file name: Traits)

structure(list(species = structure(c(1L, 12L, 23L, 34L, 45L, 56L, 60L, 61L, 62L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 57L, 58L, 59L), .Label = c("sp1", "sp10", "sp11", "sp12", "sp13", "sp14", "sp15", "sp16", "sp17", "sp18", "sp19", "sp2", "sp20", "sp21", "sp22", "sp23", "sp24", "sp25", "sp26", "sp27", "sp28", "sp29", "sp3", "sp30", "sp31", "sp32", "sp33", "sp34", "sp35", "sp36", "sp37", "sp38", "sp39", "sp4", "sp40", "sp41", "sp42", "sp43", "sp44", "sp45", "sp46", "sp47", "sp48", "sp49", "sp5", "sp50", "sp51", "sp52", "sp53", "sp54", "sp55", "sp56", "sp57", "sp58", "sp59", "sp6", "sp60", "sp61", "sp62", "sp7", "sp8", "sp9"), class = "factor"), trait1 = structure(c(3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 1L, 3L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 3L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 3L, 2L, 2L, 3L, 3L, 3L, 3L, 2L, 3L, 2L, 2L, 2L, 3L, 3L, 2L, 3L), .Label = c("<NA>", "a", "b"), class = "factor"), trait2 = structure(c(6L, 6L, 4L, 3L, 8L, 4L, 4L, 4L, 3L, 1L, 1L, 6L, 6L, 6L, 6L, 1L, 2L, 2L, 7L, 7L, 6L, 3L, 6L, 6L, 2L, 1L, 2L, 4L, 4L, 3L, 5L, 4L, 4L, 1L, 4L, 2L, 4L, 6L, 3L, 8L, 2L, 4L, 4L, 1L, 1L, 4L, 6L, 1L, 4L, 6L, 6L, 6L, 6L, 5L, 3L, 2L, 4L, 1L, 4L, 1L, 4L, 6L), .Label = c("<NA>", "c", "d", "e", "f", "g", "h", "i"), class = "factor"), trait3 = structure(c(4L, 4L, 4L, 6L, 4L, 1L, 6L, 2L, 2L, 2L, 2L, 1L, 6L, 4L, 4L, 6L, 3L, 6L, 3L, 3L, 6L, 2L, 2L, 2L, 6L, 2L, 6L, 4L, 4L, 4L, 6L, 2L, 6L, 2L, 6L, 2L, 2L, 2L, 6L, 6L, 6L, 2L, 6L, 6L, 2L, 4L, 6L, 2L, 4L, 6L, 6L, 4L, 5L, 4L, 4L, 2L, 2L, 2L, 6L, 2L, 6L, 6L), .Label = c("", "<NA>", "j", "k", "l", "m"), class = "factor"), trait4 = structure(c(4L, 4L, 2L, 2L, 5L, 3L, 2L, 3L, 1L, 5L, 1L, 3L, 3L, 3L, 3L, 5L, 1L, 4L, 2L, 2L, 1L, 5L, 1L, 3L, 3L, 1L, 5L, 3L, 1L, 5L, 5L, 2L, 2L, 1L, 1L, 3L, 1L, 3L, 5L, 2L, 4L, 3L, 1L, 1L, 1L, 4L, 3L, 5L, 3L, 2L, 2L, 2L, 3L, 5L, 5L, 4L, 3L, 1L, 2L, 1L, 2L, 4L), .Label = c("<NA>", "1", "2", "3", "4"), class = "factor"), trait5 = structure(c(12L, 14L, 9L, 2L, 21L, 9L, 1L, 14L, 9L, 22L, 2L, 12L, 12L, 4L, 12L, 14L, 6L, 14L, 12L, 11L, 12L, 3L, 1L, 1L, 3L, 1L, 2L, 12L, 5L, 12L, 2L, 17L, 2L, 1L, 16L, 19L, 15L, 19L, 12L, 10L, 3L, 12L, 18L, 12L, 19L, 8L, 12L, 20L, 16L, 12L, 7L, 12L, 13L, 9L, 11L, 5L, 2L, 1L, 14L, 11L, 12L, 16L), .Label = c("<NA>", "10", "12", "14", "15", "16", "17", "18", "20", "23", "25", "30", "32", "35", "4", "40", "50", "6", "7", "8", "9", "small"), class = "factor"), 
trait6 = structure(c(3L, 3L, 3L, 2L, 3L, 3L, 3L, 3L, 3L, 
1L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 
3L, 1L, 3L, 3L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 1L, 3L, 3L, 
2L, 2L, 3L, 2L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 1L, 3L, 3L, 3L, 3L), .Label = c("<NA>", "p", 
"q"), class = "factor"), trait7 = c(0.61, 0.64, 0.61, 0.62, 
0.64, 0.61, 0.39, 0.13, 0.73, 0.75, 0.51, 0.77, 0.75, 0.6, 
0.88, 0.55, 0.61, 0.72, 0.64, 0.72, 0.65, 0.7, 0.85, 0.74, 
0.55, 0.76, 0.76, 0.43, 0.67, 0.48, 0.73, 0.88, 0.4, 0.46, 
0.48, 0.67, 0.67, 0.88, 0.53, 0.48, 0.93, 0.64, 0.62, 0.47, 
0.64, 0.64, 0.7, 0.64, 0.9, 0.85, 0.7, 0.86, 0.72, 0.43, 
0.33, 0.74, 0.52, 0.69, 0.75, 0.45, 0.88, 0.68)), .Names = c("species", "trait1", "trait2", "trait3", "trait4", "trait5", "trait6", "trait7"), class = "data.frame", row.names = c(NA, -62L))

3) treatment data (file name: Treatment)

    structure(list(plot = structure(1:139, .Label = c("A01", "A02", "A03", "A04", "A05", "A06", "A07", "A08", "A09", "A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18", "A19", "A20", "A21", "A22", "A23", "A24", "A25", "A26", "A27", "A28", "A29", "A30", "A31", "A32", "A33", "A34", "A35", "A36", "A37", "A38", "A39", "A40", "A41", "A42", "A43", "A44", "A45", "A46", "A47", "A48", "A49", "B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B09", "B10", "B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20", "B21", "B22", "B23", "B24", "B25", "B26", "B27", "B28", "B29", "B30", "B31", "B32", "B33", "B34", "B35", "B36", "B37", "B38", "B39", "B40", "B41", "B42", "B43", "B44", "B45", "C01", "C02", "C03", "C04", "C05", "C06", "C07", "C08", "C09", "C10", "C11", "C12", "C13", "C14", "C15", "C16", "C17", "C18", "C19", "C20", "C21", "C22", "C23", "C24", "C25", "C26", "C27", "C28", "C29", "C30", "C31", "C32", "C33", "C34", "C35", "C36", "C37", "C38", "C39", "C40", "C41", "C42", "C43", "C44", "C45"), class = "factor"), site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", "C"), class = "factor")), .Names = c("plot", "site"), class = "data.frame", row.names = c(NA, -139L))

In addition, these are subsets of the 'grazedPlants.RData' file:

dput(head(grazedPlants$abundMat))

    structure(list(X = structure(1:6, .Label = c("FC01", "FC02", 
    "FC03", "FC04", "FC05", "FC06", "FC07", "FC08", "NC01", "NC02", 
    "NC03", "NC04", "NC05", "NC06", "NC07", "NC08", "ND01", "ND02", 
    "ND03", "ND04", "ND05", "ND06", "ND07", "ND08", "TD01", "TD02", 
    "TD03", "TD04", "TD05", "TD06", "TD07", "TD08"), class = "factor"), 
        CAPSBURS = c(0L, 210L, 195L, 1L, 617L, 0L), CREPVESI.HAE = c(0L, 
        0L, 0L, 0L, 0L, 0L), GALICORR = c(0L, 0L, 0L, 0L, 0L, 0L), 
        MUSCNEGL = c(0L, 0L, 0L, 0L, 0L, 0L), BROMEREC = c(1L, 0L, 
        0L, 0L, 0L, 0L), ARENSERP = c(0L, 0L, 1L, 0L, 4L, 0L), VULPMYUR = c(0L, 
        28L, 3L, 0L, 0L, 8L), CERAPUMI = c(0L, 2L, 0L, 0L, 0L, 0L
        ), ORNICOLL = c(0L, 1L, 0L, 0L, 0L, 1L), VEROARVE = c(0L, 
        0L, 0L, 0L, 0L, 0L), GERAROTU = c(0L, 0L, 0L, 80L, 0L, 41L
        ), ERODCICU = c(0L, 0L, 0L, 1L, 0L, 3L), BROMHORD = c(0L, 
        735L, 0L, 72L, 0L, 4L), POABULB = c(0L, 209L, 110L, 2L, 0L, 
        0L), HORDMURI = c(5L, 11L, 0L, 0L, 0L, 0L), TARALAEV = c(0L, 
        0L, 0L, 0L, 0L, 0L), SANGMINO = c(0L, 0L, 0L, 0L, 1L, 0L), 
        TRIFSCAB = c(0L, 1L, 25L, 0L, 0L, 0L), CREPSANC.SAN = c(0L, 
        0L, 0L, 0L, 0L, 0L), ERYNCAMP = c(1L, 0L, 0L, 0L, 0L, 0L), 
        BROMDIAN = c(0L, 0L, 0L, 0L, 0L, 1L), STELMEDI = c(0L, 0L, 
        0L, 0L, 0L, 0L), ALYSALYS = c(0L, 0L, 0L, 0L, 0L, 0L), CERAGLOM = c(17L, 
        5L, 2L, 0L, 21L, 0L), EROPVERN = c(0L, 0L, 0L, 0L, 0L, 0L
        ), FESTRUBR = c(30L, 0L, 14L, 0L, 69L, 0L), FILAPYRA = c(0L, 
        0L, 9L, 0L, 0L, 0L), SCLEANNU = c(0L, 0L, 0L, 0L, 0L, 0L), 
        CAREHALL = c(80L, 21L, 0L, 2L, 0L, 368L), POTENEUM = c(2L, 
        1L, 0L, 0L, 0L, 0L), VICITETR.GRA = c(0L, 0L, 0L, 1L, 0L, 
        0L), SHERARVE = c(0L, 0L, 0L, 0L, 0L, 0L), POLYCALC = c(0L, 
        0L, 0L, 2L, 0L, 0L), VICISATI.SAT = c(0L, 0L, 0L, 0L, 0L, 
        0L), FILIVULG = c(0L, 0L, 0L, 126L, 16L, 0L), MYOSRAMO.RAM = c(0L, 
        56L, 0L, 0L, 0L, 27L), CAREHUMI = c(0L, 3L, 9L, 0L, 0L, 7L
        ), TARAOFFI = c(2L, 0L, 0L, 0L, 0L, 0L), HIPPCOMO = c(1L, 
        1L, 0L, 3L, 1L, 0L), FESTCHRI = c(0L, 0L, 0L, 0L, 13L, 0L
        ), KOELVALL = c(0L, 0L, 0L, 0L, 0L, 0L), ONONSTRI = c(0L, 
        0L, 0L, 0L, 0L, 0L), GENIHISP = c(0L, 0L, 0L, 0L, 0L, 0L), 
        CARTMITI = c(0L, 0L, 0L, 0L, 0L, 0L), ANTHVULN = c(0L, 0L, 
        0L, 0L, 0L, 0L), CENTPECT.SUP = c(0L, 0L, 0L, 0L, 0L, 0L), 
        LOTUCORN = c(0L, 0L, 0L, 0L, 0L, 1L), PLANLANC = c(0L, 0L, 
        0L, 1L, 0L, 0L), RANUBULB = c(0L, 0L, 0L, 0L, 15L, 0L), CHAMSAGI = c(0L, 
        0L, 0L, 0L, 0L, 0L), TEUCPOLI = c(0L, 0L, 0L, 0L, 0L, 0L), 
        THYMSERP = c(0L, 0L, 0L, 0L, 0L, 0L), COROMINI = c(0L, 0L, 
        0L, 0L, 0L, 0L), SALVPRAT = c(0L, 0L, 0L, 2L, 0L, 0L), ASTRMONS = c(0L, 
        0L, 0L, 0L, 0L, 0L), HIERPILO = c(1L, 0L, 0L, 0L, 0L, 0L), 
        LINUTENU.TEN = c(0L, 0L, 0L, 0L, 0L, 0L), APHYMONS = c(0L, 
        0L, 2L, 0L, 0L, 0L), TEUCCHAM = c(0L, 0L, 1L, 0L, 0L, 0L), 
        HELICANU = c(0L, 0L, 1L, 0L, 16L, 6L), HELINUMM = c(2L, 0L, 
        0L, 0L, 0L, 0L), CAREFLAC = c(18L, 0L, 0L, 57L, 0L, 8L), 
        STIPPENN = c(125L, 26L, 0L, 0L, 0L, 0L), INULMONT = c(0L, 
        0L, 1L, 0L, 0L, 0L), TEUCMONT = c(0L, 0L, 0L, 0L, 0L, 0L), 
        GLOBVULG = c(0L, 0L, 0L, 0L, 1L, 0L), HELIAPEN = c(0L, 1L, 
        0L, 0L, 0L, 0L), RANUGRAM = c(0L, 0L, 0L, 2L, 0L, 0L)), .Names = c("X", 
    "CAPSBURS", "CREPVESI.HAE", "GALICORR", "MUSCNEGL", "BROMEREC", 
    "ARENSERP", "VULPMYUR", "CERAPUMI", "ORNICOLL", "VEROARVE", "GERAROTU", 
    "ERODCICU", "BROMHORD", "POABULB", "HORDMURI", "TARALAEV", "SANGMINO", 
    "TRIFSCAB", "CREPSANC.SAN", "ERYNCAMP", "BROMDIAN", "STELMEDI", 
    "ALYSALYS", "CERAGLOM", "EROPVERN", "FESTRUBR", "FILAPYRA", "SCLEANNU", 
    "CAREHALL", "POTENEUM", "VICITETR.GRA", "SHERARVE", "POLYCALC", 
    "VICISATI.SAT", "FILIVULG", "MYOSRAMO.RAM", "CAREHUMI", "TARAOFFI", 
    "HIPPCOMO", "FESTCHRI", "KOELVALL", "ONONSTRI", "GENIHISP", "CARTMITI", 
    "ANTHVULN", "CENTPECT.SUP", "LOTUCORN", "PLANLANC", "RANUBULB", 
    "CHAMSAGI", "TEUCPOLI", "THYMSERP", "COROMINI", "SALVPRAT", "ASTRMONS", 
    "HIERPILO", "LINUTENU.TEN", "APHYMONS", "TEUCCHAM", "HELICANU", 
    "HELINUMM", "CAREFLAC", "STIPPENN", "INULMONT", "TEUCMONT", "GLOBVULG", 
    "HELIAPEN", "RANUGRAM"), row.names = c(NA, 6L), class = "data.frame")

dput(head(grazedPlants$traits))

    structure(list(X = structure(c(9L, 20L, 28L, 42L, 7L, 4L), .Label = c("ALYSALYS", "ANTHVULN", "APHYMONS", "ARENSERP", "ASTRMONS", "BROMDIAN", "BROMEREC", "BROMHORD", "CAPSBURS", "CAREFLAC", "CAREHALL", "CAREHUMI", "CARTMITI", "CENTPECT.SUP", "CERAGLOM", "CERAPUMI", "CHAMSAGI", "COROMINI", "CREPSANC.SAN", "CREPVESI.HAE", "ERODCICU", "EROPVERN", "ERYNCAMP", "FESTCHRI", "FESTRUBR", "FILAPYRA", "FILIVULG", "GALICORR", "GENIHISP", "GERAROTU", "GLOBVULG", "HELIAPEN", "HELICANU", "HELINUMM", "HIERPILO", "HIPPCOMO", "HORDMURI", "INULMONT", "KOELVALL", "LINUTENU.TEN", "LOTUCORN", "MUSCNEGL", "MYOSRAMO.RAM", "ONONSTRI", "ORNICOLL", "PLANLANC", "POABULB", "POLYCALC", "POTENEUM", "RANUBULB", "RANUGRAM", "SALVPRAT", "SANGMINO", "SCLEANNU", "SHERARVE", "STELMEDI", "STIPPENN", "TARALAEV", "TARAOFFI", "TEUCCHAM", "TEUCMONT", "TEUCPOLI", "THYMSERP", "TRIFSCAB", "VEROARVE", "VICISATI.SAT", "VICITETR.GRA", "VULPMYUR"), class = "factor"), SM = c(0.587, -0.81, 0.446, 0.195, 0.357, -0.128), OFL = c(136.1359, 131.9466, 148.8626, 148.5717, 156.9126, 149.2011), RPH = c(1.4092797, 1.0612845, 0.729711, 0.7335391, 1.7941586, 1.4624131), VPH = c(1.12952674, 0.86161103, 0.54418607, 0.44006594, 1.49088504, 1.01091915), SLA = c(1.0997531, 1.0612576, 1.4866071, 1.0666007, 0.9769541, 1.1077584), LDMC = c(2.491384, 2.446196, 2.227827, 2.476604, 2.717145, 2.347065), LA = c(0.58658981, 0.29217792, 0.53604022, 0.07529993, 0.28781845, 0.68114674), 
    LNC = c(1.46497, 1.31874, 1.417831, 1.277482, 1.302021, 1.363778
    ), LCC = c(2.632749, 2.621413, 2.6114, 2.653969, 2.637665, 
    2.630081), LPC = c(0.40060595, 0.305709643, 0.543668171, 
    0.225178098, 0.135509475, 0.354213077)), .Names = c("X", "SM", "OFL", "RPH", "VPH", "SLA", "LDMC", "LA", "LNC", "LCC", "LPC"), row.names = c(NA, 6L), class = "data.frame"))

dput(head(grazedPlants$treatment))

structure(c(1L, 1L, 1L, 1L, 1L, 1L), .Names = c("FC01", "FC02", 
"FC03", "FC04", "FC05", "FC06"), .Label = c("Fertilised (Limestone)", 
"Pasture (Limestone)", "Pasture (Dolomite)", "Ungrazed (Dolomite)"
), class = "factor")

You seem to be following structure to the letter. It must be input data is not consistent. Can you post the dput or str of the package's grazedPlants.RData, traits, and treatment? We need to compare with yours. — Parfait
I also realised that ftP <- manyglm(AbundVec ~ as.matrix(traitVec[1])*trV, family = "poisson") works for my dataset, so only one trait column rather than all seven. Would this suggest that there is a dataframe issue? — tabtimm
Can you dump your (not author's) full (not head()) versions of Abund, Traits, and Treatment? By cutting off with head() I cannot reproduce your error. And do not send Dropbox links, try dumping in pastebin.com. Even better if you can dump into a Github Gist like I did for author's objects. — Parfait
Ok. I have dumped the dput() outputs for Traits and Treatment. dput() of Abund is a very large file that exceeds the maximum number of characters allowed. I have dumped it here: pastebin.com/07xkP2UG . Does this work? — tabtimm

Parfait Parfait · Accepted Answer · 2017-10-03T16:00:25

Using your recently posted objects here (Abund, Traits, and Treatment) and comparing them with posted link's objects (abund, traits, and treatment), your fundamental issue is Traits does not contain numeric columns. Other issues include rownames() and $ assignment.

Author's traits (numeric columns)

head(traits)
#                      SM      OFL       RPH       VPH       SLA     LDMC         LA      LNC      LCC       LPC
# CAPSBURS      0.5872967 136.1359 1.4092797 1.1295267 1.0997531 2.491384 0.58658981 1.464970 2.632749 0.4006060
# CREPVESI.HAE -0.8103945 131.9466 1.0612845 0.8616110 1.0612576 2.446196 0.29217792 1.318740 2.621413 0.3057096
# GALICORR      0.4463238 148.8626 0.7297110 0.5441861 1.4866071 2.227827 0.53604022 1.417831 2.611400 0.5436682
# MUSCNEGL      0.1951452 148.5717 0.7335391 0.4400659 1.0666007 2.476604 0.07529993 1.277482 2.653969 0.2251781
# BROMEREC      0.3571080 156.9126 1.7941586 1.4908850 0.9769541 2.717145 0.28781845 1.302021 2.637665 0.1355095
# ARENSERP     -0.1283229 149.2011 1.4624131 1.0109192 1.1077584 2.347065 0.68114674 1.363778 2.630081 0.3542131

Your Traits (factor columns + species which should be rownames)

head(Traits)
#   species trait1 trait2 trait3 trait4 trait5 trait6 trait7
# 1     sp1      b      g      k      3     30      q   0.61
# 2     sp2      b      g      k      3     35      q   0.64
# 3     sp3      a      e      k      1     20      q   0.61
# 4     sp4      a      d      m      1     10      p   0.62
# 5     sp5      a      i      k      4      9      q   0.64
# 6     sp6      a      e             2     20      q   0.61

Adjusting for rownames and converting the factor types to numeric, you should then be able to call manyglm(). If such conversion does not make sense as used here, be sure to adjust Traits to reflect numeric figures.

row.names(Traits) <- Traits$species          # ASSIGNING ROW NAMES TO SPECIES
Traits <- Traits[-1]                         # REMOVING FIRST COLUMN  
# CONVERT COLUMNS TO NUMERIC
Traits[, c(1:ncol(Traits))] <- sapply(Traits[, c(1:ncol(Traits))], as.numeric)

AbundVec <- as.vector(as.numeric(AbundMat))
N.rows <- NROW(AbundMat)
N.vars <- NCOL(AbundMat)
SppVec <- rep(row.names(Traits), each=N.rows)
TraitVec <- Traits[rep(1:N.vars, each=N.rows),]
SiteVec <- rep(row.names(AbundMat), N.vars)
Treatment = as.factor(Treatment$site)        # SELECT SPECIFIC COLUMN
levels(Treatment) <- c("A", "B", "C")       
TrV <- rep(Treatment, N.vars)

FtP <- manyglm(AbundVec ~ as.matrix(TraitVec)*TrV, family = "poisson")

Function manyglm() on species abundance and trait data

============ Additional note: ==========================

====================================================

1 Answers