With an potentially working anova

DESCRIPTION

@@ -20,4 +20,6 @@ Collate:
     'omit.action.R'
     'model.procmod.R'
     'mprocuste.R'
+    'anova.pm.R'
+    'mcov.R'
     'plot.pm.R'

NAMESPACE

@@ -23,6 +23,7 @@ S3method(subset,procmod.frame)
 export(as.procmod.frame)
 export(is.procmod.frame)
 export(mcor)
+export(mcor.partial)
 export(model.procmod.default)
 export(mvar)
 export(pm)

R/anova.pm.R

+#' @include mprocuste.R
+
+#' Anova function for `pm` objects.
+#'
+#' @author Eric Coissac
+#' @author Christelle Gonindard-Melodelima
+#' @export
+anova.pm = function(m,...) {
+
+  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 = mcor(mf)
+
+  print(reponse)
+  print(explicatives)
+
+  cors.Y = cors[reponse,explicatives]
+
+  SCT  = m$SCT
+  SCR = sum(m$residuals^2)
+
+  SCE = SCT * coef.norm * cors.Y
+
+  names(SCE) = explicatives
+
+  SCET = sum(SCE)
+
+  SCIT = SCT - SCET - SCR
+
+  SumSq  = c(SCE,Residuals=SCR)
+  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)
+}

R/mcov.R

