# Data from table 2.6 in Agresti
# we compare logistic regression and Poisson regression and show that they
# are the same

victim <- factor(rep(c("white", "black"), each = 2))
defendant <- factor(rep(c("white", "black"), times = 2))
deathpenalty <- matrix(c(53, 11, 0, 4, 414, 37, 16, 139),
        ncol = 2, dimnames = list(NULL, c("yes", "no")))
n <- rowSums(deathpenalty)
n

data.frame(victim, defendant, deathpenalty)

# logistic regression (binomial glm with logit link)
# same analysis as was done in mcmc-too.R

bout <- glm(deathpenalty ~ victim + defendant, family = binomial)
summary(bout)

# Now try to do the same thing assuming cells are independent poisson
# have to stretch data out into a vector

resp <- as.vector(deathpenalty)

# then have to make all the predictors twice as long
victim <- rep(victim, times = 2)
defendant <- rep(defendant, times = 2)
deathpenalty <- factor(rep(colnames(deathpenalty), each = 4))

data.frame(victim, defendant, deathpenalty, resp)

# have to include in model the marginal(s) that the product multinomial
# assumes.  Those are victim * defendant.  Also what were simple effects
# before now have to be "crossed" with deathpenalty to be the same in
# this model.

pout <- glm(resp ~ victim * defendant + deathpenalty * (victim + defendant),
    family = poisson)
summary(pout)

# compare mean value parameters

bout.mu <- predict(bout, type = "response")
pout.mu <- predict(pout, type = "response")
bout.mu
pout.mu

# clearly the result for binomial is "success" probabilities

bout.mu.too <- c(n * bout.mu, n * (1 - bout.mu))
bout.mu.too
pout.mu
all.equal(as.vector(bout.mu.too), as.vector(pout.mu))

# compare saturated model canonical parameters

bout.theta <- predict(bout)
pout.theta <- predict(pout)
bout.theta
pout.theta

logit <- function(x) log(x / (1 - x))
all.equal(bout.theta, logit(bout.mu))
all.equal(pout.theta, log(pout.mu))

# the result for binomial has canonical parameter zero for "failure" cells

pout.theta.too <- pout.theta[1:4] - pout.theta[5:8]
all.equal(bout.theta, pout.theta.too)