## Section 4.7 Examples:

## Example 1: 

vice.data <-
array(c(911,44,538,456,
        3,2,43,279),
      dim = c(2, 2, 2),
      dimnames = list(
          Cigarette = c("Yes", "No"),
          Marijuana = c("Yes", "No"),
          Alcohol = c("Yes", "No")))
vice.data <- as.table(vice.data)
print(vice.data)


library(MASS) # load the "MASS" package

# No-interaction model:

loglm(~ Cigarette + Marijuana + Alcohol, data=vice.data)  

# Model with only two-way interactions:

loglm(~ Cigarette + Marijuana + Alcohol + Cigarette:Marijuana + Cigarette:Alcohol + Marijuana:Alcohol, data=vice.data)


# The 'step' function picks the "best" model according to the Akaike Information Criterion (AIC):

# Start with the "saturated model" containing ALL possible interactions and work backward:

saturated.model <- loglm(~ Cigarette * Marijuana * Alcohol, data=vice.data)

step(saturated.model, direction="backward")

# The best model?

best.model <-loglm(~ Cigarette + Marijuana + Alcohol + Cigarette:Marijuana + Cigarette:Alcohol + Marijuana:Alcohol, data=vice.data)

best.model

print(fitted(best.model))

print(coef(best.model))

## Some odds ratios:

# Does the odds of a cigarette smoker using marijuana (relative to a non-cigarette smoker) depend on alcohol use?

# Odds Ratio[Cig. Smoker Odds / Cig. Non-Smoker Odds], for Alcohol=Yes:
 
(910.38159 / 538.6184) / (44.61841 / 455.3816)

# Odds Ratio[Cig. Smoker Odds / Cig. Non-Smoker Odds], for Alcohol=No:

(3.616797 / 42.38354) / (1.383203 / 279.61646)

####################################################

## Example 2:

cuttings.data <-
array(c(156,84,84,156,
        107,133,31,209),
      dim = c(2, 2, 2),
      dimnames = list(
          Survival = c("Alive", "Dead"),
          Planting.Time = c("Early", "Late"),
          Length = c("Long", "Short")))
cuttings.data <- as.table(cuttings.data)
print(cuttings.data)

library(MASS) # load the "MASS" package

# The 'step' function picks the "best" model according to the Akaike Information Criterion (AIC):

# Start with the "saturated model" containing ALL possible interactions and work backward:

saturated.model <- loglm(~ Survival * Planting.Time * Length, data=cuttings.data)

step(saturated.model, direction="backward")

# Slightly simpler model:

simp.model <- loglm(~ Survival + Planting.Time + Length + Survival:Planting.Time + Survival:Length + Planting.Time:Length, data=cuttings.data)

simp.model

# No-interaction model:

v.simp.model <- loglm(~ Survival + Planting.Time + Length, data=cuttings.data)

v.simp.model

print(fitted(simp.model))

print(coef(simp.model))

## Some odds ratios:

# Does the odds of a early plant surviving (relative to a late plant) depend on cutting length?

# Odds Ratio[Early Plant Odds / Late Plant Odds], for Length=Long:
 
(161.09 / 78.91) / (78.90 / 161.09)

# Odds Ratio[Early Plant Odds / Late Plant Odds], for Length=Short:

(101.90 / 138.10) / (36.10 / 203.90)


####################################################

## Example 3:

data(Titanic)  # loading the built-in data set

Titanic3 <- margin.table(Titanic, c(4,2,1))[2:1,,]

print(Titanic3)  # 2-by-2-by-4 table

library(MASS) # load the "MASS" package

# The 'step' function picks the "best" model according to the Akaike Information Criterion (AIC):

# Start with the "saturated model" containing ALL possible interactions and work backward:

saturated.model <- loglm(~ Class * Sex * Survived, data=Titanic3)

step(saturated.model, direction="backward")

# Is it possible to leave off the second-order interaction?

loglm(~ Class + Sex + Survived + Class:Sex + Class:Survived + Sex:Survived, 
    data = Titanic3)


best.model <-  loglm(~ Class * Sex * Survived, data = Titanic3)

best.model

print(fitted(best.model))

print(coef(best.model))

## Some odds ratios:

# Does the odds of a female surviving (relative to a male) depend on Class?

# Odds Ratio[Female Odds / Male Odds], Class=1st:

(141/4)/(62/118)

# Odds Ratio[Female Odds / Male Odds], Class=2nd:

(93/13)/(25/154)

# Odds Ratio[Female Odds / Male Odds], Class=3rd:

(90/106)/(88/422)

# Odds Ratio[Female Odds / Male Odds], Class=Crew:

(20/3)/(192/670)