# Some more sophisticated play
# - load data file
# (please look at the data file first)

ions <- read.table("/home/units/r.data/potato-R.dat", header=TRUE)
ions
names(ions)
summary(ions)

Duration                        # An error
ions$Duration                   # Correct

attach(ions)                    # Learn about attaching
Duration                        # Now its OK

# - plot
plot(Duration, Absorption)
plot(Absorption ~ Duration )
plot(Absorption ~ Ions)         # Boxplot for different types of variables
plot(Absorption ~ Duration, pch=19 )
plot(Absorption ~ Duration, pch="&" )
lines(Absorption ~ Duration)

?par



# - bust up by ions type
Absorption[Ions=="B"]
plot( Absorption[Ions=="B"] ~ Duration[Ions=="B"], pch="B")
plot( Absorption[Ions=="R"] ~ Duration[Ions=="R"], pch="R")         # Wrong!

plot( Absorption[Ions=="B"] ~ Duration[Ions=="B"], pch="B")
points( Absorption[Ions=="R"] ~ Duration[Ions=="R"], pch="R")         # Better!

plot( Absorption ~ Duration, type="n")
points( Absorption[Ions=="B"] ~ Duration[Ions=="B"], pch="B")         # Works
points( Absorption[Ions=="R"] ~ Duration[Ions=="R"], pch="R")

# - other types?
plot(Absorption ~ Duration, type="l")
plot(Absorption ~ Duration, type="b")

# - annotations, continuations and named function calls
plot( Absorption ~ Duration,
    type="n",
    xlab="Duration of immersion (in hours)",
    ylab="Ion absorption (mg/100g of water in the tissue)",
    main="Durations against Absorption",
    sub="(Source:  Hand et al;, dataset 305)")
points( Absorption[Ions=="B"] ~ Duration[Ions=="B"], pch="B", col=1)
points( Absorption[Ions=="R"] ~ Duration[Ions=="R"], pch="R", col=2)

# - legend
legend(60, 7, legend=c("Bromide", "Rubidium"), pch=c("B" , "R"), col=c(1,2) )
# - beautiful

# Regression
# - basic
lm( Absorption ~ Duration )

# more info
ions.lm <- lm( Absorption ~ Duration )
summary(ions.lm)

names(ions.lm)

win.graph()                         # New graphics window
plot( ions.lm$residuals)
plot( resid(ions.lm) )              # Same

dev.off()

# - interactions and factors
ions.lm2 <- lm( Absorption ~ Duration + Ions )
summary(ions.lm2)

# - full interaction
ions.lm3 <- lm( Absorption ~ Duration * Ions )
summary(ions.lm3)

# - anova
anova(ions.lm3)

plot( resid(ions.lm3) )

qqnorm( resid(ions.lm3) )
qqline( resid(ions.lm3) )

detach(ions)