Code Monkey home page Code Monkey logo

mvabund's Introduction

mvabund

License CRAN Downloads Codecov test coverage

The goal of mvabund is to provide tools for a model-based approach to the analysis of multivariate abundance data in ecology (Yi Wang et al. 2012), in particular, testing hypotheses about the community-environment association. Abundance measures include counts, presence/absence data, ordinal or biomass data.

This package includes functions for visualising data, fitting predictive models, checking model assumptions, as well as testing hypotheses about the community–environment association.

Installation

mvabund is available on CRAN and can be installed directly in R:

install.packages("mvabund")

library(mvabund)

This is an archived version -- you can access the development version of mvabund at ecostats/mvabund.

mvabund's People

Contributors

dwarton avatar aliceyiwang avatar eddelbuettel avatar johnwilshire avatar julianfbyrnes avatar fontikar avatar gordy2x avatar olivroy avatar

Stargazers

Kyle Austin Gervers avatar Victor Tsang avatar avatar Skip Woolley avatar Joshua avatar Andreas Scharmüller avatar avatar Russell Dinnage avatar Matthias Grenié avatar Eduard Szöcs avatar

Watchers

avatar avatar Russell Dinnage avatar avatar avatar avatar avatar

Forkers

eveslavich julianfbyrnes johnwilshire niveaaaa yu-yang2-unsw-edu-au mattocci27 humburg joshualiuxu nico-mbio ecologysarah fontikar gerverska mmcgillycuddy fhui28

mvabund's Issues

Error message in summary.manyglm: length of 'dimnames' [1] not equal to array extent

I get this error message when I use summary() following model fitting with manyglm. When p.uni="unadjusted" is not included, summary() works fine, but whenever p.uni= is included, error message appears. I am using R version 3.2.4 and just installed mvabund via install.packages function. If you could find out what I am doing wrong, that would be very helpful. Below are my data, code, and the error message. Thank you very much.
Sincerely,
Kiyoshi

mydata <- structure(list(Area = structure(1:12, .Label = c("CTV", "Falconbridge",
"Franco", "Hwy144", "Kukagami", "Landfill", "Laurentian", "Nepewassi Lake",
"ONeil", "Sand Bay Rd", "Tilton Lake", "Trout Lake"), class = "factor"),
Rank = structure(c(4L, 4L, 2L, 3L, 2L, 4L, 3L, 1L, 3L, 1L,
2L, 1L), .Label = c("Reference", "Moderate", "Semibarren",
"Barren"), class = "factor"), ANAM = c(1L, 4L, 0L, 1L, 19L,
1L, 4L, 88L, 9L, 41L, 6L, 21L), HYVE = c(10L, 7L, 12L, 2L,
3L, 7L, 6L, 26L, 4L, 9L, 5L, 11L), LICA = c(0L, 2L, 2L, 37L,
24L, 16L, 25L, 19L, 0L, 21L, 49L, 85L), LICL = c(71L, 32L,
52L, 33L, 77L, 45L, 40L, 110L, 91L, 80L, 56L, 93L), LIPI = c(16L,
61L, 69L, 14L, 4L, 34L, 51L, 96L, 52L, 38L, 11L, 20L), LISE = c(28L,
56L, 32L, 5L, 13L, 1L, 0L, 41L, 72L, 24L, 1L, 46L), LISY = c(0L,
4L, 13L, 2L, 9L, 11L, 1L, 72L, 0L, 14L, 5L, 6L), PSCR = c(7L,
13L, 18L, 9L, 16L, 13L, 9L, 36L, 9L, 17L, 13L, 29L)), .Names = c("Area",
"Rank", "ANAM", "HYVE", "LICA", "LICL", "LIPI", "LISE", "LISY",
"PSCR"), class = "data.frame", row.names = c(NA, -12L))

mv <- mvabund(mydata[,3:10])
fit <- manyglm(mv ~ Rank, data=mydata, family="negative.binomial")
summary(fit, resamp="monte.carlo", p.uni="unadjusted", test="LR")

