###################################################
### chunk number 1: 
###################################################
options(SweaveHooks=list(fig=function(){par(cex.main=1.1,
                         mar=c(4.1,4.1,2.6,2.1), pty="s",
                         mgp=c(2.25,0.5,0), tck=-0.02)}))
graphics.off()
x11(width=8, height=8)


###################################################
### chunk number 2: place-files
###################################################
library(DAAGxtras)
dataFile(c("molclock1", "molclock2"))  # NB dataFile, not datafile


###################################################
### chunk number 3: paste
###################################################
paste("Leo", "the", "lion")
paste("a", "b")
paste("a", "b", sep="")
webpage <- "http://www.maths.anu.edu.au/~johnm/datasets/text/"
paste(webpage, "molclock.txt", sep="")


###################################################
### chunk number 4: webread
###################################################
webpage <- "http://www.maths.anu.edu.au/~johnm/datasets/text/"
molclock <- read.table(url(paste(webpage, "molclock.txt", sep="")))


###################################################
### chunk number 5: bostonc eval=FALSE
###################################################
## datafile("bostonc")


###################################################
### chunk number 6: Acmena
###################################################
library(DAAG)
Acmena <- subset(rainforest, species=="Acmena smithii")


###################################################
### chunk number 7: Acmena-2species
###################################################
AcSpecies <- subset(rainforest, species %in% c("Acacia mabellae", 
                                               "Acmena smithii"))


###################################################
### chunk number 8: table-NA
###################################################
with(rainforest, table(complete.cases(root), species))


###################################################
### chunk number 9: Acmena-scatterplot
###################################################
Acmena <- subset(rainforest, species=="Acmena smithii")
plot(wood ~ dbh, data=Acmena)
plot(wood ~ dbh, data=Acmena, log="xy")


###################################################
### chunk number 10: Acmena-lmline
###################################################
plot(wood~dbh, data=Acmena, log="xy")
## The lm() command will fit a line; more details later
## abline() then plots this line.
abline(lm(log10(wood) ~ log10(dbh), data=Acmena))
## Now print the coefficents, using a log10 scale
coef(lm(log10(wood) ~ log10(dbh), data=Acmena))
## Now check the coefficents, but fitting on a natural log scale
coef(lm(log(wood) ~ log(dbh), data=Acmena))


###################################################
### chunk number 11: ex2
###################################################
orings$Included <- logical(23)  # orings has 23 rows
orings$Included[c(1,2,4,11,13,18)] <- TRUE


###################################################
### chunk number 12: ex2-i
###################################################
plot(Total ~ Temperature, data=orings, pch=orings$included+1)
plot(Total ~ Temperature, data=orings, col=orings$Included+1)


###################################################
### chunk number 13: sapply-class eval=FALSE
###################################################
## sapply(ant111b, class)


###################################################
### chunk number 14: sapply-class eval=FALSE
###################################################
## unlist(sapply(ant111b, class))


###################################################
### chunk number 15: factor-levels
###################################################
par(mfrow=c(1,2), pty="s")
plot(weight ~ volume, pch=unclass(cover), data=allbacks)
plot(weight ~ volume, data=allbacks, type="n")
with(allbacks, text(weight ~ volume, labels=as.character(cover)))
par(mfrow=c(1,1))


###################################################
### chunk number 16: factors
###################################################
gender <- factor(c(rep("female", 91), rep("male", 92)))
table(gender)
gender <- factor(gender, levels=c("male", "female"))
table(gender)
gender <- factor(gender, levels=c("Male", "female")) # Note the mistake
                            # The level was "male", not "Male"
table(gender)
rm(gender)                  # Remove gender


###################################################
### chunk number 17: ex5
###################################################
gender <- factor(c(rep("female", 91), rep("male", 92)))
table(gender)
gender <- factor(gender, levels=c("male", "female"))


