R/summary.manylm.R In mvabund: Statistical Methods for Analysing Multivariate Abundance Data

Documented in summary.manylm

```###############################################################################
# R user interface to summary test of a multivarivate linear model
# Author: Yi Wang (yi dot wang at computer dot org)
# 05-Jan-2010
###############################################################################

summary.manylm <- function(object, nBoot=1000,resamp="residual", test="F", cor.type=object\$cor.type, shrink.param=NULL, p.uni="none", studentize=TRUE, R2="h", show.cor = FALSE, show.est=FALSE, show.residuals=FALSE, symbolic.cor = FALSE, tol=1.0e-6, ... )
{
allargs <- match.call(expand.dots = FALSE)
dots <- allargs\$...
if ("rep.seed" %in% names(dots)) rep.seed <- eval(parse(text=dots\$rep.seed))
else rep.seed <- FALSE
if ("ld.perm" %in% names(dots)) ld.perm <- eval(parse(text=dots\$ld.perm))
else ld.perm <- FALSE
if ("bootID" %in% names(dots)) bootID <- eval(parse(text=dots\$bootID))
else bootID <- NULL

if(!any(class(object)=="manylm"))
stop("The function 'summary.manylm' can only be used for a manylm object.")

######## Construct Summary Input Arguments ########
nRows = nrow(object\$y)
nVars = ncol(object\$y)
nParam = ncol(object\$x)
Y <- as.matrix(object\$y, nrow=nRows, ncol=nVars)
X <- as.matrix(object\$x, "numeric")
if(substr(p.uni,1,1) == "n"){
pu <- 0
calc.pj <- FALSE
} else if(substr(p.uni,1,1) == "u"){
pu <- 1
calc.pj <- TRUE
} else if(substr(p.uni,1,1)=="a"){
pu <- 2
calc.pj <- TRUE
} else if(substr(p.uni,1,1) == "s"){
pu <- 3
calc.pj <- TRUE
} else

# the following values need to be converted to integer types
if (cor.type == "I") corr <- 1
else if (cor.type == "R") corr <- 0
else if (cor.type == "shrink") corr <- 2
else stop("No such correlation type.")

# To exclude case resampling
#if (resamp=="case") stop("Sorry, case resampling is not yet available.")
if (resamp=="case") resam <- 0
else if (resamp == "residual") resam <- 1
else if (resamp == "score") resam <- 2
else if (resamp == "perm.resid") resam <- 3
else stop("No such resampling method.")

if (test=="LR") testype <- 0
else if (test == "F") testype <- 1
else stop("No such test method.")

# estimate ridge parameter if cor.type is not "shrink" when fitting the model
w <- object\$weights
if (is.null(w)){
### Fit the multivariate LM.
w   <- rep(1, times=nRows)
}
else {
if (!is.numeric(w))  stop("'weights' must be a numeric vector")
if (any(w < 0)) stop("negative 'weights' not allowed")
}

if (cor.type=="shrink" & is.null(shrink.param)) {
shrink.param <- ridgeParamEst(dat=Y, X=X, weights=w, only.ridge=TRUE, doPlot=FALSE, tol=tol)\$ridgeParameter
# to simplify later computation
if (shrink.param == 0) cor.type <- "I"
if (shrink.param == 1) cor.type <- "R"
if (abs(shrink.param)>1)
stop("the absolute 'shrink.param' should be between 0 and 1")
}
else if (cor.type == "I")
shrink.param <- 1
else if (cor.type == "R")
shrink.param <- 0

if (!is.null(bootID)) {
if (max(bootID)>nRows) {
bootID <- as.null()
cat(paste("Invalid bootID -- sample id larger than no. of observations. Switch to generating bootID matrix on the fly (default nBoot=999).","\n"))
}
else {
if (is.matrix(bootID)) nBoot <- dim(bootID)[1]
else nBoot <- as.integer(length(bootID)/nRows)

if ((resamp == "score")) {
if (is.numeric(bootID)) {
cat(paste("Using <double> bootID matrix from input for 'score' resampling.","\n"))
bootID <- matrix(as.numeric(bootID), nrow=nBoot, ncol=nRows)
}
else {
cat(paste("Invalid bootID -- 'score' resampling should use <double> matrix. Switch to generating bootID matrix on the fly.","\n"))
bootID <- as.null()

}
}
else{
if (is.integer(bootID)){
cat(paste("Using <int> bootID matrix from input.","\n"))
bootID <- matrix(as.integer(bootID-1), nrow=nBoot, ncol=nRows)
}
else {
cat(paste("Invalid bootID -- sample id for methods other than 'score' resampling should be integer numbers up to the no. of observations. Switch to generating bootID matrix on the fly.","\n"))
bootID <- as.null()
}
}
}
}

if (studentize) st <- 1
else st <- 0
if (substr(R2,1,1) == "h") {
rsq <- 0
R2name <- c("Hooper's R-squared")
}
else if (substr(R2, 1, 1) == "v") {
rsq <- 1
R2name <- c("Vector's R-squared")
}
else
stop ("No such R2 method.")

# construct for param list
params <- list(tol=tol, nboot=nBoot, cor_type=corr, shrink_param=shrink.param, test_type=testype, resamp=resam, reprand=rep.seed, studentize=st, punit=pu, rsquare=rsq)

######## Call Summary Rcpp #########
val <- RtoSmryCpp(params, Y, X, bootID)

######## Collect Summary Values ########
smry  <- list()
rankX <- object\$rank

# residual statistics  (from the previous R codes)
r <- as.matrix(object\$residuals)
rdf <- nRows - rankX   # residual rdf
genVar <- det(resvar)

if (is.null(object\$terms))
stop("invalid 'lm' object:  no 'terms'")
Qr <- object\$qr
if (is.null(Qr)) Qr <- qr(object\$x)
p1 <- 1:rankX
R  <- chol2inv(Qr\$qr[p1, p1, drop = FALSE])   # = inv X'X   pxp matrix
genVarB <- genVar*(c(R)[1+0:(rankX-1)*(rankX+1)])
est <- object\$coefficients[Qr\$pivot[p1], , drop=FALSE]
est <- cbind(genVarB, est)
dimnames(est)[[2]][1] <-  c("Gen. Var")

# residual correlations
if (show.cor)  {
nrowR <- nrow(R)
correlation <- matrix(NA,nrow=nrowR*nrow(resvar),ncol= nrowR*nrow(resvar))
se <- c(outer(sqrt(c(R)[1+0:(rankX-1)*(rankX+1)]),sqrt(c(resvar)[1+0:(nVars-1)*(nVars+1)] ) ))
for (i in 1:nrow(resvar))
for (j in 1:nrow(resvar))
correlation[((1:nrowR)+nrowR*(i-1)),((1:nrowR) + nrowR*(j-1))]<-(R*resvar[i,j])
correlation <- correlation / outer(se, se)
corrnames <- rep(paste("b_", 1:rankX, sep=""), times=nVars)
corrnames <- paste(corrnames, rep(1:nVars, each=rankX), sep="")
colnames(correlation) <- rownames(correlation) <- corrnames
}
else correlation <- NULL
rownames(R) <- rownames(X)[p1]
colnames(R) <- colnames(X)[p1]

# test statistics
if(!is.null(test)) {
testname <- paste(test,"value")
pname <- paste("Pr(>",test,")", sep="")
} else {
testname <- "no test"
pname    <- ""
}
responseNames <- colnames(object\$y)
explanNames <- colnames(X)
# significance
significance <- cbind(val\$signific, val\$Psignific)
dimnames(significance) <- list(NULL, c(testname, pname))
xnames <- colnames(X)
dimnames(significance)[[1]] <- xnames
# unitvariate tests
if (calc.pj & !is.null(test)){
unit_signic <- t(val\$unitsign)
unit_signic.p <- t(val\$Punitsign)
rownames(unit_signic) <- rownames(unit_signic.p) <- c(responseNames)
colnames(unit_signic) <- colnames(unit_signic.p) <- c(explanNames)
}
else {
unit_signic <- NULL
unit_signic.p <- NULL
}

# overall statistics
df.int  <- if (attr(object\$terms, "intercept")) 1 else 0
overaltest <- c(val\$multstat, val\$Pmultstat, numdf = rankX - df.int, dendf = rdf)
names(overaltest) <- c(testname, pname, "num df", "den df")
if (calc.pj) {
unit_overaltest <- cbind(val\$unitmult, val\$Punitmult)
dimnames(unit_overaltest) <- list(responseNames, c(testname, pname))
}
else unit_overaltest <- NULL

# resampling flags
smry\$p.uni <- p.uni
smry\$test  <- test
smry\$cor.type <- cor.type
smry\$resamp <- resamp
# display flags
smry\$show.est <- show.est
smry\$show.cor <- show.cor
smry\$show.residuals <- show.residuals
smry\$symbolic.cor <- symbolic.cor
# parameter values
smry\$shrink.param <- shrink.param
smry\$nBoot <- nBoot
smry\$n.bootsdone <- val\$nSamp
smry\$n.iter.sign <- nBoot - val\$nSamp
smry\$aliased <- c(is.na(object\$coef)[,1])
# residual stats
smry\$residuals <- r
smry\$df <- c(rankX, rdf, NCOL(Qr\$qr))
smry\$genVar <- genVar
smry\$est <- est
smry\$cov.unscaled <- R
smry\$correlation <- correlation
# test stats
smry\$coefficients <- significance
smry\$uni.test <- unit_signic
smry\$uni.p <- unit_signic.p
smry\$R2 <- R2name
smry\$r.squared <- c(val\$R2)
smry\$statistic <- overaltest
smry\$statistic.j <- unit_overaltest

class(smry) <- "summary.manylm"
return(smry)
}
setGeneric("summary")
setMethod("summary", "manylm", summary.manylm)
```

Try the mvabund package in your browser

Any scripts or data that you put into this service are public.

mvabund documentation built on Jan. 20, 2018, 9:19 a.m.