TV <- read.csv("http://www.webpages.uidaho.edu/~chrisw/stat407507/tvshowsoct19.csv",na.strings = ".",header=T)


TV$network  <-  as.factor(TV$network)

TV$type  <-  as.factor(TV$type)

TV$year  <-  as.factor(TV$year)


  # change contrast options before fitting models, this way Type III SS 
  # from the car package will agree with output from SAS
options(contrasts=c(unordered="contr.sum", ordered="contr.poly"))

TV.lm  <- lm(rating ~ network*type*year, data = TV)

anova(TV.lm)


library(car)

 # this produces the results using Type II SS 
Anova(lm(rating ~ network*type*year, data = TV),type="II")

Anova(lm(rating ~ network*type*year, data = TV),type="III")



library(afex)

library(emmeans)

 # the mixed function in the afex package produces a mixed-model analysis like that
 # produced by Proc Mixed in SAS

TV.afex1  <- mixed(rating ~ network + type + network:type + (1|year) +
                     (1|network:year) + (1|type:year) + (1|network:type:year), data = TV)

summary(TV.afex1)

nice(TV.afex1)

emmeans(TV.afex1, ~network)

emmeans(TV.afex1, ~type)