knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "man/figures/README-", out.width = "100%" )
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.