Problem Description

The objective of this problem is to determine how sedation and pain management at IU Health Ball Memorial Hospital (BMH) compares to other IU Health hospitals in respect to the usage of the agent ‘propofol’ for sedation in the ICU. At IU BMH, propofol continuous infusion rates are significantly higher than other hospitals within the IU Health system. A detailed exploratory data analysis will be conducted and general linear models will be fit to solve this problem. In addition, variable selection using deviance tables and testing for evidence of interaction will occur.

List of column names with labels for the propofol dataset

#load datasets and libraries
library(rsample)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.4
## ✔ ggplot2   3.4.4     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(naniar)
library(corrplot)
## corrplot 0.92 loaded
Propdata <- readRDS("JBlain_Rfile copy.rds")
head(Propdata)
## # A tibble: 6 × 44
##   Org_Name   Age   Sex   Height Weight IBW   COVID IVDA_Hx Home_Opiod Home_Benzo
##   <labelled> <lab> <lab> <labe> <labe> <lab> <lab> <label> <labelled> <labelled>
## 1 IU Health… 43    0     154.9  117.7  47.7… 0     0       1          0         
## 2 IU Health… 65    0     149.9   73.9  43.2… 0     0       0          0         
## 3 IU Health… 64    0     180.3  128.3  70.7… 1     0       0          0         
## 4 IU Health… 67    0     157.5  155.3  50.1… 0     0       0          0         
## 5 IU Health… 63    0     167.6   79.4  59.2… 1     0       0          0         
## 6 IU Health… 37    0     172.7  164.0  63.8… 0     0       1          1         
## # ℹ 34 more variables: Opiate_Tol <labelled>, Triglyc_Max <labelled>,
## #   Creatinine <labelled>, CrCl <labelled>, Map_Range <labelled>,
## #   HR_Range_low <labelled>, Total_days_prop <labelled>,
## #   Dur_prop_minutes <labelled>, Per_prop_over50 <labelled>,
## #   Max_Dose_Prop <labelled>, Min_Rass <labelled>, Max_Rass <labelled>,
## #   Dex_ord <labelled>, Prop_dose_when_Dex <labelled>,
## #   Dex_rate_when_prop <labelled>, Prop_dec_afterDex <labelled>, …

Creating Analysis Dataset

IBW, COVID, IVDA_Hx, Triglyc_Max, Creatinine, CrCl, Map_Range, and HR_Range_low will all be dropped from the analysis dataset since they are not relevant to the analysis being completed.

#creating analysis dataset
df_eda <- select(Propdata, -c(IBW, COVID, IVDA_Hx,Triglyc_Max,Creatinine,CrCl,Map_Range,HR_Range_low))
df_eda[df_eda=="x"] <- NA
df_eda[df_eda=="X"] <- NA

x/X that represent NA/missing values is replaced with NA to be able to visualize and replace NA values.

#visualizing and appropriately replacing missing variables
vis_miss(df_eda[,1:12])

df_eda$Height <- ifelse(is.na(df_eda$Height), "Missing", df_eda$Height)
vis_miss(df_eda[,13:25])

df_eda$Min_Rass <- ifelse(is.na(df_eda$Min_Rass), "NotApp", df_eda$Min_Rass)
df_eda$Max_Rass <- ifelse(is.na(df_eda$Max_Rass), "NotApp", df_eda$Max_Rass)
df_eda$Dex_ord <- ifelse(is.na(df_eda$Dex_ord), "NotApp", df_eda$Dex_ord)
df_eda$Prop_dose_when_Dex <- ifelse(is.na(df_eda$Prop_dose_when_Dex), "NotApp", df_eda$Prop_dose_when_Dex)
df_eda$Dex_rate_when_prop <- ifelse(is.na(df_eda$Dex_rate_when_prop), "NotApp", df_eda$Dex_rate_when_prop)
df_eda$Prop_dec_afterDex <- ifelse(is.na(df_eda$Prop_dec_afterDex), "NotApp", df_eda$Prop_dec_afterDex)
df_eda$Fent_before_prop <- ifelse(is.na(df_eda$Fent_before_prop), "NotApp", df_eda$Fent_before_prop)
df_eda$Prop_dose_when_fent <- ifelse(is.na(df_eda$Prop_dose_when_fent), "NotApp", df_eda$Prop_dose_when_fent)
df_eda$Fent_max_when_prop <- ifelse(is.na(df_eda$Fent_max_when_prop), "NotApp", df_eda$Fent_max_when_prop)
df_eda$Prop_dec_after_fent <- ifelse(is.na(df_eda$Prop_dec_after_fent), "NotApp", df_eda$Prop_dec_after_fent)
df_eda$Numb_fent_bolus <- ifelse(is.na(df_eda$Numb_fent_bolus), "NotApp", df_eda$Numb_fent_bolus)
df_eda$Mode_CPOT <- ifelse(is.na(df_eda$Mode_CPOT), "NotApp", df_eda$Mode_CPOT)
df_eda$Was_ket_ord <- ifelse(is.na(df_eda$Was_ket_ord), "NotApp", df_eda$Was_ket_ord)
vis_miss(df_eda[,26:36])

