I ran an ordered logit difference-in-difference estimation and got stuck in terms of interpretation. The dependent variable takes on values 1 2 and 3 and I am trying to obtain coefficients for the variables for each of the three categories as I see in papers but I am not sure how to go about it.
Solved – ordered logit-difference-in-difference
difference-in-differenceeconometricsinterpretationordered-logit
Related Solutions
To manually verify the predictions derived from using polr()
from package MASS
, assume a situation with a categorical dependent variable $Y$ with ordered categories $1, \ldots, g, \ldots, k$ and predictors $X_{1}, \ldots, X_{j}, \ldots, X_{p}$. polr()
assumes the proportional odds model
$$ \text{logit}(p(Y \leqslant g)) = \ln \frac{p(Y \leqslant g)}{p(Y > g)} = \beta_{0_g} - (\beta_{1} X_{1} + \dots + \beta_{p} X_{p}) $$
For possible choices implemented in other functions, see this answer. The logistic function is the inverse of the logit-function, so the predicted probabilities $\hat{p}(Y \leqslant g)$ are
$$ \hat{p}(Y \leqslant g) = \frac{e^{\hat{\beta}_{0_{g}} - (\hat{\beta}_{1} X_{1} + \dots + \hat{\beta}_{p} X_{p})}}{1 + e^{\hat{\beta}_{0_{g}} - (\hat{\beta}_{1} X_{1} + \dots + \hat{\beta}_{p} X_{p})}} $$
The predicted category probabilities are $\hat{P}(Y=g) = \hat{P}(Y \leq g) - \hat{P}(Y \leq g-1)$. Here is a reproducible example in R with two predictors $X_{1}, X_{2}$. For an ordinal $Y$ variable, I cut a simulated continuous variable into 4 categories.
set.seed(1.234)
N <- 100 # number of observations
X1 <- rnorm(N, 5, 7) # predictor 1
X2 <- rnorm(N, 0, 8) # predictor 2
Ycont <- 0.5*X1 - 0.3*X2 + 10 + rnorm(N, 0, 6) # continuous dependent variable
Yord <- cut(Ycont, breaks=quantile(Ycont), include.lowest=TRUE,
labels=c("--", "-", "+", "++"), ordered=TRUE) # ordered factor
Now fit the proportional odds model using polr()
and get the matrix of predicted category probabilities using predict(polr(), type="probs")
.
> library(MASS) # for polr()
> polrFit <- polr(Yord ~ X1 + X2) # ordinal regression fit
> Phat <- predict(polrFit, type="probs") # predicted category probabilities
> head(Phat, n=3)
-- - + ++
1 0.2088456 0.3134391 0.2976183 0.1800969
2 0.1967331 0.3068310 0.3050066 0.1914293
3 0.1938263 0.3051134 0.3067515 0.1943088
To manually verify these results, we need to extract the parameter estimates, from these calculate the predicted logits, from these logits calculate the predicted probabilities $\hat{p}(Y \leqslant g)$, and then bind the predicted category probabilities to a matrix.
ce <- polrFit$coefficients # coefficients b1, b2
ic <- polrFit$zeta # intercepts b0.1, b0.2, b0.3
logit1 <- ic[1] - (ce[1]*X1 + ce[2]*X2)
logit2 <- ic[2] - (ce[1]*X1 + ce[2]*X2)
logit3 <- ic[3] - (ce[1]*X1 + ce[2]*X2)
pLeq1 <- 1 / (1 + exp(-logit1)) # p(Y <= 1)
pLeq2 <- 1 / (1 + exp(-logit2)) # p(Y <= 2)
pLeq3 <- 1 / (1 + exp(-logit3)) # p(Y <= 3)
pMat <- cbind(p1=pLeq1, p2=pLeq2-pLeq1, p3=pLeq3-pLeq2, p4=1-pLeq3) # matrix p(Y = g)
Compare to the result from polr()
.
> all.equal(pMat, Phat, check.attributes=FALSE)
[1] TRUE
For the predicted categories, predict(polr(), type="class")
just picks - for each observation - the category with the highest probability.
> categHat <- levels(Yord)[max.col(Phat)] # category with highest probability
> head(categHat)
[1] "-" "-" "+" "++" "+" "--"
Compare to result from polr()
.
> facHat <- predict(polrFit, type="class") # predicted categories
> head(facHat)
[1] - - + ++ + --
Levels: -- - + ++
> all.equal(factor(categHat), facHat, check.attributes=FALSE) # manual verification
[1] TRUE
That is not a robustness check because the ordinary linear model is guaranteed not to fit. It will yield probabilities estimates outside $[0,1]$. A better approach to checking the assumptions of an ordinal regression model are:
- First, relax the assumptions allowing for nonlinear effects using regression splines
- Then, check the equal slopes (parallelism) assumption.
For the logistic ordinal model (proportional odds model) the equal slopes (proportional odds) assumption can be checked in several ways, including:
- Fit a series of binary logistic models for different cutoffs of $Y$ and plot the regression coefficients vs. cutoff and check for constancy
- Construct partial residual plots for all cutoffs of $Y$. One often has to collapse infrequent $Y$ categories to carry this out.
- Plot the logit of the empirical cumulative distribution of $Y$ stratified by very important predictors and check for parallelism.
Best Answer
If you have specified the difference in differences as the interaction of your time and group variable then you need to be aware that interactions in non-linear models (like ordered logit) have a different interpretation. This relates to a similar question here which asked how difference in differences works in the context of logit models. In this reply there are alternative models and links to articles that explain the intuition and interpretation of difference in differences in the context of nonlinear models. Have a look and then see whether ordered logit is actually what you want to do. Particularly relevant for your case would be Karaca-Mandic et al (2012). Note though that what you are trying to do is not very well explored in the econometric literature which is mostly for the fact that difference in differences in nonlinear models is far from trivial (for an overview of the problem see Athey and Imbens, 2006).