library(car)

attach(Duncan)

boxplot(prestige ~ type)
identify(type, prestige, rownames(Duncan))

boxplot(logit(prestige) ~ type)
identify(type, logit(prestige), rownames(Duncan))
prob.axis()

# one-way ANOVA

  # contrast coding in R

    # contr.treatment = dummy coding

type
contrasts(type)

options("contrasts")

mod.duncan.1 <- lm(logit(prestige) ~ type)
mod.duncan.1
anova(mod.duncan.1)

    # contr.sum = deviation coding

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

mod.duncan.2 <- lm(logit(prestige) ~ type)
mod.duncan.2
anova(mod.duncan.2)

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

contrasts(type) <- "contr.Sum"
contrasts(type)
summary(lm(logit(prestige) ~ type))

# two-way ANOVA

detach(Duncan)
data(Moore)
attach(Moore)
some(Moore)
dim(Moore)

  # means, standard deviations, counts

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

    # reordering factor levels

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

    # tables of cell means, standard deviations, counts

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

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

tapply(conformity,
       list(Status=partner.status, Authoritarianism=fcategory), length)
       
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"))

mod.moore <- lm(conformity ~ partner.status*fcategory)
mod.moore

Anova(mod.moore)

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

  # using dummy-coding

contrasts(partner.status) <- contrasts(fcategory) <- "contr.Treatment"
contrasts(partner.status)
contrasts(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!