df_eda$Prop_dose_when_ket <- ifelse(is.na(df_eda$Prop_dose_when_ket), "NotApp", df_eda$Prop_dose_when_ket)
df_eda$Ket_max_when_prop <- ifelse(is.na(df_eda$Ket_max_when_prop), "NotApp", df_eda$Ket_max_when_prop)
df_eda$Was_prop_dec <- ifelse(is.na(df_eda$Was_prop_dec), "NotApp", df_eda$Was_prop_dec)
df_eda$Con_Prop_NMBA <- ifelse(is.na(df_eda$Con_Prop_NMBA), "NotApp", df_eda$Con_Prop_NMBA)
vis_miss(df_eda)

View(df_eda)

Exploratory Data Analysis

Visualizing categorical variables

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=Org_Name))

Most patients are from the IU Health Ball Memorial Hospital. This is helpful since this is where we are concentrating our analysis.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Sex), fill=factor(Sex)))

There are more males than females.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Home_Opiod), fill=factor(Home_Opiod)))

Most patients do not have home opiod usage.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Home_Benzo), fill=factor(Home_Benzo)))

Most patients do not have home benzo usage.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Opiate_Tol), fill=factor(Opiate_Tol)))

Most patients have no tolerance for Opiate.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Min_Rass), fill=factor(Min_Rass)))

Min RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients minimum RASS when propofol rate >50 mcg/kg/min was that they were unarousable. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Max_Rass), fill=factor(Max_Rass)))

Max RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients maximum RASS when propofol rate >50 mcg/kg/min was that they would pull at or attempt to remove tubes alongside generally aggressive behavior. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Dex_ord), fill=factor(Dex_ord)))

Dexmedetomidine was ordered before propofol rate reached >50 mcg/kg/min at any point. Most patients were not ordered dexmedetomidine.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Prop_dec_afterDex), fill=factor(Prop_dec_afterDex)))

Propofol was decreased to <50 mcg/kg/min after dexmedetomidine was added (did not go back above 50). Since most patients were not ordered dexmedetomidine, this is not applicable to most patients.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Fent_before_prop), fill=factor(Fent_before_prop)))

Fentanyl was ordered before propofol rate was >50 mcg/kg/min. Most patients did have fentanyl ordered for them before the propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Prop_dec_after_fent), fill=factor(Prop_dec_after_fent)))

Propofol was decreased to <50 mcg/kg/min after fentanyl was added. This was not applicable to most patients.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Mode_CPOT), fill=factor(Mode_CPOT)))

Mode CPOT when propofol rate >50 mcg/kg/min . “The CPOT is a behavioural assessment pain scale for patients unable to verbalise pain.” (2) Most patients had a mode CPOT of 0, which means that overall, the patient was relaxed, had no muscle tension, and had an absence of movement. Link to learn more about the scale: https://ccs-sth.org/resources/Documents/Sedation%20Analgesia%20Delirium%20in%20CC/CCS%20STH%20Guideline%20template%20CPOT.pdf

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Was_ket_ord), fill=factor(Was_ket_ord)))

Most patients did not have ketamine prescribed to them before propofol rate was >50 mcg/kg/min and out of the ones who did, most did not have it ordered.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Was_prop_dec), fill=factor(Was_prop_dec)))

Since most patients did not have ketamine prescribed to them, it was not applicable to most patients if propofol decreased to <50 mcg/kg/min after ketamine was added or the dose was increased.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Rec_pain_med), fill=factor(Rec_pain_med)))

Most patients did not receive any scheduled pain medications when propofol >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(On_benzo), fill=factor(On_benzo)))

Most patients were not on scheduled benzo.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Vas), fill=factor(Con_Prop_Vas)))

More patients did not receive concurrent propofol and vasopressor when propofol >50 mcg/kg/min than those who did.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Dex), fill=factor(Con_Prop_Dex)))

More patients did not receive concurrent propofol and dexmedetomidine than those who did.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Fen), fill=factor(Con_Prop_Fen)))

More patients did receive concurrent propofol and fentanyl infusion than those who did not.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Ket), fill=factor(Con_Prop_Ket)))

Most patients did not receive concurrent propofol and ketamine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_Ben), fill=factor(Con_Prop_Ben)))

Most patients did not receive concurrent propofol and benzodiazepine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=factor(Con_Prop_NMBA), fill=factor(Con_Prop_NMBA)))

More patients did not receive concurrent propofol and NMBA infusion than those who did.

Visualizing covariation between two categorical variables

Creating initialisms of Org_Name values to make bar graphs easier to interpret.

df_eda$Org_Name[df_eda$Org_Name=="IU Health Arnett Hospital"] <- 'IU HAH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Ball Memorial Hospital"] <- 'IU HBMH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Bloomington Hospital"] <- 'IU HBloomH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health Methodist Hospital"] <- 'IU HMH'
df_eda$Org_Name[df_eda$Org_Name=="IU Health University Hospital"] <- 'IU HUH'

