## The latticist and playwith packages
## The latticist package
library(DAAGxtras)
library(latticist)
latticist(nihills)
## The playwith package
library(playwith)
playwith(xyplot(time ~ dist, data=nihills))

## Regression with 2 explanatory variables - the nihills data
splom(nihills)  # splom(~ nihills) has the same effect
lognihills <- log(nihills) 
names(lognihills) <- c("ldist", "lclimb", "ltime", "ltimef") 
splom(lognihills)
nihills.lm <- lm(ltime ~ lclimb + ldist, data=lognihills)
summary(nihills.lm)
# May need x11() to open graphics window
par(mfrow=c(2,2))
plot(nihills.lm)  # By default, which=c(1:3,5)
par(mfrow=c(1,1))
plot(nihills.lm, which=5)  
plot(nihills.lm, which=4)
plot(nihills.lm, which=6)  # Takes a while to understand
## In search of moremeaningful coefficients
lognihills$logGrad <- with(nihills, log(climb/dist))
nihillsG.lm <- lm(logtime ~ logdist + logGrad, data=lognihills)
nihillsG.lm <- lm(ltime ~ ldist + logGrad, data=lognihills)
summary(nihillsG.lm)
par(mfrow=c(1,2))
termplot(nihillsG.lm, partial=TRUE, smooth=panel.smooth)
par(mfrow=c(1,1))

## Save all or part of the workspace to an "image" file
save.image()
args(save)
save(list=ls(), file="f.RData")

## Attach an image file.  NB: .RData is now preferred to .rda
attach("d:/rda/housetemps.rda")
ls()
search()
ls(pos=6)   # Assumes file was attached in position 6
plotT()     # Run function from attached image file
detach("file:d:/rda/housetemps.rda")
  # detach(6) is an alternative if image file was in position 6

## Time series
LakeHuron
## help(LakeHuron)
plot(LakeHuron)
## Create multivariate time series
library(DAAG)
jobs
jobts <- ts(jobs[,1:6], start=1995, frequency=12)
help(jobs)
plot(jobs)
x11()    # Open another graphics window
plot(jobts)

library(DAAGxtras)
head(bomregions)
bomdts <- ts(bomregions[, -1], start=1900, frequency=1)
args(window)
## help(window)
window(bomdts, start=1900, end=1905)   ## NB: head(bomdts) does not work
window(bomdts, start=2005, end=2008)
plot(bomdts)
## plot(bomdts[, -(1:8)])
## plot(bomdts[, (1:8)])
## plot(bomdts[, 7])

plot(bomdts[, "mdbAVt"], panel=panel.smooth)
plot(bomdts[, "mdbAVt"], panel=points, type="p")

## For greater flexibility, plot using the data frame
plot(bomregions[, "Year"], bomregions[, "mdbAVt"], type="b")
plot(bomregions[, "Year"], bomregions[, "mdbAVt"], type="l")
lines(lowess(bomregions[-(1:10), "Year"], bomregions[-(1:10), "mdbAVt"]))
lines(lowess(bomregions[-(1:10), "Year"], bomregions[-(1:10), "mdbAVt"],
             f=0.1), col="red")
lines(lowess(bomregions[-(1:10), "Year"], bomregions[-(1:10), "mdbAVt"],
             f=0.2), col="blue")

plot(bomregions[, "Year"], bomregions[, "mdbAVt"], type="l")
with(bomregions[-(1:10), ], lines(lowess(mdbAVt ~ Year)))
with(bomregions[-(1:10), ], lines(lowess(mdbAVt ~ Year, f=0.1), 
                                  col="red"))
with(bomregions[-(1:10), ], lines(lowess(mdbAVt ~ Year, f=0.2), 
                                  col="blue"))

## Use of xyplot()  (lattice package)
library(lattice)
xyplot(mdbAVt ~ Year, data=bomregions, type=c("l","smooth"),
        span=0.2)
xyplot(mdbAVt ~ Year, data=bomregions, type=c("l","smooth"),
       span=0.2, par.settings=simpleTheme(col.line="red"))