###################################################
### chunk number 18: 
###################################################
table(gender)
gender <- factor(gender, levels=c("Male", "female")) # Note the mistake
                            # The level was "male", not "Male"
table(gender)
rm(gender)                  # Remove gender


###################################################
### chunk number 19: subset
###################################################
Acmena <- subset(rainforest, species=="Acmena smithii")
ord <- order(Acmena$dbh)
acm <- Acmena[ord, ]


###################################################
### chunk number 1: prelim2
###################################################
library(DAAG)


###################################################
### chunk number 2: loop-states
###################################################
oldpar <- par(mfrow=c(2,4))
for (i in 2:9){
plot(austpop[,1], log(austpop[, i]), xlab="Year", ylab=names(austpop)[i], 
     pch=16, ylim=c(0,10))}
par(oldpar)


###################################################
### chunk number 3: rnorm500
###################################################
y <- rnorm(500)


###################################################
### chunk number 4: 200x4
###################################################
av <- numeric(500)
for(i in 1:500){
  av[i] <- mean(rnorm(4))
}


###################################################
### chunk number 5: 500x4-matrices
###################################################
mat <- matrix(rnorm(500*4), nrow=500)
av <- apply(mat, 2, mean)


###################################################
### chunk number 6: meanANDsd
###################################################
meanANDsd <- function(x){
    av <- mean(x)
    sdev <- sd(x)
    c(mean=av, sd = sdev) # The function returns this vector
}


###################################################
### chunk number 7: Acmena-distn
###################################################
library(DAAG)      # The data frame rainforest is from DAAG
Acmena <- subset(rainforest, species=="Acmena smithii")
hist(Acmena$dbh)
hist(Acmena$dbh, prob=TRUE)   # Use a density scale


###################################################
### chunk number 8: overlay
###################################################
hist(Acmena$dbh, prob=TRUE, xlim=c(0,50))   # Use a density scale
lines(density(Acmena$dbh, from=0))


###################################################
### chunk number 9: countNA
###################################################
library(MASS)
## Create a function that counts NAs
count.na <- function(x)sum(is.na(x))
## Check that the function works as expected
count.na(c(1, 5, NA, 5, NA,8))
## Apply the function to each column of the Pima.tr2
sapply(Pima.tr2, count.na)


###################################################
### chunk number 10: anymiss
###################################################
missIND <- complete.cases(Pima.tr2)
Pima.tr2$anymiss <- c("some","none")[missIND+1]


###################################################
### chunk number 11: Pima-plus
###################################################
stripplot(anymiss ~ npreg + bmi | type, data=Pima.tr2, outer=TRUE,
          scales=list(relation="free"), xlab="Measure")


###################################################
### chunk number 12: densityplot
###################################################
library(lattice)
densityplot( ~ npreg + bmi | type, groups=anymiss, data=Pima.tr2,
            auto.key=list(columns=2))


###################################################
### chunk number 13: qq-Pima
###################################################
qq(anymiss ~ npreg | type, data=Pima.tr2)


###################################################
### chunk number 14: miss-vs-nomiss
###################################################
densityplot(~ bmi, groups=anymiss, data=Pima.tr2, auto.key=list(columns=2))


###################################################
### chunk number 15: boot-miss
###################################################
complete <- complete.cases(Pima.tr2)
rownum <- seq(along=complete)  # generate row numbers
allpresSample <- sample(rownum[complete], replace=TRUE)
  # By default, sample size is #values of the first argument
NApresSample <- sample(rownum[!complete], replace=TRUE)
densityplot(~bmi, groups=completeF, data=Pima.tr2, 
            auto.key=list(columns=2), subset=c(allpresSample, 
            NApresSample))