Initialisms for hospital names:

  • IU Health Arnett Hospital <- ‘IU HAH’
  • IU Health Ball Memorial Hospital <- ‘IU HBMH’
  • IU Health Bloomington Hospital <- ‘IU HBloomH’
  • IU Health Methodist Hospital <- ‘IU HMH’
  • IU Health University Hospital <- ‘IU HUH’
ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Opiod)), position = "dodge")

The organization with the most patients that have home opiod usage is IU HBMH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Benzo)), position = "dodge")

The organization that has the most patients with home benzo usage is IU HBloomH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Opiate_Tol)), position = "dodge")

IU HBMH has the most patients with high tolerance of opiates.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_afterDex)), position = "dodge")

Only IU HBMH and IU HBloomH had propofol decreases to <50 mcg/kg/min after dexmedetomidine was added.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Fent_before_prop)), position = "dodge")

IU HBMH had the highest number of fentanyl orders before propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_after_fent)), position = "dodge")

Only IU HBloomH and IU HMH had propofol decreases to <50 mcg/kg/min after fentanyl was added.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_ket_ord)), position = "dodge")

IU HMH had the highest number of ketamine orders before propofol rate was >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_prop_dec)), position = "dodge")

Propofol was only decreased to <50 mcg/kg/min after ketamine was added or dose increased at IU HBMH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Rec_pain_med)), position = "dodge")

Only patients at IU HMH and IU HUH received scheduled pain medications when propofol >50.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Vas)), position = "dodge")

Patients received concurrent propofol and vasopressor when propofol >50 mcg/kg/min at IU HBMH, IU HMH, and IU HBloomH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Dex)), position = "dodge")

Only IU HUH did not have patients who received concurrent propofol and dexmedetomidine.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Fen)), position = "dodge")

Concurrent propofol and fentanyl infusion is a popular choice at all organizations except IU HUH.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ket)), position = "dodge")

On the contrary, concurrent propofol and ketamine infusion was not a popular choice at the organizations.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ben)), position = "dodge")

IU HMH had the most patients that received concurrent propofol and benzodiazepine infusion.

ggplot(data=df_eda)+
  geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_NMBA)), position = "dodge")

IU HBMH had the most patients that received concurrent propofol and NMBA infusion.

Visualizing continuous variables

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Age), binwidth = 3)

Most patients are between 50 and 70 years old.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Height)), binwidth = 10)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Most patients are around 180cm in height (roughly 5 ft 11in).

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Weight), binwidth = 5)

Weights are pretty evenly distributed with a few obvious outliers.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Total_days_prop), binwidth = 1)

Most patients were on propofol for 15 days or less, with a few outliers towards 20 days and 30 days.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Dur_prop_minutes), binwidth = 2000)

Most patients were on propofol for less than 10000 minutes (which is roughly less than a week) with a few higher outliers.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Per_prop_over50), binwidth = 5)

The % of time of propofol above 50 mcg/kg/min is relatively spread across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = Max_Dose_Prop), binwidth = 3)

The max dose of propofol in mcg/kg/min greatly varies across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_Dex)), binwidth = 3)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).

This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when dexmedetomidine was added varied across patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Dex_rate_when_prop)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the dexmedetomidine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was usually around 0.5 mcg/kg/hr.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_fent)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 59 rows containing non-finite values (`stat_bin()`).

This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was varied.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Fent_max_when_prop)), binwidth = 5)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 22 rows containing non-finite values (`stat_bin()`).

The fentanyl max rate in mcg/hr when propofol >50 mcg/kg/min was between 100-150 mcg/hr for most patients.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Numb_fent_bolus)), binwidth = 5)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 13 rows containing non-finite values (`stat_bin()`).

Most patients did not receive fentanyl bolus when propofol rate >50 mcg/kg/min.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_ket)), binwidth = 3)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 65 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was usually 50 mcg/kg/min or above.

ggplot(data=df_eda)+
  geom_histogram(mapping = aes(x = as.numeric(Ket_max_when_prop)), binwidth = 0.25)
## Warning in FUN(X[[i]], ...): NAs introduced by coercion

## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 61 rows containing non-finite values (`stat_bin()`).

This was not applicable for most patients, but for the patients where it was applicable, the ketamine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was spread between 0.0 to 2.0 mcg/kg/hr.

Visualizing covariation between a categorical variable and a continuous variable

The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. Based off evaluating the EDA results, the selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA. Below are some additional visualizations showing covariation between per_prop_over50 and the selected covariates.

ggplot(data = df_eda, mapping = aes(x = factor(Org_Name), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Home_Opiod), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Home_Benzo), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Opiate_Tol), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Min_Rass), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Max_Rass), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Dex_ord), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Prop_dec_afterDex), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Fent_before_prop), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Mode_CPOT), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Was_ket_ord), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Was_prop_dec), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Rec_pain_med), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(On_benzo), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Vas), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Dex), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Fen), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ket), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ben), y = Per_prop_over50)) + geom_boxplot()

ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_NMBA), y = Per_prop_over50)) + geom_boxplot()

Fitting GLM by selecting a quantitative response variable column and a set of covariates

The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. The selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA.