xyplot(mdbAVt ~ Year, data=bomregions, type="l", span=0.2)
trellis.focus()
with(bomregions[-(1:10), ], panel.loess(Year, mdbAVt, span=2/3,
     col="red"))
trellis.unfocus()

names(ais)
xyplot(mdbAVt ~ Year, data=bomregions, type="l", span=0.2)
trellis.focus()
with(bomregions[-(1:10), ], panel.loess(Year, mdbAVt, span=2/3,
     col="red"))
with(bomregions[-(1:10), ], panel.identify(Year, mdbAVt, 
       labels=paste(1910:2008)))
## Click, click, ... etc.
trellis.unfocus()

## Plots of weather data, again
xyplot(mdbAVt ~ Year, data=bomregions, type="l", span=0.2)
trellis.focus()
with(bomregions[-(1:10), ], panel.loess(Year, mdbAVt, span=2/3,
     col="red"))
with(bomregions[-(1:10), ], panel.identify(Year, mdbAVt, 
       labels=paste(1910:2008)))
trellis.unfocus()

## OR (probably better)
xyplot(mdbAVt ~ Year, data=bomregions, type="l", span=0.2)
trellis.focus()
xy <- trellis.panelArgs()
with(xy, panel.loess(x[-(1:10)], y[-(1:10)], span=2/3,
     col="red"))
with(xy, panel.identify(x[-(1:10)], y[-(1:10)], 
       labels=paste(1910:2008)))
## Click, click, click, ..., then right clicklibrary
trellis.unfocus()

## Time series analysis
plot(bomdts[,"northRain"])
with(bomregions, lines(lowess(northRain ~ Year)))

## Lag plot
with(bomregions, lag.plot(na.omit(mdbAVt), lags=7, 
                          do.lines=FALSE))
## Autocorrelation function
with(bomregions, acf(na.omit(mdbAVt)))
with(bomregions, pacf(na.omit(mdbAVt)))

with(bomregions, lag.plot(na.omit(mdbAVt), lags=1, 
                           do.lines=FALSE))
with(bomregions, lag.plot(na.omit(mdbAVt), lags=7, 
                           do.lines=FALSE))
with(bomregions, acf(na.omit(mdbAVt)))
with(bomregions, pacf(na.omit(mdbAVt)))
with(bomregions, plot(northRain ~ SOI))

## Time series regression
library(forecast)
north.arima <- with(bomregions, 
                    auto.arima(northRain^(1/3), xreg=SOI))
plot(resid(north.arima) ~ bomregions$Year)
lines(lowess(resid(north.arima) ~ bomregions$Year))
northCO2.arima <- with(bomregions, 
                    auto.arima(northRain^(1/3), 
                               xreg=cbind(SOI,CO2))) 
hat <- fitted(northCO2.arima) 
plot(hat~bomregions$Year)   
plot(unclass(hat)~bomregions$Year)
lines(unclass(hat)~bomregions$Year, col="blue")

lines(bomregions$northRain^(1/3)~bomregions$Year, col=2)

library(forecast)
north.arima <- with(bomregions, auto.arima(northRain, xreg=SOI))                         
north.arima
plot(resid(north.arima) ~ bomregions$Year)
north.arima <- with(bomregions, 
                    auto.arima(northRain^(1/3), xreg=SOI)) 
north.arima

## Incorporate CO2 as explanatory variable
plot(resid(north.arima) ~ bomregions$Year)
plot(resid(north.arima) ~ bomregions$Year)
lines(lowess(resid(north.arima) ~ bomregions$Year))
northCO2.arima <- with(bomregions, 
                    auto.arima(northRain^(1/3), 
                               xreg=cbind(SOI,CO2))) 
northCO2.arima

hat <- fitted(northCO2.arima) 
plot(hat~bomregions$Year)  
lines(hat~bomregions$Year, col="blue")
plot(unclass(hat) ~ bomregions$Year)  
lines(unclass(hat) ~ bomregions$Year, col="blue")