###################################################
### chunk number 16: boot-diff-of-means
###################################################
twot <- function(n=99){
  complete <- complete.cases(Pima.tr2)
  rownum <- seq(along=complete)  # generate row numbers
  d2 <- numeric(n+1)
  d2[1] <- with(Pima.tr2, mean(bmi[complete], na.rm=TRUE)-
                                 mean(bmi[!complete], na.rm=TRUE))
  for(i in 1:n){
    allpresSample <- sample(rownum[complete], replace=TRUE)
    NApresSample <- sample(rownum[!complete], replace=TRUE)
    d2[i+1] <- with(Pima.tr2, mean(bmi[allpresSample], na.rm=TRUE)-
                                   mean(bmi[NApresSample], na.rm=TRUE))
}
  d2
}
  # NB: There are n+1 values in all; the mean difference for the
  # actual sample values, plus differences for n bootstrap samples
## Now estimate and plot the density function
d2 <- twot(n=999)
dens <- density(d2)
plot(dens)
## Check the proportion of values of d2 that are less than or equal to 0
sum(d2<0)/length(d2)


###################################################
### chunk number 1: prelim2
###################################################
library(DAAG)


###################################################
### chunk number 2: strip-plot
###################################################
## First, use regular graphics
with(ant111b, stripchart(harvwt ~ site))
## Next, use lattice graphics
library(lattice)
stripplot(site ~ harvwt, data=ant111b) 


###################################################
### chunk number 3: groupDis
###################################################
library(DAAGxtras)
groupDis <- as.character(rockArt$District)
tab <- table(rockArt$District)
le4 <- rockArt$District %in% names(tab)[tab <= 4]
groupDis[le4] <- "other"
groupDis <- factor(groupDis)


###################################################
### chunk number 4: smooths
###################################################
Acmena <- subset(rainforest, species=="Acmena smithii")
## Use lowess()
plot(wood ~ dbh, data=Acmena)
ord <- order(Acmena$dbh)
with(Acmena[ord, ], lines(predict(loess(wood ~ dbh)) ~ dbh))
## Now use panel.smooth()
plot(wood ~ dbh, data=Acmena)
with(Acmena, panel.smooth(dbh, wood))
  # Observe that panel.smooth() does not at this time allow either 
  # a graphics formula as argument, or the data parameter.


###################################################
### chunk number 5: sapply-code
###################################################
oldpar <- par(mfrow=c(2,4))
invisible(
sapply(2:9, function(i, df)
       plot(df[,1], log(df[, i]), 
            xlab="Year", ylab=names(df)[i], pch=16, ylim=c(0,10)),
            df=austpop)
)
par(oldpar)


###################################################
### chunk number 6: create-mat
###################################################
xxMAT <- matrix(runif(480000), ncol=50)
xxDF <- as.data.frame(xxMAT)


###################################################
### chunk number 7: add-one
###################################################
system.time(invisible(xxMAT+1))[1:3]
system.time(invisible(xxDF+1))[1:3]


###################################################
### chunk number 8: apply-funs eval=FALSE
###################################################
## ## Use apply() [matrix argument], or sapply() [data frame argument]
## system.time(av1 <- apply(xxMAT, 2, mean))[1:3]
## system.time(av1 <- sapply(xxDF, mean))[1:3]
## ## Use a loop that does the calculation for each column separately
## system.time({av2 <- numeric(50);
##             for(i in 1:50)av[i] <- mean(xxMAT[,i])
##             })[1:3]
## system.time({av2 <- numeric(50);
##             for(i in 1:50)av[i] <- mean(xxDF[,i])
##             })[1:3]
## ## Matrix multiplication
## system.time({colOFones <- rep(1, dim(xxMAT)[2])
##              av3 <- xxMAT %*% colOFones / dim(xxMAT)[2]
##             })[1:3]


