version 0.1-18

cran · Feb 22, 2011 · 6646899 · 6646899
1 parent cf7f8da
commit 6646899
Show file tree

Hide file tree

Showing 20 changed files with 295 additions and 154 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,23 +1,22 @@
 Package: abd
 Type: Package
 Title: The Analysis of Biological Data
-Version: 0.1-12
-Date: 2011-1-12
+Version: 0.1-18
+Date: 2011-02-22
 Author: Kevin M. Middleton <[email protected]>, Randall Pruim
         <[email protected]>
-Maintainer: Kevin M. Middleton <[email protected]>, Randall Pruim
-        <[email protected]>
+Maintainer: Kevin M. Middleton <[email protected]>
 Depends: R (>= 2.10), nlme, lattice, grid
 Suggests: boot, car, ggplot2, HH, plyr, vcd, Hmisc
-Description: The abd package contains data sets and sample code for the
-        book, The Analysis of Biological Data by Michael Whitlock and
-        Dolph Schluter (2009; Roberts & Company Publishers).
+Description: The abd package contains data sets and sample code for The
+        Analysis of Biological Data by Michael Whitlock and Dolph
+        Schluter (2009; Roberts & Company Publishers).
 License: GPL-2
 LazyLoad: yes
 LazyData: yes
 Encoding: UTF-8
 Repository: CRAN
 Repository/R-Forge/Project: abd
-Repository/R-Forge/Revision: 142
-Date/Publication: 2011-01-26 08:24:30
-Packaged: 2011-01-12 21:46:40 UTC; rforge
+Repository/R-Forge/Revision: 150
+Date/Publication: 2011-04-05 07:14:39
+Packaged: 2011-02-23 22:30:37 UTC; rforge
diff --git a/R/as.xtabs.R b/R/as.xtabs.R
@@ -0,0 +1,37 @@
+as.xtabs <- function(x, ...) { UseMethod('as.xtabs') }
+
+as.xtabs.data.frame <- function(x, rowvar=NULL, colvar=NULL, labels=1, ...) {
+
+  if (labels >= 1) {
+  	cnames <- names(x)[-1]
+  	m <- as.matrix(x[,-c(labels)])
+  } else {
+  	cnames <- names(x)
+	m <- as.matrix(x)
+  }
+
+  rnames <- x[,labels]
+  rownames(m) <- rnames
+  if (! is.character(rowvar) ) { rowvar <- "variable.1" }
+  if (! is.character(colvar) ) { colvar <- "variable.2" }
+
+  dn <- list( rnames, cnames)
+  names(dn) <- c(rowvar, colvar)
+  attr(m,'dimnames') <- dn
+  class(m) <- c('xtabs', 'table')
+  return(m)
+}
+
+as.xtabs.matrix <- function(x, rowvar=NULL, colvar=NULL, ...)  {
+  rnames <- rownames(x)
+  cnames <- colnames(x)
+  rownames(m) <- rnames
+  if (! is.character(rowvar) ) { rowvar <- "variable.1" }
+  if (! is.character(colvar) ) { colvar <- "variable.2" }
+
+  dn <- list( rnames, cnames)
+  names(dn) <- c(rowvar, colvar)
+  attr(m,'dimnames') <- dn
+  class(m) <- c('xtabs', 'table')
+  return(m)
+}
diff --git a/R/repeatability.R b/R/repeatability.R
@@ -3,25 +3,45 @@ repeatability <- function(x){
     var.among <- x[["var.among"]]
     var.within <- x[["var.within"]] 
     R <- var.among / (var.among + var.within)
+    class(R) <- c("repeatability", "rep.varcomps")
+  }
+
+  ## Note: could use ICC1 from multilevel package, but need to code
+  ## dv as factor. This value matches ICC from repeatability.lme
+  ## and from ICC1.lme, but not repeatability.varcomps above.
+  if(inherits(x, "aov")){
+    smry <- summary(x)
+    n <- smry[[1]]$Df[1] + 1
+    MSw <- smry[[1]]$"Mean Sq"[2]
+    MSa <- smry[[1]]$"Mean Sq"[1]
+    s2a <- (MSa - MSw)/n					
+    R <- s2a/(MSw + s2a)
     class(R) <- c("repeatability", "rep.aov")
   }