## The playwith package
library(DAAGxtras)
library(playwith)
playwith(xyplot(age ~ distance, data=hotspots), labels=hotspots$name)

## Sweave
Sweave("r-exercises", stylepath=TRUE)
help(Sweave)
dir()
file.show("r-exercises.tex")
dir(paste(R.home(),"share/texmf/", sep="/"))

## Distinguish between x and "x"
x <- 2
x
xtxt <- "x"
ls()
rm(list=ls())
rm(list="x")
ls()
rm(list=ls())

## Data frames as lists
library(DAAGxtras)
nihills["climb"]
nihills[["climb"]]
as.matrix(nihills["climb"])
class(as.matrix(nihills["climb"]))

## Factors and ordered factors
fac1 <- c("claudio", "kate", "Pramad")
fac2 <- c("John", "Mark", "claudio")
c(1:4, 6:9)
fac1
fac2
c(fac1, fac2)
fac1 <- factor(fac1)
fac2 <- factor(fac2)
fac1

## Use c() only to concatenate factors with the same levels vector
c(fac1, fac2)

## Ordering issues
fac1 <- factor(c("high", "medium", "low", "medium"),
               levels=c("low","medium","high"))
fac1
fac1 > "low"  # Returns an error

## Ordered factors
fac1 <- ordered(c("high", "medium", "low", "medium"),
                levels=c("low","medium","high"))
fac1
fac1 > "low" # This now makes sense

## Lengths and number of characters
chr<-character(0)
length(chr)
chr3 <- character(3)
chr3
character()
w <- "word"
length(w)
s <- "sex"
nchar(w)
nam <- c("claudio", "kate", "Pramad")
nchar(nam)
length(nam)

## Missing values
z <- c(1,NA,4, 3, NA)
na.omit(z)
zo <- na.omit(z)
length(zo)
is.na(z)
1==1
1==2
z == NA
z

## Use of any()
nn <- 3:9
nn
any(nn>10)
any(nn>8)

## User functions
meanANDsd <- function(x){
  av <- mean(x)
  sdev <- sd(x)
  c(mean = av, sd=sdev)
}

meanANDsd
meanANDsd(1:3)
meanANDsd <- edit(meanANDsd)  
  # Edit to give x a default, e.g. x=rnorm(20)
meanANDsd
meanANDsd()
meanANDsd()


## Thursday afternoon

## S3 model objects are lists
names(nihills)
library(DAAGxtras)
nihills0.lm <- lm(time ~ dist + climb, data=nihills)
nihills0.lm
length(nihills0.lm)
names(nihills0.lm)
print(nihills0.lm)
nihills0.lm$call
nihills0.lm$rank
nihills0.lm$coefficients
coefficients(nihills0.lm)  # coefficients is an extractor function
head(nihills0.lm$model)
names(nihills0.lm)
summary(nihills0.lm)

## Tabulation
library(DAAGxtras)
## help(nassCDS)
dim(nassCDS)
with(nassCDS, table(sex,dead))
with(nassCDS, xtabs(~ sex + dead))
with(nassCDS, xtabs(weight ~ sex + dead))
## help(excessRisk)
tab <- with(nassCDS, xtabs(weight ~ sex + dead))
tab
class(tab)
z <- "tab"
class(z)
## class("tab")  # Tells nothing that we did not know!
class(tab)
tab <- with(nassCDS, xtabs(weight ~ sex + dead +dvcat))
tab
tab <- with(nassCDS, xtabs(weight ~ dvcat + sex + dead))
tab
tab <- with(nassCDS, xtabs(weight ~ dvcat + dead + sex))
tab
## prop.table(tab, margin="dead")  # Does not work
## prop.table(tab, margin=dead)    # Does not work
prop.table(tab, margin=2)
prop.table(tab, margin=c(1,3))
prop.table(tab, margin=c(1))


## Friday, on the lab PC

## Multilevel models
library(lme4)
library(DAAG)
vfit <- lmer(yield ~ shade + (1|block/shade), kiwishade)
vfit

