##------------------------------------------------------------------------##
## Script for Soc 740  Lecture 6                                          ##
##    John Fox                                                            ##
##    Winter 2012                                                         ##
##------------------------------------------------------------------------##

library(car)

Boxplot(prestige ~ type, data=Duncan)

Boxplot(logit(prestige) ~ type, data=Duncan)
probabilityAxis()

# one-way ANOVA

  # contrast coding in R

    # contr.treatment = dummy coding

Duncan$type
contrasts(Duncan$type)

options("contrasts")

mod.duncan.1 <- lm(logit(prestige) ~ type, data=Duncan)
mod.duncan.1
anova(mod.duncan.1)
Anova(mod.duncan.1) # equivalent

    # contr.sum = deviation coding

contrasts(Duncan$type) <- "contr.sum"
contrasts(Duncan$type)

mod.duncan.2 <- lm(logit(prestige) ~ type, data=Duncan)
mod.duncan.2
anova(mod.duncan.2) # same as before

    # Note: contr.Treatment and contr.Sum from the car package

contrasts(Duncan$type) <- "contr.Sum"
contrasts(Duncan$type)
summary(lm(logit(prestige) ~ type, data=Duncan))

# two-way ANOVA

some(Moore) # Mooore and Krupat conformity data
dim(Moore)

  # means, standard deviations, counts

with(Moore, tapply(conformity,
       list(Status=partner.status, Authoritarianism=fcategory), mean))

    # reordering factor levels

Moore$fcategory <- factor(Moore$fcategory, levels=c("low", "medium", "high"))
Moore$partner.status <- factor(Moore$partner.status, levels=c("low", "high"))

    # tables of cell means, standard deviations, counts

with(Moore, tapply(conformity,
       list(Status=partner.status, Authoritarianism=fcategory), mean))

with(Moore, tapply(conformity,
       list(Status=partner.status, Authoritarianism=fcategory), sd))

xtabs(~ partner.status + fcategory, data=Moore) # cell counts

with(Moore, {       
    interaction.plot(fcategory, partner.status, conformity, type="b", pch=c(1,16),
        cex=2, ylim=range(conformity))
    points(jitter(as.numeric(fcategory), factor=0.5), conformity, 
        pch=ifelse(partner.status == "low", "L", "H"))
    identify(fcategory, conformity)
})

  # Types II and III tests

options("contrasts")
options(contrasts=c("contr.Sum", "contr.poly"))
contrasts(Moore$partner.status)
contrasts(Moore$fcategory)

mod.moore <- lm(conformity ~ partner.status*fcategory, data=Moore)
mod.moore

Anova(mod.moore)

Anova(mod.moore, type="III")

  # using dummy-coding

contrasts(Moore$partner.status) <- contrasts(Moore$fcategory) <- "contr.Treatment"
contrasts(Moore$partner.status)
contrasts(Moore$fcategory)

mod.moore.1 <- update(mod.moore)
mod.moore.1

Anova(mod.moore.1) # Type-II tests OK

Anova(mod.moore.1, type="III")  # Type-III tests WRONG!