glm_mod1 <- glm(Per_prop_over50~ Home_Opiod+Opiate_Tol+Dex_ord+Prop_dose_when_Dex+Fent_before_prop+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+On_benzo+Con_Prop_Vas+Con_Prop_Dex+Con_Prop_Ket+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod1)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Home_Opiod + Opiate_Tol + Dex_ord + 
##     Prop_dose_when_Dex + Fent_before_prop + Prop_dec_after_fent + 
##     Mode_CPOT + Was_ket_ord + Rec_pain_med + On_benzo + Con_Prop_Vas + 
##     Con_Prop_Dex + Con_Prop_Ket + Con_Prop_Ben + Con_Prop_NMBA, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (3 not defined because of singularities)
##                                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                           16.9747    25.8444   0.657   0.5163  
## Home_Opiod                             9.5381    11.4020   0.837   0.4095  
## Opiate_Tol                             2.6929     6.9633   0.387   0.7017  
## Dex_ord1                              18.2757    41.8090   0.437   0.6652  
## Dex_ordNotApp                         25.9547    41.0530   0.632   0.5320  
## Prop_dose_when_Dex10                  37.5382    34.6597   1.083   0.2874  
## Prop_dose_when_Dex100                 28.5396    23.8664   1.196   0.2411  
## Prop_dose_when_Dex21.240400000000001  29.9235    29.2526   1.023   0.3145  
## Prop_dose_when_Dex40                 -10.7162    23.9682  -0.447   0.6580  
## Prop_dose_when_Dex50                  13.2324    24.4940   0.540   0.5930  
## Prop_dose_when_Dex50.0002              8.1402    23.1175   0.352   0.7272  
## Prop_dose_when_Dex51.2821             62.7275    60.2609   1.041   0.3062  
## Prop_dose_when_Dex52                  31.8747    37.0702   0.860   0.3967  
## Prop_dose_when_Dex53.864199999999997  44.2100    30.4603   1.451   0.1570  
## Prop_dose_when_Dex54.146299999999997  34.0348    38.0389   0.895   0.3780  
## Prop_dose_when_Dex59.991              23.1032    24.6688   0.937   0.3565  
## Prop_dose_when_Dex60                  25.8142    29.7748   0.867   0.3928  
## Prop_dose_when_Dex65                   9.1500    23.8664   0.383   0.7041  
## Prop_dose_when_Dex70                   0.8011    21.3461   0.038   0.9703  
## Prop_dose_when_Dex80                  12.2866    35.3059   0.348   0.7303  
## Prop_dose_when_Dex89.826099999999997  33.7013    28.4573   1.184   0.2456  
## Prop_dose_when_DexNotApp             -15.2214    44.8005  -0.340   0.7364  
## Fent_before_prop1                      4.9861    17.8279   0.280   0.7816  
## Fent_before_propNotApp                20.8589    18.1684   1.148   0.2600  
## Prop_dec_after_fent1                 -24.2186    16.3259  -1.483   0.1484  
## Prop_dec_after_fentNotApp             -7.1677    21.9951  -0.326   0.7468  
## Mode_CPOT1                             4.9364    24.4242   0.202   0.8412  
## Mode_CPOT2                            11.2787    12.1800   0.926   0.3618  
## Mode_CPOT3                            15.2220    11.0428   1.378   0.1783  
## Mode_CPOT4                             4.7574    12.4450   0.382   0.7050  
## Mode_CPOT5                            23.8243     9.7429   2.445   0.0206 *
## Mode_CPOT6                            37.0329    13.7906   2.685   0.0117 *
## Mode_CPOT7                                 NA         NA      NA       NA  
## Mode_CPOT8                           -11.6850    26.7368  -0.437   0.6652  
## Mode_CPOTNotApp                       -4.4039    17.8279  -0.247   0.8066  
## Was_ket_ord1                         -11.4571    24.0958  -0.475   0.6379  
## Was_ket_ordNotApp                    -22.9554    20.2429  -1.134   0.2658  
## Rec_pain_med                         -15.3959    13.6534  -1.128   0.2684  
## On_benzo                              -8.9283    18.1684  -0.491   0.6267  
## Con_Prop_Vas                           3.4379     8.1570   0.421   0.6764  
## Con_Prop_Dex                               NA         NA      NA       NA  
## Con_Prop_Ket                               NA         NA      NA       NA  
## Con_Prop_Ben                         -18.0266     9.6409  -1.870   0.0713 .
## Con_Prop_NMBA1                        11.0285     6.6217   1.666   0.1062  
## Con_Prop_NMBANotApp                    1.0988    23.5232   0.047   0.9631  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 287.3556)
## 
##     Null deviance: 31538.8  on 71  degrees of freedom
## Residual deviance:  8620.7  on 30  degrees of freedom
## AIC: 634.87
## 
## Number of Fisher Scoring iterations: 2

Creating new model with the removal of home_opiod, opiate_tol, dex_ord, fent_before_prop, on_benzo, con_prop_vas, con_prop_dex, and con_prop_ket due to insignificance.