Error in rownames<-(*tmp*, value = c("ANAM", "HYVE", "LICA", "LICL", :
length of 'dimnames' [1] not equal to array extent

univariate P's for anova.manylm

Originally in an email from David to Julian, will close if fixed

This is a weird one in the C code.

For the below design, there are 3 reps in 2 treatments, so when using resampling the smallest P-value you should be able to get is 1/10 (10 possible permutations of 6 obs in 2 treatments). But the univariate P-value for the second response is 0.002.

My guess is that there is some floating point error in the code – we should be using statements like

t*>t-eps

rather than t*>=t

but perhaps in the code for univariate P-values this isn’t happening?

library(mvabund)
load("kos_bug.RData")
kos2

kosFac
# fitting a linear model to log-transformed proportions:
prop.mva = kos2.mva/kosFac$sampleTotals
kos2.lm = manylm(log(prop.mva) ~ condition, data = kosFac)
anova(kos2.lm,p.uni="unadjusted",nBoot=999)

Analysis of Variance Table

Model: manylm(formula = log(prop.mva) ~ condition, data = kosFac)

Overall test for all response variables
Test statistics:
            Res.Df Df.diff val(F) Pr(>F) 
(Intercept)      5                       
condition        4       1    282  0.094 .
---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Univariate Tests

Test statistics:

             K02014         K00059         K02033         K02034         K07316      

            F value Pr(>F) F value Pr(>F) F value Pr(>F) F value Pr(>F) F value Pr(>F)

(Intercept)                                                                          

condition    51.331  0.094 148.638  0.002  44.148  0.094  37.196  0.094   0.643  0.591

Arguments: with 999 resampling iterations using residual (without replacement) resampling and response assumed to be uncorrelated

kos_bug.rdata

fitted values truncated to 10^6 or 10^-6

fitted values are typically truncated to a max of 10^6 or a min of 10^-6. Can control of this value be given to the user, but maybe defaulting to maxFit=10^6? In case they have a dataset with some huge counts?

BootID with case resampling- bug in implementation of offset

When doing an anova using a randomly generated (independent and identically distributed) BootID matrix, with resamp= "case", the p-values are all around 0.5. If no BootID is used, p-values on the same data are 0.001. This suggests that when BootID argument is supplied, the offset is not being implemented correctly.

CODE TO REPLICATE BELOW:

rm(list=ls())

set.seed(42)

simulate some rough and ready play data

co-ords

x = runif(500)
y = runif(500)

covariates

library(RandomFields)
temperature = RandomFields::RFsimulate(model = RandomFields::RMexp(var=1, scale=0.1),
x=x, y=y,grid=F,n=1)$variable1
#20 species

nspec = 20
responseMultiSpecies = matrix(nrow=500, ncol=nspec)
Beta_multi = rnorm(nspec,10,1)

for (j in 1:nspec){
spatial_residuals4 = RandomFields::RFsimulate(model = RandomFields::RMexp(var=1, scale=0.05),
x=x, y=y,grid=F, n=1)$variable1
responseMultiSpecies[,j] = temperature*Beta_multi[j] + spatial_residuals4 + rnorm(500)
responseMultiSpecies[,j]=rbinom(n=500, size=1, prob=exp( responseMultiSpecies[,j])/c(1+ exp(responseMultiSpecies[,j])))
}

multispecies_dat=data.frame(responseMultiSpecies, temperature)

library(mvabund)

responseMultiSpecies=mvabund(multispecies_dat[,1:20]) #20 species multivariate reponse

models

mod.1 = manyglm(responseMultiSpecies1, data = multispecies_dat,family="binomial")
mod.2 = manyglm(responseMultiSpeciestemperature, data = multispecies_dat,family="binomial")

IID boot ID matrix

bootID_IID = matrix(sample(1:500,500*1000,replace=T),nrow=1000,ncol=500)

anova's

mvabund.case.bootID.0 = anova(mod.1, mod.2, resamp="case",bootID=bootID_IID)
mvabund.case.bootID.none = anova(mod.1, mod.2, resamp="case")

With a boot ID matrix, the p value is close of 0.5

mvabund.case.bootID.0$table

Without a boot ID matrix, the p value is close to 0.001

mvabund.case.bootID.none$table

try it a number of times (warning ~30 minute running time)

mvabund.case.bootID.0_ = c()
mvabund.case.bootID.none_ = c()
for (i in 1:50){
bootID_IID = matrix(sample(1:500,500*1000,replace=T),nrow=1000,ncol=500)

anova's

mvabund.case.bootID.0_[i] = anova(mod.1, mod.2, resamp="case",bootID=bootID_IID)$table[2,4]
mvabund.case.bootID.none_[i] = anova(mod.1, mod.2, resamp="case")$table[2,4]
}

Summary : with bootID argument p - values from 0.5025 - 0.5714, mean 0.5357.

Summary: without bootID argument p - values all 0.001

Gamma family in manyglm

Add support for the gamma family using the inverse link

  • update manyglm.R

  • update resampTest.h

  • update glm.cpp

  • update glmtest.cpp

  • update Rinterface.cpp

  • update manyglm.Rd

  • write some tests for manyglm

plot.manyglm ignoring small values

plot.manyglm is ignoring small predicted values, rather than truncating them to keep them on the plot. In default.plot.manyglm, lines 410-4:
yhtmp <- c(yh)
yh.is.zero <- yhtmp < (-6)

yh.is.zero <- yhtmp < max(-6,(-max(yhtmp)))#this line is wrong - it kicks out any value more negative than max(yh)

   yh0 <- yhtmp[!yh.is.zero]

should probably do something like yh0=pmax(-6,yh) instead

unexpected output during re-sampling with method="pit-trap" and offset

Hello - I am getting many lines of output printed to the console when running anova.manyglm on a model fit as follows:
mv.full <- manyglm(mvDat ~ Treatment*Genotype, offset=effort, family="negative_binomial")
Where effort = log(rowSums(mvDat))
When I run a test with:
anova(mv.full, nBoot=199, cor.type="R", test="wald")
The expected test results are returned, but only after printing many lines to the console that look like this:
l=nan, theta=1000000.0000, yi=44.0000, mu=-nan
l=nan, theta=1000000.0000, yi=24.0000, mu=-nan
l=nan, theta=1000000.0000, yi=43.0000, mu=-nan
l=nan, theta=1000000.0000, yi=41.0000, mu=-nan
l=nan, theta=1000000.0000, yi=45.0000, mu=-nan
l=nan, theta=1000000.0000, yi=8.0000, mu=-nan
There is one line for every sample*nBoot and occasionally yi=inf. This output does not appear if I remove the offset or use to a different resampling technique. Plotting the model diagnostics does not reveal any obvious problems and I get nearly identical statistics setting resamp = "perm.resid". However, this behavior does not occur when I add a similar offset to any of the example data provided with the package, so I am not sure how to provide a simple reproducible example. I am currently using the newest mvabund v4.1.1. Could this be some type of minor bug? Or perhaps my model is set up incorrectly in some way that is not obvious to me?

Predict not working with poly formulas

Hi Alice,
predict.manyglm is not working when the formula contains polynomial terms. This isn't urgent for me as I can work around.
Cheers,
Eve

#Minimal example below
#Simulate toy data
Y=rpois(100, 2)
X=data.frame(x1=rnorm(100), x2=rnorm(100))

#fit model and predict response from original data
mod= manyglm(Y~poly(x1,2)+poly(x2,2), data=X)
predict(mod)
predict(mod, newdata=X)
predict(mod, newdata=mod$x)

#This error is received for all 3 attempts
#Error in poly(x1, 2) : object 'x1' not found

family as vector input

can users specify a vector of families as input, so that different variance/link functions are used on different columns of data?

error with manyglm using $ in the formula and composition = T

Hello,

I'm getting the error ‘variable lengths differ’ when using the a $ in the manyglm formula, when composition =T

Herbivores <- read.csv(file = "Herbivore_specialisation.csv", header = TRUE)
Herb_spp <- mvabund(Herbivores[,5:11])
mod3 <- manyglm(Herb_spp ~ Herbivores$Habitat*Herbivores$DayNight, composition = T, family="negative_binomial")

Error in model.frame.default(formula = formLong, data = dat, subset = subset, :
variable lengths differ (found for 'Herbivores$Habitat')

I am able to get the model to work using:

  1. manyglm(Herb_spp~Habitat*DayNight, data = Herbivores))
  2. hab = Herbivores$Habitat
    DN = Herbivores$DayNight
    manyglm(Herb_spp~hab*DN)

Can this be corrected?
Thank you

Update .travis.yml

I noticed that the .travis.yml file is a little stale, and would like to suggest two changes:

  • update to what I currently use (see any of my repos, eg here for RcppGSL)
  • install all required CRAN packags as r-cran-* binaries (and I have checked they all exist in Michael's PPA).

[ Now, I don't mean to impose "my" style here. I am happy with this older / existing / maintained (by me) setup but I am of course aware that most of the world moved to another format (see the official Travis docs). I don't use that, but you are of course free to switch. I just can't help. ]

With the disclaimer out of the way, happy to support/update the existing one. Or not if you feel what is there now is good enough.

anova.manyglm fails with cor.type='shrink' 


## Visualise the effect of treatment on copepod abundance
tasm.cop <- mvabund(Tasmania$copepods)
#treatment <- Tasmania$treatment
# create >2 treatment variables variables
treatment = c(rep(c("Undis", "Dist"), 4),rep(c("Undis2", "Dist2"), 4))
block <- Tasmania$block
#plot(tasm.cop ~ treatment, col=as.numeric(block))

## Fitting predictive models using a negative binomial model for counts:
tasm.cop.nb <- manyglm(tasm.cop ~ treatment, family="negative.binomial")
# works fine...
anova_tasm.cop.nb = anova.manyglm(tasm.cop.nb, pairwise.comp = treatment, nBoot = 3, test="wald")

## Fitting predictive models using a negative binomial model for counts, and shrink:
tasm.cop.nb_shrink <- manyglm(tasm.cop ~ treatment, family="negative.binomial", cor.type = "shrink", test="wald")
# doesn't work
anova_pairwise_tasm.cop.nb = anova.manyglm(tasm.cop.nb_shrink, pairwise.comp = treatment, nBoot = 3, test="wald")

Time elapsed: 0 hr 0 min 0 sec
Error in dimnames(x) <- dn :
length of 'dimnames' [2] not equal to array extent
In addition: Warning messages:
1: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
2: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
3: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
4: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
5: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.
6: In anova.manyglm(m, show.time = "none", keep.boot = TRUE, nBoot = anova_obj$nBoot, :
The likelihood ratio test can only be used if correlation matrix of the abundances is is assumed to be the Identity matrix. The Wald Test will be used.

I've done my best to find where it's going wrong and as far as I can tell it's in do_pairwise_comp ()

am <- anova.manyglm(m,
                          show.time = 'none',
                          keep.boot = TRUE,
                          nBoot = anova_obj$nBoot,resamp = resamp, cor.type = cor.type)

which gives a table like
Analysis of Variance Table

Model: subY ~ subWhat

Multivariate test:
Res.Df Df.diff wald Pr(>wald)
(Intercept) 5
subWhat 4 2 0.25
Arguments:
Test statistics calculated assuming correlated response via ridge regularization
P-value calculated using 3 iterations via PIT-trap resampling.

and am$table[2, 3] is NaN

whereas running it without cor.type = 'shrink' gives am$table[2, 3] a numerical value (165.5)
Analysis of Deviance Table

Model: subY ~ subWhat

Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 5
subWhat 4 2 165.5 0.25
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 3 iterations via PIT-trap resampling.

Is there a fix for this? I'm using the latest version of mvabund from github

Error with GAMMA manyglm var.est

Currently manyglm is struggling to handle majority of gamma distributed data that I have been attempting to model.

The RtoGlm function easily results in NA's for the z$var.estimator when the gamma family is used (below code from manyglm)....

      z <- RtoGlm(modelParam, Y, X, O)
      if (any(z$var.est == 0)) {
        z$var.estimator = pmax(z$var.est, 1e-06) 
        z$residuals = (Y - z$fit)/sqrt(z$var.est)

Which stops the model estimation at the above step because any(z$var.est == 0) returns an NA. Not sure why RtoGlm is returning NA's here quite often for gamma data which glm has no issues fitting. Perhaps best to replace NA's with 0 and then run the model replacing the var.estimator appropriately for example below (using tidyverse replace_na).

      z <- RtoGlm(modelParam, Y, X, O)
      z$var.est[is.na(z$var.est)]<-0
      if (any(replace_na(z$var.est,0) == 0)) {
        z$var.estimator = pmax(z$var.est, 1e-06)
        z$residuals = (Y - z$fit)/sqrt(z$var.est)

Happy to provide data and further discussion if need be. This issue makes it really hard for gamma function to be taken up and applied widely!

cheers

Ben
[email protected]

pairwise.comp factor environment usage

Thanks Ben for catching this issue
pairwise.comp wasn't looking for the factor in the manyglm fit environment, however when it is put the in the global environment it works.

library(mvabund)

#import data
tb_dat <- read.csv("tb_dat.csv",header=T)
fishspp  <- mvabund(tb_dat[,8:13])
#Round off fishspp to integer values
fishspp <- round(fishspp)

#fit a manyglm object
pref.mod <- manyglm(fishspp ~ region + spp.zone + spp + offset(log(time_viewed)), data = tb_dat)

#pairwise comp doesn't work as it can't find spp within the model
pref.aov <- anova(pref.mod, pairwise.comp = spp,nBoot = 199)
pref.aov

#now if we put species in the global environment it works!
spp <- as.factor(tb_dat$spp)
pref.aov <- anova(pref.mod, pairwise.comp = spp, nBoot = 199)
pref.aov

a solution is maybe to deparse the name or to make the documentation and usage examples more clear

summary P-value calculation

bootstrap P-values in summary.manyglm and summary.manylm are computed as (# resampled stats >= observed stat) / nBoot rather than (# resampled stats >= observed stat + 1) / ( nBoot +1). The latter is preferred, but is currently only implemented in summary for perm.resid, please fix.

Feature request: Please consider adding the univariate test and p-value extraction as one table

Thank you for creating and maintaining mvabund. I want to extract the Univariate section, which comes after the Multivariate ANOVA table as one table, but the current attributes have the p-value and the test table separated. We could extract the 2 tables and merge them, but I think being able to extract the original format would be great.
https://stackoverflow.com/questions/65316449/extract-mvabundanova-manyglms-p-values-and-tests-as-one-table

remove set seed on anova.manylm

there is a set seed somewhere in anova.manylm, such that it always return identical P-values. e.g.
example(manylm)
anova(fit, fit0, fit1, fit2, nBoot=5000, test="F", p.uni="adjust")
anova(fit, fit0, fit1, fit2, nBoot=5000, test="F", p.uni="adjust")
We should have different P-values due to MC error but don't.

This does not seem to be an issue in summary.manylm nor in the manyglm functions.

Var.est error with gamma manyglm

Hi,

I am trying to run a manyglm with a gamma distribution but am getting a strange error which I don't understand.

My error can be replicated with (shortened version):

mydata <- read.csv("https://raw.githubusercontent.com/HaydenSchilling/MGLMs-Otoliths/master/Data/Otolith_data_mmol_mol_Ca.csv")
E_data <- mydata[,c(2:8,10:14)]
elements <- mvabund(E_data) # select only element data
gamma_glm <- manyglm(elements ~ 1, family="gamma")

The error says:

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

It is possible to run a univariate glm with gamma distribution on all of variables I'm trying to use in the manyglm

glm1 <- glm(elements[,1] ~ 1, family =Gamma(link="log"))

Any ideas would be very much appreciated,
Thanks!

bug in anova.manyany for comp vs no comp

if calling anova.manyany on two objects, where one has comp=TRUE and the other has COMP=FALSE, nonsensical answers are returned (because only the first variable from COMP=FALSE is used)

summary.manyglm fails if family is "binomial(link=logit)"

Fails:
[from the demo code]

glm.spid.bin <- manyglm(pres.abs~soil.dry+bare.sand+moss, data=X, family="binomial")
summary(glm.spid.bin)

Error in summary.manyglm(glm.spid.bin) :
'family' not defined. Choose one of 'poisson', 'negative.binomial', 'binomial' for an manyglm object

The family is:

glm.spid.bin$family
[1] "binomial(link=logit)"

But summary.manyglm looks for an exact string match on the family (object$family == "binomial"), so it fails. Probably this test should be more like the one used in anova.manyglm (substr(object$family, 1, 1) == "b")

sessionInfo()
R version 3.2.2 (2015-08-14)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] mvabund_3.11.5

anova.manyglm overriding cloglog fits with binomial("logit")

For manyglm objects fitted using binomial("cloglog") model, anova.manyglm currently ignores the cloglog bit and tests significance using a logit link. Alice is working on this one as we speak...
e.g. (using code at pull request #15 which fixes a previous bug with family calls)

example(manyglm)
glm2.ft=manyglm(pres.abssoil.dry,data=X,family=binomial("cloglog"))
glm1.ft=manyglm(pres.abs1,data=X,family=binomial("cloglog"))
sum(deviance(glm1.ft)-deviance(glm2.ft)) # returns [1] 129.3183
anova(glm2.ft,nBoot=9)$table[2,3] #returns the wrong answer, [1] 131.0454

glm1.ftb=manyglm(pres.abs1,data=X,family=binomial())
glm2.ftb=manyglm(pres.abssoil.dry,data=X,family=binomial())
sum(deviance(glm1.ftb)-deviance(glm2.ftb)) #returns [1] 131.0454
anova(glm2.ftb,nBoot=9)$table[2,3] #returns [1] 131.0454

inf residuals

when fitting a Poisson model to a count dataset with outliers, outlying residuals come back as Inf which gives an error when you try to plot. Maybe truncate Inf to some big value and return a warning?

manylm not accepting input bootID

example(anova.manylm)
bootid = matrix(NA,4,28)
for(i in 1:4)
{
bootid[i,1] = i
bootid[i,2:28] = (1:28)[-i]
}
anova(fit,bootID=bootid)
anova(fit,bootID=as.numeric(bootid))

note in contrast that manyglm has the appropriate behaviour

block function in anova/summary

From a forwarded email of Davids to Julian:
"
There is a block function in manyglm. I would like to:

  • Extend it to manylm
  • Remove the error “needs to be balanced” if resamp=”case”
    "

pairwise.comp error with anova.manyglm

Hi, I have two questions. First, when I use the pairwise.comp option with the manyglm function I am receiving an error that I can not seem to figure out. All the rest of the functions work with my data up to this point. Attached is a simplified data file and R-script that reproduced the error, and below is the script again for easy viewing.

####script
Data<-read.csv("ComMat.csv")
Community<-(Data[,3:5]+0.1) #Creating matrix, and adding tiny constant
Factor<-Data$Factor

Mv_Matrix<-mvabund::mvabund(Community) #works
Mod_Mv_Matrix<-mvabund::manyglm(Mv_Matrix ~ Factor, family="gamma") #works
mvabund::anova.manyglm(Mod_Mv_Matrix, nBoot = 1000, test="LR") #works
mvabund::anova.manyglm(Mod_Mv_Matrix,nBoot = 1000, test="LR", p.uni="adjusted") #works
mvabund::anova.manyglm(Mod_Mv_Matrix,nBoot = 1000, pairwise.comp = ~ Factor, test="LR", p.uni="adjusted")

Everything works until the last line with the pairwise comparisons, when I get this error:
Error in apply(resampled_stats, 1, per_row) :
dim(X) must have a positive length

I would be grateful for any help
Thanks
Jamie

Archive.zip

unable to define booID within summary.manyglm

Hi,
Thanks so much for developing this great package.

I seem to have an issue with passing user-defined permutation to summary.manyglm by bootID. The version of mvabund is v3.11.9.

I only have 6 rows in a multivariate abundance data, so I was trying to perform an exact permutation by 'bootID = allPerm(6)'. Then, I get the error message below.
'Error in max(bootID) : invalid 'type' (symbol) of argument'. Please see the example.

require(mvabund)
require(permute)
data(spider)
spiddat <- mvabund(spider$abund)
X <- spider$x
spiddat2 <- spiddat[1:6, ]
X2 <- X[1:6, ]
glm.spid <- manyglm(spiddat2[,1:3]~X2, family="negative.binomial")
summary(glm.spid, bootID = allPerms(6))

Error in max(bootID) : invalid 'type' (symbol) of argument

I am wondering if I'm not defining it right. Any advice would be extremely helpful.

Regards,

Shun

Support for the tweedie family

The tweedie family is a very flexible family of distributions, We are only interested in tweedie distributions with power between 1 and 2, this means that they have mass at 0 but otherwise are similar to a gamma distribution,

Useful Links

File src/mvabund_types.h now creates problems for Rcpp

mvabund uses what was once a "pre-release" of the Rcpp::Nullable<> class in src/mvabund_types.h.

Rcpp itself now has a slightly newer version, and including both now leads to errors.

In order to sort this out, I would suggest to simply rely on Rcpp alone. I have verified that mvabund then build fine.

We can either just remove src/mvabund_types.h, or keep it around (if we every need it again) but keep it empty. I would be happy to prepare a corresponding pull request for either option.

With best wishes for 2017, Dirk

binstruc + binomial manyglm matrix size error

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?

CI updates

The Travis CI support via file .travis.yml could be updated to using Ubuntu 20.04 and slightlt nicer internal logic.

However, a bigger / more important change may be to also switch to / enable GitHub Actions. I am doing that for my repos now one-by-one because I have run out of so-called "free" credits at Travis. This repo runs under Alice's handle and allocation so we get by with the 10,000 credits from Travis.

More details (and examples) are at https://eddelbuettel.github.io/r-ci/ -- if it helps I would be happy to simplify and improve the setup. We can simply turn both Travis CI and GitHub Actions on (with a small update to Travis).

summary.manyglm() fails with p.uni = "adjusted" or "unadjusted"

This fails in two cases / two reasons (with the mvabund version here in github):

data(spider)
spiddat <- mvabund(spider$abund)
X <- spider$x
glm.spid <- manyglm(spiddat[,1:3]~X, family="negative.binomial")
summary(glm.spid, nBoot = 100, p.uni = 'none') 
summary(glm.spid, nBoot = 100, p.uni = 'unadjusted') 
# Error in `colnames<-`(`*tmp*`, value = c("(Intercept)", "Xsoil.dry", "Xbare.sand",  : 
#  length of 'dimnames' [2] not equal to array extent
summary(glm.spid, nBoot = 100, p.uni = 'adjusted') 
# Error in `colnames<-`(`*tmp*`, value = c("(Intercept)", "Xsoil.dry", "Xbare.sand",  : 
#  length of 'dimnames' [2] not equal to array extent

This fails in line 255 because there are not enough columns in unit_signic.

require(mvabund)
require(vegan)
# data
data(pyrifos)
take <- c('binitent', 'olchaeta', 'caenhora', 'cloedipt',
          'chaoobsc', 'gammpule', 'libellae', 'agdasphr') 
abu <- pyrifos[ , names(pyrifos) %in% take] 
abu <- round((exp(abu) - 1)/10) 
env <- data.frame(time = rep(c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24), each = 12),
                  treatment = rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11),
                  replicate = gl(12, 1, length=132))
abu <- mvabund(abu)
# subset to only one time point
take_abu <- abu[env$time == '-4', ]
take_env <- env[env$time == '-4', ]

# model + summary
mod <- manyglm(take_abu ~ factor(treatment), data = take_env)
summary(mod, nBoot = 100,  p.uni = 'adjusted')

This fails because unit_signic is empty!

In both cases I suspect the error coming from RtoGlmSmry() in Line 199.
However, I haven't tracked it further down (I not a C++ Ninja and not experienced in debugging C++), though I expect the problem in Rinterface.cpp or further down in glmtest.cpp.

manyglm.summary fails on intercept model

for say poisson data:

library(mvabund)
Y <- mvabund(sapply(1:4, function(x)  rpois(30, 5)))
m <- manyglm(Y ~ 1 , family = 'poisson') 
summary(m)

Causes a segmentation fault

Weights in manyglm

a good thing to add to manyglm would be observation weights (like the weights argument on glm).
Notes from David:

  • These could be a vector (applied across rows of data) or a matrix (different weight for each entry in the response matrix)
  • The main potential use for this is fitting logistic regression models to binomial counts (using number of trials as the weights and sample proportion as the response), so I would test if it was working against glm applied in this way.
  • Another possible application is fitting multivariate presence-only data but the details need to be fleshed out

Changes needed to comply with R-devel

Running RD CMD check --as-cran where RD is an alias for a very recent R-devel build leads to a lot of line noise. I can fix that -- as I have fixed it for other packages:

* using log directory ‘/home/edd/git/mvabund.Rcheck’
* using R Under development (unstable) (2015-08-31 r69241)
* using platform: x86_64-pc-linux-gnu (64-bit)
* using session charset: UTF-8
* using option ‘--as-cran’
* checking for file ‘mvabund/DESCRIPTION’ ... OK
* this is package ‘mvabund’ version ‘3.11.3.2’
* checking CRAN incoming feasibility ... NOTE
Possibly mis-spelled words in DESCRIPTION:
  Rcpp (12:48)
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for executable files ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package ‘mvabund’ can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking use of S3 registration ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... NOTE
RcppVersion: no visible global function definition for
  ‘installed.packages’
absspp: no visible global function definition for ‘var’
anova.manyany: no visible global function definition for ‘logLik’
anova.manyglm : <anonymous>: no visible global function definition for
  ‘formula’
anova.manylm : <anonymous>: no visible global function definition for
  ‘formula’
as.mvformula: no visible global function definition for ‘as.formula’
as.mvformula: no visible global function definition for ‘terms’
as.mvformula: no visible global function definition for ‘model.frame’
best.r.sq: no visible global function definition for ‘terms’
best.r.sq: no visible global function definition for ‘reformulate’
best.r.sq: no visible global function definition for ‘update’
betaest: Error while checking: invalid 'envir' argument
bootAnova: no visible global function definition for ‘logLik’
case.names.manyglm: no visible global function definition for ‘weights’
case.names.manyglm: no visible global function definition for
  ‘residuals’
case.names.manylm: no visible global function definition for ‘weights’
case.names.manylm: no visible global function definition for
  ‘residuals’
cooks.distance.manylm: no visible global function definition for
  ‘residuals’
cooks.distance.manylm: no visible global function definition for
  ‘deviance’
cooks.distance.manylm: no visible global function definition for
  ‘df.residual’
covratio.manylm: no visible global function definition for
  ‘weighted.residuals’
cv.glm1path: no visible global function definition for ‘flush.console’
cv.glm1path: no visible binding for global variable ‘sd’
default.boxplot.mvformula: no visible global function definition for
  ‘terms’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘dev.list’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘postscript’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘pdf’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘jpeg’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘bmp’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘png’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘dev.off’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘var’
default.meanvar.plot.mvabund: no visible global function definition for
  ‘lm’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘dev.list’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘terms’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘model.response’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘model.frame’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘var’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘palette’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘rainbow’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘postscript’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘pdf’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘jpeg’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘bmp’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘png’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘dev.off’
default.meanvar.plot.mvformula: no visible global function definition
  for ‘lm’
default.plot.manyglm: no visible global function definition for
  ‘dev.list’
default.plot.manyglm: no visible global function definition for
  ‘postscript’
default.plot.manyglm: no visible global function definition for ‘pdf’
default.plot.manyglm: no visible global function definition for ‘jpeg’
default.plot.manyglm: no visible global function definition for ‘bmp’
default.plot.manyglm: no visible global function definition for ‘png’
default.plot.manyglm: no visible global function definition for
  ‘dev.off’
default.plot.manyglm: no visible global function definition for
  ‘residuals’
default.plot.manyglm: no visible global function definition for
  ‘rainbow’
default.plot.manyglm: no visible global function definition for
  ‘df.residual’
default.plot.manyglm: no visible global function definition for
  ‘dev.interactive’
default.plot.manylm: no visible global function definition for
  ‘residuals’
default.plot.manylm: no visible global function definition for
  ‘dev.list’
default.plot.manylm: no visible global function definition for
  ‘postscript’
default.plot.manylm: no visible global function definition for ‘pdf’
default.plot.manylm: no visible global function definition for ‘jpeg’
default.plot.manylm: no visible global function definition for ‘bmp’
default.plot.manylm: no visible global function definition for ‘png’
default.plot.manylm: no visible global function definition for
  ‘dev.off’
default.plot.manylm: no visible global function definition for
  ‘predict’
default.plot.manylm: no visible global function definition for
  ‘weights’
default.plot.manylm: no visible global function definition for
  ‘rainbow’
default.plot.manylm: no visible global function definition for
  ‘df.residual’
default.plot.manylm: no visible global function definition for
  ‘na.omit’
default.plot.manylm: no visible global function definition for
  ‘weighted.residuals’
default.plot.manylm: no visible global function definition for
  ‘deviance’
default.plot.manylm: no visible global function definition for
  ‘dev.interactive’
default.plot.mvabund: no visible global function definition for
  ‘dev.list’
default.plot.mvabund: no visible global function definition for
  ‘as.formula’
default.plot.mvabund: no visible global function definition for
  ‘update’
default.plot.mvabund: no visible global function definition for
  ‘postscript’
default.plot.mvabund: no visible global function definition for ‘pdf’
default.plot.mvabund: no visible global function definition for ‘jpeg’
default.plot.mvabund: no visible global function definition for ‘bmp’
default.plot.mvabund: no visible global function definition for ‘png’
default.plot.mvabund: no visible global function definition for
  ‘dev.cur’
default.plot.mvabund: no visible global function definition for
  ‘dev.off’
default.plot.mvabund: no visible global function definition for
  ‘dev.set’
default.plot.mvformula: no visible global function definition for
  ‘dev.list’
default.plot.mvformula: no visible global function definition for
  ‘as.formula’
default.plot.mvformula: no visible global function definition for
  ‘terms’
default.plot.mvformula: no visible global function definition for
  ‘model.response’
default.plot.mvformula: no visible global function definition for
  ‘model.frame’
default.plot.mvformula: no visible global function definition for
  ‘postscript’
default.plot.mvformula: no visible global function definition for ‘pdf’
default.plot.mvformula: no visible global function definition for
  ‘jpeg’
default.plot.mvformula: no visible global function definition for ‘bmp’
default.plot.mvformula: no visible global function definition for
  ‘dev.cur’
default.plot.mvformula: no visible global function definition for
  ‘dev.off’
default.plot.mvformula: no visible global function definition for
  ‘dev.set’
default.plotMvaFactor: no visible global function definition for
  ‘dev.list’
default.plotMvaFactor: no visible global function definition for
  ‘postscript’
default.plotMvaFactor: no visible global function definition for ‘pdf’
default.plotMvaFactor: no visible global function definition for ‘jpeg’
default.plotMvaFactor: no visible global function definition for ‘bmp’
default.plotMvaFactor: no visible global function definition for ‘png’
default.plotMvaFactor: no visible global function definition for
  ‘dev.cur’
default.plotMvaFactor: no visible global function definition for
  ‘dev.off’
default.plotMvaFactor: no visible global function definition for
  ‘dev.set’
default.plotMvaFactor: no visible global function definition for
  ‘rainbow’
default.print.anova.manyglm: no visible global function definition for
  ‘printCoefmat’
default.print.summary.manyglm: no visible global function definition
  for ‘symnum’
default.print.summary.manyglm: no visible binding for global variable
  ‘quantile’
default.print.summary.manyglm: no visible global function definition
  for ‘quantile’
deviance.manylm: no visible global function definition for ‘na.fail’
deviance.manylm: no visible global function definition for ‘na.omit’
deviance.manylm: no visible global function definition for ‘na.exclude’
dfbeta.manyglm: no visible global function definition for
  ‘variable.names’
dfbeta.manylm: no visible global function definition for
  ‘variable.names’
dfbetas.manylm: no visible global function definition for ‘dfbeta’
dffits.manylm: no visible global function definition for
  ‘weighted.residuals’
effic.make.link : mu.eta.fct: no visible global function definition for
  ‘dnorm’
effic.make.link : linkinv.fct: no visible global function definition
  for ‘qnorm’
effic.make.link : linkinv.fct: no visible global function definition
  for ‘pnorm’
effic.make.link : mu.eta.fct: no visible global function definition for
  ‘dcauchy’
effic.make.link : linkinv.fct: no visible global function definition
  for ‘qcauchy’
effic.make.link : linkinv.fct: no visible global function definition
  for ‘pcauchy’
effic.make.link: no visible global function definition for ‘make.link’
extend.x.formula: no visible global function definition for ‘terms’
extend.x.formula: no visible global function definition for
  ‘reformulate’
extend.x.formula: no visible global function definition for ‘update’
extractAIC.manyglm: no visible global function definition for ‘logLik’
formulaMva: no visible global function definition for ‘terms’
formulaMva: no visible global function definition for ‘model.frame’
formulaUnimva: no visible global function definition for ‘update’
formulaUnimva: no visible global function definition for ‘terms’
formulaUnimva: no visible global function definition for ‘as.formula’
get.design: no visible global function definition for ‘as.formula’
get.design: no visible global function definition for ‘model.matrix’
get.design: no visible global function definition for ‘contrasts’
get.design: no visible global function definition for ‘terms’
get.polys: no visible global function definition for ‘contrasts’
get.polys: no visible global function definition for ‘contrasts<-’
get.polys: no visible global function definition for ‘poly’
make.vslink: no visible global function definition for ‘make.link’
manyany: Error while checking: invalid 'envir' argument
manyglm: no visible global function definition for ‘model.offset’
manyglm: no visible global function definition for ‘model.response’
manyglm: no visible global function definition for ‘model.matrix’
manyglm: no visible global function definition for ‘.getXlevels’
manylm: no visible global function definition for ‘model.response’
manylm: no visible global function definition for ‘model.weights’
manylm: no visible global function definition for ‘model.offset’
manylm: no visible global function definition for ‘is.empty.model’
manylm: no visible global function definition for ‘model.matrix’
manylm: no visible global function definition for ‘.getXlevels’
manylm.fit: no visible binding for global variable ‘data’
manylm.influence: no visible global function definition for
  ‘weighted.residuals’
manylm.influence.measures : is.influential: no visible global function
  definition for ‘pf’
manylm.influence.measures: no visible global function definition for
  ‘weighted.residuals’
manylm.influence.measures: no visible global function definition for
  ‘df.residual’
manylm.influence.measures: no visible global function definition for
  ‘variable.names’
mu.update: Error while checking: invalid 'envir' argument
mvformula: no visible global function definition for ‘formula’
mvformula: no visible global function definition for ‘terms’
plot.glm1path: no visible global function definition for ‘loess.smooth’
plot.manyany: no visible global function definition for ‘qnorm’
plot.manyany: no visible global function definition for ‘rainbow’
plotFormulafeature: no visible global function definition for
  ‘model.matrix’
plotFormulafeature: no visible global function definition for
  ‘model.frame’
plotFormulafeature: no visible global function definition for ‘terms’
predict.manyglm: no visible binding for global variable ‘na.pass’
predict.manyglm: no visible global function definition for ‘formula’
predict.manyglm: Error while checking: invalid 'envir' argument
predict.manylm: no visible binding for global variable ‘na.pass’
predict.manylm: no visible global function definition for ‘formula’
predict.manylm: no visible global function definition for ‘as.formula’
predict.manylm: no visible global function definition for ‘model.frame’
predict.manylm: no visible global function definition for ‘lm’
predict.manylm: no visible global function definition for ‘predict.lm’
predict.traitglm: no visible global function definition for ‘coef’
print.anova.manyany: no visible global function definition for
  ‘printCoefmat’
print.anova.manylm: no visible global function definition for
  ‘printCoefmat’
print.manyany: no visible global function definition for ‘qnorm’
print.manyglm: no visible global function definition for ‘coef’
print.manyglm: no visible global function definition for ‘naprint’
print.manylm: no visible global function definition for ‘coef’
print.manylm: no visible global function definition for ‘naprint’
print.summary.manylm: no visible global function definition for
  ‘symnum’
print.summary.manylm: no visible binding for global variable ‘quantile’
print.summary.manylm: no visible global function definition for
  ‘quantile’
quasibinomialMany: no visible global function definition for
  ‘make.link’
quasipoisson: no visible global function definition for ‘make.link’
rcalc: no visible global function definition for ‘var’
residuals.glm1path: no visible global function definition for ‘runif’
residuals.glm1path: no visible global function definition for ‘qnorm’
residuals.manyany: no visible global function definition for ‘runif’
residuals.manyglm: no visible global function definition for ‘coef’
residuals.manyglm: no visible global function definition for ‘runif’
residuals.manyglm: no visible global function definition for ‘qnorm’
ridgeParamEst: no visible global function definition for ‘runif’
ridgeParamEst: no visible global function definition for ‘cov’
rstandard.manylm: no visible global function definition for ‘deviance’
rstandard.manylm: no visible global function definition for
  ‘df.residual’
rstandard.manylm: no visible global function definition for
  ‘weighted.residuals’
shiftpoints: no visible global function definition for ‘na.omit’
traitglm: no visible global function definition for ‘coef’
transAxis: no visible global function definition for ‘terms’
transAxis: no visible global function definition for ‘dev.off’
weighted.residuals.manylm: no visible global function definition for
  ‘weights’
weighted.residuals.manylm: no visible global function definition for
  ‘residuals’
weights.manyglm: no visible global function definition for ‘naresid’
Undefined global functions or variables:
  .getXlevels as.formula bmp coef contrasts contrasts<- cov data
  dcauchy dev.cur dev.interactive dev.list dev.off dev.set deviance
  df.residual dfbeta dnorm flush.console formula installed.packages
  is.empty.model jpeg lm loess.smooth logLik make.link model.frame
  model.matrix model.offset model.response model.weights na.exclude
  na.fail na.omit na.pass naprint naresid palette pcauchy pdf pf png
  pnorm poly postscript predict predict.lm printCoefmat qcauchy qnorm
  quantile rainbow reformulate residuals runif sd symnum terms update
  var variable.names weighted.residuals weights
Consider adding
  importFrom("grDevices", "bmp", "dev.cur", "dev.interactive",
             "dev.list", "dev.off", "dev.set", "jpeg", "palette", "pdf",
             "png", "postscript", "rainbow")
  importFrom("stats", ".getXlevels", "as.formula", "coef", "contrasts",
             "contrasts<-", "cov", "dcauchy", "deviance", "df.residual",
             "dfbeta", "dnorm", "formula", "is.empty.model", "lm",
             "loess.smooth", "logLik", "make.link", "model.frame",
             "model.matrix", "model.offset", "model.response",
             "model.weights", "na.exclude", "na.fail", "na.omit",
             "na.pass", "naprint", "naresid", "pcauchy", "pf", "pnorm",
             "poly", "predict", "predict.lm", "printCoefmat", "qcauchy",
             "qnorm", "quantile", "reformulate", "residuals", "runif",
             "sd", "symnum", "terms", "update", "var", "variable.names",
             "weighted.residuals", "weights")
  importFrom("utils", "data", "flush.console", "installed.packages")
to your NAMESPACE.
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking contents of ‘data’ directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking line endings in C/C++/Fortran sources/headers ... OK
* checking line endings in Makefiles ... OK
* checking compilation flags in Makevars ... OK
* checking for GNU extensions in Makefiles ... OK
* checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
* checking compiled code ... OK
* checking examples ... [53s/54s] OK
Examples with CPU or elapsed time > 5s
                 user system elapsed
anova.manyany   7.868  0.056   7.939
traitglm        6.812  0.000   6.821
summary.manyglm 4.996  0.004   5.016
* checking PDF version of manual ... OK
* DONE
Status: 2 NOTEs

It looks worse than it is. The suggested change to NAMESPACE will do most of this. But it helps to have a version of R-devel nearby ...

Obtain univariate p-values for all pairwise comparisons of factor levels in anova.manyglm

Dear,

I was wondering if there is a way to combine the pairwise comparison option with the p.uni="adjusted" option in anova.manyglm to obtain adjusted p-values for every response variable for every possible pairwise combination of factor levels.

While p.uni="adjusted" or pairwise comparisons work fine when specified separately, I get the following error message when I use both arguments together:
Error in dimnames(x) <- dn :
length of 'dimnames' [2] not equal to array extent

Many thanks in advance for your time and consideration.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.