Skip to content

P
ProcMod
Commit 31958291 by Eric Coissac
Options

New version of pm more similar to lm

parent eda9e6c1
Showing with 391 additions and 30 deletions
DESCRIPTION
......@@ -16,5 +16,8 @@ VignetteBuilder: knitr
Collate:
'internals.R'
'formula.procmod.frame.R'
'mprocuste.R'
'procmod.frame.R'
'omit.action.R'
'model.procmod.R'
'mprocuste.R'
'plot.pm.R'
NAMESPACE
......@@ -14,6 +14,8 @@ S3method(deviance,pm)
S3method(dim,procmod.frame)
S3method(extractAIC,pm)
S3method(formula,procmod.frame)
S3method(na.exclude,procmod.frame)
S3method(na.omit,procmod.frame)
S3method(plot,pm)
S3method(print,pm)
S3method(residuals,pm)
......@@ -21,7 +23,9 @@ S3method(subset,procmod.frame)
export(as.procmod.frame)
export(is.procmod.frame)
export(mcor)
export(model.procmod.default)
export(mvar)
export(pm)
export(procmod.frame)
export(vars.procmod)
export(weighted.residuals)
R/model.procmod.R 0 → 100644
#' @include internals.R
#' @include omit.action.R
#' @include procmod.frame.R
#' @include formula.procmod.frame.R
#'
NULL
.modelprocmodframe=function(term, row.names,
variables,varnames,
dots,dotnames,
subset, na.action) {
stopifnot(is.list(variables))
stopifnot(is.character(varnames))
nvars = length(variables)
stopifnot(nvars == length(varnames))
stopifnot(is.list(dots))
ndots = length(dots)
stopifnot(ndots==length(dotnames))
nactualdots = 0;
for (i in seq_len(ndots))
if (! is.null(dots[[i]])) nactualdots=nactualdots+1
data = vector(mode="list", nvars + nactualdots)
names= character(nvars + nactualdots)
for (i in seq_len(nvars)) {
data[[i]]=variables[[i]]
names[i]=varnames[i]
}
j=1
for (i in seq_len(ndots))
if (! is.null(dots[[i]])) {
if (nchar(dotnames[i]) + 3 > 256)
stop(sprintf("overlong names in '%s'",dotnames[i]))
buf=paste('(',dotnames[i],')',sep="")
data[nvars + j]=dots[[i]]
names[nvars + j]=buf
j=j+1
}
names(data)=names
data$row.names=row.names
data = do.call(ProcMod::procmod.frame,data)
if (! is.null(subset) ) {
if (is.numeric(subset))
subset = which(seq_len(nrow(data)) %in% subset)
if (is.character(subset))
subset = which(rownames(data) %in% subset)
data = subset(data,subset=subset)
}
if ( ! is.null(na.action)) {
# some na.actions need this to distinguish responses from
# explanatory variables
attr(data,"terms",terms)
if ( is.character(na.action) && length(na.action) > 0)
na.action = match.fun(na.action[1]);
ans = do.call(na.action,list(data))
if ( !is.procmod.frame(ans) ||
length(ans) != length(data))
stop("invalid result from na.action")
}
else
ans = data;
return(ans)
}
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @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) &&
nargs() == 1 &&
is.list(formula) &&
!is.null(m <- formula$model))
return(m)
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(procmod::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)
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)
if (!inherits(formula, "terms"))
formula <- terms(formula, data = data)
env <- environment(formula)
rownames <- .row_names_info(data, 0L)
vars <- attr(formula, "variables")
predvars <- attr(formula, "predvars")
if (is.null(predvars))
predvars <- vars
varnames <- sapply(vars, function(x) paste(deparse(x, width.cutoff = 500),
collapse = " "))[-1L]
variables <- eval(predvars, data, env)
resp <- attr(formula, "response")
if (is.null(rownames) && resp > 0L) {
lhs <- variables[[resp]]
rownames <- rownames(lhs)
}
if (possible_newdata &&
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),
call. = FALSE, domain = NA)
}
if (is.null(attr(formula, "predvars"))) {
for (i in seq_along(varnames))
predvars[[i + 1L]] <- makepredictcall(variables[[i]],
vars[[i + 1L]])
attr(formula, "predvars") <- predvars
}
extras <- substitute(list(...))
extranames <- names(extras[-1L])
extras <- eval(extras, data, env)
subset <- eval(substitute(subset), data, env)
data = .modelprocmodframe(formula, rownames,
variables,varnames,
extras, extranames,
subset, na.action)
# if (length(xlev)) {
# for (nm in names(xlev)) if (!is.null(xl <- xlev[[nm]])) {
# xi <- data[[nm]]
# if (is.character(xi))
# xi <- as.factor(xi)
# if (!is.factor(xi) || is.null(nxl <- levels(xi)))
# warning(gettextf("variable '%s' is not a factor",
# nm), domain = NA)
# else {
# ctr <- attr(xi, "contrasts")
# xi <- xi[, drop = TRUE]
# nxl <- levels(xi)
# if (any(m <- is.na(match(nxl, xl))))
# stop(sprintf(ngettext(length(m), "factor %s has new level %s",
# "factor %s has new levels %s"), nm, paste(nxl[m],
# collapse = ", ")), domain = NA)
# data[[nm]] <- factor(xi, levels = xl, exclude = NULL)
# if (!identical(attr(data[[nm]], "contrasts"),
# ctr))
# warning(gettext(sprintf("contrasts dropped from factor %s",
# nm), domain = NA), call. = FALSE)
# }
# }
# }
# else if (drop.unused.levels) {
# for (nm in names(data)) {
# x <- data[[nm]]
# if (is.factor(x) && length(unique(x[!is.na(x)])) <
# length(levels(x))) {
# ctr <- attr(x, "contrasts")
# data[[nm]] <- x[, drop = TRUE]
# if (!identical(attr(data[[nm]], "contrasts"),
# ctr))
# warning(gettext(sprintf("contrasts dropped from factor %s due to missing levels",
# nm), domain = NA), call. = FALSE)
# }
# }
# }
attr(formula, "dataClasses") <- vapply(data, .MFclass, "")
attr(data, "terms") <- formula
data
}
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
vars.procmod = function(object, data = environment(object),
xlev = NULL, ...) {
t <- if (missing(data))
terms(object)
else terms(object, data = data)
if (is.null(attr(data, "terms")))
data <- model.procmod.default(object, data, xlev = xlev)
else {
deparse2 <- function(x) paste(deparse(x, width.cutoff = 500L),
collapse = " ")
reorder <- match(vapply(attr(t, "variables"), deparse2,
"")[-1L], names(data))
if (anyNA(reorder))
stop("model frame and formula mismatch in model.matrix()")
if (!identical(reorder, seq_len(ncol(data))))
data <- data[, reorder, drop = FALSE]
}
attr(data,"response")=attr(t, "response")
return(data)
}
R/plot.pm.R 0 → 100644
#' Default plot function for `pm` objects.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
plot.pm = function(m,
points.cex=0.6,
labels.cex=0.6,
arrow.length=0.05,
asp=1,
...) {
if (is.null(m$y))
stop("Cannot plot procrust model without the y values")
layout(matrix(c(1,1,1,1,2,3),nrow = 2))
xlim = c(min(m$y[,1],m$fitted.values[,1]),
max(m$y[,1],m$fitted.values[,1])
)
ylim = c(min(m$y[,2],m$fitted.values[,2]),
max(m$y[,2],m$fitted.values[,2])
)
centers = (m$y + m$fitted.values)/2
plot(m$y[,1:2],
xlim=xlim,ylim=ylim,
cex=points.cex,
col="blue",
asp=asp)
points(m$fitted.values[,1:2],
xlim=xlim,ylim=ylim,
cex=points.cex,
col="red")
arrows(m$fitted.values[,1],m$fitted.values[,2],
m$y[,1],m$y[,2],
length = arrow.length,
col="black")
text(centers[,1],centers[,2],
rownames(m$fitted.values),
xlim=xlim,ylim=ylim,
cex=labels.cex,
col="black",
pos=3)
anova = anova(m)
v=anova[,"Sum Sq"]
names(v)=rownames(anova)
barplot(v/sum(v),las=2)
}
R/procmod.frame.R
......@@ -86,6 +86,10 @@ procmod.frame = function(...,
x <- list(...)
n <- length(x)
# message(x)
# message(n)
# message(row.names)
if ((! has.row.names || is.null(row.names)) && n >= 1)
row.names = rownames(x[[1]])
......@@ -113,7 +117,7 @@ procmod.frame = function(...,
}
stopifnot(all(nrows[i]==nrows))
message(row.names," : ",length(row.names),",",nrows[i])
# message(row.names," : ",length(row.names),",",nrows[i])
if (length(row.names)==nrows[i]) {
attr(value,"row.names")=row.names
if (check.rows)
......@@ -157,7 +161,8 @@ as.procmod.frame = function(data,...) {
#' @author Christelle Gonindard-Melodelima
#' @export
as.procmod.frame.list = function(data,...) {
do.call(procmod.frame,data,...)
data=c(data,list(...))
do.call(procmod.frame,data)
}
#' Coerce to a ProcMod Frame.
......@@ -184,7 +189,9 @@ as.procmod.frame.array = function(data,...) {
if (length(attr(data,"dimnames"))==3)
names(l)=attr(data,"dimnames")[[3]]
do.call(procmod.frame,l,...)
data=c(l,list(...))
do.call(procmod.frame,l)
}
#' Dimensions of a Matrix Frame.
......@@ -249,7 +256,7 @@ dim.procmod.frame = function(x)
has.drop = !missing(drop)
Narg = nargs() - 2 + (has.i | has.j | has.drop) - has.drop
message("Nargs = ",Narg," i:",has.i," j:",has.j," drop:",has.drop)
# message("Nargs = ",Narg," i:",has.i," j:",has.j," drop:",has.drop)
if (!all(names(sys.call()) %in% c("", "drop")))
warning("named arguments other than 'drop' are discouraged")
......@@ -260,21 +267,21 @@ dim.procmod.frame = function(x)
return(x)
else if (! has.i && has.j) {
# Case 2 : X[,j] ou x[i]
message('Case 2 : X[,j]')
# message('Case 2 : X[,j]')
y <- as.list(x)[j]
}
else if ( has.i && Narg==1) {
# Case 3 : X[,j] ou x[i]
message('Case 3 : X[i]')
# message('Case 3 : X[i]')
y <- as.list(x)[i]
}
else if ( has.i && !has.j && Narg>1) {
# Case 4 : X[i,]
message('Case 4 : X[i,]')
# message('Case 4 : X[i,]')
y = lapply(x, function(m) m[i,,drop=FALSE])
}
else if ( has.i && has.j ) {
message('Case 5 : X[i,j]')
# message('Case 5 : X[i,j]')
y = x[j,drop=FALSE]
y = y[i,,drop=FALSE]
}
......
vignettes/.gitignore
ProcMod.html
ProcMod.R
ProcMod_files
vignettes/ProcMod.Rmd
......@@ -427,8 +427,8 @@ mat.scale = function(mat) {
Coordinates of the sites are extracted from the PCoA analysis
```{r}
euk.pco.li = mat.scale(euk.pco$li)
bac.pco.li = mat.scale(bac.pco$li)
euk.pco.li = euk.pco$li
bac.pco.li = bac.pco$li
```
```{r, fig.show='hold'}
......@@ -449,7 +449,7 @@ The environmental variable are transformed using a Principal Component Analysis.
```{r}
env.pca = dudi.pca(env,scannf = FALSE,nf=nrow(env)-1)
env.pca.li = mat.scale(env.pca$li)
env.pca.li = env.pca$li
plot(env.pca.li[,1:2],cex=0,
main="Environmental data")
text(env.pca.li[,1:2],
......@@ -457,13 +457,13 @@ text(env.pca.li[,1:2],
cex=0.4)
soil.pca = dudi.pca(soil,scannf = FALSE,nf=nrow(soil)-1)
soil.pca.li = mat.scale(soil.pca$li)
soil.pca.li = soil.pca$li
climat.pca = dudi.pca(climat,scannf = FALSE,nf=nrow(climat)-1)
climat.pca.li = mat.scale(climat.pca$li)
climat.pca.li = climat.pca$li
geo.pco = dudi.pco(dist(geo),full = TRUE)
geo.pco.li = mat.scale(geo.pco$li)
geo.pco.li = geo.pco$li
```
......@@ -473,7 +473,7 @@ dh = matrix(1,nrow = 62,ncol=62)
diag(dh)=0
dh=as.dist(dh)
dh.pco =dudi.pco(dh,full = TRUE)
dh.pco.li = mat.scale(dh.pco$li)
dh.pco.li = dh.pco$li
```
### Using the package to analyse relationship among the tables
......@@ -487,9 +487,11 @@ library(ProcMod)
#### Computing the variance/covariance matix
```{r}
vars = mvar(euk=euk.pco.li,bac=bac.pco.li,
climat=climat.pca.li,soil=soil.pca.li,
geo=geo.pco.li,hist=dh.pco.li)
data = procmod.frame(euk=euk.pco.li,bac=bac.pco.li,
climat=climat.pca.li,soil=soil.pca.li,
geo=geo.pco.li,hist=dh.pco.li)
vars = mvar(data)
```
```{r echo=FALSE}
......@@ -500,13 +502,7 @@ knitr::kable(vars)
#### Computing the correlation matix
```{r}
cors = mcor(euk=euk.pco.li,bac=bac.pco.li,
climat=climat.pca.li,soil=soil.pca.li,
geo=geo.pco.li,hist=dh.pco.li)
```
```{r echo=FALSE}
knitr::kable(cors)
cors = mcor(data)
```
```{r echo=FALSE}
......@@ -516,7 +512,7 @@ knitr::kable(cors)
#### Building the multiprocruste model
```{r}
euk.pm = pm(euk.pco.li ~ soil.pca.li + climat.pca.li + geo.pco.li + dh.pco.li)
euk.pm = pm(euk ~ soil + climat + geo + hist,data=data)
euk.pm
```
......@@ -543,8 +539,7 @@ partition
```{r}
bac.pm = pm(bac.pco.li ~ soil.pca.li + climat.pca.li + geo.pco.li + dh.pco.li)
#bac.pm = pm(bac.pco.li ~ euk.pco.li + dh.pco.li )
bac.pm = pm(bac ~ soil + climat + geo + hist,data=data)
bac.pm
```
......@@ -564,7 +559,7 @@ partition
```
```{r}
soil.pm = pm( dh.pco.li ~ bac.pco.li + euk.pco.li + climat.pca.li + geo.pco.li + soil.pca.li)
soil.pm = pm( hist ~ bac + euk + climat + geo + soil,data=data)
soil.pm
```
......
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