glm_mod2 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod2)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT + Was_ket_ord + Rec_pain_med + Con_Prop_Ben + Con_Prop_NMBA, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                            6.8742    18.8962   0.364  0.71803   
## Prop_dose_when_Dex10                  38.0836    31.1457   1.223  0.22895   
## Prop_dose_when_Dex100                 29.5253    22.4311   1.316  0.19597   
## Prop_dose_when_Dex21.240400000000001  38.0810    23.1727   1.643  0.10856   
## Prop_dose_when_Dex40                   9.5802    22.4311   0.427  0.67172   
## Prop_dose_when_Dex50                   3.0937    22.2819   0.139  0.89031   
## Prop_dose_when_Dex50.0002             14.8410    22.2819   0.666  0.50940   
## Prop_dose_when_Dex51.2821             60.4206    36.5719   1.652  0.10675   
## Prop_dose_when_Dex52                  52.4912    31.9293   1.644  0.10843   
## Prop_dose_when_Dex53.864199999999997  44.0642    25.5473   1.725  0.09269 . 
## Prop_dose_when_Dex54.146299999999997  59.0718    26.6983   2.213  0.03301 * 
## Prop_dose_when_Dex59.991              31.5426    21.8944   1.441  0.15787   
## Prop_dose_when_Dex60                  37.9108    23.8074   1.592  0.11958   
## Prop_dose_when_Dex65                  10.1357    22.4311   0.452  0.65394   
## Prop_dose_when_Dex70                  16.2055    19.0705   0.850  0.40077   
## Prop_dose_when_Dex80                  32.2903    29.7630   1.085  0.28479   
## Prop_dose_when_Dex89.826099999999997  24.7197    21.8944   1.129  0.26596   
## Prop_dose_when_DexNotApp              18.5013    14.0309   1.319  0.19519   
## Prop_dec_after_fent1                 -25.0874    14.5919  -1.719  0.09370 . 
## Prop_dec_after_fentNotApp             -4.8870     9.2117  -0.531  0.59884   
## Mode_CPOT1                            -4.3380    14.7124  -0.295  0.76971   
## Mode_CPOT2                             4.0649    11.3549   0.358  0.72234   
## Mode_CPOT3                            10.2141     9.2816   1.100  0.27805   
## Mode_CPOT4                            13.1764    10.2097   1.291  0.20465   
## Mode_CPOT5                            28.3823     8.7954   3.227  0.00258 **
## Mode_CPOT6                            32.5280    13.0564   2.491  0.01721 * 
## Mode_CPOT7                                 NA         NA      NA       NA   
## Mode_CPOT8                            -3.2120    23.6690  -0.136  0.89277   
## Mode_CPOTNotApp                      -11.1861    17.5033  -0.639  0.52660   
## Was_ket_ord1                           6.4647    19.0272   0.340  0.73591   
## Was_ket_ordNotApp                     -8.4424    13.8822  -0.608  0.54671   
## Rec_pain_med                          -6.3054    12.6832  -0.497  0.62195   
## Con_Prop_Ben                         -13.3830     7.9786  -1.677  0.10168   
## Con_Prop_NMBA1                         8.7512     5.3445   1.637  0.10980   
## Con_Prop_NMBANotApp                    0.4982    22.2819   0.022  0.98228   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 291.0054)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 11058  on 38  degrees of freedom
## AIC: 636.79
## 
## Number of Fisher Scoring iterations: 2

Creating a new model with the removal of was_ket_ord and rec_pain_med due to insignificance.

glm_mod3 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod3)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT + Con_Prop_Ben + Con_Prop_NMBA, family = "gaussian", 
##     data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                           -1.2829    13.6982  -0.094  0.92584   
## Prop_dose_when_Dex10                  36.7878    30.3507   1.212  0.23242   
## Prop_dose_when_Dex100                 30.0312    22.1141   1.358  0.18189   
## Prop_dose_when_Dex21.240400000000001  35.3591    22.7713   1.553  0.12816   
## Prop_dose_when_Dex40                  10.0860    22.1141   0.456  0.65073   
## Prop_dose_when_Dex50                   3.2226    21.9997   0.146  0.88426   
## Prop_dose_when_Dex50.0002             14.7122    21.9997   0.669  0.50741   
## Prop_dose_when_Dex51.2821             61.0174    27.5644   2.214  0.03248 * 
## Prop_dose_when_Dex52                  47.0792    28.5646   1.648  0.10696   
## Prop_dose_when_Dex53.864199999999997  58.8424    21.9997   2.675  0.01070 * 
## Prop_dose_when_Dex54.146299999999997  68.0200    22.1141   3.076  0.00373 **
## Prop_dose_when_Dex59.991              32.4626    21.5769   1.505  0.14012   
## Prop_dose_when_Dex60                  37.5762    23.4872   1.600  0.11731   
## Prop_dose_when_Dex65                  10.6415    22.1141   0.481  0.63293   
## Prop_dose_when_Dex70                  19.7994    17.9888   1.101  0.27747   
## Prop_dose_when_Dex80                  41.8048    25.7058   1.626  0.11155   
## Prop_dose_when_Dex89.826099999999997  25.6397    21.5769   1.188  0.24155   
## Prop_dose_when_DexNotApp              18.9386    13.8551   1.367  0.17910   
## Prop_dec_after_fent1                 -24.1454    14.0111  -1.723  0.09237 . 
## Prop_dec_after_fentNotApp             -5.6781     8.9085  -0.637  0.52742   
## Mode_CPOT1                            -1.6622    12.6487  -0.131  0.89609   
## Mode_CPOT2                             1.9816    11.0306   0.180  0.85832   
## Mode_CPOT3                             9.3170     9.0612   1.028  0.30987   
## Mode_CPOT4                            12.7623    10.0047   1.276  0.20927   
## Mode_CPOT5                            29.2227     8.4547   3.456  0.00129 **
## Mode_CPOT6                            31.9617    12.8696   2.484  0.01719 * 
## Mode_CPOT7                                 NA         NA      NA       NA   
## Mode_CPOT8                             2.7059    18.6796   0.145  0.88553   
## Mode_CPOTNotApp                      -11.1176    17.2299  -0.645  0.52236   
## Con_Prop_Ben                         -10.7900     7.6380  -1.413  0.16530   
## Con_Prop_NMBA1                         9.3859     5.2143   1.800  0.07922 . 
## Con_Prop_NMBANotApp                    0.6271    21.9997   0.029  0.97740   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 283.9154)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 11641  on 41  degrees of freedom
## AIC: 634.49
## 
## Number of Fisher Scoring iterations: 2

