# probability in mixed multinomial logistic regression

I am using MCMCglmm package in R for my multilevel multinomial logistic regression model. I have a level-1 binary outcome 'Sex', which was coded as 1,2, and a level-1 three category unordered multinomial outcome 'mathach', which was coded as 0,1,2. I also have level-1 continous predictor 'SES'. The ID variable is 'School'. The reproducible code using R dataset MathAchieve and MathAchSchool from nlme package is as below:

library(MCMCglmm)
library(nlme)
data(MathAchieve,package='nlme')
data(MathAchSchool,package='nlme')
dat=merge(MathAchSchool,MathAchieve,by='School')
dat$$mathach[dat$$MathAch<5]=0
dat$$mathach[dat$$MathAch>=5 & dat$$MathAch<15]=1 dat$$mathach[dat$$MathAch>15]=2 dat$$mathach=as.factor(dat\$mathach)
str(dat)

set.seed(9689724)
m1=MCMCglmm(mathach~SES,random=~School+SES,data=dat,rcov=~us(trait):units,family='categorical',verbose=F)
summary(m1)

Iterations = 3001:12991
Thinning interval  = 10
Sample size  = 1000

DIC: 9944.844

G-structure:  ~School

post.mean l-95% CI u-95% CI eff.samp
School    0.9069   0.6566    1.219    25.86

~SES

post.mean  l-95% CI u-95% CI eff.samp
SES 0.0007462 2.728e-06 0.004624    13.33

R-structure:  ~us(trait):units

post.mean l-95% CI u-95% CI eff.samp
traitmathach.1:traitmathach.1.units   36.6564 12.44909  63.0278    1.016
traitmathach.2:traitmathach.1.units    0.2105 -0.79940   1.1571    3.582
traitmathach.1:traitmathach.2.units    0.2105 -0.79940   1.1571    3.582
traitmathach.2:traitmathach.2.units    0.1275  0.07205   0.2024   11.354

Location effects: mathach ~ SES

post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)    1.1074   0.8905   1.3358    12.12 <0.001 ***
SES            1.0324   0.8863   1.1513    12.21 <0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


My questions are:

1. The regression coeffecient of SES in the multinomial logistic regression is 1.0324, but which category of 'mathach' does the coefficient refer to? how to interpret it?

2. If I want to make predictions with individual having SES=0.5 and come from School 1224, how do to do that?

The regression coeffecient of SES in the multinomial logistic regression is 1.0324, but which category of 'mathach' does the coefficient refer to? how to interpret it?

Both, because your model has SES as both a random and fixed effect.

If I want to make predictions with individual having SES=0.5 and come from School 1224, how do to do that?

You would call the predict() function, passing the fitted model object and a data frame with SES of 0.5 and SchoolID of 1224

• Thank you very much for the quick response. But mathach is a 3-category multinomial variable, would there be two equations with two different SES regression coefficients? – user11806155 Jun 28 '20 at 12:29
• I tried this, I get a strange error > predict(m1,data.frame(SES=0.5,School='1224'),marginal=NULL,type='response') Error in FUN(X[[i]], ...) : object 'mathach' not found – user11806155 Jun 28 '20 at 12:31
• @user11806155 That looks strange. You might want to ask on the r-sig-me mailing list because this site is for statistical questions, not for programming questions or questions about particular software. – Robert Long Jun 28 '20 at 13:11
• So maybe I can manually calculate the predictive probability myself. Statistically, what would be the correct way of the calculation using the value of SES and ID? – user11806155 Jun 28 '20 at 13:15
• As you said, the outcome is multinomial but there appears to be only one estimate, so there doesn't appear to be enough information to make a prediction. We would need to look into the workings of mcmcglmm. – Robert Long Jun 28 '20 at 13:22