@@ -24,32 +24,28 @@ NULL
 mvar = function(...) {
 
   Xs <- list(...)
-  if (length(Xs)==1)
+  if (length(Xs)==1) {
     if (is.list(Xs[[1]]))
       Xs=as.procmod.frame(Xs[[1]])
     else if (is.procmod.frame(Xs[[1]]))
       Xs=Xs[[1]]
     else
-      Xs=procmod.frame(Xs[[1]])
-    else
+      Xs=procmod.frame(Xs[[1]]) }
+  else
       Xs=as.procmod.frame(Xs)
 
     Xnames=names(Xs)
 
-    Xs <- lapply(Xs,scale,scale = FALSE)
+    Xs <- as.procmod.frame(lapply(Xs,scale,scale = FALSE))
 
     nX = length(Xs)
 
     Xx <- rep(1:nX,nX)
     Xy <- rep(1:nX,rep(nX,nX))
 
-    XXs <- mapply(function(x,y) crossprod(Xs[[x]], Xs[[y]]),
+    CovXXs <- mapply(function(x,y) sum(svd(crossprod(Xs[[x]], Xs[[y]]))$d),
                   Xx,Xy,
-                  SIMPLIFY = FALSE)
-
-    sol_xxs <- lapply(XXs,svd)
-
-    CovXXs  = sapply(sol_xxs, function(sol) sum(sol$d))
+                  SIMPLIFY = TRUE) / (nrow(Xs)-1)
 
 
     dim(CovXXs)=c(nX,nX)

R/mprocuste.R

@@ -8,92 +8,27 @@
 NULL
 
 
-#' Compute the variance, covariance matrix of K coordinate matrices.
-#'
-#' Covariance between two matrices is defined as the sum of the
-#' sigular values of the X'Y matrix. All the matrices must have
-#' the same number of rows.
-#'
-#' @param ... the set of matrices
-#'
-#' @examples
-#'    # Build Three matrices of 3 rows.
-#'    A <- matrix(1:9,nrow=3)
-#'    B <- matrix(10:15,nrow=3)
-#'    C <- matrix(20:31,nrow=3)
-#'    # compute the variance covariance matrix
-#'    mvar(A,B,C)
-#'    mvar(A=A,B=B,C=C)
-#'
-#' @author Eric Coissac & Christelle Gonindard-Melodelima
-#' @export
-mvar = function(...) {
-
-  Xs <- list(...)
-  if (length(Xs)==1)
-    if (is.list(Xs[[1]]))
-      Xs=as.procmod.frame(Xs[[1]])
-    else if (is.procmod.frame(Xs[[1]]))
-      Xs=Xs[[1]]
-    else
-      Xs=procmod.frame(Xs[[1]])
-  else
-    Xs=as.procmod.frame(Xs)
-
-  Xnames=names(Xs)
-
-  Xs <- lapply(Xs,scale,scale = FALSE)
-
-  nX = length(Xs)
-
-  Xx <- rep(1:nX,nX)
-  Xy <- rep(1:nX,rep(nX,nX))
-
-  XXs <- mapply(function(x,y) crossprod(Xs[[x]], Xs[[y]]),
-                Xx,Xy,
-                SIMPLIFY = FALSE)
-
-  sol_xxs <- lapply(XXs,svd)
-
-  CovXXs  = sapply(sol_xxs, function(sol) sum(sol$d))
-
-
-  dim(CovXXs)=c(nX,nX)
-  colnames(CovXXs)=Xnames
-  rownames(CovXXs)=Xnames
-
-  return(CovXXs)
-}
-
-#' Compute the person correlation matrix of K coordinate matrices
-#'
-#' @author Eric Coissac
-#' @author Christelle Gonindard-Melodelima
-#' @export
-mcor = function(...) {
-  cov = mvar(...)
-  s = sqrt(diag(cov))
-  vv= outer(s,s)
-  return(cov/vv)
-}
-
 pm.fit = function(covmat,y,xs,
                   singular.ok = singular.ok,
                   tol) {
 
-  xy.cov = covmat[xs,y]
+  xy.cov = covmat[xs,y,drop=FALSE]
   xx.cov = covmat[xs,xs,drop=FALSE]
 
-  xx.qr  = qr(xx.cov,tol=tol)
+  #xx.qr  = qr(xx.cov,tol=tol)
 
-  if (!singular.ok && xx.qr$rank < ncol(xx.cov))
-    stop("singular fit encountered")
+  #if (!singular.ok && xx.qr$rank < ncol(xx.cov))
+  #  stop("singular fit encountered")
 
-  coef   = qr.coef(xx.qr,xy.cov)
+  coef = as.numeric(solve(xx.cov) %*% xy.cov)
+  names(coef)=rownames(xx.cov)
+  #coef   = qr.coef(xx.qr,xy.cov)
+  #effects= qr.qty(xx.qr,xy.cov)
 
 
-  return(list(coefficients=coef,
-              rank=xx.qr$rank + 1
+  return(list(coefficients=coef
+              #rank=xx.qr$rank + 1
+              #qr=xx.qr,
              )
         )
 }
@@ -167,7 +102,9 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
     variances= diag(data.cov)
     std.dev = sqrt(variances)
 
-    vars.norm = as.procmod.frame(mapply(function(x) scale(x,scale = FALSE), vars))
+    vars.norm = as.procmod.frame(mapply(function(x) scale(x,scale = FALSE),
+                                        vars,
+                                        SIMPLIFY = FALSE))
 
     if (is.null(w)) {
       subset.w=rep(TRUE,nvars)
@@ -188,6 +125,8 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
                singular.ok = singular.ok,
                tol = tol)
 
+    z$cov = data.cov
+
     y.sd = std.dev[irep]
     xs.sd= std.dev[seq_len(nvars)!=irep]
 
@@ -195,7 +134,7 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
 
 
     Y = vars.norm[[irep]]
-    Xs= vars.norm[seq_len(nvars)!=irep]
+    Xs= vars.norm[seq_len(nvars)!=irep,drop=FALSE]
 
     Ymean <- colMeans(Y)
     Xmeans<- lapply(Xs,colMeans)
@@ -239,14 +178,14 @@ pm = function (formula,data, subset, weights, na.action, method = "qr",
     z$residuals = Y - yhat
     z$fitted.values = yhat
     z$weights = w
-    z$df.residual = if (!is.null(w))
+    z$df.residual = (if (!is.null(w))
                       sum(w != 0)
-                    else nrow(Y)
+                     else nrow(Y)) - nvars + 1
 
     if (ret.x)
       z$x <- vars.norm[-irep]
 
-    if (!A)
+##    if (!A)
       z$A <- A_xys
   }
 
@@ -283,53 +222,6 @@ print.pm = function(m,...) {
 }
 
 
-#' Anova function for `pm` objects.
-#'
-#' @author Eric Coissac
-#' @author Christelle Gonindard-Melodelima
-#' @export
-anova.pm = function(m,...) {
-
-  VarY  = m$SCT
-  SCE = sum(m$residuals^2)
-
-
-  VarYhats = VarY * m$varpart
-  names(VarYhats) = names(m$coefficients)
-  SChat = sum(VarYhats)
-
-
-  if (length(Xs)>1) {
-    SCI = VarY - SCE - sum(SChat)
-    SumSq  = c(VarYhats,Interactions=SCI,Residuals=SCE)
-    Df     = c(rep(1,length(VarYhats)+1),m$df.residual)
-  }
-  else {
-    SumSq  = c(VarYhats,Residuals=SCE)
-    Df     = c(1,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)
-}
 
 summary.pm = function(object,
                       correlation = FALSE,

vignettes/ProcMod.Rmd

@@ -512,7 +512,7 @@ knitr::kable(cors)
 #### Building the multiprocruste model
 
 ```{r}
-euk.pm = pm(euk ~ soil + climat +  geo + hist,data=data)
+euk.pm = pm(euk ~ soil + climat +  geo + hist ,data=data)
 euk.pm
 ```
 
@@ -525,7 +525,7 @@ plot(euk.pm)
 ```{r}
 W=1/rowSums(euk.pm$residuals^2)
 W=W/max(W)
-euk.pm.w = pm(euk ~ soil + climat +  geo + hist,data=data,weights = W)
+euk.pm.w = pm(euk ~ soil + climat +  geo + hist ,data=data,weights = W)
 euk.pm.w
 
 ```
@@ -580,7 +580,7 @@ partition
 ```
 
 ```{r}
-soil.pm = pm( hist ~ bac + euk + climat +  geo + soil,data=data)
+soil.pm = pm( soil ~ bac + euk + climat +  geo,data=data)
 soil.pm
 ```
 
@@ -591,6 +591,7 @@ plot(soil.pm)
 ```{r}
 soil.anova = anova(soil.pm)
 soil.anova
+soil.anova$`Sum Sq`/sum(soil.anova$`Sum Sq`)
 ```
 
 
@@ -602,3 +603,35 @@ pm(eukariote ~ bacterie)
 cor(bac.dist,euk.dist)
 ```
 
+```{r}
+X1 = rnorm(10,0,2)
+X2 = rnorm(10,0,2)
+X3 = rnorm(10,0,2)
+if (cor(X1,X2)< 0)
+  X2=-X2
+if (cor(X1,X3)< 0)
+  X3=-X3
+X1C = scale(X1,scale = TRUE)
+X2C = scale(X2,scale = TRUE)
+X3C = scale(X3,scale = TRUE)
+Y = 1 *X1C+ 2 *X2C + 3 * X3C + rnorm(10)
+YC = scale(Y,scale = TRUE)
+XXXX = procmod.frame(YC,X1C,X2C,X3C)
+
+L = lm(YC~X1C+X2C+X3C,data = XXXX)
+P = pm(YC~X1C+X2C+X3C,data = XXXX)
+mcor(XXXX)
+cor(do.call(data.frame,XXXX))
+anova(L)
+anova(P)
+L
+P
+```
+
+```{r}
+y = c(4,2,3,4,5,5,7,5,4,9,7,6)
+x1= c(1,1,1,1,1,1,2,2,2,2,2,2)
+x2= c(8,7,7,9,5,4,3,6,7,2,3,2) 
+scherrer = procmod.frame(y,x1,x2)
+```
+