########################################################################
##
## Gj: 19-Apr-2000
##
## Splus commands used for Short Course on Splus, part 2.
##
########################################################################

### ON SEARCH ###
search()  ## where S-Plus looks for objects
args(ls)  ## the arguments to ls(), and def.
attach(tox)  ## attaching the tox data.frame
search()
ls(pos=2)
detach(2)  ## then detach the data.frame
search()

########################################################################
### DATA ###

## the tox data:
tox <- read.table("/home/gardar/pub/toxicity.dat",
                  header=T)
names(tox) <- c('d','n','r')  ## change names
tox$p <- tox$r/tox$n * 100    ## percentage killed

## the fuel.frame:
class(fuel.frame)
dim(fuel.frame)
summary(fuel.frame)

########################################################################
### TRADITIONAL GRAPHICS ###

## the plot function:
attach(tox)
plot(d,p)
plot(d,p,main="Toxicity Data",ylab="% killed",
     xlab="Concentration",sub="(what ever!)")
plot(d,p,type="b")
plot(d,p,type="b",lty=3,lwd=3,pch=16,cex=2)
plot(d,p,type='n',axes=F,xlab='',ylab='')

## updating tox:
detach('tox')
tox$log.d <- log(tox$d)
tox$logit.p <- log(tox$p/(100-tox$p)) 
tox
attach(tox)
ls(pos=2)

## plotting transformed variables -- big example (skip)
par(mar=c(5,5,6,5))  ## change the margins
plot(log.d,logit.p, type='n',  ## plot nothing
     axes=F,                   ## no axes
     xlab='',ylab='')          ## no labels
box()  ## make a box, def. is bty='o', try 'l'
axis(1)  ## add axis to side 1
axis(2)  ## and to side 2
lines(log.d,logit.p,type='b',pch=16)
axis(4,at=logit.p,  ## put ticks at logit.p 
     labels=round(p,1),  ## the labels used
     srt=90,  ## rotate the labels
     tck=-par()$tck,  ## ticks inside
     cex=0.8,  ## reduce the character size
     mgp=c(2,0.5,0))  ## lines for (Ti,La,Ax)
axis(3,at=log.d,labels=d,
     tck=-par()$tck,cex=0.8,mgp=c(2,0.5,0))
title(xlab='log concentration',
      ylab='logit of % killed')
mtext('The toxicity data',  ## margin text
      adj=0,  ## adjust to left (1 is right)
      side=3,line=4,cex=1.5)
mtext('% killed',side=4,line=2.5,cex=0.8)
mtext('concentration',side=3,line=2.5,cex=0.8)
title(sub=date(),cex=0.5,adj=1)
abline(h=0,lty=2)  ## horizontal line
legend(locator(1),  ## use mouse to pick loc.
       '50% killed',bty='n')  ## bty='n', no box

## use text:
plot(log.d,logit.p,type='n')
lines(log.d,logit.p,type='b',pch='    ',cex=3)
text(log.d,logit.p,paste(round(p),'%',sep=''))

detach('tox')

## other functions:
attach(fuel.frame)
boxplot(split(Mileage,Type),varwidth=T)
qqnorm(Mileage)
qqline(Mileage)
hist(Mileage)

########################################################################
### TRELLIS ###

xyplot(logit.p ~ log.d, data=tox)

xyplot(Mileage ~ Weight | Type, data=fuel.frame)
bwplot(Type ~ Mileage, data=fuel.frame)
histogram(~ Mileage | Type, data=fuel.frame)
qqmath(~ Mileage | Type, data=fuel.frame)

## to export to ps-file:
trellis.device(postscript,file='myplot.ps',
               horizontal=T)
xyplot(logit.p ~ log.d, data=tox)
dev.off()


########################################################################
### UNIVARIATE TESTS ###

attach(fuel.frame)
t.test(Mileage[Type=='Medium'],
       Mileage[Type=='Sporty'],
       var.equal=F)
detach('fuel.frame')


########################################################################
### MODELS ###

tox$logit.p <- log(tox$p/(100-tox$p))  ## don't add 5
ula <- lm(logit.p ~ log.d, data=tox)
ula
summary(ula)  
par(mfrow=c(2,3))  ## 2 rows and 2 col of plots
plot(ula) ## get 6 plots on one page.

## LC50
co <- coef(ula)
co   ## just to show what it does...
lc50 <- exp(-co[1]/co[2])
lc50

## plot:
par(mfrow=c(1,1))  ## one plot on device
plot(tox$log.d, tox$logit.p,
     xlab="log Con.",ylab="logit P")
abline(coef(ula))  ## or lines(log.d,fitted(ula))
abline(h=0,lty=2)  ## horizontal 50% line
abline(v=log(lc50),lty=2)  ## the LC50 line


## glm:
la.glm <- glm(cbind(r,n-r) ~ log.d, data=tox,
               family=binomial(link=logit))
la.glm
summary(la.glm)

## LC50
co <- coef(la.glm)
co   ## just to show what it does...
lc50 <- exp(-co[1]/co[2])
lc50

########################################################################