##------------------------------------------------------------------------##
## Script for Soc 740  Lecture 5                                          ##
##    John Fox                                                            ##
##    Winter 2014                                                         ##
##------------------------------------------------------------------------##

# dichtomous factor

library(car)

Davis[12, c(2, 3)] <- Davis[12, c(3, 2)]  # correct data

Davis$sex
class(Davis$sex)
contrasts(Davis$sex)  # dummy regressor for sex

scatterplot(weight ~ repwt | sex, smooth=FALSE, data=Davis)

mod.davis.1 <- lm(weight ~ repwt + sex, data=Davis)  # additive model
summary(mod.davis.1)

mod.davis.2 <- lm(weight ~ repwt*sex, data=Davis)  # model with interaction
summary(mod.davis.2)

linearHypothesis(mod.davis.2,      # slopes AND intercepts
  c("sexM = 0", "repwt:sexM = 0")) # the same for M and F
                                                        
# polytomous factor

View(Prestige)
Prestige <- na.omit(Prestige)
head(Prestige)

    # treatment of factors in R

Prestige$type
class(Prestige$type)
contrasts(Prestige$type)  # dummy regressors

Prestige$type <- factor(Prestige$type, 
    levels=c("bc", "wc", "prof"))  # change order of levels
contrasts(Prestige$type)

    # additive dummy-regression model
    
scatter3d(prestige ~ income + education | type, id.method="none", 
    data=Prestige, parallel=TRUE)

mod.prestige.1 <- lm(prestige ~ income + education + type, data=Prestige)
summary(mod.prestige.1)

mod.prestige.0 <- update(mod.prestige.1, . ~ . - type) # remove type
anova(mod.prestige.0, mod.prestige.1)  # incremental F-test for type

Anova(mod.prestige.1)  # "type II" tests

linearHypothesis(mod.prestige.1, c("typewc = 0", "typeprof = 0")) # equivalent

    # dummy regression with interactions
    
scatter3d(prestige ~ income + education | type, id.method="none", 
    data=Prestige, parallel=FALSE)

mod.prestige.2 <- lm(prestige ~ income + education +
                   type + income:type + education:type, data=Prestige)

mod.prestige.2 <- lm(prestige ~ income*type + education*type, 
    data=Prestige) # equivalent

mod.prestige.2 <- lm(prestige ~ type*(income + education), 
    data=Prestige) # equivalent

mod.prestige.2 <- update(mod.prestige.1, 
    . ~ . + income:type + education:type) # equivalent
    
summary(mod.prestige.2)

Anova(mod.prestige.2)  # type-II tests

    # "effect" plots

library(effects)  # effects package must be installed
plot(allEffects(mod.prestige.2))
plot(allEffects(mod.prestige.2), multiline=TRUE)

# one-way ANOVA

Boxplot(prestige ~ type, data=Duncan)

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

Duncan$type
contrasts(Duncan$type)

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

# two-way ANOVA

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

    # 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)
})

    # ANOVA

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

Anova(mod.moore) # Type II tests
anova(mod.moore) # Type I tests -- do not use!