Interpreting categorical variables in regression

Question

When running a regression with a categorical independent variable, we get results for each level of the variable except for the base, which we can choose.

Now I've always had a hard time on how to interpret these results.

Say we have a study of aneurysm locations. They can be located in, say, 10 different areas.

We want to see if smokers develop aneurysms in other areas than non-smokers. We have our dependent variable (smoker, no/yes) and our independent variable of location with 10 levels.

If we run the regression we might get a significant hit on 3 locations. But this is compared to the base location which let's say is level "5".

So yes, smokers are significantly more likely to get aneurysms in location 1, 2 and 3 compared to location 5. But this doesn't answer my research question of "which areas are smokers more likely to develop aneurysms in?".

What I would like to do is to maybe make a "dummy level" to my categorical variable in which half the patients have that location and half do not and then use that as base to see if ANY of the 10 true levels have a significantly higher risk for smokers. I don't want to compare the levels to each other and I want every level included and not for one of them to be used as a base.

Is there some way to do this or am I using the wrong model to answer my research question?

I assume splitting the categorical variable into 10 dummy variables is probably not so smart.

Answer 1

Given the development of aneurysm in one area mutually excludes developments of aneurysm in another it is possible to let the dependent variable be $Y_{ij}=1$ if individual $i$ has developed aneurysm in the area $j$ and $Y_{ij}=0$ otherwise. For each individual you then have a vector of measurements $Y_i =(Y_{i1},...,Y_{iJ})$ as a multivariate output. The vector follows a multinomial distribution the probabilities of which are simply estimated by sample relative frecquencies.

You can fit a multinomial logit model and do testing comparing the relative frecquencies of smoker with non-smokers.

$$H_o : \pi(smoker) = \pi(non-smoker)$$

where $\pi(smoker)$ is a vector of probabilities $$\pi_1,...,\pi_J$$ with $\pi_j$ being the probability that a smoker develops aneurysm in area $j$.

Here is some R code simulating such data and running a test

    library(data.table)
library(mnlogit)

true_difference <- 0
N <- 1000
J <- 10
smoker <- as.numeric(runif(N)<0.5)

# select area effects (parameters to generate probabilities)
AE <- 0.5*rnorm(J)
AE[1] <- AE[1] + true_difference
AE_smoker <- rep(0,J)
AE_smoker[3] <- AE_smoker[3]+true_difference
AE_smoker[4] <- AE_smoker[4]+true_difference
AE_smoker[5] <- AE_smoker[5]+true_difference


AE
AE_smoker


p <- exp(AE)
p <- p/sum(p)
p_smoker <- exp(AE + AE_smoker)
p_smoker <- p_smoker/sum(p_smoker)
p
p_smoker


area <- rep(NA,N)
for (i in 1:N)
    {
        i_probs <- smoker[i]*p_smoker + (1-smoker[i])*p
        area[i] <- sample(1:10,1,prob=i_probs)
    } 


dt <- data.table(id=1:N,area=area,smoker=smoker)
agg_dt <- dt[, .(count = .N), by = .(area,smoker)]
setkey(agg_dt,area,smoker)
agg_dt[,.(difference=diff(count)),by=area]


# Transform data to long format
choice <- rep(0,J*N)
ii <- 1
for (i in 1:N)
    {
        choice[ii:(ii+J-1)][dt$area[i]] <- 1
        ii <- ii + J    
    }

mydata <- data.frame(id=sort(rep(dt$id,J)),
choice=as.logical(choice),
smoker=rep(dt$smoker,each=J),
alt=rep(1:J,N))

# Specify model
# https://arxiv.org/pdf/1404.3177.pdf
fm <- formula(choice ~  1 | smoker | 1)
fit <- mnlogit(fm,data=mydata,choiceVar="alt")


fm.c <- formula(choice ~  1 | 1 | 1)
fit.c <- mnlogit(fm.c,data=mydata,choiceVar="alt")
lrtest(fit, fit.c)
waldtest(fit, fit.c)
# Scoretest do not seem to be working ...
scoretest(fit, fit.c)

Answer 2

You're not using the wrong model; you're just not looking at an output that tells you what you want know....

My suspicion is that you are using R, just because this question comes up a lot with newer R users, because so often examples in R rely on the output from the summary function, which provides model coefficients, and usually t tests or other relevant tests for the coefficients.

Instead, what it sounds like what you want, is 1) an anova-like table of effects, and b) (pairwise) comparisons among group means or estimated marginal means. b) will tell you if group 1 is different than group 2, and so on.

In R, a) is shown with e.g. library(car); Anova(model), where 'model' is your model object. b) is shown using the emmeans package. Something like library(emmeans); marginal = emmeans(model, ~ Group); pairs(marginal, adjust="tukey").

Answer 3

A way I use to help me interpret is to create a table with the two variables I'm looking at - the dependent and the independent I'm focusing at that moment. In a regression, the other independent variables may change this relation a bit, but usually I can get a good idea regarding what's going on.

Following your example, I could have a table with conditional frequencies by row like: $$ \begin{array}{c|cc} Area \backslash Smoke & Yes\ (1) & No\ (0)\\ \hline 1 & 60\% & 40\% \\ 2 & 50\% & 50\% \\ 3 & 30\% & 70\% \\ \ldots & \ldots & \ldots \\ 10 & 55\% & 45% \\ \end{array} $$ Of course that we could also use conditional frequencies by column.

By looking at the table, we could say that the smokers were more likely to have an aneurysm in areas $1$ and $10$ and that the non-smokers were more likely to have an aneurysm in area $3$.

---

Edit:

If you have a logistic regression, you could look at the odds ratio ($OR$) or at the coefficients ($B$) to reach a similar conclusion from the reverse path. By the way, $OR = e^B$.

For instance, if the 2nd category is the reference category for the area variable (I admit that I chose the $50\% - 50\%$ category to make my life easier), and non-smoking is the reference category for the smoker variable, we could observe that $OR_1=1.4$ for area $1$, $OR_3=0.44$ for area $3$, and $OR_{10}=1.1$ for area $10$ (I've made up some numbers). This would mean that the smoking individuals, when compared to the non-smokers, were more likely to have aneurysms in areas $1$ and $10$ and less likely to have aneurysms in area $3$. By looking at the $p$-values for each coefficient, you could also say in which cases the $OR$ is significantly different from $1$.

We compare the value of the $OR$ with $1$ because $OR_j=1$ would mean that the odds of finding a smoker in category $j$ would be the same as the odds of finding a smoker in category $2$, i.e., being a smoker would have the same odds of having an aneurysm in area $j$ or in area $2$.