Creating a new model with the removal of con_prop_ben and con_prop_NMBA due to insignificance.

glm_mod4 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT ,data=df_eda,family="gaussian")
summary(glm_mod4)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent + 
##     Mode_CPOT, family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                                       Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                            2.87253   11.91249   0.241  0.81057   
## Prop_dose_when_Dex10                  34.25348   30.27421   1.131  0.26400   
## Prop_dose_when_Dex100                 26.32687   21.79393   1.208  0.23350   
## Prop_dose_when_Dex21.240400000000001  30.25075   21.79393   1.388  0.17211   
## Prop_dose_when_Dex40                   6.38170   21.79393   0.293  0.77104   
## Prop_dose_when_Dex50                   0.06822   20.27029   0.003  0.99733   
## Prop_dose_when_Dex50.0002             20.39379   21.79393   0.936  0.35451   
## Prop_dose_when_Dex51.2821             59.39517   27.55450   2.156  0.03663 * 
## Prop_dose_when_Dex52                  56.24703   27.55450   2.041  0.04725 * 
## Prop_dose_when_Dex53.864199999999997  64.52406   21.79393   2.961  0.00493 **
## Prop_dose_when_Dex54.146299999999997  64.31573   21.79393   2.951  0.00506 **
## Prop_dose_when_Dex59.991              29.75937   20.70055   1.438  0.15762   
## Prop_dose_when_Dex60                  34.72992   23.33553   1.488  0.14381   
## Prop_dose_when_Dex65                   6.93725   21.79393   0.318  0.75175   
## Prop_dose_when_Dex70                  20.99152   17.45961   1.202  0.23568   
## Prop_dose_when_Dex80                  30.79671   24.76247   1.244  0.22020   
## Prop_dose_when_Dex89.826099999999997  22.93645   20.70055   1.108  0.27388   
## Prop_dose_when_DexNotApp              17.31647   13.22757   1.309  0.19729   
## Prop_dec_after_fent1                 -23.68429   13.55431  -1.747  0.08755 . 
## Prop_dec_after_fentNotApp             -6.12923    8.72825  -0.702  0.48624   
## Mode_CPOT1                            -4.44635   12.49053  -0.356  0.72356   
## Mode_CPOT2                            -1.03662   10.35968  -0.100  0.92075   
## Mode_CPOT3                             9.23034    9.10922   1.013  0.31646   
## Mode_CPOT4                            11.76122    9.63274   1.221  0.22860   
## Mode_CPOT5                            28.36474    8.50821   3.334  0.00175 **
## Mode_CPOT6                            39.26555   12.49053   3.144  0.00299 **
## Mode_CPOT7                                  NA         NA      NA       NA   
## Mode_CPOT8                           -10.16625   17.38056  -0.585  0.56159   
## Mode_CPOTNotApp                      -13.19973   17.38056  -0.759  0.45163   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 292.141)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 12854  on 44  degrees of freedom
## AIC: 635.63
## 
## Number of Fisher Scoring iterations: 2

Performing variable selection using deviance tables

anova(glm_mod4,test="Chisq")
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Per_prop_over50
## 
## Terms added sequentially (first to last)
## 
## 
##                     Df Deviance Resid. Df Resid. Dev Pr(>Chi)   
## NULL                                   71      31539            
## Prop_dose_when_Dex  17  10746.1        54      20793 0.003599 **
## Prop_dec_after_fent  2   1076.0        52      19717 0.158555   
## Mode_CPOT            8   6862.4        44      12854 0.002789 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The final model for the % of time of propofol above 50 mcg/kg/min (Per_prop_over50) based on the deviance table can be written as follows:

Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT

