I am trying to tie the odds ratio from a 2x2 cross classification table to the intercepts of a logistic regression on those 2 variables. I have a cross classification table that produces 2 odds ratios and the results of a logistic regression of PLACE3 ~ VIOL should produce intecepts should match the odds ratio of the contingency table. i.e. Odds ratio = exp(intercepts) BUT the POLR package is not producing the correct intercepts.
Here is the data. In the logistic regression PLACE3 is the outcome and VIOl is the independent variable. You can see the PLACE3 vs. VIOL contingency table below and the logistic regression of PLACE3 ~ VIOL. The odds ratios in the contingency table 1.79 and 3.1 are correct but the polr function seems off. Any thoughts on why exp(summary(m)$zeta) does not produce 1.79 and 3.1?
For reference this is from Lemeshow's Applied Logisitic Regression book page 274.
library(data.table)
aps <- fread('http://www.umass.edu/statdata/statdata/data/aps.dat')
colnames(aps) = c("ID","PLACE","PLACE3","AGE","RACE","GENDER","NEURO","EMOT","DANGER","ELOPE","LOS","BEHAV","CUSTD",
"VIOL")
head(aps)
Here is a cross classification table of PLACE3 vs. VIOl variables
table(aps$PLACE3,aps$VIOL)
0 1
0 80 179
1 26 104
2 15 104
using PLACE3 = 0 as the reference the 2 odds ratios from the contingency table are
(104*80)/(179*26) #1.79
(104*80)/(179*15) #3.10
These odds ratios should be the same as exponentiating the slope coefficients from a logistic model PLACE3 ~ VIOL which is below
aps$constant = rep(1,dim(aps)[1])
m <- polr(as.factor(PLACE3) ~ constant + as.factor(VIOL), data = aps, Hess=TRUE,model=TRUE,method = c("logistic"))
summary(m)
> summary(m)
Call:
polr(formula = as.factor(PLACE3) ~ constant + as.factor(VIOL),
data = aps, Hess = TRUE, model = TRUE, method = c("logistic"))
Coefficients:
Value Std. Error t value
as.factor(VIOL)1 0.8454 0.2112 4.003
Intercepts:
Value Std. Error t value
0|1 0.6869 0.1884 3.6464
1|2 1.8608 0.2032 9.1557
Residual Deviance: 1031.75
AIC: 1037.75
But you can see the exponentiation of the zeta vector is not 1.79 and 3.10
exp(summary(m)$zeta)
> exp(summary(m)$zeta)
0|1 1|2
1.987495 6.429049