+
+  ## Note: could use ICC1.lme from psychometric package. The
+  ## results here is numerically equal.
   if(inherits(x, "lme")){
     varcomps <- VarCorr(x)
     var.among <- as.numeric(varcomps[1, 1])
     var.within <- as.numeric(varcomps[2, 1])
     R <- var.among / (var.among + var.within)
     class(R) <- c("repeatability", "rep.lme")
   }
-  R
+  return(R)
 }
 
 print.repeatability <- function(x, ...){
-  if (inherits(x, "rep.aov")){
+  if (inherits(x, "rep.varcomps")){
     cat("Repeatability measured by random effects ANOVA.\n")
-    cat("\tRepeatability is", round(x, digits = 3))
+    cat("\tRepeatability is", round(x, digits = 3), "\n")
+  } 
+  if (inherits(x, "rep.aov")){
+    cat("Repeatability measured by one-way ANOVA.\n")
+    cat("\tRepeatability is", round(x, digits = 3), "\n")
   } 
   if (inherits(x, "rep.lme")){
     cat("Repeatability measured by linear mixed-effects model.\n")
-    cat("\tRepeatability is", round(x, digits = 3))
+    cat("\tRepeatability is", round(x, digits = 3), "\n")
   }
-}
+}
diff --git a/R/varcomps.R b/R/varcomps.R
@@ -19,16 +19,16 @@ varcomps <- function(fm, n){
                       "var.within" = var.within,
                       "var.among" = var.among)
   class(varcomp.obj) <- "varcomps"
-  varcomp.obj
+  return(varcomp.obj)
 }
 
-print.varcomps <- function(x, ...){
+print.varcomps <- function(x, digits = 3, ...){
   cat("Mean Squares\n")
   cat("\tGroups\t", x$MS.groups, "\n")
   cat("\tError\t", x$MS.error, "\n")
   cat("Variance Components\n")
   cat("\tWithin\t", x$var.within, "\n")
   cat("\tAmong\t", x$var.among, "\n")
-  cat("F = ", x$Fstat, "\n")
-  cat("p = ", x$p, "on", x$df1, "and", x$df2, "degrees of freedom.\n")
-}
+  cat("F = ", signif(x$Fstat, digits = digits), "\n")
+  cat("p = ", signif(x$p, digits = digits), "on", x$df1, "and", x$df2, "degrees of freedom.\n")
+}
diff --git a/abd-Ex.R b/abd-Ex.R
@@ -55,13 +55,7 @@ flush(stderr()); flush(stdout())
 
 ### ** Examples
 
-data(Aspirin)
-Aspirin
-Aspirin.expanded <- expand.dft(Aspirin, "count")
-xtabs(~ cancer + treatment, Aspirin.expanded)
-if (require(vcd)) {
-  mosaic(~cancer + treatment, Aspirin.expanded)
-}
+demo(sec9.2)
 
 
 
@@ -2485,10 +2479,7 @@ flush(stderr()); flush(stdout())
 
 ### ** Examples
 
-data(Trematodes)
-xtabs(~ infection.status + eaten, Trematodes)
-chisq.test( xtabs(~ infection.status + eaten, Trematodes) )
-summary(chisq.test( xtabs(~ infection.status + eaten, Trematodes) ) )
+demo(sec9.3)
 
 
 
@@ -2600,20 +2591,7 @@ flush(stderr()); flush(stdout())
 
 ### ** Examples
 
-data(VampireBites)
-VampireBites
-
-xtabs(count ~ estrous + bitten, data = VampireBites)
-fisher.test(xtabs(count ~ estrous + bitten, data = VampireBites))
-
-# With G-test
-# Source from http://www.psych.ualberta.ca/~phurd/cruft/
-try({
-  source("http://www.psych.ualberta.ca/~phurd/cruft/g.test.r");
-  g.test(xtabs(count ~ estrous + bitten, data = VampireBites));
-  g.test(xtabs(count ~ estrous + bitten, data = VampireBites))$expected
-  }
-)
+demo(sec9.4)
 
 
 