glm_mod5 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Mode_CPOT,data=df_eda,family="gaussian")
summary(glm_mod5)
## 
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Mode_CPOT, 
##     family = "gaussian", data = df_eda)
## 
## Coefficients: (1 not defined because of singularities)
##                                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                            4.0461    11.7992   0.343  0.73323   
## Prop_dose_when_Dex10                   8.7601    27.1855   0.322  0.74874   
## Prop_dose_when_Dex100                 19.0241    20.9654   0.907  0.36892   
## Prop_dose_when_Dex21.240400000000001  22.9480    20.9654   1.095  0.27941   
## Prop_dose_when_Dex40                  -0.9211    20.9654  -0.044  0.96515   
## Prop_dose_when_Dex50                   3.7196    20.2538   0.184  0.85510   
## Prop_dose_when_Dex50.0002             13.0910    20.9654   0.624  0.53544   
## Prop_dose_when_Dex51.2821             51.4568    27.1855   1.893  0.06469 . 
## Prop_dose_when_Dex52                  48.3087    27.1855   1.777  0.08218 . 
## Prop_dose_when_Dex53.864199999999997  57.2213    20.9654   2.729  0.00896 **
## Prop_dose_when_Dex54.146299999999997  57.0129    20.9654   2.719  0.00920 **
## Prop_dose_when_Dex59.991              27.2815    20.2538   1.347  0.18458   
## Prop_dose_when_Dex60                  26.1812    22.4645   1.165  0.24985   
## Prop_dose_when_Dex65                  -0.3655    20.9654  -0.017  0.98616   
## Prop_dose_when_Dex70                  16.1304    17.4200   0.926  0.35929   
## Prop_dose_when_Dex80                  22.8584    24.2670   0.942  0.35114   
## Prop_dose_when_Dex89.826099999999997  20.4586    20.2538   1.010  0.31773   
## Prop_dose_when_DexNotApp               9.3781    11.7641   0.797  0.42944   
## Mode_CPOT1                            -3.8108    12.6201  -0.302  0.76404   
## Mode_CPOT2                            -0.4011    10.4505  -0.038  0.96955   
## Mode_CPOT3                             9.8659     9.1753   1.075  0.28786   
## Mode_CPOT4                             6.9363     9.3809   0.739  0.46342   
## Mode_CPOT5                            29.6107     8.5984   3.444  0.00123 **
## Mode_CPOT6                            39.9011    12.6201   3.162  0.00278 **
## Mode_CPOT7                                 NA         NA      NA       NA   
## Mode_CPOT8                            -9.5307    17.5907  -0.542  0.59057   
## Mode_CPOTNotApp                      -12.5642    17.5907  -0.714  0.47868   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 300.3285)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 13815  on 46  degrees of freedom
## AIC: 636.82
## 
## Number of Fisher Scoring iterations: 2

Running diagnostics and checking model adequacy

par(mfrow=c(2,2))
plot(glm_mod5)
## Warning: not plotting observations with leverage one:
##   1, 2, 4, 7, 8, 9, 12, 17, 18, 19, 20, 21, 22, 26, 71
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Evaluating evidence of interaction between a quantitative and a qualitative variable

There is evidence of interaction between a quantitative predictor and a qualitative predictor based off the EDA. There seems to be some evidence of interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative). An ANCOVA will be run.

#Starting with ANOVA
Propdata_anova = select(df_eda, Per_prop_over50,Org_Name, Total_days_prop)
glm_anova <- glm(Per_prop_over50~ factor(Org_Name),data=Propdata_anova,family="gaussian")
summary(glm_anova)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name), family = "gaussian", 
##     data = Propdata_anova)
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                  11.048      9.357   1.181   0.2419  
## factor(Org_Name)IU HBloomH    3.072     11.246   0.273   0.7856  
## factor(Org_Name)IU HBMH      19.356      9.814   1.972   0.0527 .
## factor(Org_Name)IU HMH       -4.315     10.462  -0.412   0.6813  
## factor(Org_Name)IU HUH       -2.403     14.293  -0.168   0.8670  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 350.2303)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 23465  on 67  degrees of freedom
## AIC: 632.96
## 
## Number of Fisher Scoring iterations: 2

There is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.

#Running ANCOVA
glm_ancova <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop, data=Propdata_anova,family="gaussian")
summary(glm_ancova)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop, 
##     family = "gaussian", data = Propdata_anova)
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)  
## (Intercept)                11.64445   10.04539   1.159   0.2506  
## factor(Org_Name)IU HBloomH  2.98745   11.33897   0.263   0.7930  
## factor(Org_Name)IU HBMH    19.21253    9.92115   1.937   0.0571 .
## factor(Org_Name)IU HMH     -4.36936   10.54302  -0.414   0.6799  
## factor(Org_Name)IU HUH     -2.58352   14.43647  -0.179   0.8585  
## Total_days_prop            -0.05424    0.31577  -0.172   0.8641  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 355.3779)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 23455  on 66  degrees of freedom
## AIC: 634.93
## 
## Number of Fisher Scoring iterations: 2
anova(glm_ancova)
## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: Per_prop_over50
## 
## Terms added sequentially (first to last)
## 
## 
##                  Df Deviance Resid. Df Resid. Dev
## NULL                                71      31539
## factor(Org_Name)  4   8073.4        67      23465
## Total_days_prop   1     10.5        66      23455