## Aside: the function datafile() that can be used to place selected
## files in the working directory
datafile("oneBadRow")
dir()
file.show("oneBadRow")
file.show("oneBadRow.txt")

## Detection of data input problems
read.table("oneBadRow.txt")
count.fields("oneBadRow.txt")
tab <- table(count.fields("oneBadRow.txt"))
tab

## Outputting tables as LaTeX code
library(xtable)
xtable(tab)

## Sink output to a file
args(sink)
sink("tab.txt")
xtable(tab)
sink()    # Stop sinking

dir()
## help(read.table)
read.table("oneBadRow.txt", fill=TRUE, header=FALSE)  # Read regardless
datafile("bostonc")
file.show("bostonc.txt")
readLines("bostonc.txt", n=10)
## Challenge.  Now read these data into a data frame  

## The quickplot() function
library(DAAGxtras)
library(ggplot2)
quickplot(Year, mdbRain, data=bomregions, geom=c("point","smooth"))
quickplot(Year, mdbRain, data=bomregions, geom=c("point","smooth"), span=0.1)
## help(quickplot)

library(DAAG)
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
      facets=sex ~ .)
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
      facets=. ~ sex)

quickplot(sport, ht, data=ais, geom=c("boxplot"), facets=. ~ sex) 
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
      color=sex)
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
      facets=. ~ sex, size=2)  # Can replace size by cex
## The following throws an error; this is a bug (It should increase line
## and text size by a factor of 2, without generating a label on the side)
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
       facets=. ~ sex, size=I(2))  # Can replace size by cex
## Mappings (similar effect to groups in lattice)
quickplot(wt, ht, data=ais, geom=c("point", "density2d"), 
      shape=sex, color=sport)
quickplot(wt, ht, data=ais, geom=c("point"), 
      shape=sex, color=sport)
quickplot(sport, ht, data=ais, geom=c("boxplot"), 
      shape=sex, color=sport)
quickplot(sport, ht, data=ais, geom=c("boxplot"), 
       color=sex)
quickplot(wt, ht, data=ais, geom=c("boxplot"))

## 20%, 50% & 80% quantiles 
##    5 d.f. normal splines 
library(splines)
quickplot(Year, seRain, data=bomsoi, geom=c("point", "quantile"),  
          formula = y ~ ns(x,5), quantiles=c(0.2,0.5,0.8) ) 

## Methods (Generic functions)
print          # Display generic function print
plot
summary
summary.lm     # Summary method for lm objects
head(summary.lm)
methods("lm")
methods("summary")
getAnywhere("summary.rlm")

## Classes
x  <- 1:10
class(x) <- "lm"   # Silly; but R does not complain!
x
ls(pattern=".lm$")
print.ordered
ordfac <- factor(1:6)
class(ordfac)
ordfac <- ordered(1:6)
class(ordfac)
ordered
factor

## S3 Classes
## Create a function that generates objects of a new class "address"
## Create a print method for objects of this class.
addfun<- function(nam="", address=""){
  zz <- list(name=nam, address=address)
  class(zz) <- "address"
  zz
}

print.address <- function(add=list(NULL,NULL)){
  cat("\n")
  cat(add$name, "\n")
 cat("Address is:\n", add$address, "\n")
  invisible()
}

add1 <- addfun(nam="John",address="Lyneham") 
print(add1)

addclaude <- addfun("Claudio", "Belconnen")
print(addclaude)

## S4 classes
library(lme4)
head(humanpower)
library(DAAG)
head(humanpower1)
hp.lmList <- lmList(o2 ~ wattsPerKg | id, data=humanpower1)
  # lmList() returns an S4 object
names(hp.lmList)
hp.lmList[["1"]]
slotNames(hp.lmList)
hp.lmList@call

## Interface to SQL database (an aside)
help(sqlQuery)
library(RODBC)
help(sqlQuery)
channel <- odbcConnect("test")
sqlSave(channel, USArrests, rownames = "State", verbose = TRUE)