FEprovideR

Fixed effects logistic model with high-dimensional parameters

A stuctured profile likelihood algorithm for the logistic fixed effects model and an approximate expectation maximization (EM) algorithm for the logistic mixed effects model.

Installation

You can install the released version of FEprovideR from Github with:

install.packages("devtools") # you need devtools to install packages from Github
devtools::install_github("umich-biostatistics/FEprovideR")

You can install directly from CRAN with:

install.packages("FEprovideR")

Example

This tutorial simulates a data set to demonstrate the functions provided by FRprovideR.

# load the package
library(FEprovideR)

# other imports
library(Matrix)
library(poibin)
library(ggplot2)

To simulate a data set, use the following code chunk:

# Simulate a data set
m <- 500
prov.size <- pmax(round(rnorm(m, 50, 15)),11)
gamma <- rnorm(m, log(3/7), 0.4)
beta <- c(1,0.5,-1)
Y.char <- 'Y'
prov.char <- 'prov.ID'
Z.char <- paste0('z', 1:length(beta))
sim.fe.prov <- function(m, prov.size, gamma, beta, Y.char, Z.char, prov.char) {
  N <- sum(prov.size) # total number of discharges
  gamma.dis <- rep(gamma, times=prov.size)
  prov <- rep(1:m, times=prov.size) # provider IDs
  Z <- matrix(rnorm(N*length(beta)), ncol=length(beta))
  Y <- rbinom(N, 1, plogis(gamma.dis+Z%*%beta))
  data <- as.data.frame(cbind(Y, prov, Z))
  colnames(data) <- c(Y.char, prov.char, Z.char) 
  return(data)
}
data <- sim.fe.prov(m, prov.size, gamma, beta, Y.char, Z.char, prov.char)

This data is also available in the included data sets that come with the package. To use the included data, run:

data(hospital)            # raw data
data(hospital_prepared)   # processed data

Now, set relevant parameters and fit a model to the prepared data:

# a small positive number specifying stopping criterion of Newton-Raphson algorithm
tol <- 1e-5  
# Name input variables and other parameters
Y.char <- 'Y'
prov.char <- 'prov.ID'
Z.char <- paste0('z', 1:3)
data(hospital_prepared) # build in data set
fe.ls <- fe.prov(hospital_prepared, Y.char, Z.char, prov.char, tol) # model fitting

Conduct hypothesis tests on the estimated standardized readmission ratios (SSRs):

# hypothesis testing
null <- "median"
n <- 10000
alpha <- 0.05
score.fe <- test.fe.prov(hospital_prepared, fe.ls, Y.char, Z.char, prov.char, test="score", null, alpha)
exact.pb <- test.fe.prov(hospital_prepared, fe.ls, Y.char, Z.char, prov.char, test="exact.poisbinom", null, alpha)
exact.bs <- test.fe.prov(hospital_prepared, fe.ls, Y.char, Z.char, prov.char, test="exact.bootstrap", null, alpha, n)
exact.binom <- test.fe.prov(hospital_prepared, fe.ls, Y.char, Z.char, prov.char, test="exact.binom", null="median", alpha)

Compute confidence intervals for the estimated SSRs:

# confidence intervals
confint.df <- confint.fe.prov(fe.ls, parm = "all", level = 0.88, hospital_prepared, Y.char, Z.char, prov.char)
confint.df <- confint.fe.prov(fe.ls, parm = "all", level = 0.90, hospital_prepared, Y.char, Z.char, prov.char)
confint.df <- confint.fe.prov(fe.ls, level = 0.90, data = hospital_prepared, Y.char = Y.char, Z.char = Z.char, prov.char = prov.char)

# CIs for a subset of providers
confint.df3 <- confint.fe.prov(fe.ls, hospital_prepared, Y.char, Z.char, prov.char, parm=c(1,2,50), level=0.95) 

Funnel plots for SRRs (Standardized readmission ratios)

# format input data for funnel plot
input.dis <- data.frame(ID=hospital_prepared[hospital_prepared$included==1, prov.char],
                        prob=fe.ls$Exp)
input.prov <- data.frame(SRR=fe.ls$df.prov$SRR, flag=score.fe$flag)

Score test based funnel plot:

target <- c(1)
alphas <- c(0.1, 0.5, 0.01)
input.prov <- data.frame(SRR=fe.ls$df.prov$SRR, flag=score.fe$flag)
funnel.SRR(input.dis, input.prov, target, alphas, type="FE.score")