Once again, there is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.

#interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative)
par(mar=c(1,1,1,1))
par(mfrow=c(1,1))
ggplot(data = Propdata_anova, aes(x=Total_days_prop, y=Per_prop_over50, color=factor(Org_Name)))+
  geom_point() + geom_smooth() + xlab("total days on propofol infusion (Total_days_prop)") + ylab("% of time of propofol above 50 mcg/kg/min (Per_prop_over50)")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small.  fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4.93
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 8.07
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 145.68
## Warning in sqrt(sum.squares/one.delta): NaNs produced
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.  fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 4.93
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 8.07
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 145.68
## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1.1299e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 121
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 4
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 1.1299e-17
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 121
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small.  fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 3.97
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 5.03
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1.0609
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.  fewer
## data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 3.97
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 5.03
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0
## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1.0609
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

#testing for interaction
glm_ancova1 <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop +factor(Org_Name)*Total_days_prop, data=Propdata_anova, family="gaussian")
summary(glm_ancova1)
## 
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop + 
##     factor(Org_Name) * Total_days_prop, family = "gaussian", 
##     data = Propdata_anova)
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 20.9456    21.3392   0.982    0.330
## factor(Org_Name)IU HBloomH                   2.7760    23.3274   0.119    0.906
## factor(Org_Name)IU HBMH                      6.7984    21.8265   0.311    0.756
## factor(Org_Name)IU HMH                     -12.6966    23.0936  -0.550    0.584
## factor(Org_Name)IU HUH                     -19.1727    40.0650  -0.479    0.634
## Total_days_prop                             -0.8998     1.7366  -0.518    0.606
## factor(Org_Name)IU HBloomH:Total_days_prop  -0.1169     1.8869  -0.062    0.951
## factor(Org_Name)IU HBMH:Total_days_prop      1.2184     1.7854   0.682    0.498
## factor(Org_Name)IU HMH:Total_days_prop       0.7482     1.8892   0.396    0.693
## factor(Org_Name)IU HUH:Total_days_prop       1.7962     4.5305   0.396    0.693
## 
## (Dispersion parameter for gaussian family taken to be 361.8778)
## 
##     Null deviance: 31539  on 71  degrees of freedom
## Residual deviance: 22436  on 62  degrees of freedom
## AIC: 639.73
## 
## Number of Fisher Scoring iterations: 2
#running model diagnostics
anova(glm_ancova, glm_ancova1)
## Analysis of Deviance Table
## 
## Model 1: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop
## Model 2: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop + factor(Org_Name) * 
##     Total_days_prop
##   Resid. Df Resid. Dev Df Deviance
## 1        66      23455            
## 2        62      22436  4   1018.5

The interaction term is significant as the difference in deviance between the larger model and the smaller model is statistically significant.

Final model with interaction term is as follows:

Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT - 0.518 * Total_days_prop

Summarizing findings and making recommendations

In summary, the final explanatory variables for the % of time of propofol above 50 mcg/kg/min are Prop_dose_when_dex, Mode_CPOT, and Total_days_prop. My recommendation would be for the Ball Memorial Hospital to evaluate their usage of dexmedetomidine, how nurses use the Critical Care Pain Observation Tool, and the number of days patients are kept on propofol.

Reflection on project completion

Completing the section of creating the analysis dataset helped me better understand when it is more convenient to create an analysis dataset by “deleting” columns I do not want to analyze or “selecting” the columns I do want to analyze. For this analysis, I realized it would be better to subtract the columns I do not want to analyze since many of the columns do have some relation to propofol.

Completing the section of conducting an exploratory data analysis (EDA) helped me better understand that the EDA should actually be where you spend a majority of your time when it comes to doing any kind of data modeling. With the number of variables in this dataset, it helped me practice interpreting different kinds of graphs and determine how they related to each other in a more efficient manner.

Completing the section of selecting quantitative information regarding ‘propofol’ and a set of covariates to fit an appropriate general linear model (GLM) helped me to understand the fundamental concepts of GLMs better. This dataset required additional research into what each of the variables mean to know how a GLM would be helpful in this case. Completing the section of performing variable selection helped me to explore deviance tables. It was proving difficult to perform LASSO here as doing a correlation plot was not possible due to the number of categorical/qualitative variables present.

Completing the section of running diagnostics helped me to better understand the different statistics of Residuals vs. Fitted, Q-Q Residuals, Scale-Location, and Residuals vs. Leverage on real-world data. Completing the section of running tests to see if there is evidence of interaction between a quantitative exploratory variable and a qualitative explanatory variable helped me to better understand how to interpret interaction.

Bibliography

References

  1. Physiopedia. (n.d.). Richmond agitation-sedation scale (RASS). Physiopedia. https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS) Description of Richmond Agitation-Sedation Scale (RASS).

  2. South Tees Hospitals NHS FT. (2023, January 20). CPOT = Critical-Care Pain Observation Tool. Critical Care Services - The James Cook University Hospital. https://ccs-sth.org/resources/Documents/Sedation Analgesia Delirium in CC/CCS STH Guideline template CPOT.pdf Description of CPOT = Critical-Care Pain Observation Tool.