Commit 433144c8 by Eric Coissac

Add a first version of an anova.pm function and clean a little the plot.pm function

parent 73fdb570
......@@ -162,6 +162,10 @@ print.pm = function(m,...) {
#' @author Christelle Gonindard-Melodelima
#' @export
plot.pm = function(m,...) {
points.cex=0.6
labels.cex=0.6
xlim = c(min(m$model[[1]][,1],m$fitted.values[,1]),
max(m$model[[1]][,1],m$fitted.values[,1])
)
......@@ -169,36 +173,74 @@ plot.pm = function(m,...) {
max(m$model[[1]][,2],m$fitted.values[,2])
)
plot(m$model[[1]][,1:2],
xlim=xlim,ylim=ylim,
cex=0)
centers = (m$model[[1]] + m$fitted.values)/2
text(m$model[[1]][,1],m$model[[1]][,2],
rownames(m$model[[1]]),
plot(m$model[[1]][,1:2],
xlim=xlim,ylim=ylim,
cex=0.8,
cex=points.cex,
col="blue")
text(m$fitted.values[,1],m$fitted.values[,2],
rownames(m$fitted.values),
points(m$fitted.values[,1:2],
xlim=xlim,ylim=ylim,
cex=0.8,
cex=points.cex,
col="red")
arrows(m$fitted.values[,1],m$fitted.values[,2],
m$model[[1]][,1],m$model[[1]][,2],
length = 0.05,
col="black")
# text(m$model[[1]][,1],m$model[[1]][,2],
# rownames(m$model[[1]]),
# xlim=xlim,ylim=ylim,
# cex=labels.cex,
# col="blue",
# pos=3)
text(centers[,1],centers[,2],
rownames(m$fitted.values),
xlim=xlim,ylim=ylim,
cex=labels.cex,
col="black",
pos=3)
}
anova.pm = function(m,...) {
ctrace <- function(MAT) sum(MAT^2)
Xs <- lapply(m$model,scale,scale = FALSE)
Y <- Xs[[1]]
Xs <- Xs[-1]
Ymean <- colMeans(Y)
Xrots <- mapply(function(x,a) x %*% a,
Xs,m$A,
SIMPLIFY = FALSE)
Chat = mapply(function(coef,xrot) sweep(coef * xrot,
2,
Ymean,
'-'),
m$coefficients, Xrots,
SIMPLIFY = FALSE)
VarY = ctrace(Y)
VarYhats = sapply(Chat,ctrace)
names(VarYhats) = names(m$coefficients)
SChat = sum(VarYhats)
SCE = sum(m$residuals^2)
SCI = VarY - SCE - sum(SChat)
return(list(SCI=SCI,
SCX=VarYhats,
SCE=SCE))
}
#' Estimate a procuste regression of several table on a reference one.
#'
#' @author Eric Coissac
#' @author Christelle Gonindard-Melodelima
#' @export
mprocuste = function (Y, ...)
{
ctrace <- function(MAT) sum(MAT^2)
......
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