@@ -2959,6 +2937,26 @@ findData('Finch')                 # look for data sets with 'finch' in name
 
 
 
+cleanEx()
+nameEx("as.xtabs")
+### * as.xtabs
+
+flush(stderr()); flush(stdout())
+
+### Name: as.xtabs
+### Title: Convert objects to xtabs format
+### Aliases: as.xtabs as.xtabs.data.frame as.xtabs.matrix
+### Keywords: manip
+
+### ** Examples
+
+# example from example(fisher.test)
+df <- data.frame( X=c('Tea','Milk'), Tea=c(3,1), Milk=c(1,3) )
+xt <- as.xtabs(df, rowvar="Guess", colvar="Truth"); xt
+if (require(vcd)) { mosaic(xt) }
+
+
+
 cleanEx()
 nameEx("col.abd")
 ### * col.abd
@@ -3212,11 +3210,20 @@ flush(stderr()); flush(stdout())
 ### ** Examples
 
 data(WalkingStickFemurs)
-aovfit <- aov(femur.length ~ 1 + Error(specimen), data = WalkingStickFemurs)
-vc <- varcomps(aovfit, n = 2)
+# With aov() and Error()
+Error.fit <- aov(femur.length ~ 1 + Error(specimen), data = WalkingStickFemurs)
+vc <- varcomps(Error.fit, n = 2)
 vc
-R.varcomps <- repeatability(vc)
-R.varcomps
+repeatability(vc)
+
+# With aov()
+aov.fit <- aov(femur.length ~ specimen, data = WalkingStickFemurs)
+repeatability(aov.fit)
+
+# With lme()
+lme.fit <- lme(femur.length ~ 1, random = ~ 1 | specimen, 
+               data = WalkingStickFemurs)
+repeatability(lme.fit)
 
 
 

diff --git a/demo/00Index b/demo/00Index
@@ -13,4 +13,9 @@ sec6.4	Examples from Section 6.4
 sec6.5	Examples from Section 6.5
 sec7.1	Examples from Section 7.1
 sec7.2	Examples from Section 7.2
+sec9.2	Examples from Section 9.2
+sec9.3	Examples from Section 9.3
+sec9.4	Examples from Section 9.4
+sec9.5	Examples from Section 9.5
+sec17.8	Examples from Section 17.8
 examples Examples yet to be put into sections
diff --git a/demo/examples.R b/demo/examples.R
@@ -1,44 +1,9 @@
-
-### Aspirin
-
-Aspirin.expanded <- expand.dft(Aspirin, "count")
-str(Aspirin.expanded)
-
-# Plot 2X2 Contingency tables
-plot( ~ treatment + cancer, data = Aspirin.expanded)
-plot(table(Aspirin.expanded), main = "")
-
-# Calculate odds
-(Pr.asp <- 18496 / (18496 + 1438))
-(Odds.asp <- Pr.asp / (1 - Pr.asp))
-(Pr.no.asp <- 18515 / (18515 + 1427))
-(Odds.no.asp <- Pr.no.asp / (1 - Pr.no.asp))
-(Odds <- Odds.asp / Odds.no.asp)
-ln.Odds <- log(Odds)
-
-(SE.Odds <- sqrt(sum(1/Aspirin$count)))
-Z <- 1.96
-(CI.low <- ln.Odds - Z * SE.Odds)
-(CI.high <- ln.Odds + Z * SE.Odds)
-
-exp(CI.low)
-exp(CI.high)
-
-# Using odds.ratio()
-# First reformat the data so that "No cancer" is in column 1
-# and "Aspirin" is in row 2.
-x <- matrix(c(18515, 18496, 1427, 1438), nrow = 2)
-x
-odds.ratio(x)
-
 ### Blackbirds