###################################################
### chunk number 9: calcCT
###################################################
"calcCT" <-
  function(df=fumig, times=c(5,10,30,60,90,120), ctcols=3:8){
    usualfac <- c(7.5,12.5,25,30,30,15)
    modfac <- c(20,25,30,30,15)
    modtimes <- c(15,30,60,120)
    require(splines)
    m <- dim(df)[1]
    x1 <- times[-1]
    conc15 <- numeric(m)
    usualct <- numeric(m)
    modct <- numeric(m)
    for(i in 1:m){
      y <- unlist(df[i, ctcols])
      y1 <- y[-1]
      ct.lm <- lm(y1 ~  ns(x1,4))
      xy = data.frame(x1=c(15,30,60,120))
      hat <- predict(ct.lm, newdata=xy)
      conc15[i] <- hat[1]
      usualct[i] <- sum(usualfac*y)/60
      modct[i] <- sum(modfac*y1)/60
    }
    df <- cbind(usualct=usualct, modct=modct, df[,-ctcols],
                estconc15=conc15)
    df
  }


###################################################
### chunk number 10: possum1
###################################################
possum1 <- rnorm(10)
object.size(possum1)


###################################################
### chunk number 11: 
###################################################
nam <- ls(pattern="^possum")


###################################################
### chunk number 12: sizes eval=FALSE
###################################################
## sapply(nam, function(x)object.size(get(x)))


###################################################
### chunk number 13: betterls
###################################################
workls <- function()ls(name=".GlobalEnv")
workls()


###################################################
### chunk number 14: anymiss
###################################################
anymiss <- complete.cases(Pima.tr2)
Pima.nomiss <- Pima.tr2[!anymiss, ]
Pima.miss <- Pima.tr2[anymiss, ]


###################################################
### chunk number 15: densityplot
###################################################
library(lattice)
densityplot( ~ npreg, groups=anymiss, data=Pima.tr2)


###################################################
### chunk number 16: binops-arith
###################################################
"+"(2,5)
"-"(10,3)
"/"(2,5)
"*"("+"(5,2), "-"(3,7))


###################################################
### chunk number 17: binops-arith2
###################################################
(0:25) %/% 5
(0:25) %% 5


###################################################
### chunk number 18: subset notation
###################################################
z <- c(2, 6, -3, NA, 14, 19)
"["(z, 5)
heights <- c(Andreas=178, John=185, Jeff=183) 
"["(heights, c("Jeff", "John"))


###################################################
### chunk number 1: 
###################################################
options(SweaveHooks=list(fig=function(){par(cex.main=1.1,
                         mar=c(4.1,4.1,2.6,2.1), pty="s",
                         mgp=c(2.25,0.5,0), tck=-0.02)}))


###################################################
### chunk number 2: density-cdf
###################################################
curve(dnorm(x, mean=2.5, sd=1.5), from = 2.5-3*1.5, to = 2.5+3*1.5)
curve(pnorm(x, mean=2.5, sd=1.5), from = 2.5-3*1.5, to = 2.5+3*1.5)


###################################################
### chunk number 3: area
###################################################
pnorm(4, mean=2.5, sd=1.5) - pnorm(0.5, mean=2.5, sd=1.5)


###################################################
### chunk number 4: rnorm-hist
###################################################
par(mfrow=c(1,3), pty="s")
x <- rnorm(50)
hist(x, probability=TRUE)
lines(density(x))
xval <- pretty(c(-3,3), 50)
lines(xval, dnorm(xval), col="red")


###################################################
### chunk number 5: hist-body eval=FALSE
###################################################
## library(MASS)
## plot(density(Animals$body))
## logbody <- log(Animals$body)
## plot(density(logbody))
## av <- mean(logbody)
## sdev <- sd(logbody)
## xval <- pretty(c(av-3*sdev, av+3*sdev), 50)
## lines(xval, dnorm(xval, mean=av, sd=sdev))


###################################################
### chunk number 6: rnorm-density
###################################################
plot(density(rnorm(50)), type="l")


###################################################
### chunk number 7: ex-bw
###################################################
plot(density(rnorm(500)), type="l")


###################################################
### chunk number 8: ex-bw
###################################################
plot(density(rnorm(50), bw=0.2), type="l")
plot(density(rnorm(50), bw=0.6), type="l")


