Removing intercept from GLM for multiple factorial predictors only works for first factor in model I am running a binomial logistic regression with a logit link function in R. My response is factorial [0/1] and I have two multilevel factorial predictors - let's call them $a$ and $b$ where $a$ has 4 factor levels  $(a_1,a_2,a_3,a_4)$ and $b$ has 9 factor levels $(b_1,b_2,\dotsc,b_9)$. Therefore:
mod <- glm(y ~ a + b, family = binomial(logit), data=pretend)
summary(mod)

The model output would then show all the information about the model as well as the coefficients.
There is a factor level for both a and b missing (a1 and b1) from the summary output. I understand that it is constrained in the "intercept" of the model.
I have read that if I want to remove the intercept term and see the estimates for those factor levels I can just add -1 to the model formula, i.e.:
mod2 <- glm(y ~ a + b - 1, family=binomial(logit), data=pretend)
summary(mod2)

In the new model (mod2) the intercept term is then gone and variable a's factor-level a1 is given amongst the list of coefficients. But, variable b's factor-level b1 is still missing and given that there is no intercept term anymore, how can I interpret the odds-ratio for that factor level then?
Could someone please explain to me how to get the coefficient for b1 too and why this is happening?        
 A: @kjetil b halvorsen has done a good job outlining the main ideas here.  Let me add a couple supplementary points.  
With a categorical variable, suppressing the intercept results in level means coding, instead of the default reference level coding.  I explain this in greater detail here:  How can logistic regression have a factorial predictor and no intercept? 
You can use level means coding with multiple categorical variables, but in essence you have to fit the full interaction.  In your case, you only wanted to fit the additive model (y~a+b); that is what you cannot do, as previously explained.  
Should you be committed to using level means coding, the procedure is fairly straightforward.  You first create a new, single variable as the Cartesian product (the combinations) of all possible levels of your various categorical variables.  For example, in place of your original two categorical variables (a, with 4 levels, and b, with 9), you would have a single variable with 36 levels (a1b1, a1b2, a1b3, a1b4, a1b5, a1b6, a1b7, a1b8, a1b9, a2b1, ..., a4b9).  Then you fit your model using level means coding (i.e., suppressing the intercept) with the new variable:  
mod <- glm(y~0+ab, family=binomial(logit), data=pretend)
summary(mod)

Note again that this is equivalent to glm(y~a*b, ...); it is only that the output will be presented differently.  
