0
$\begingroup$

I'm working with the pneumo dataset from the VGAM package. Here is the data:

glimpse(pneumo)
## Rows: 8
## Columns: 4
## $ exposure.time <dbl> 5.8, 15.0, 21.5, 27.5, 33.5, 39.5, 46.0, 51.5
## $ normal        <dbl> 98, 51, 34, 35, 32, 23, 12, 4
## $ mild          <dbl> 0, 2, 6, 5, 10, 7, 6, 2
## $ severe        <dbl> 0, 1, 3, 8, 9, 8, 10, 5

I want to model the severity categories (mild, normal, severe) using exposure time as a predictor. First, I use pivot_longer to reformat the data and make sure the outcome variable is ordered properly. I don't think there are any issues here:

pneumo_long <- pneumo %>% pivot_longer(!exposure.time,
                                       names_to = "severity",
                                       values_to = "count")

pneumo_long <- pneumo_long %>% mutate(severity=factor(severity,
                                                      levels=c("mild",
                                                               "normal",
                                                               "severe"),
                                    ordered=TRUE))

Next I fit a proportional odds model:

mod21_vglm <- vglm(formula = severity ~ exposure.time, 
                   family = cumulative(parallel=TRUE), weights = count, 
                   data=pneumo_long[pneumo_long$count != 0,])

summary(mod21_vglm)
## 
## Call:
## vglm(formula = severity ~ exposure.time, family = cumulative(parallel = TRUE), 
##     data = pneumo_long[pneumo_long$count != 0, ], weights = count)
## 
## Coefficients: 
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1 -1.844887   0.256770  -7.185 6.72e-13 ***
## (Intercept):2  2.367550   0.277947   8.518  < 2e-16 ***
## exposure.time -0.015662   0.008913  -1.757   0.0789 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2])
## 
## Residual deviance: 502.1551 on 41 degrees of freedom
## 
## Log-likelihood: -251.0775 on 41 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
## exposure.time 
##     0.9844597

Looks fine. But then I try a non-proportional odds model and get a warning about the Hauck-Donner effect being detected in the intercepts:

mod21_vglm_2 <- vglm(formula = severity ~ exposure.time, 
                   family = cumulative(parallel=FALSE), weights = count, 
                   data=pneumo_long[pneumo_long$count != 0,])

summary(mod21_vglm_2)
## 
## Call:
## vglm(formula = severity ~ exposure.time, family = cumulative(parallel = FALSE), 
##     data = pneumo_long[pneumo_long$count != 0, ], weights = count)
## 
## Coefficients: 
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1   -3.83843    0.47160  -8.139 3.98e-16 ***
## (Intercept):2    4.56682    0.53583   8.523  < 2e-16 ***
## exposure.time:1  0.05893    0.01352   4.358 1.31e-05 ***
## exposure.time:2 -0.08578    0.01458  -5.883 4.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2])
## 
## Residual deviance: 421.9943 on 40 degrees of freedom
## 
## Log-likelihood: -210.9971 on 40 degrees of freedom
## 
## Number of Fisher scoring iterations: 7 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1', '(Intercept):2'
## 
## 
## Exponentiated coefficients:
## exposure.time:1 exposure.time:2 
##       1.0607046       0.9177929

I can't wrap my head around why this is happening. Surely separation can't be an issue here? Or am I missing something?

A couple additional notes:

Likelihood ratio test strongly prefers the non-proportional odds model:

dev_diff_21 <- deviance(mod21_vglm) - deviance(mod21_vglm_2)

df_21 <- mod21_vglm@df.residual - mod21_vglm_2@df.residual

p.value_21 <- 1 - pchisq(q = dev_diff_21, df = df_21)

data.frame(dev_diff_21, df_21, p.value_21)
##   dev_diff_21 df_21 p.value_21
## 1    80.16081     1          0

But the proportional odds model gives more intuitive results:

exposure <- data.frame("exposure.time" =
                         c(5,10,15,20,25))

pred_21_p <- predict(object = mod21_vglm, newdata = exposure, 
                     type="response")

pred_21_nop <- predict(object = mod21_vglm_2, newdata = exposure, 
                     type="response")

pred_21_p
##         mild    normal     severe
## 1 0.12750534 0.7804766 0.09201811
## 2 0.11904450 0.7821820 0.09877347
## 3 0.11107362 0.7829595 0.10596689
## 4 0.10357371 0.7828082 0.11361813
## 5 0.09652521 0.7817282 0.12174660
pred_21_nop
##         mild    normal     severe
## 1 0.02809222 0.9562021 0.01570567
## 2 0.03735920 0.9387244 0.02391636
## 3 0.04952735 0.9142113 0.03626134
## 4 0.06538954 0.8799886 0.05462189
## 5 0.08587295 0.8326341 0.08149299

My hunch is that the proportional odds model is better here and that as exposure time increases, the probability for mild should go down and the probability for severe should go up. But what's going on with the non-proportional odds model?

New contributor
mbrookz is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
1
  • 1
    $\begingroup$ Two things occur to me. The proportional odds model imposes a constraint which the non-p model does not. The fact that your zeroes are in one extreme of the covariate may be relevant. However I am sure others more expert than me will be able to give a definitive answer. $\endgroup$ – mdewey Jun 8 at 15:50

Your Answer

mbrookz is a new contributor. Be nice, and check out our Code of Conduct.

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.