###################################################
### chunk number 9: sample-with-replacement eval=FALSE
###################################################
## sample(1:5, replace=TRUE)   # With replacement sample, size 5
## for(i in 1:10)print(sample(1:5, replace=TRUE))  # 10 such samples
## plot(density(log10(Animals$body)))
## lines(density(sample(log10(Animals$body), replace=TRUE)), col="red")


###################################################
### chunk number 10: rnorm-qq
###################################################
qqnorm(rnorm(10))
qqnorm(rnorm(50))
qqnorm(rnorm(200))


###################################################
### chunk number 11: qref
###################################################
library(DAAG)
qreference(m=10)
qreference(m=50)
qreference(m=200)


###################################################
### chunk number 12: qreference-ex
###################################################
qreference(possum$totlngth)


###################################################
### chunk number 1: 
###################################################
options(SweaveHooks=list(fig=function(){par(cex.main=1.1,
                         mar=c(4.1,4.1,2.6,2.1),
                         mgp=c(2.25,0.5,0), tck=-0.02)}))


###################################################
### chunk number 2: cuckoos
###################################################
par(mfrow=c(1,2))
library(DAAG)
library(lattice) 
table(cuckoos$species)   # Check numbers of each species
names(cuckoos)[1] <- "length"
stripplot(species ~ length, data=cuckoos)
## Look also at the relationship between length and breadth;
## is it the same for all species?
stripplot(breadth ~ length, groups=species, 
          data=cuckoos, auto.key=list(columns=3))
par(mfrow=c(1,1))


###################################################
### chunk number 3: cuckoo-lm
###################################################
lm(length ~ -1 + species, data=cuckoos) 


###################################################
### chunk number 4: cuckoo-lm2
###################################################
lm(length ~ species, data=cuckoos)


###################################################
### chunk number 5: cuckoos-sd
###################################################
with(cuckoos, sapply(split(length, species), sd))


###################################################
### chunk number 6: cuckoos-se
###################################################
sdev <- with(cuckoos, sapply(split(length, species), sd))
n <- with(cuckoos, sapply(split(length, species), length))
sdev/sqrt(n)


###################################################
### chunk number 7: cuckoos-sefun
###################################################
se <- function(x)sd(x)/sqrt(length(x))
with(cuckoos, sapply(split(length, species), se))


###################################################
### chunk number 8: cuckoos-anon
###################################################
with(cuckoos, 
sapply(split(length, species), function(x)sd(x)/sqrt(length(x)))
)


###################################################
### chunk number 9: reg-spline
###################################################
library(splines)
plot(ohms ~ juice, data=fruitohms)
hat <- with(fruitohms, fitted(lm(ohms ~ ns(juice, 4))))
with(fruitohms, lines(juice,hat, col=3))
## Check the locations of the internal knots.
attributes(ns(fruitohms$juice, 4))$knots 
##
## Now compare with lowess
with(fruitohms, lines(lowess(juice, ohms), col="red"))
with(fruitohms, lines(lowess(juice, ohms, f=0.2), col="cyan"))


###################################################
### chunk number 10: hills2k
###################################################
lhills2k <- log(subset(races2000[, 7:9], subset=races2000$type=="hill"))
names(lhills2k) <- c("ldist",  "lclimb", "ltime")
lhills2k.lm <- lm(ltime ~ ldist + lclimb, data = lhills2k)


###################################################
### chunk number 11: spline-ns34
###################################################
lhills2k.ns <- lm(ltime ~ ns(ldist,4) + ns(lclimb,3), data = lhills2k)


###################################################
### chunk number 12: rlm-fit
###################################################
library(MASS)
lhills2k.rlm <- rlm(ltime ~ ns(ldist,4) + ns(lclimb,3), data = lhills2k)
par(mfrow=c(2,2))
plot(lhills2k.rlm)
par(mfrow=c(1,1))