The required dataset is already uploaded to the project directory in the format of a CSV (comma separated value)
R can also import another document format such as xls, dta (from Stata), or dat (from SPSS)
You can use this code to import the dataset and label the variable
# Importing dataset from csv filehospital <-read_csv("hospital.csv") %>%as_tibble()# Labeling the datasethospital <- hospital %>%set_variable_labels(id_rs ="Hospital ID",gp_FTE ="Number of General Practitioners",spec_FTE ="Number of Medical Specialists",nurse_total ="Number of Nurses",other_prof ="Number of other employees",beds ="Number of beds",outpatients ="Number of outpatient services",bed_days ="Total Bed Days",admissions1 ="Number of admissions",bed_occ ="Bed Occupancy",throughput ="Throughput",alos ="Average Length of Stay",totalcost_op ="Total cost of outpatient services",total_cost_bd ="Total cost of beddays",totalcost_ip ="Total cost of admission/inpatient",publichospital ="Public Hospital",poorins ="Proportion with Social Health Insurance",classCD ="Class C/D",non_JavaBali ="Non Java-Bali") %>%set_value_labels(publichospital =c("Public"=1, "Private"=0),classCD =c("Class C/D"=1, "Class A/B"=0),non_JavaBali =c("Non Java-Bali"=1, "Java-Bali"=0))
Part 2. Descriptive analysis
Summary statistics table
In this section, we use the gtsummary package that can create a publication-ready table.
As an example, we create a table describing several variables namely: number of GPs, number of specialists, and ownership.
hospital %>%select(gp_FTE,spec_FTE,publichospital) %>%# Select variableto_factor() %>%# Set categorical variable to "factor"tbl_summary(digits =all_continuous() ~2, # Number of digits for cont. variablesstatistic =all_continuous() ~"{mean} ({sd})") %>%# What type of summary statisticsas_flex_table() # Presented as flex_table, easy to copy to Word processing software
Characteristic
N = 198
Number of General Practitioners
15.41 (9.95)
Number of Medical Specialists
19.86 (20.52)
Public Hospital
Public
120 (61%)
Private
78 (39%)
Mean (SD); n (%)
Exercise
Create a descriptive table for number of outpatient services (outpatients), number of admissions (admissions1), and location of hospital (non_JavaBali)
Some basic plot like histogram, scatter plot, and box-plot will be used as an example.
# Load packagelibrary(ggplot2)
Histogram of number of GPs using base function
# Using base functionhist(hospital$gp_FTE)
Histogram of number of GPs using ggplot
# Using ggplothospital %>%ggplot(aes(x=gp_FTE)) +geom_histogram(color="black", fill="darkseagreen")
Scatter plot of number of medical specialist vs number of outpatients using base function
# Using base functionplot(hospital$spec_FTE,hospital$outpatients)
Scatter plot of number of medical specialist vs number of outpatients using ggplot
# Using ggplothospital %>%ggplot(aes(x=spec_FTE,y=outpatients)) +geom_point(color="darkseagreen") +labs(x="Number of specialists",y="Number of outpatient services")
Using ggplot + color code based on hospital class [EXTRA]
Box Plot of total cost of outpatient based on hospital class using base function
# Using base functionboxplot(hospital$classCD,hospital$totalcost_op)
Box Plot of total cost of outpatient based on hospital class using ggplot
# Using ggplothospital %>%to_factor() %>%ggplot(aes(x=classCD,y=totalcost_op)) +geom_boxplot(color="darkseagreen") +labs(x="Hospital Class",y="Total cost of outpatient")
Exercise
Histogram of number of specialists
Scatter plot of number of nurses (nurse_total) vs number of admissions (admissions1)
Boxplot of total cost of inpatients based on hospital ownership (publichospital)
Part 3. Concept of Efficiency
Part 4. Efficiency Analysis
Cost analysis
Simple calculation and creating new variables using mutate function from the dplyr package.
# Calculating unit cost per servicehospital <- hospital %>%mutate(unitcost_ip = totalcost_ip/admissions1,unitcost_bd = total_cost_bd/bed_days)
After creating new variables, we may want to label them.
## Adding variable label to the newly created variableshospital <- hospital %>%set_variable_labels(unitcost_ip ="Inpatient Unit Cost",unitcost_bd ="Bed days Unit Cost")
Again, we will create a table using the tbl_summary function from the gtsummary package.
# Create summary table of the unit costhospital %>%select(unitcost_ip, unitcost_bd) %>%tbl_summary(digits =all_continuous() ~2,statistic =all_continuous() ~"{mean} ({sd})") %>%as_flex_table()
This is the result.
Characteristic
N = 198
Inpatient Unit Cost
297.58 (194.62)
Bed days Unit Cost
82.35 (54.61)
Mean (SD)
Exercise
Please calculate the unit cost of outpatient service (totalcost_op divided by outpatients) and present the table
Pabon-Lasso Model/Diagram
It is basically a scatter plot comparing bed turnover ratio (throughput) and bed occupancy ratio (bed-occ).
Hospitals are then divided into 4 sectors based on the mean of bed turnover and bed occupancy ratios.
First, let’s calculate them
# Calculating the sample mean of bed turnover and bed occupancy ratiomean_btr <-mean(hospital$throughput)mean_bor <-mean(hospital$bed_occ)
This is a simple PL diagram using ggplot for all hospital
Based on the Pabon-Lasso model, we will define whether hospitals are considered efficient or not.
The efficient hospitals are those who have both higher than average of bed turnover and bed occupancy rate.
# Construct a binary indicator of hospitals efficiency based on the PL modelhospital <- hospital %>%mutate(eff_pl =if_else ( (throughput > mean_btr) & (bed_occ > mean_bor),1,0))
This code is used to give label to the new variable eff_pl.
# Labelling variable using the labelled packagevar_label(hospital$eff_pl) <-"Efficiency (Pabon-Lasso)"val_labels(hospital$eff_pl) <-c("Efficient"=1,"Not efficient"=0)
Exercise
Create a simple Pabon-Lasso diagram for private hospitals.
# Combining results from DEA dea_id <-data.frame(cbind(id_rs=hospital$id_rs, EI=c(dea_input$eff), EI_bc=c(dea_input$eff.bc), EI_CI_low=c(dea_input$conf.int[,1]), EI_CI_hi=c(dea_input$conf.int[,2]),EO=c(dea_output$eff), EO_bc=c(dea_output$eff.bc), EO_CI_low=c(dea_output$conf.int[,1]), EO_CI_hi=c(dea_output$conf.int[,2])))# Merge with hospital datahospital_dea <-merge( hospital, dea_id, by="id_rs", all =TRUE)
Part 5. Bivariate analysis
Introduction
In efficiency studies, bivariate analysis can be used to explore factors associated with some efficiency measures or scores.
For example, we can test whether public hospitals and private hospitals have different unit cost for their services. In this case, cost variable might be converted to a log format.
Converting variables to log form
First, we may want to convert the cost variable to a log form for further analyses
## Adding variable label to the newly created variableshospital <- hospital %>%set_variable_labels(lnuc_ip ="Log (Inpatient Unit Cost)",lnuc_bd ="Log (Bed days Unit Cost)")
Visualise the cost difference
A boxplot can be used to visualise the difference in unit cost
# A boxplot comparing inpatient unit cost between public and private hospitalhospital %>%to_factor() %>%ggplot(aes(x=publichospital, y=unitcost_ip)) +xlab("Hospital type") +ylab("Inpatient unit cost") +geom_boxplot()
We can also present the cost variable in a log form
# Cost in a log formhospital %>%to_factor() %>%ggplot(aes(x=publichospital, y=lnuc_ip)) +geom_boxplot() +labs(x="Hospital type",y="Inpatient unit cost")
Statistical test
Is the difference in means statistically significant?
In this case, we can use a two-sample t-test and produce the table at the same time using tbl_summary function.
Do a statistical test comparing the difference between inpatient unit cost of hospitals located in Java/Bali and in Non-Java/Bali
Part 6. Multivariable analysis
Linear regressions
When the outcome (y) variable is continues and we want to explore multiple factors (x) that are associated with the outcome variable, we might use an Ordinary Least Square (OLS) or linear regression.
For example, we want to analyse the relationship between some contextual factors and the number of GPs.
# Estimating a linear modelols_gp <- hospital %>%lm(gp_FTE ~ publichospital + classCD + non_JavaBali + poorins, data=.)
Logistic regression
If the outcome is binary (yes or no, 1 or 0, efficient or not), another approach that can be used is a logistic regression.
# Linear regression table using modelsummary# Labelling the covariates or regressorsols_gp_label <- ols_gpnames(ols_gp_label$coefficients) <-c("Constant", "Class C/D", "Public", "Non Java-Bali", "Proportion with SHI")# Produce tablemodelsummary(list("Number of GP"= ols_gp_label), # List of the model(s) and their labelfmt=3, # Number of digitsestimate ="{estimate}{stars}", # Format of coefficientsgof_omit =".*") # Goodness of fit statistics are excluded
Linear regression table using modelsummary
Number of GP
Constant
17.747***
(1.731)
Class C/D
7.180***
(1.330)
Public
-4.245**
(1.517)
Non Java-Bali
-0.860
(1.330)
Proportion with SHI
-14.336***
(4.286)
Logistic regression table using modelsummary (code)
# Labelling the covariates or regressorslogit_pl_label <- logit_plnames(logit_pl_label$coefficients) <-c("Constant", "Class C/D", "Public", "Non Java-Bali", "Proportion with SHI")# Produce tablemodelsummary(list("Efficient Hospital (PL)"= logit_pl_label), # List of the model(s) and their labelfmt=NULL, # Number of digitsestimate ="{round(exp(estimate), 3)}{stars}", # Format of coefficientsstatistic ="({round(exp(estimate) * std.error, 3)})",gof_omit =".*") # Goodness of fit statistics are excluded
Logistic regression table using modelsummary (result)
Efficient Hospital (PL)
Constant
0.57
(0.242)
Class C/D
1.688
(0.564)
Public
0.731
(0.272)
Non Java-Bali
0.422**
(0.138)
Proportion with SHI
8.957*
(9.299)
Exercise
Run a linear regression with number of specialists as the outcome variable (y) and hospital ownership, location, class, and district-level proportion of population with social health insurance as the explanatory variables (x).
Present the result/coefficient using either base function or other functions.
