This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0/).
The R markdown source for this document is https://www.stat.umn.edu/geyer/5421/notes/ch8.Rmd.
The version of R used to make this document is 4.5.1.
The version of the bookdown package used to make
this document is 0.45.
The version of the rmarkdown package used to make
this document is 2.30.
The version of the CatDataAnalysis package used to
make this document is 0.1.5.
The version of the MASS package used to make this
document is 7.3.65.
The version of the glmbb package used to make this
document is 0.5.1.
library(CatDataAnalysis)
library(MASS)
library(glmbb)To say we are treating one dimension of a contingency table as a “response” is the same as saying we want a product multinomial model with all of the multinomials being for the “response” variable. This means two things.
In order for product multinomial and Poisson to be the same, we must have the interaction of all other factors in the model formula,
In order to allow different levels of the response to have different means, we must have the “response” in the model.
In order to have any other term on the right-hand side of the formula, it must be “interacted with” the “response” variable (because otherwise it is one of the terms from the previous items).
The latter is because we want to know about the effect of other terms on the “response”.
We do the example on alligator food choice (Section 8.1.2 in Agresti)
data(table_8.1)
names(table_8.1)## [1] "lake" "gender" "size" "food" "count"sapply(table_8.1, class)## lake gender size food count
## "integer" "integer" "integer" "integer" "integer"lapply(table_8.1, unique)## $lake
## [1] 1 2 3 4
##
## $gender
## [1] 1 2
##
## $size
## [1] 1 2
##
## $food
## [1] 1 2 3 4 5
##
## $count
## [1] 7 1 0 5 4 2 16 3 13 6 9 8 10These data are given in really inconvenient coding. Let’s recode using the variable names in Table 8.1, which we hope are given in the correct order. (We’ll check that.)
foo <- transform(table_8.1,
lake = factor(lake,
labels = c("Hancock", "Oklawaha", "Trafford", "George")),
gender = factor(gender, labels = c("Male", "Female")),
size = factor(size, labels = c("<=2.3", ">2.3")),
food = factor(food,
labels = c("Fish", "Invertebrate", "Reptile", "Bird", "Other")))
bar <- xtabs(count ~ size + food + gender + lake, data = foo)
bar## , , gender = Male, lake = Hancock
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 7 1 0 0 5
## >2.3 4 0 0 1 2
##
## , , gender = Female, lake = Hancock
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 16 3 2 2 3
## >2.3 3 0 1 2 3
##
## , , gender = Male, lake = Oklawaha
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 2 2 0 0 1
## >2.3 13 7 6 0 0
##
## , , gender = Female, lake = Oklawaha
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 3 9 1 0 2
## >2.3 0 1 0 1 0
##
## , , gender = Male, lake = Trafford
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 3 7 1 0 1
## >2.3 8 6 6 3 5
##
## , , gender = Female, lake = Trafford
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 2 4 1 1 4
## >2.3 0 1 0 0 0
##
## , , gender = Male, lake = George
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 13 10 0 2 2
## >2.3 9 0 0 1 2
##
## , , gender = Female, lake = George
##
## food
## size Fish Invertebrate Reptile Bird Other
## <=2.3 3 9 1 0 1
## >2.3 8 1 0 0 1That matches Table 8.1 in Agresti. Note that we had to write the
formula count ~ size + food + gender + lake for R function
xtabs with the terms in that order to get our printout to
match the table in Agresti.
Now we want to consider food to be the “response”
variable. So we fit all of the models in the first three columns of
Table 8.2 in Agresti using R function loglm in R package
MASS because it gives both Pearson Chi-Square and
likelihood ratio test statistics. (We could also use R function
glm but it does not give the Pearson statistics.)
names(dimnames(bar))## [1] "size" "food" "gender" "lake"lout <- list()
lout[["empty"]] <- loglm(~ size * gender * lake + food, bar)
lout[["gender"]] <- loglm(~ size * gender * lake + food * gender, bar)
lout[["size"]] <- loglm(~ size * gender * lake + food * size, bar)
lout[["lake"]] <- loglm(~ size * gender * lake + food * lake, bar)
lout[["size.lake"]] <- loglm(~ size * gender * lake + food * (size + lake), bar)
lout[["size.gender.lake"]] <- loglm(~ size * gender * lake + food * (size + gender + lake), bar)
woof <- lapply(lout, function(x) c(lrt = x$lrt, pearson = x$pearson, df = x$df))
woof <- as.data.frame(woof)
woof <- as.matrix(woof)
woof <- t(woof)
round(woof, 1)## lrt pearson df
## empty 116.8 106.5 60
## gender 114.7 101.2 56
## size 101.6 86.9 56
## lake 73.6 79.6 48
## size.lake 52.5 58.0 44
## size.gender.lake 50.3 52.5 40This agrees with the left half of Table 8.2 in Agresti. Agresti says LRT is more reliable than Pearson, but since the sample size is not “large” with many small cell counts, neither is accurate. If we really wanted accurate \(P\)-values, we would need to parametric bootstrap, as discussed in the notes on Chapter 9.
R function glmbb in R package glmbb fits
models using R function glm. It considers all models in a
range of models. There will be a handout for
that. But for now, we just show it without much explanation.
bout <- glmbb(big = count ~ lake * gender * size * food,
little = ~ lake * gender * size + food,
family = poisson, data = foo)## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurredsummary(bout)##
## Results of search for hierarchical models with lowest AIC.
## Search was for all models with AIC no larger than min(AIC) + 10
## These are shown below.
##
## criterion weight formula
## 288.0 0.903304 count ~ lake*food + size*food + lake*gender*size
## 293.7 0.050072 count ~ lake*food + gender*food + size*food + lake*gender*size
## 294.9 0.028388 count ~ lake*gender*size + lake*size*food
## 296.8 0.010643 count ~ size*food + lake*gender*size + lake*gender*food
## 297.5 0.007593 count ~ gender*food + lake*gender*size + lake*size*foodThis says that by far the best model is the one Agresti calls \(L + S\). As we have seen, the formula for
this model, considered as a GLM or as a log-linear model is
count ~ lake*food + size*food + lake*gender*size.
These models having AIC within 10 of the minimum AIC are in Agresti’s notation \[\begin{gather*} L + S \\ G + L + S \\ L * S \\ G * L + S \\ G + L * S \end{gather*}\]
In the notation of R function glmbb these are the models
between
bout$little## count ~ lake * gender * size + food
## <environment: 0x5a55ca232b20>bout$big## count ~ lake * gender * size * food(inclusive) that have AIC no larger than min(AIC) + 10.
This is all very confusing (because). There are three sampling models, Poisson, multinomial, and product multinomial which
can all have the same maximum likelihood estimates of mean values and hence can all be fit by the same software,
have the same likelihood ratio test statistics and hence all have the same likelihood ratio test \(P\)-values (regardless of software used)
and have, for large sample sizes, nearly the same Rao (score) and Wald test statistics and \(P\)-values.
So we can ignore (mostly, except when we are using exact sampling distributions via the parametric bootstrap) the sampling distribution.
So there are three algorithms (something we have mentioned but not fussed about before)
R function glm does maximum likelihood on the exponential family canonical parameter scale
which it calls “coefficients” using the iteratively
reweighted least squares algorithm,
R function loglin in core R and R function
loglm in R recommended package MASS do maximum
likelihood on the exponential family mean
value parameter scale using the iterative
proportional fitting algorithm and may or may not use a formula
(loglm uses formulas, loglin just specifies marginals to match),
and
R function multinom in R recommended package
nnet uses neural
nets to fit multinomial response regression.
So there are three different algorithms and two different parameterizations being used but the same models being fit so all work just as well (except neural nets is sloppy so inaccurate).
So there are three different ways of specifying exactly the same models (or nine = three times three if you think about both sampling distributions and algorithms) that are seemingly unrelated. So it only makes sense if you take the sophisticated view that we have not fully covered yet.
For simplicity, we just covered one way of thinking about these models.
Later, we will learn abouti another way.
This topic is Section 8.2 in Agresti.
We are going to use R function polr in R package
MASS to do what it calls proportional odds logistic
regression (POLR). Because its definitions may differ from
Agresti’s, we will take its help as authoritative. In the “details”
section of the help for this function, it says
This model is what Agresti [they cite the second edition of our textbook, whereas we are using the third edition] calls a cumulative link model. The basic interpretation is as a coarsened version of a latent variable \(Y_i\) which has a logistic or normal or extreme-value or Cauchy distribution with scale parameter one and a linear model for the mean. The ordered factor which is observed is which bin \(Y_i\) falls into with breakpoints \[ \zeta_0 = -\infty < \zeta_1 < \cdots < \zeta_K = \infty \] This leads to the model \[ \mbox{logit} P(Y \le k | x) = \zeta_k - \eta \] with logit replaced by probit for a normal latent variable, and \(\eta\) being the linear predictor, a linear function of the explanatory variables (with no intercept). Note that it is quite common for other software to use the opposite sign for \(\eta\) (and hence the coefficients
beta).In the logistic case, the left-hand side of the last display is the log odds of category \(k\) or less, and since these are log odds which differ only by a constant for different \(k\), the odds are proportional. Hence the term proportional odds logistic regression.
Comparing with equation (8.5) in the third edition of Agresti (our
textbook) we see that R function polr does indeed use the
opposite sign for regression coefficients from what Agresti uses (and
apparently also from what the software he uses for his examples uses).
So we had better stick with what R does.
But we also see that other than a possible change of sign of
regression coefficients, R function polr agrees with the
model discussed in Section 8.2 of Agresti.
As an example, we do example 8.2.4 in Agresti.
# clean up R global environment
rm(list = ls())
data(table_8.5)
names(table_8.5)## [1] "race" "trauma" "happy"sapply(table_8.5, class)## race trauma happy
## "integer" "integer" "integer"lapply(table_8.5, unique)## $race
## [1] 0 1
##
## $trauma
## [1] 0 1 2 3 4 5
##
## $happy
## [1] 1 2 3It is annoying that Agresti has coded these data numerically when they should be categorical. Let’s fix that like we did for the alligator data.
foo <- transform(table_8.5,
race = factor(race, labels = c("White", "Black")),
happy = ordered(happy,
labels = c("Very.happy", "Pretty.happy", "Not.too.happy")))Now we are ready to invoke R function polr
pout <- polr(happy ~ race + trauma, data = foo)
summary(pout)##
## Re-fitting to get Hessian## Call:
## polr(formula = happy ~ race + trauma, data = foo)
##
## Coefficients:
## Value Std. Error t value
## raceBlack 2.0361 0.6859 2.968
## trauma 0.4056 0.1830 2.216
##
## Intercepts:
## Value Std. Error t value
## Very.happy|Pretty.happy -0.5181 0.3400 -1.5238
## Pretty.happy|Not.too.happy 3.4006 0.5680 5.9872
##
## Residual Deviance: 148.407
## AIC: 156.407We see that the “coefficients” agree with those in the SAS output
shown in Agresti Table 8.6 except for sign, which may be due to R
dropping a different dummy variable for race than SAS does
or due to R using a different sign convention (discussed in the quote
from the help page above). So we must be fitting the same model.
If we want to be sure what the signs of the coefficients mean, we can look at predictions.
predict(pout, type = "probs",
newdata = data.frame(race = c("White", "Black"),
trauma = rep(2, 2)))## Very.happy Pretty.happy Not.too.happy
## 1 0.20928237 0.7208972 0.06982045
## 2 0.03339528 0.6015130 0.36509170predict(pout, type = "probs",
newdata = data.frame(race = rep("White", 6), trauma = 0:5))## Very.happy Pretty.happy Not.too.happy
## 1 0.37329388 0.5944293 0.03227686
## 2 0.28420782 0.6681410 0.04765122
## 3 0.20928237 0.7208972 0.06982045
## 4 0.14997095 0.7488215 0.10120759
## 5 0.10523160 0.7502560 0.14451244
## 6 0.07269745 0.7251240 0.20217855One says keeping trauma fixed at 2, there are more
Very happy whites and more Not too happy
blacks. And the other says keeping race fixed at
"White", more trauma makes for fewer
Very happy and more Not too happy.
So in both cases we have positive “coefficient” means worse outcome in a problem like this where later in the order of the ordered categorical response means worse.