R/varcomp.R

   1 ## varcomp.R (2004-10-29)
   2
   3 ##   Variance Component of Mixed-Effect Linear Model
   4
   5 ## Copyright 2004 Julien Dutheil
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 varcomp <- function(x, scale = FALSE, cum = FALSE)
  11 {
  12   if (!("lme" %in% class(x))) stop("Object \"x\" is not of class \"lme\"")
  13   res <- seq(along = x$modelStruct$reStruct)
  14   var <- vector(length = length(res) + 1)
  15   for(i in res) {
  16     var[length(var) - i] <- attr(summary(x$modelStruct$reStruct[[i]]),"stdDev")[1]*x$sigma
  17   }
  18   var[length(var)] <- x$sigma
  19   var <- var^2
  20   if(scale) var <- var/sum(var)
  21   if(cum) var <- cumsum(var)
  22   names(var) <- c(rev(names(x$modelStruct$reStruct)), "Within")
  23   class(var) <- "varcomp"
  24   return(var)
  25 }
  26
  27 plot.varcomp <- function(x, xlab = "Levels", ylab = "Variance", type = "b", ...) {
  28   if (!("varcomp" %in% class(x))) stop("Object \"x\" is not of class \"varcomp\"")
  29   return(lattice::xyplot(x ~ ordered(names(x), levels=rev(names(x))), xlab=xlab, ylab=ylab, type=type, ...))
  30 }
  31
  32 # For debuging:
  33 #data(carnivora)
  34 #m <- lme(log10(SW) ~ 1, random = ~ 1|Order/SuperFamily/Family/Genus, data=carnivora)
  35 #v <- varcomp(m,T,T)
  36 #plot(v)
  37