Hospital visits– Log-odds analysis

Load packages and the data

library(vcd)
library(vcdExtra)
library(ggplot2)

data(HospVisits, package="vcdExtra")
HospVisits

##             stay
## visit        2-9 10-19 20+
##   Regular     43    16   3
##   Infrequent   6    11  10
##   Never        9    18  16

Calculate log odds

vcd::lodds calculates log odds for comparisons of adjacent categories, for one response variable

lodds(HospVisits, response="stay") 

## log odds for visit by stay 
## 
##            visit
## stay        Regular Infrequent  Never
##   2-9:10-19   0.989    -0.6061 -0.693
##   10-19:20+   1.674     0.0953  0.118

Make it into a data.frame, including the ASE

hosp.lodds <- as.data.frame(lodds(HospVisits, response="stay"))

Fit models

Here, I use 1/ASE^2 as the weights, in a weighted least squares

mod.null  <- lm(logodds ~ -1, weights=1/ASE^2, data=hosp.lodds)     # null
mod.const <- lm(logodds ~ 1, weights=1/ASE^2, data=hosp.lodds)      # constant 
mod.unif  <- lm(logodds ~ visit, weights=1/ASE^2, data=hosp.lodds)  # uniform

Parallel model

A parallel log odds model has a main effect for stay

mod.par   <- lm(logodds ~ visit + stay, data=hosp.lodds)

Compare model fits

anova(mod.null, mod.const, mod.unif, mod.par)

## Analysis of Variance Table
## 
## Model 1: logodds ~ -1
## Model 2: logodds ~ 1
## Model 3: logodds ~ visit
## Model 4: logodds ~ visit + stay
##   Res.Df   RSS Df Sum of Sq    F  Pr(>F)    
## 1      6 22.95                              
## 2      5 19.90  1      3.05 1308 0.00076 ***
## 3      3  4.38  2     15.51 3323 0.00030 ***
## 4      2  0.00  1      4.38 1875 0.00053 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LRstats(mod.null, mod.const, mod.unif, mod.par)

## Likelihood summary table:
##             AIC   BIC LR Chisq Df Pr(>Chisq)    
## mod.null   16.8  16.5    22.95  6    0.00081 ***
## mod.const  17.9  17.5    19.90  5    0.00131 ** 
## mod.unif   12.8  12.0     4.38  3    0.22305    
## mod.par   -15.9 -17.0     0.00  2    0.99767    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fit models with visit as an ordinal factor

mod1a <- lm(logodds ~ as.numeric(visit), weights=1/ASE^2, data=hosp.lodds)
mod2a <- lm(logodds ~ as.numeric(visit) + stay, weights=1/ASE^2, data=hosp.lodds)
anova(mod.const, mod2a, mod.par)

## Analysis of Variance Table
## 
## Model 1: logodds ~ 1
## Model 2: logodds ~ as.numeric(visit) + stay
## Model 3: logodds ~ visit + stay
##   Res.Df   RSS Df Sum of Sq    F  Pr(>F)    
## 1      5 19.90                              
## 2      3  4.24  2     15.66 3353 0.00030 ***
## 3      2  0.00  1      4.24 1815 0.00055 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plotting

## Basic plot of points

gg <- ggplot(hosp.lodds, aes(x=visit, y=logodds, 
                             group=stay, color=stay)) + 
  geom_point(size=5) +
  ylab("log odds of shorter stay\n") + 
  xlab("Visit frequency") + theme_bw() + 
  theme(axis.title = element_text(size = rel(1.5))) + 
  theme(axis.text = element_text(size = rel(1.2))) +
  theme(legend.position = c(.85, .85), 
        legend.background = element_rect(colour = "black")) +
  theme(legend.title = element_text(size=rel(1.3))) +
  theme(legend.text = element_text(size=rel(1.2))) +
  scale_colour_manual(values=c("blue", "red"), 
                      name = "Length of stay")

gg + geom_line(size=1.2, linetype="dotted") 

showing model fits

grid <- hosp.lodds[,1:2]

Make a function to plots lines from a given model

gg_lines <- function(grid, mod, size=1.2, color=NULL, ...) {
  grid$logodds <- stats::predict(mod, grid)
  if(is.null(color)) geom_line(data=grid, size=size, ...) 
  else geom_line(data=grid, size=size, color=color, ...)
}

Add fitted lines to the plot

gg2 <-  
  gg + gg_lines(grid, mod.null, color="gray", size=1, linetype="dashed") +
  gg_lines(grid, mod.const, color=gray(.5), size=1) +
  gg_lines(grid, mod.unif, color="black", size=1) +
  gg_lines(grid, mod.par)

gg2 + annotate("text", label="null",      x=3.3, y=0, color="gray") +
  annotate("text", label="constant",  x=3.3, y=predict(mod.const, grid)[6], color=gray(.5)) + 
  annotate("text", label="uniform",   x=3.3, y=predict(mod.unif, grid)[6], color="black") +
  annotate("text", label="parallel",  x=3.3, y=hosp.lodds[5:6, 3], color=c("blue", "red")) +
  annotate("text", label="Log odds\nModel",     x=3.3, y=0.6, color="black", size=6)