##------------------------------------------------------------------------##
## Script for Chapter 7, An R and S-PLUS Companion to Applied Regression  ##
##    John Fox                                                            ##
##    Sage Publications, 2002                                             ##
##------------------------------------------------------------------------##

# preliminaries

options(width=65)
options(digits=5)

# Basics

args(plot.default)

plot(c(0,1), c(0,1), type='n', xlab='', ylab='')

par('col')

par()

plot(1:25, xlab='Symbol Number', ylab='', type='n')
for (pch in 1:25) points(pch, pch, pch=pch)
lines(1:25, type='h', lty=2)

plot(1:26, xlab='letters', ylab='', type='n', axes=F)
box()
for (letter in 1:26) 
    points(letter, 27 - letter, pch=letters[letter])

  ## R alternative
    
    plot(1:25, pch=1:25, xlab='Symbol Number', ylab='')
    lines(1:25, type='h', lty='dashed')

    plot(26:1, xlab='letters', ylab='', 
        pch=letters, axes=F, frame=T)

  
plot(c(1,7), c(0,1), type='n', axes=F, 
    xlab='Line Type (lty)', ylab='')
box()
axis(1, at=1:6)
for (lty in 1:6) lines(c(lty, lty, lty + 1), 
                    c(0, 0.5, 1), lty=lty)


  ## R alternative:
  
    plot(c(1,7), c(0,1), type='n', axes=F, frame=T,
        xlab='Line Type (lty)', ylab='')
    axis(1, at=1:6)
    for (lty in 1:6) lines(c(lty, lty, lty + 1), 
                        c(0, 0.5, 1), lty=lty)

# text 

par(mfrow=c(1,2))
    
plot(c(0,1), c(0,1), axes=F, type='n', xlab='', ylab='')
box()
text(x=c(.2, .5), y=c(.2, .7), 
    c('example text', 'another string'))
title('(a)')

plot(c(0,1), c(0,1), axes=F, type='n', xlab='', ylab='')
box()
text(locator(3), c('one','two','three'))
title('(b)')

locator()

# arrows

plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='')
arrows(x0=1:5, y0=rep(0.1, 5), 
    x1=1:5, y1=seq(0.3, 0.9, len=5), code=3)
title('(a) arrows')


plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='')
segments(x0=1:5, y0=rep(0.1, 5), 
    x1=1:5, y1=seq(0.3, 0.9, len=5))
title('(b) segments')

    ## S-PLUS alternative:
    
    plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='')
    arrows(x1=1:5, y1=rep(0.1, 5), 
        x2=1:5, y2=seq(0.3, 0.9, len=5))
    title('(a) arrows')

    plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='')
    segments(x1=1:5, y1=rep(0.1, 5), 
        x2=1:5, y2=seq(0.3, 0.9, len=5))
    title('(b) segments')

# polygons

plot(c(0,1), c(0,1), type='n', xlab='', ylab='')
polygon(c(.2,.8,.8), c(.2,.2,.8), col=1)
polygon(c(.2,.2,.8), c(.2,.8,.8))

# legend

plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='', 
    frame.plot=T)
legend(locator(1), legend=c('group A', 'group B', 'group C'),
    lty=1:3, pch=1:3)
    
    # S-Plus alternative:
    
    plot(c(1,5), c(0,1), axes=F, type='n', xlab='', ylab='')
    legend(locator(1), legend=c('group A', 'group B', 'group C'),
        lty=1:3, marks=1:3)

# colors

piechart(rep(1, length(palette())), col=palette())
    
palette()

rainbow(10)

gray(0:8/8)

piechart(rep(1,100), col=rainbow(100), labels=rep('',100))

piechart(rep(1,100), col=gray(0:100/100), labels=rep('',100))

    # S-PLUS alternative:
  
    piechart(rep(1,16), col=1:16)


# Extended example: effect plots

library(car)
data(Cowles)
dim(Cowles)
Cowles[sort(sample(1421, 10)),]

mod.cowles <- glm(volunteer ~ neuroticism * extraversion + sex, 
    data=Cowles, family=binomial)
summary(mod.cowles)

    # show full response surface
    
neuroticism <- 0:24
extraversion <- 0:24
sex <- c('male', 'female')
graph.data <- expand.grid(neuroticism=neuroticism, extraversion=extraversion, sex=sex)
graph.data$fit <- predict(mod.cowles, newdata=graph.data, type='response')
graph.data

par(mfrow=c(2,1))    

prob <- matrix(graph.data$fit[graph.data$sex=='male'], 25, 25)
persp(neuroticism, extraversion, prob,
     phi=30, theta=45, expand=0.65, d=2, shade=0.75, 
     ticktype='detailed', zlab='Probability(Volunteer)', main='Males')
    
prob <- matrix(graph.data$fit[graph.data$sex=='female'], 25, 25)
persp(neuroticism, extraversion, prob,
     phi=30, theta=45, expand=0.65, d=2, shade=0.75, 
     ticktype='detailed', zlab='Probability(Volunteer)', main='Females')

    # make effect plot    
    
extraversion <- 0:24
neuroticism <- seq(0, 24, by=6)
graph.data <- expand.grid(neuroticism=neuroticism, extraversion=extraversion)
X <- cbind(constant=1, as.matrix(graph.data), sex=0.5, 
    neuro.extra=graph.data$neuroticism * graph.data$extraversion)
X

logit <- X %*% coefficients(mod.cowles)
se <- sqrt(diag(X %*% Var(mod.cowles) %*% t(X)))
prob <- matrix(1/(1+exp(-logit)), 5, 25,)
low <- matrix(1/(1+exp(-(logit - se))), 5, 25)
high <- matrix(1/(1+exp(-(logit + se))), 5, 25)
prob

plot(range(extraversion), range(c(low, high + 0.05)), type='n', 
    xlab='Introversion-Extraversion', 
    ylab='Probability of Volunteering',
    xaxt='n', xlim=c(0,30))
    
axis(1, at=seq(0, 24, by=6))

for (neuro in 1:5){
    lines(extraversion, prob[neuro,])
    text(25, prob[neuro,25], paste('N = ', neuroticism[neuro]), adj=0)
    }

extra <- seq(1, 25, by=6)
for (neuro in c(1, 3, 5)){
        arrows(extraversion[extra], low[neuro, extra],
            extraversion[extra], high[neuro, extra], 
            angle=90, code=3, lty=2, length=0.05)
            }

    ## S-PLUS alternative:
                
    for (neuro in c(1, 3, 5)){
            segments(extraversion[extra], low[neuro, extra],
                extraversion[extra], high[neuro, extra])
                }

text(locator(1), 'Stability-Neuroticism', adj=1, cex=0.75)