0
$\begingroup$

I am trying to understand why for the following example I do not get exactly the coefficient estimates and furthermore why the standard errors for the coefficients are so large in the case of vglm.

Example setup:

library(reshape2)
library(VGAM)
library(nnet)

group <- c("ctl", "trt")
ldl.begin <- c("<=3.4","3.4-4.1","4.1-4.9",">4.9")
ldl.end <- c("<=3.4","3.4-4.1","4.1-4.9",">4.9")
count.ctl <- c(18,8,0,0,
           16,30,13,2,
           0,14,28,7,
           0,2,15,22)
count.trt <- c(21, 4, 2, 0, 
           17,25,6,0,
           11,35,36,6,
           1,5,14,12)

df.4 <- expand.grid(ldl.end = ldl.end, ldl.begin = ldl.begin, group = group)
df.4$count <- c(count.ctl, count.trt)
df.4$ldl.begin <- factor(df.4$ldl.begin, levels = c("<=3.4","3.4-4.1","4.1-4.9",">4.9"))
df.4$ldl.end <- factor(df.4$ldl.end, levels = c("<=3.4","3.4-4.1","4.1-4.9",">4.9"))
df.4$group <- factor(df.4$group, levels = c("ctl", "trt"))

df.5 <- dcast(df.4, ldl.begin + group ~ ldl.end, value.var = "count")
df.5$ldl.begin <- factor(df.5$ldl.begin, levels = c("<=3.4","3.4-4.1","4.1-4.9",">4.9"))
df.5$group <- factor(df.5$group, levels = c("ctl", "trt"))
names(df.5) <- c("ldl.begin", "group", "end.ldl.1", "end.ldl.2", "end.ldl.3", "end.ldl.4") 

model.3.1 <- multinom(ldl.end ~ group + ldl.begin, weights = count, df.4)
summary(model.3.1)

# The placement of end.ldl.1 makes it the baseline.

model.3.2 <- vglm(cbind(end.ldl.2,end.ldl.3,end.ldl.4,end.ldl.1)~ group + ldl.begin, data=df.5, family=multinomial)
summary(model.3.2)

A segment of the summary for multinom below:

Coefficients:
(Intercept)   grouptrt ldl.begin3.4-4.1 ldl.begin4.1-4.9 ldl.begin>4.9
3.4-4.1  -0.9131126 -0.5496199         1.697319         2.835446      3.248406
4.1-4.9  -2.4848193 -1.1989666         2.435842         5.103357      6.619781
>4.9    -12.0279955 -1.8157173         9.865304        13.366615     16.585710

and the summary from vglm:

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept):1        -0.9134     0.3641  -2.509 0.012107 *  
(Intercept):2        -2.4846     0.7395  -3.360 0.000780 ***
(Intercept):3       -19.5738  2355.2418  -0.008 0.993369    
grouptrt:1           -0.5497     0.3318  -1.657 0.097572 .  
grouptrt:2           -1.1990     0.3888  -3.084 0.002045 ** 
grouptrt:3           -1.8156     0.4976  -3.649 0.000264 ***
ldl.begin3.4-4.1:1    1.6978     0.4002   4.242 2.21e-05 ***
ldl.begin3.4-4.1:2    2.4355     0.7872   3.094 0.001975 ** 
ldl.begin3.4-4.1:3   17.4126  2355.2419   0.007 0.994101    
ldl.begin4.1-4.9:1    2.8360     0.4862   5.833 5.44e-09 ***
ldl.begin4.1-4.9:2    5.1033     0.8199   6.224 4.84e-10 ***
ldl.begin4.1-4.9:3   20.9129  2355.2418   0.009 0.992915    
ldl.begin>4.9:1       3.2537     1.1246   2.893 0.003815 ** 
ldl.begin>4.9:2       6.6247     1.2648   5.238 1.62e-07 ***
ldl.begin>4.9:3      24.1364  2355.2420   0.010 0.991823    

Note in particular the discrepency for the >4.9 estimates.

A comparison between the results for multinorm and vglm with another dataset (gator.txt from https://onlinecourses.science.psu.edu/stat504/sites/onlinecourses.science.psu.edu.stat504/files/lesson08/gator.txt) produces identical coefficient and standard errors.

$\endgroup$
  • $\begingroup$ Apparently, there is some sort of complete separation in the data, so that the coefficient estimates for the >4.9 group diverge to +/- infinity. The fact that the two packages stop at different finite values is just due to different numerical optimizations. Note also the huge standard errors for these coefficients. $\endgroup$ – Achim Zeileis Oct 5 '15 at 22:00
  • $\begingroup$ Thanks Achim. Yes, I did notice the standard errors. I get a similar result with stata (mlogit). I will have a look into separation and pseudo-separation. I think there are options for dealing with this in the context of logistic regression but not multinomial. $\endgroup$ – t-student Oct 6 '15 at 12:05
  • $\begingroup$ I would probably use an orderer logit model here. The proportional odds assumption seems to be reasonable for this data. You could probably even use a linear coding (or something similar) for ldl.begin rather than separate effects for each level. Both of these improve the AIC of the model. $\endgroup$ – Achim Zeileis Oct 6 '15 at 22:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.