##------------------------------------------------------------##
##       Script for Session 2: Statistical Models in R        ##
##                       John Fox                             ##
##  Introduction to the R Statistical Computing Environment   ##
##                ICPSR Summer Program                        ##
##                         2010                               ##
##------------------------------------------------------------##

# Linear models in R
    
    # multiple regression (Prestige data)
    
library(car)
names(Prestige)
some(Prestige)  # 10 randomly sampling observations

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

confint(prestige.mod)

    # dummy regression
    
Prestige$type # a factor
class(Prestige$type) 
str(Prestige$type) # structure

Prestige.2 <- na.omit(Prestige) # filter out missing data
nrow(Prestige)
nrow(Prestige.2)
levels(Prestige.2$type)
Prestige.2$type <- with(Prestige.2, factor(type, levels=c("bc", "wc", "prof"))) # reorder levels
Prestige.2$type

        # generating contrasts from factors
        
getOption("contrasts")
contrasts(Prestige.2$type)
model.matrix(~ type, data=Prestige.2)

contrasts(Prestige.2$type) <- contr.treatment(levels(Prestige.2$type), base=2)  # changing baseline category
contrasts(Prestige.2$type)

contrasts(Prestige.2$type) <- "contr.helmert"  # Helmert contrasts
contrasts(Prestige.2$type)

contrasts(Prestige.2$type) <- "contr.sum"  # "deviation" contrasts
contrasts(Prestige.2$type)

contrasts(Prestige.2$type) <- NULL  # back to default

Prestige.2$type.ord <- ordered(Prestige.2$type, levels=c("bc", "wc", "prof")) # ordered factor
Prestige.2$type.ord
round(contrasts(Prestige.2$type.ord), 3)   # orthogonal polynomial contrasts

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

anova(prestige.mod.1)   # sequential ("type-I") tests

prestige.mod.0 <- lm(prestige ~ income + education, data=Prestige.2) # note: NA's filtered!
summary(prestige.mod.0)

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

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

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

lm(prestige ~ log2(income*type) + education*type, data=Prestige.2) # equivalent specifications
lm(prestige ~ (log2(income) + education)*type, data=Prestige.2)

        # effect displays
    
library(effects)
plot(allEffects(prestige.mod.3), ask=FALSE)

    # Anova Models

some(Moore)

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

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

with(Moore, tapply(conformity,
    list(Authoritarianism=fcategory, "Partner's Status"=partner.status),
    mean))
with(Moore, tapply(conformity,
    list(Authoritarianism=fcategory, "Partner's Status"=partner.status),
    sd))

            # graph of means:
            
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 tables

contr <- options(contrasts=c("contr.sum", "contr.poly"))   # contr.sum = deviation contrasts
moore.mod <- lm(conformity ~ fcategory*partner.status, data=Moore)
summary(moore.mod)

Anova(moore.mod)    # type II sums of squares
Anova(moore.mod, type="III")    # type III sums of squares

options(contr) # restore defaults


    # more on lm

args(lm)

some(Davis)
lm(weight ~ repwt, data=Davis, subset=sex == "F")  # observation selection (women only)
lm(weight ~ repwt, data=Davis, subset=1:100)

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

lm(conformity ~ partner.status*fcategory,  # specifying contrasts
    contrasts=list(partner.status=contr.sum, fcategory=contr.poly),
    data=Moore)

lm(100*conformity/40 ~ partner.status*fcategory, data=Moore)  # data argument; note computation of y

lm(prestige~I(income + education), data=Duncan)  # "protecting" expresssion on RHS of the model


# Generalized linear models
    
    # binary logit model
    
some(Mroz)

mroz.mod <- glm(lfp ~ k5 + k618 + age + wc + hc + lwg + inc,
   data=Mroz, family=binomial)
summary(mroz.mod)

round(exp(cbind(Estimate=coef(mroz.mod), confint(mroz.mod))), 2) # odds ratios
    
mroz.mod.2 <- update(mroz.mod, . ~ . - k5 - k618)
anova(mroz.mod.2, mroz.mod, test="Chisq") # likelihood-ratio test

Anova(mroz.mod)  # analysis-of-deviance table

plot(allEffects(mroz.mod), ask=FALSE)

    # Poisson regression

some(Ornstein)
nrow(Ornstein)
(tab <- xtabs(~interlocks, data=Ornstein))

x <- as.numeric(names(tab)) # the names are the distinct values of interlocks
plot(x, tab, type="h", xlab="Number of Interlocks", ylab="Frequency")
points(x, tab, pch=16)

mod.ornstein <- glm(interlocks ~ log2(assets) + nation + sector,
   family=poisson, data=Ornstein)
summary(mod.ornstein)
Anova(mod.ornstein)

        # quasi-Poisson model, allowing for overdispersion
        
mod.ornstein.q <- update(mod.ornstein, family=quasipoisson)
summary(mod.ornstein.q)

plot(allEffects(mod.ornstein.q, default.levels=50), ask=FALSE)

# repeated-measures ANOVA and MANOVA 

some(OBrienKaiser)
?OBrienKaiser
contrasts(OBrienKaiser$treatment)
contrasts(OBrienKaiser$gender)

    # defining the within-subjects design
    
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
    levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
idata

    # fitting the multivariate linear model
    
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, 
                     post.1, post.2, post.3, post.4, post.5, 
                     fup.1, fup.2, fup.3, fup.4, fup.5) ~  treatment*gender, 
                data=OBrienKaiser)
mod.ok
    
    # multivariate and univariate tests

(av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour)) 

summary(av.ok)

    # graphing the means
    
        # reshape the data from "wide" to "long"
    
OBrien.long <- reshape(OBrienKaiser, 
    varying=c("pre.1", "pre.2", "pre.3", "pre.4", "pre.5", 
        "post.1", "post.2", "post.3", "post.4", "post.5", 
        "fup.1", "fup.2", "fup.3",  "fup.4", "fup.5"),
    v.names="score",
    timevar="phase.hour", direction="long")
OBrien.long$phase <- ordered(c("pre", "post", "fup")[1 + ((OBrien.long$phase.hour - 1) %/% 5)], 
    levels=c("pre", "post", "fup"))
OBrien.long$hour <- ordered(1 + ((OBrien.long$phase.hour - 1) %% 5))
dim(OBrien.long)
head(OBrien.long, 25) # first 25 rows

        # compute means

Means <- as.data.frame(ftable(with(OBrien.long, 
    tapply(score, list(treatment=treatment, gender=gender, phase=phase, hour=hour), mean))))
names(Means)[5] <- "score"
dim(Means)
head(Means, 25)

        # graph of means
        
library(lattice)
xyplot(score ~ hour | phase + treatment, groups=gender, type="b", 
    strip=function(...) strip.default(strip.names=c(TRUE, TRUE), ...), 
    ylab="Mean Score", data=Means, auto.key=list(title="Gender", cex.title=1))