A stuctured profile likelihood algorithm for the logistic fixed effects model and an approximate expectation maximization (EM) algorithm for the logistic mixed effects model.
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")
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)
# 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")