################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Appendix") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=65) ## palette(gray((0:8)/8)) #library(lattice) trellis.par.set(canonical.theme(color=FALSE)) library(primer) ################################################### ### chunk number 2: eval=FALSE ################################################### ## help.start() ################################################### ### chunk number 3: ################################################### # Here I Get my Working Directory; that is, # I find out which folder R is currently operating from. getwd() ################################################### ### chunk number 4: eval=FALSE ################################################### ## setwd("~/Documents/Rwork") ################################################### ### chunk number 5: eval=FALSE ################################################### ## setwd("C:/Documents and Settings/Users/Jane/My Documents/Rwork") ## ################################################### ### chunk number 6: ################################################### # This will calculate the mean of 10 random standard normal variables. mean( rnorm( 10 ) ) ################################################### ### chunk number 7: ################################################### ?mean ################################################### ### chunk number 8: help ################################################### help(mean) ################################################### ### chunk number 9: eval=FALSE ################################################### ## help.search("mean") ## apropos("mean") ################################################### ### chunk number 10: eval=FALSE ################################################### ## RSiteSearch("violin") ## RSiteSearch("violin",restrict = c("functions")) ################################################### ### chunk number 11: eval=FALSE ################################################### ## help(RSiteSearch) ################################################### ### chunk number 12: ass ################################################### a <- 2+3 a ################################################### ### chunk number 13: a2 ################################################### b <- a+a ################################################### ### chunk number 14: semi ################################################### a+a; a+b ################################################### ### chunk number 15: ################################################### Y <- c(8.3,8.6,10.7,10.8,11,11,11.1,11.2,11.3,11.4) ################################################### ### chunk number 16: ################################################### Y = c(8.3,8.6,10.7,10.8,11,11,11.1,11.2,11.3,11.4) ################################################### ### chunk number 17: ################################################### 1:4 4:1 -1:3 -(1:3) ################################################### ### chunk number 18: ################################################### seq(from=1, to=3, by = .2) seq(1, 3, by=.2) seq(1, 3, length=7) ################################################### ### chunk number 19: ################################################### rep(1, 3) rep(1:3, 2) rep(1:3, each=2) ################################################### ### chunk number 20: ################################################### sum(Y) mean(Y) max(Y) length(Y) summary(Y) ################################################### ### chunk number 21: ################################################### Names <- c("Sarah", "Yunluan") Names b <- c(TRUE, FALSE) b ################################################### ### chunk number 22: ################################################### class(Y) class(b) ################################################### ### chunk number 23: ################################################### Y > 10 Y > mean(Y) ################################################### ### chunk number 24: ################################################### Y == 11 ################################################### ### chunk number 25: ################################################### Y != 11 ################################################### ### chunk number 26: ################################################### a <- 1:3 b <- 4:6 a + b ################################################### ### chunk number 27: ################################################### a*b a/b ################################################### ### chunk number 28: ################################################### a + 1 a * 2 1/a ################################################### ### chunk number 29: ################################################### a*1:2 ################################################### ### chunk number 30: ################################################### a*c(1,2,1) ################################################### ### chunk number 31: ################################################### 1:4 * 1:2 ################################################### ### chunk number 32: ################################################### 1:4 * c(1,2,1,2) ################################################### ### chunk number 33: ################################################### Y[1] Y[1:3] ################################################### ### chunk number 34: ################################################### Y > mean(Y) Y[ Y > mean(Y) ] ################################################### ### chunk number 35: ################################################### a <- c(5,3,6, NA) a is.na(a) !is.na(a) a[ !is.na(a) ] na.exclude(a) ################################################### ### chunk number 36: ################################################### mean(a) mean(a, na.rm=TRUE) d <- na.exclude(a) mean(d) ################################################### ### chunk number 37: ################################################### matrix(letters[1:4], ncol=2) ################################################### ### chunk number 38: ################################################### M <- matrix( 1:4, nrow=2 ) M ################################################### ### chunk number 39: ################################################### M2 <- matrix( 1:4, nrow=2, byrow=TRUE ) M2 ################################################### ### chunk number 40: ################################################### I <- diag(1, nrow=2) I ################################################### ### chunk number 41: solve ################################################### Minv <- solve(M) M %*% Minv ################################################### ### chunk number 42: ################################################### M[1,2] M[1, 1:2] ################################################### ### chunk number 43: ################################################### M[,2] M[,] ################################################### ### chunk number 44: ################################################### N <- matrix(0:3, nrow=2) N ################################################### ### chunk number 45: ################################################### M * N ################################################### ### chunk number 46: ################################################### M%*%N ################################################### ### chunk number 47: commutative ################################################### N%*%M ################################################### ### chunk number 48: ################################################### 1:2%*%M M%*%1:2 ################################################### ### chunk number 49: ################################################### V <- matrix(1:2, ncol=1) ################################################### ### chunk number 50: ################################################### M %*% V try(V %*% M) ################################################### ### chunk number 51: ################################################### M + N M + 2 ################################################### ### chunk number 52: ################################################### t(M) ################################################### ### chunk number 53: ################################################### dat <- data.frame(species=c("S.altissima", "S.rugosa", "E.graminifolia", "A. pilosus"), treatment=factor( c("Control", "Water", "Control", "Water") ), height = c(1.1, 0.8, 0.9, 1.0), width =c(1.0, 1.7, 0.6, 0.2) ) dat ################################################### ### chunk number 54: ################################################### dat[2,] dat[3,4] ################################################### ### chunk number 55: ################################################### dat[,2]=="Water" dat[ dat[,2]=="Water", ] ################################################### ### chunk number 56: ################################################### subset(dat, treatment == "Water") ################################################### ### chunk number 57: ################################################### c("Control", "Medium", "High") rep( c("Control", "Medium", "High"), each=3 ) Treatment <- factor( rep( c("Control", "Medium", "High"), each=3 ) ) Treatment ################################################### ### chunk number 58: ################################################### quartz(,4,4) ################################################### ### chunk number 59: ################################################### levels(Treatment) stripchart(1:9 ~ Treatment) ################################################### ### chunk number 60: ################################################### dev.print(pdf, "strip.pdf") ################################################### ### chunk number 61: ################################################### quartz(,4,4) ################################################### ### chunk number 62: ################################################### Treatment <- factor( rep( c("Control", "Medium", "High"), each=3 ), levels=c("Control", "Medium", "High") ) levels(Treatment) stripchart(1:9 ~ Treatment) ################################################### ### chunk number 63: ################################################### dev.print(pdf, "strip2.pdf") ################################################### ### chunk number 64: ################################################### my.list <- list(My.Y = Y, b = b, Names, Weed.data=dat, My.matrix = M2, my.no = 4 ) my.list ################################################### ### chunk number 65: ################################################### my.list[["b"]] my.list[[2]] ################################################### ### chunk number 66: ################################################### my.list[["b"]] my.list$b ################################################### ### chunk number 67: ################################################### my.list[1:2] ################################################### ### chunk number 68: ################################################### my.list[["b"]][1] ################################################### ### chunk number 69: ################################################### mean( dat$height ) ################################################### ### chunk number 70: eval=FALSE ################################################### ## help(mean) ################################################### ### chunk number 71: ################################################### mean(1:4) mean(1:4, trim=0) ################################################### ### chunk number 72: ################################################### class(1:10) class(warpbreaks) summary( 1:10) summary(warpbreaks) ################################################### ### chunk number 73: ################################################### summary( lm(breaks ~ wool, data=warpbreaks) ) ################################################### ### chunk number 74: ################################################### MyBogusMean <- function(x, cheat=0.05) { SumOfX <- sum(x) n <- length(x) trueMean <- SumOfX/n (1 + cheat) * trueMean} RealSales <- c(100, 200, 300) MyBogusMean(RealSales) ################################################### ### chunk number 75: ################################################### MyBogusMean(RealSales, cheat=0.1) MyBogusMean(RealSales, cheat=0) ################################################### ### chunk number 76: ################################################### e <- c(5, 4, 2, 1, 3) e sort(e) sort(e, decreasing=TRUE) ################################################### ### chunk number 77: ################################################### e order(e) e[ order(e) ] ################################################### ### chunk number 78: ################################################### dat order.nos <- order(dat$height) order.nos ################################################### ### chunk number 79: ################################################### dat[order.nos, ] ################################################### ### chunk number 80: ################################################### dat[ rev(order.nos), ] ################################################### ### chunk number 81: ################################################### m <- matrix( 1:10, nrow=2) m apply(m, MARGIN=1, mean) apply(m, MARGIN=2, sum) ################################################### ### chunk number 82: ################################################### sapply(1:10, function(i) mean( rnorm(5) ) ) ################################################### ### chunk number 83: ################################################### gens <- 10 output <- numeric(gens + 1) output[1] <- 25 for( t in 1:gens ) output[t + 1] <- output[t] + round( rnorm(n=1, mean=0, sd=2), 0) output ################################################### ### chunk number 84: ################################################### summary(CO2) CO2.wide <- reshape(CO2, v.names="uptake", idvar="Plant", timevar="conc", direction="wide") names(CO2.wide) ################################################### ### chunk number 85: ################################################### CO2.long <- reshape(CO2.wide, v.names="Uptake", varying=list(4:10), timevar="Concentration", times=c(95, 175, 250, 350, 500, 675,1000) ) head(CO2.long) ################################################### ### chunk number 86: ################################################### CO2.long2 <- with( CO2.long, CO2.long[order(Plant, Concentration),]) head(CO2.long2) ################################################### ### chunk number 87: ################################################### tapply(CO2[["uptake"]], list(CO2[["Treatment"]]), mean) ################################################### ### chunk number 88: ################################################### tapply(CO2[["uptake"]], list(CO2[["Treatment"]], CO2[["Type"]]), sd ) ################################################### ### chunk number 89: ################################################### tapply(CO2[["uptake"]], list(CO2[["Treatment"]], CO2[["Type"]]), function(x) c(mean(x), sd(x)) ) ################################################### ### chunk number 90: ################################################### aggregate(CO2[,4:5], list(Plant = CO2[["Plant"]]), mean) ################################################### ### chunk number 91: ################################################### dat <- data.frame(Name=rep(c("Control","Treatment"), each=5), First=runif(10), Second=rnorm(1)) write.table(dat, file="dat.txt") write.csv(dat, file="dat.csv") ################################################### ### chunk number 92: ################################################### dat.new <- read.csv("dat.csv") dat.new2 <- read.table("dat.txt", header=TRUE) ################################################### ### chunk number 93: ################################################### mod.out <- summary( aov(First ~ Name, data=dat)) mod.out[[1]] write.csv(mod.out[[1]], "ModelANOVA.csv") ################################################### ### chunk number 94: ################################################### qnorm(p=c(0.025,0.975)) ################################################### ### chunk number 95: myhist ################################################### myplot <- hist(rnorm(20, m=11, sd=6), probability=TRUE) myplot lines(myplot$mids, dnorm(myplot$mids, m=11, sd=6) ) ################################################### ### chunk number 96: ################################################### logGrowth <- function(t, y, p){ N <- y[1] with(as.list(p), { dN.dt <- r * N * (1 - a * N) return( list( dN.dt ) ) } ) } ################################################### ### chunk number 97: ################################################### logGrowth <- function(t, y, p){ dN.dt <- p[1] * y[1] * (1 - p[2] * y[1]) return( list( dN.dt ) ) } ################################################### ### chunk number 98: ################################################### p <- c(r=1, a = 0.001) y0 <- c(N=10) t <- 1:20 ################################################### ### chunk number 99: ################################################### library(deSolve) out <- ode(y=y0, times=t, func=logGrowth, parms=p) out[1:5,] ################################################### ### chunk number 100: ################################################### LVComp <- function(t, y, p){ N <- y with(as.list(p), { dN1.dt <- r[1] * N[1] * (1 - a[1,1]*N[1] - a[1,2]*N[2]) dN2.dt <- r[2] * N[2] * (1 - a[2,1]*N[1] - a[2,2]*N[2]) return( list( c(dN1.dt, dN2.dt) ) ) } ) } ################################################### ### chunk number 101: ################################################### a <- matrix(c(0.02, 0.01, 0.01, 0.03), nrow=2) r <- c(1,1) p2 <- list(r, a) N0 <- c(10,10) t2 <- c(1,5,10,20) out <- ode(y=N0, times=t2, func=LVComp, parms=p2) out[1:4,] ################################################### ### chunk number 102: ################################################### EV <- function(t, y, p){ with(as.list(p), { dv.dt <- b * y[1]*(1-.005*y[1]) - a*y[1]*y[2] de.dt <-a*e*y[1]*y[2] - s*y[2] return( list( c(dv.dt, de.dt) ) ) } ) } ################################################### ### chunk number 103: ################################################### rootfun <- function (t,y,p) { dstate <- unlist(EV(t,y,p)) #rate of change vector return(sum(abs(dstate))-1e-10) } ################################################### ### chunk number 104: ################################################### p <- c(b = 0.5, a = 0.02, e=0.1, s=0.2) t <- c(0,1e10) ################################################### ### chunk number 105: ################################################### out <- ode(y=c(45,200), t, EV, parms=p, rootfun=rootfun, method="lsodar") out[,] ################################################### ### chunk number 106: ################################################### y <- c(1, 0:10) f <- function(y, mu) { sum( (y-mu)^2 ) } ################################################### ### chunk number 107: LS ################################################### guesses <- seq(4, 6, by=.05) LS.criterion <- sapply(guesses, function(mu) f(mu=mu, y=y)) plot(guesses, LS.criterion, type='l') ################################################### ### chunk number 108: optimize ################################################### (results <- optimize(f, c(0,10), y=y)) ################################################### ### chunk number 109: ################################################### LL <- function(mu, SD){ -sum( dnorm(y, mean=mu, sd=SD, log=TRUE) ) } ################################################### ### chunk number 110: mle2 ################################################### library(bbmle) (fit <- mle2(LL, start=list(mu=5, SD=1), control=list(maxit=10^5) ) ) ################################################### ### chunk number 111: ################################################### mle2(y ~ dnorm(mu, sd=SD), start=list(mu=1, SD=2)) ################################################### ### chunk number 112: ################################################### summary(fit) pr <- profile(fit) ################################################### ### chunk number 113: profile ################################################### par(mar=c(5,4,3,2)) plot(pr) ################################################### ### chunk number 114: ################################################### host1 <- expression(R*H*(1+a*P)^-k) D(host1, "H") ################################################### ### chunk number 115: Myplot ################################################### data(trees); attach(trees) plot(Girth, Height) ################################################### ### chunk number 116: eval=FALSE ################################################### ## par(mar=c(5,4,3,2)) ## plot(Girth, Volume, type='n', main = "My Trees") ## points(Girth, Volume, type='h' , col="lightgrey", pch=19) ################################################### ### chunk number 117: eval=FALSE ################################################### ## hts <- (Height-min(Height))/max(Height-min(Height)) ## my.colors <- hcl(h=30 + 270 * hts, alpha=0.90) ## text(Girth, Volume, Height, col=my.colors, cex = 0.5 + hts ) ################################################### ### chunk number 118: MyTrees1 ################################################### par(mar=c(5,4,3,2)) plot(Girth, Volume, type='n', main = "My Trees") points(Girth, Volume, type='h' , col="lightgrey", pch=19) hts <- (Height-min(Height))/max(Height-min(Height)) my.colors <- hcl(h=30 + 270 * hts, l=50, alpha=0.90) text(Girth, Volume, Height, col=my.colors, cex = 0.5 + hts ) ################################################### ### chunk number 119: Mymatplot ################################################### trees.sort <- trees[order(trees$Girth, trees$Height),] matplot(trees.sort$Girth, trees.sort[,2:3], type='b') text(18, 40, "Volume", col="darkred") text(10, 58, "Height") ################################################### ### chunk number 120: ################################################### quartz(,4,4) par(mar=c(5,4,2,4)) plot(Girth, Volume, main = "My Trees") ################################################### ### chunk number 121: ################################################### par(new=TRUE) plot(Girth, Height, axes=FALSE, bty="n", xlab="", ylab="", pch=3) ################################################### ### chunk number 122: ################################################### axis(4) mtext("Height", side=4, line=3) ################################################### ### chunk number 123: MyTrees2Y ################################################### par(mar=c(5,4,2,4)) plot(Girth, Volume, main = "My Trees") par(new=TRUE) plot(Girth, Height, axes=FALSE, bty="n", xlab="", ylab="", pch=3) axis(4) mtext("Height", side=4, line=3) ################################################### ### chunk number 124: ################################################### quartz(width=5, height=3) ################################################### ### chunk number 125: eval=FALSE ################################################### ## windows(width=5, height=3) ################################################### ### chunk number 126: ################################################### quartz(,5,5) layout( matrix( c(1,2,3,3), nrow=2, byrow=TRUE) ) plot(Girth, Height) ################################################### ### chunk number 127: ################################################### par(mar=c(3,3,1,1), mgp=c(1.6, .2, 0), tcl=.2) plot(Girth, Height) par(mar=c(3,3,2,1), mgp=c(1.6, .2, 0), tcl=.2) plot(Girth, Height, axes=FALSE, xlim=c(8,22) ) axis(1, tcl=-.3) axis(2, tick=F) rug(Height, side=2, col=2) title("A Third, Very Wide, Plot") ################################################### ### chunk number 128: ################################################### dev.print(pdf, 'MyTrees4.pdf') ################################################### ### chunk number 129: ################################################### getwd() quartz(,4,4) plot(Height, Volume, main="Tree Data") dev.print(pdf, "MyTree.pdf") ################################################### ### chunk number 130: ################################################### summary(Girth) stem(Girth) ################################################### ### chunk number 131: ################################################### quartz(,6,7) ################################################### ### chunk number 132: ################################################### layout( matrix( c(1,2,2,3,4,4), nrow=2, byrow=TRUE ) ) plot(1:length(Girth), Girth, xlab="Order of Sample Collection?") hist(Girth, prob=TRUE) rug(Girth) lines(density(Girth)) boxplot(Girth, main="Boxplot of Girth") points(jitter( rep(1,length(Girth)) ), Girth) qqnorm(log(Girth)) qqline(log(Girth)) title(sub="Log transformed data") ################################################### ### chunk number 133: ################################################### dev.print(pdf, "Distributions.pdf") ################################################### ### chunk number 134: ################################################### A <- matrix(c(0, .1, 10, .5), nrow=2) eig.A <- eigen(A) str(eig.A)