-
-
 plot(log.after ~ log.before, data = Blackbirds,
-  xlim = c(3.9, 5.1), ylim = c(3.9, 5.1),
-  pch = 16, col = "red",
-  ylab = "log Antibody production after implant",
-  xlab = "log Antibody production before implant")
+     xlim = c(3.9, 5.1), ylim = c(3.9, 5.1),
+     pch = 16, col = "red",
+     ylab = "log Antibody production after implant",
+     xlab = "log Antibody production before implant")
 abline(b = 1, a = 0)
 
 hist(Blackbirds$diff.in.logs,
@@ -224,11 +189,12 @@ plot(count.sort, count.relfreq,
   ylab = "Cumulative relative frequency")
 
 \dontrun{
-p <- ggplot(data.frame(count.sort, count.relfreq), 
-  aes(count.sort, count.relfreq))
-p + geom_step(direction = "vh") +
-  scale_x_continuous("Species abundance") +
-  scale_y_continuous("Cumulative relative frequency")}
+  p <- ggplot(data.frame(count.sort, count.relfreq), 
+              aes(count.sort, count.relfreq))
+  p + geom_step(direction = "vh") +
+    scale_x_continuous("Species abundance") +
+    scale_y_continuous("Cumulative relative frequency")
+}
 
 ### FingerRatio
 plot(FingerRatio$CAGrepeats,
@@ -300,14 +266,14 @@ cv(GlidingSnakes$undulation.rate)
 
 ### GreatTitMalaria
 \dontrun{
-# Fig. 2.3-1
-require(ggplot2)
-bar <- ggplot(GreatTitMalaria, 
-  aes(x = Treatment, y = count, fill = Response))
-bar + geom_bar(stat = "identity", position = "dodge")
-
-# Fig. 2.3-2
-bar + geom_bar(stat = "identity", position = "fill")
+  ## Fig. 2.3-1
+  require(ggplot2)
+  bar <- ggplot(GreatTitMalaria, 
+                aes(x = Treatment, y = count, fill = Response))
+  bar + geom_bar(stat = "identity", position = "dodge")
+
+  ## Fig. 2.3-2
+  bar + geom_bar(stat = "identity", position = "fill")
 }
 
 ### Guppies
@@ -619,18 +585,19 @@ legend("bottomright", c("Confidence Bands", "Prediction Interval"),
 
 ### Lynx
 \dontrun{
-# Alternate form converting to Date class.
-Year <- as.Date(paste("01jan", Lynx$date, sep = ""),
-  "\%d\%b\%Y")
-Lynx <- cbind(Lynx, Year)
-
-require(ggplot2)
-p <- ggplot(Lynx, aes(Year, no.pelts))
-p + geom_line() + 
-  geom_point(color = "red") +
-  scale_y_continuous("Lynx fur returns") +
-  opts(panel.grid.minor = theme_blank()) +
-  opts(panel.grid.major = theme_line(size = 0.25, colour = "white"))}
+  ## Alternate form converting to Date class.
+  Year <- as.Date(paste("01jan", Lynx$date, sep = ""),
+                  "\%d\%b\%Y")
+  Lynx <- cbind(Lynx, Year)
+
+  require(ggplot2)
+  p <- ggplot(Lynx, aes(Year, no.pelts))
+  p + geom_line() + 
+    geom_point(color = "red") +
+    scale_y_continuous("Lynx fur returns") +
+    opts(panel.grid.minor = theme_blank()) +
+    opts(panel.grid.major = theme_line(size = 0.25, colour = "white"))
+}
 
 ### MarineReserve
 hist(MarineReserve$biomass.ratio)

diff --git a/demo/sec17.8.R b/demo/sec17.8.R
@@ -0,0 +1,5 @@
+
+# Figure 17.8-3 but using LOESS rather than splines
+xyplot(length ~ age, ShrinkingSeals, cex=.5, alpha=.8, col.line='black',
+	type=c('p','smooth'), span=.1, degree=2, eval=500
+	)
diff --git a/demo/sec2.5.R b/demo/sec2.5.R
@@ -4,7 +4,7 @@ xyplot(son.attract ~ father.ornament, data=Guppies,
 	ylab="Son's attractiveness")
 
 # Figure 2.5-2
-xyplot(pelts ~ year, LynxPopulationCycles, 
+xyplot(pelts ~ year, Lynx,
 	type='b',
 	ylab="Lynx fur returns",
 	xlab="Year")