R: How to fit a partial proportional odds model with vglm

Question

I am trying to fit an ordinal regression model in R with the vglm function from the VGAM package. My dependent variable is Fnreports which is a factor with the levels 0 to 5, the independent variable RA01 indicates the 4 treatment groups from my between-group experiment (also a factor). I assume that the proportional odds assumption does not hold (I fitted the model with the polr function (package MASS) and performed the Brant test (package Brant), which the model failed) and I like to test this again with the Likelihood ratio test (lrtest), and then use a partial proportional-odds model instead. When fitting the proportional odd model I get the following results:

> pom <- vglm(FNreports ~ RA01, data=ds_c, 
            family = cumulative(parallel = TRUE, reverse = TRUE))
> summary(pom)

Call:
vglm(formula = FNreports ~ RA01, family = cumulative(parallel = TRUE, 
    reverse = TRUE), data = ds_c)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)    
(Intercept):1   0.6408     0.3572   1.794   0.0728 .  
(Intercept):2  -0.2704     0.3512  -0.770   0.4414    
(Intercept):3  -0.8494     0.3608  -2.354   0.0186 *  
(Intercept):4  -1.7106     0.3934  -4.348 1.37e-05 ***
(Intercept):5  -3.0970     0.5249  -5.900 3.64e-09 ***
RA012           0.5918     0.4957   1.194   0.2325    
RA013           0.7710     0.5388   1.431   0.1524    
RA014           0.1583     0.4851   0.326   0.7441    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of linear predictors:  5 

Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]), logitlink(P[Y>=4]), logitlink(P[Y>=5]), 
logitlink(P[Y>=6])

Residual deviance: 324.8854 on 472 degrees of freedom

Log-likelihood: -162.4427 on 472 degrees of freedom

Number of Fisher scoring iterations: 5 

No Hauck-Donner effect found in any of the estimates


Exponentiated coefficients:
   RA012    RA013    RA014 
1.807249 2.161890 1.171545 

However, when fitting the model without the proportional odds assumption I get the following error:

npom <- vglm(FNreports ~ RA01, data=ds_c, 
              family = cumulative(parallel = FALSE, reverse = TRUE))
Fehler in tapplymat1(ccump, "diff") : 
  NA/NaN/Inf in externem Funktionsaufruf (arg 1)
Zusätzlich: Warnmeldungen:
1: In Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals = residuals,  :
  fitted values close to 0 or 1
2: In Deviance.categorical.data.vgam(mu = mu, y = y, w = w, residuals = residuals,  :
  fitted values close to 0 or 1

I assume it is because the RA01 variable is a factor, because when I change it to an integer I get some results, but I don't think that comparing the two models now is appropriate.

> ds_c$RA01 <- as.integer(ds_c$RA01)
> npom1 <- vglm(FNreports ~ RA01, data=ds_c, 
              family = cumulative(parallel = FALSE, reverse = TRUE))
> summary(npom1)

Call:
vglm(formula = FNreports ~ RA01, family = cumulative(parallel = FALSE, 
    reverse = TRUE), data = ds_c)

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)  
(Intercept):1  0.91241    0.52283   1.745   0.0810 .
(Intercept):2 -0.23255    0.47121  -0.494   0.6216  
(Intercept):3 -1.06322    0.50066  -2.124   0.0337 *
(Intercept):4 -0.99884    0.55081  -1.813   0.0698 .
(Intercept):5 -2.23570    0.90520  -2.470   0.0135 *
RA01:1         0.01192    0.19374   0.062   0.9509  
RA01:2         0.11362    0.17441   0.651   0.5147  
RA01:3         0.22309    0.18047   1.236   0.2164  
RA01:4        -0.13274    0.21115  -0.629   0.5296  
RA01:5        -0.19197    0.36362  -0.528   0.5975  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of linear predictors:  5 

Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]), 
     logitlink(P[Y>=4]), logitlink(P[Y>=5]), 
     logitlink(P[Y>=6])

Residual deviance: 322.756 on 470 degrees of freedom

Log-likelihood: -161.378 on 470 degrees of freedom

Number of Fisher scoring iterations: 8 

No Hauck-Donner effect found in any of the estimates


Exponentiated coefficients:
   RA01:1    RA01:2    RA01:3    RA01:4    RA01:5 
1.0119889 1.1203286 1.2499292 0.8756900 0.8253353

> lrtest(npom1, pom)
Likelihood ratio test

Model 1: FNreports ~ RA01
Model 2: FNreports ~ RA01
  #Df  LogLik Df  Chisq Pr(>Chisq)
1 470 -161.38                     
2 472 -162.44  2 2.1294     0.3448

Moreover, if it is legit to perform a model fitting with the an integer as independent variable instead of a factor, I don't know how to interpret the output, especially the Exponentiated coefficients.

Answer 1

I have had convergence problems with partial proportional odds models with vglm. Occasionally vgam will work. Since you have only one independent variable the partial PO model is equivalent to a multinomial logistic model. Try getting the log-likelihood from that.

If the Brent test is the score test for PO it can be anticonservative, i.e., make you worry too much about PO.

If the dependent variable has one level that has only a few observations you might temporarily re-run your two models after collapsing that category with one of its neighboring categories.

The partial PO model has too many parameters to be efficient, so you might try a constrained partial PO model. VGAM will fit that. For example you can fit the full covariate effect but restrict departures from PO to be linear in the dependent variable's numeric code.