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

# preliminaries

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


# simple regression (Davis data)

library(car)
data(Davis)
Davis[1:10,]

attach(Davis)
davis.mod<-lm(weight~repwt)
davis.mod
summary(davis.mod)

plot(repwt, weight)
abline(davis.mod)
abline(0, 1, lty=2)

identify(repwt, weight)

davis.mod.2<-update(davis.mod, subset=-12)
summary(davis.mod.2)
abline(davis.mod.2, lty=3) # readers: add this line

detach(Davis)


# multiple regression (Prestige data)

data(Prestige)
names(Prestige)

attach(Prestige)
prestige.mod<-lm(prestige~income+education+women)
summary(prestige.mod)


Prestige.scaled<-data.frame(
    scale(Prestige[,c("prestige","income","education","women")]))
Prestige.scaled[1:5,]

summary(lm(prestige~income+education+women, data=Prestige.scaled))
summary(lm(prestige~income+education+women-1, data=Prestige.scaled))


    
# dummy regression

type

detach(Prestige)
Prestige.2<-na.omit(Prestige)
attach(Prestige.2)
type
length(type)

class(type)
unclass(type)

options("contrasts")

contrasts(type)

contrasts(type)<-contr.treatment(levels(type), base=2)
contrasts(type)
contrasts(type)<-NULL

contrasts(type) <- 'contr.helmert'
contrasts(type)

contrasts(type) <- 'contr.sum'
contrasts(type)

type.ord<-ordered(type, levels=c("bc", "wc", "prof"))
type.ord
round(contrasts(type.ord), 3)

prestige.mod.1<-lm(prestige~income+education+type)
summary(prestige.mod.1)

anova(prestige.mod.1)

prestige.mod.0<-lm(prestige~income+education)
summary(prestige.mod.0)
anova(prestige.mod.0, prestige.mod.1)

Anova(prestige.mod.1)

prestige.mod.2<-lm(prestige~income+education+type-1)
summary(prestige.mod.2)
Anova(prestige.mod.2)

prestige.mod.3<-update(prestige.mod.1, .~. + income:type + education:type)
summary(prestige.mod.3)
Anova(prestige.mod.3)

lm(prestige~income*type+education*type)

lm(prestige ~ type + (income+education) %in% type)
lm(prestige ~ type + (income+education) %in% type - 1)

detach(Prestige.2)


# Anova Models

data(Moore)
attach(Moore)
Moore

fcategory<-factor(fcategory, levels=c("low","medium","high"))
fcategory

means<-tapply(conformity, list(fcategory, partner.status), mean)
means
tapply(conformity, list(fcategory, partner.status), 
    function(x) sqrt(var(x)))
tapply(conformity, list(fcategory, partner.status), length)

Fcat<-as.numeric(fcategory)
plot(c(0.5, 3.5), range(conformity), xlab="F category", 
    ylab="Conformity", type="n", axes=F)
axis(1, at=1:3, labels=c("low", "medium", "high"))
axis(2)
box()
points(jitter(Fcat[partner.status=="low"]), 
    conformity[partner.status=="low"], pch="L")
points(jitter(Fcat[partner.status=="high"]), 
    conformity[partner.status=="high"], pch="H")
lines(1:3, means[,1], lty=1, lwd=3, type="b", pch=19, cex=2)
lines(1:3, means[,2], lty=3, lwd=3, type="b", pch=1, cex=2)
legend(locator(1), c("high","low"), lty=c(1,3), lwd=c(3,3), pch=c(19,1))

identify(Fcat, conformity)

options(contrasts=c("contr.sum", "contr.poly"))
moore.mod<-lm(conformity~fcategory*partner.status)
summary(moore.mod)
Anova(moore.mod)
Anova(moore.mod, type="III")

    ## -- for the reader
    
options(contrasts=c("contr.helmert", "contr.poly"))
moore.mod.1<-lm(conformity~fcategory*partner.status)
summary(moore.mod.1)
Anova(moore.mod.1)
Anova(moore.mod.1, type="III")

options(contrasts=c("contr.treatment", "contr.poly"))
moore.mod.2<-lm(conformity~fcategory*partner.status)
summary(moore.mod.2)
Anova(moore.mod.2)
Anova(moore.mod.2, type="III")

options(contrasts=c("contr.sum", "contr.poly"))
moore.mod.3<-update(moore.mod, subset=-c(16,19))
Anova(moore.mod.3)

    ## -- end

detach(Moore)

# Contrasts

data(Baumann)
attach(Baumann)
Baumann

tapply(post.test.3, group, mean)
tapply(post.test.3, group, function(x) sqrt(var(x)))
plot(group, post.test.3, xlab="Group", ylab="Reading Score")

contrasts(group)<-matrix(c(1, -0.5, -0.5,  0, 1, -1), 3, 2)
contrasts(group)

summary(lm(post.test.3 ~ group))

detach(Baumann)


# linear hypotheses

data(Duncan)
attach(Duncan)
duncan.mod<-lm(prestige~income+education)
summary(duncan.mod)

linear.hypothesis(duncan.mod, c(0, 1, -1))

davis.mod
davis.mod.2

diag(2)  # order-2 identify matrix
linear.hypothesis(davis.mod, diag(2), c(0,1))
linear.hypothesis(davis.mod.2, diag(2), c(0,1))

# data and confidence ellipses

par(mfrow=c(1,2))
data.ellipse(income, education, levels=c(.5,.75,.9,.95))
identify(income, education, row.names(Duncan))
  title('(a)')
  
confidence.ellipse(duncan.mod)
  title('(b)')

# lm

args(lm)

lm(weight ~ repwt, data=Davis, subset=sex == "F")

lm(weight ~ repwt, data=Davis, subset=1:100)

lm(prestige ~ income + education, data=Duncan, subset=-c(6, 16))

detach(Duncan)
attach(Moore)
lm(conformity ~ partner.status*fcategory, 
    contrasts=list(partner.status=contr.sum, fcategory=contr.poly))
detach(Moore)

lm(100*conformity/40 ~ partner.status*fcategory, data=Moore)

lm(logit(conformity/40) ~ partner.status*fcategory, data=Moore)

lm(prestige~I(income+education), data=Duncan)