Commit b1e73f7f by Eric Coissac

Adds anova and continue pm model cleaning

parent c68f1565
......@@ -30,12 +30,14 @@ Collate:
'internals.R'
'procmod_frame.R'
'multivariate.R'
'procmod.R'
'covls.R'
'corls_test.R'
'formula_procmod_frame.R'
'omit_action.R'
'model_procmod.R'
'mprocuste.R'
'anova_pm.R'
'procmod.R'
'covls.R'
'corls_test.R'
'procuste.R'
'simulate.R'
'terms.R'
......@@ -7,6 +7,7 @@ S3method("[",procmod_frame)
S3method("[[<-",procmod_frame)
S3method(AIC,pm)
S3method(BIC,pm)
S3method(anova,pm)
S3method(as.data.frame,dist)
S3method(as.list,procmod_frame)
S3method(as_procmod_frame,array)
......@@ -19,6 +20,9 @@ S3method(dim,procmod_frame)
S3method(extractAIC,pm)
S3method(formula,procmod_frame)
S3method(model_procmod,default)
S3method(model_procmod,formula)
S3method(model_procmod,pm)
S3method(model_procmod,terms)
S3method(na.exclude,procmod_frame)
S3method(na.omit,procmod_frame)
S3method(names,procmod_corls)
......
#' @include mprocuste.R
NULL
#' Anova function for `pm` objects.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
anova.pm <- function(m, permutations = how(nperm = 999), ...) {
mf <- m$model
coef <- m$coefficients
sdxy <- sqrt(diag(m$cov))
mt <- terms(mf)
reponse <- attr(mt, "response")
explicatives <- attr(mt, "term.labels")
nXs <- length(explicatives)
coef_norm <- coef * sdxy[explicatives] / sdxy[reponse]
cors <- m$cov / outer(sdxy, sdxy)
# print(reponse)
# print(explicatives)
cors_Y <- cors[reponse, explicatives]
TSS <- m$TSS
RSS <- sum(m$residuals ^ 2)
ESS <- TSS * abs(coef_norm * cors_Y)
names(ESS) <- explicatives
SumSq <- c(ESS, Residuals = RSS)
Df <- c(rep(1, nXs), m$df.residual)
MeanSq <- SumSq / Df
Fvalue <- c(MeanSq[-length(MeanSq)] / MeanSq[length(MeanSq)], NA)
Pvalue <- 1 - pf(Fvalue,
Df[-length(MeanSq)],
Df[length(MeanSq)])
result <- data.frame(Df,
SumSq,
MeanSq,
Fvalue,
Pvalue)
colnames(result) <-
c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
attr(result, "heading") <- paste("Analysis of Variance Table\n\nResponse:",
all.vars(m$terms)[1])
class(result) <- c("anova", class(result))
return(result)
}
#' Generate permutation matrix according to a schema.
#'
#' The permutation schema is defined using the `how` function.
#' The implementation of this function is currectly extraced
#' from the VEGAN package and just copied pasted here to avoid
#' dependency on an hidden vegan function.
#'
getPermuteMatrix <- function(perm, N, strata = NULL) {
if (length(perm) == 1) {
perm <- how(nperm = perm)
}
if (!missing(strata) && !is.null(strata)) {
if (inherits(perm, "how") && is.null(getBlocks(perm))) {
setBlocks(perm) <- strata
}
}
if (inherits(perm, "how")) {
perm <- shuffleSet(N, control = perm)
} else {
if (!is.integer(perm) && !all(perm == round(perm))) {
stop("permutation matrix must be strictly integers: use round()")
}
}
if (is.null(attr(perm, "control"))) {
attr(perm, "control") <- structure(list(
within = list(type = "supplied matrix"),
nperm = nrow(perm)
), class = "how")
}
perm
}
......@@ -7,6 +7,21 @@
#'
NULL
.MFclass <- function(x) {
if (is.logical(x))
return("logical")
if (is.ordered(x))
return("ordered")
if (is.factor(x))
return("factor")
if (is.character(x))
return("character")
if (is.matrix(x) && is.numeric(x))
return(paste0("nmatrix.", ncol(x)))
if (is.numeric(x))
return("numeric")
return("other")
}
.modelprocmodframe <- function(term, row.names,
variables, varnames,
......@@ -102,6 +117,12 @@ NULL
)
}
#' Construct Design ProcMod Frame.
#'
#' \code{model_procmod} creates a design (or model) \code{\link[ProcMod]{procmod_frame}},
#' e.g., by expanding factors to a set of dummy variables (depending on the contrasts)
#' and expanding interactions similarly.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
......@@ -111,13 +132,14 @@ model_procmod <- function(formula,...) {
#' @rdname model_procmod
#' @export
model_procmod.default <- function(formula,
model_procmod.terms <- function(formula,
data = NULL,
subset = NULL,
na.action = na.fail,
drop.unused.levels = FALSE,
xlev = NULL,
...) {
possible_newdata <- !missing(data) &&
is_procmod_frame(data) &&
identical(
......@@ -127,88 +149,21 @@ model_procmod.default <- function(formula,
(nr <- nrow(data)) > 0
if (!missing(formula) &&
nargs() == 1 &&
is.list(formula) &&
!is.null(formula$model)) {
return(formula$model)
}
if (!missing(formula)
&& nargs() == 1
&& is.list(formula)
&& all(c("terms", "call") %in% names(formula))) {
fcall <- formula$call
m <- match(c("formula", "data",
"subset", "weights",
"na.action"),
names(fcall),
0)
fcall <- fcall[c(1, m)]
fcall[[1L]] <- quote(model_procmod.default)
env <- environment(formula$terms)
if (is.null(env)) {
env <- parent.frame()
}
return(eval(fcall, env))
}
if (missing(formula)) {
if (!missing(data)
&& inherits(data, "procmod_frame")
&& length(attr(data, "terms"))) {
return(data)
}
formula <- as.formula(data)
}
else if (missing(data)
&& inherits(formula, "procmod.frame")) {
if (length(attr(formula, "terms"))) {
return(formula)
}
data <- formula
formula <- as.formula(data)
}
formula <- as.formula(formula)
#
# na.action is first set to the na.action parameter if provided
# secondly set to the na.action of data if defined
# secondly to the global na.action option if defined
# lastly to the default na.fail
#
if (missing(na.action)) {
if (!is.null(naa <- attr(data, "na.action")) &
mode(naa) != "numeric") {
mode(naa) != "numeric") {
na.action <- naa
} else if (!is.null(naa <- getOption("na.action"))) {
na.action <- naa
}
}
if (missing(data)) {
data <- environment(formula)
} else if (!is_procmod_frame(data) &&
!is.environment(data) &&
!is.null(attr(data, "class"))) {
data <- as_procmod_frame(data)
} else if (is.array(data)) {
data <- as_procmod_frame(data)
}
if (!inherits(formula, "terms")) {
formula <- terms(formula, data = data)
}
env <- environment(formula)
rownames <- .row_names_info(data, 0L)
......@@ -221,6 +176,8 @@ model_procmod.default <- function(formula,
predvars <- vars
}
varnames <- sapply(vars, .deparse2)[-1L]
variables <- eval(predvars, data, env)
......@@ -233,26 +190,22 @@ model_procmod.default <- function(formula,
}
if (possible_newdata &&
length(variables)) {
length(variables)) {
nr2 <- max(sapply(variables, NROW))
if (nr2 != nr) {
warning(sprintf(
paste0(
ngettext(
nr, "'newdata' had %d row",
"'newdata' had %d rows"
),
" ",
ngettext(
nr2,
"but variable found had %d row",
"but variables found have %d rows"
)
),
nr, nr2
warning(sprintf(paste0(
ngettext(nr, "'newdata' had %d row",
"'newdata' had %d rows"),
" ",
ngettext(
nr2,
"but variable found had %d row",
"but variables found have %d rows"
)
),
call. = FALSE, domain = NA
)
nr, nr2),
call. = FALSE,
domain = NA)
}
}
......@@ -324,6 +277,145 @@ model_procmod.default <- function(formula,
attr(formula, "dataClasses") <- vapply(data, .MFclass, "")
attr(data, "terms") <- formula
data
}
#' @rdname model_procmod
#' @export
model_procmod.formula <- function(formula,
data = NULL,
subset = NULL,
na.action = na.fail,
drop.unused.levels = FALSE,
xlev = NULL,
...) {
if (missing(data)) data <- environment(formula)
if (!is.environment(data)) data <- as_procmod_frame(data)
formula <- terms(formula, data = data)
env <- environment(formula)
if (is.null(env)) env <- parent.frame()
if (is.environment(data)) {
data <- as_procmod_frame(eval(attr(formula, "variables"),
envir = env))
}
fcall <- match.call()
fcall[[1]] <- quote(model_procmod)
fcall$formula <- quote(formula)
fcall$data <- quote(data)
eval(fcall)
}
#' @rdname model_procmod
#' @export
model_procmod.pm <- function(formula,...) {
if (nargs() == 1) {
#
# looks if provided formula contains
# a model and returns it
#
if (!is.null(formula$model)) return(formula$model)
#
# looks if provided formula contains
# a term and a call attributs
#
if (all(c("terms", "call") %in% names(formula))) {
fcall <- formula$call
m <- match(c("formula", "data",
"subset", "weights",
"na.action"),
names(fcall),
0)
fcall <- fcall[c(1, m)]
fcall[[1L]] <- quote(model_procmod)
env <- environment(formula$terms)
if (is.null(env)) {
env <- parent.frame()
}
return(eval(fcall, env))
}
}
}
#' @rdname model_procmod
#' @export
model_procmod.default <- function(formula,
data = NULL,
subset = NULL,
na.action = na.fail,
drop.unused.levels = FALSE,
xlev = NULL,
...) {
possible_newdata <- !missing(data) &&
is_procmod_frame(data) &&
identical(
substitute(data),
quote(newdata)
) &&
(nr <- nrow(data)) > 0
if (missing(formula)) {
if (!missing(data)
&& inherits(data, "procmod_frame")
&& length(attr(data, "terms"))) {
return(data)
}
formula <- as.formula(data)
}
else if (missing(data)
&& inherits(formula, "procmod.frame")) {
if (length(attr(formula, "terms"))) {
return(formula)
}
data <- formula
formula <- as.formula(data)
}
formula <- as.formula(formula)
if (missing(na.action)) {
if (!is.null(naa <- attr(data, "na.action")) &
mode(naa) != "numeric") {
na.action <- naa
} else if (!is.null(naa <- getOption("na.action"))) {
na.action <- naa
}
}
if (missing(data)) {
data <- environment(formula)
} else if (!is_procmod_frame(data) &&
!is.environment(data) &&
!is.null(attr(data, "class"))) {
data <- as_procmod_frame(data)
} else if (is.array(data)) {
data <- as_procmod_frame(data)
}
}
......
......@@ -126,7 +126,7 @@ pm <- function(formula, data, subset, weights, na_action, method = "qr",
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- quote(ProcMod:::model_procmod.default)
mf[[1L]] <- quote(model_procmod)
mf <- eval(mf, parent.frame())
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/anova_pm.R
\name{anova.pm}
\alias{anova.pm}
\title{Anova function for `pm` objects.}
\usage{
\method{anova}{pm}(m, permutations = how(nperm = 999), ...)
}
\description{
Anova function for `pm` objects.
}
\author{
Eric Coissac
Christelle Gonindard-Melodelima
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/anova_pm.R
\name{getPermuteMatrix}
\alias{getPermuteMatrix}
\title{Generate permutation matrix according to a schema.}
\usage{
getPermuteMatrix(perm, N, strata = NULL)
}
\description{
The permutation schema is defined using the `how` function.
The implementation of this function is currectly extraced
from the VEGAN package and just copied pasted here to avoid
dependency on an hidden vegan function.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/model_procmod.R
\name{model_procmod}
\alias{model_procmod}
\alias{model_procmod.terms}
\alias{model_procmod.formula}
\alias{model_procmod.pm}
\alias{model_procmod.default}
\title{Construct Design ProcMod Frame.}
\usage{
model_procmod(formula, ...)
\method{model_procmod}{terms}(formula, data = NULL, subset = NULL,
na.action = na.fail, drop.unused.levels = FALSE, xlev = NULL, ...)
\method{model_procmod}{formula}(formula, data = NULL, subset = NULL,
na.action = na.fail, drop.unused.levels = FALSE, xlev = NULL, ...)
\method{model_procmod}{pm}(formula, ...)
\method{model_procmod}{default}(formula, data = NULL, subset = NULL,
na.action = na.fail, drop.unused.levels = FALSE, xlev = NULL, ...)
}
\description{
\code{model_procmod} creates a design (or model) \code{\link[ProcMod]{procmod_frame}},
e.g., by expanding factors to a set of dummy variables (depending on the contrasts)
and expanding interactions similarly.
}
\author{
Eric Coissac
Christelle Gonindard-Melodelima
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment