Highlights
Information about this assessment:
Introduction
In this assessment we pretend that you are a data analyst at the Department of Health during an outbreak of a novel viral infection (COVID 19). The manager of your unit has asked that you apply your skills to health administrative data to provide information on this novel disease. The results will be used to assist clinicians in treating patients and help public health officials set priorities.
Part 1 - Data management and missing data
The first step in any analysis is to explore the data, check for missing values and clean or reformat as necessary. Let’s start by loading the patient data and storing it in an R object called patient.data. Next, we will remove patients who did not have COVID 19 (all of this is done for you in the code below, run it in R). Don’t forget to set your working directory before running the code below.
patient.data<-read.csv("Synthea_patient_covid.csv", na.strings = "")
patient.data<-patient.data[which(patient.data$covid_status==1),]
Question 1
The first thing you want to know is how many episodes of care ended in the patient’s death (remember a row represents an episode of care not a patient)? Use the table() function to find and report this number (paste your code, output and write a statement such as “In total XX patients died”).
Question 2
Next you want to know how many missing values each variable has. Use the same method and code from week 3 to obtain these results (note the code is complex but all you need to do is change the name of the data object from the one in the week three exercise to the one used here). Underneath your output, list the variables which do have some missing data):
Question 3
One of the variables which has missing data is DEATHDATE. However, this is expected because only patients who died will have a value on this variable. But there could still be genuine missing values on this variable if any patients who died also does not have a value for DEATHDATE.
Make a cross-classification table showing death status by missingness on DEATHDATE and report how many (if any) patients who died were missing a DEATHDATE (remember is.na() will return the values TRUE or FALSE for missingness.
Question 4
Our data does not include a variable for age, but just as we did in class we can make this variable by running the three lines of code given below (run them now):
patient.data$dob<-as.Date(patient.data$BIRTHDATE, tryFormats = "%d/%m/%Y")
patient.data$enc.date<-as.Date(patient.data$DATE, tryFormats = "%d/%m/%Y")
patient.data$age<-as.numeric((patient.data$enc.date-patient.data$dob)/365.25)
Next, use the hist() and summary() functions to plot and summarise the age variable (include your code, plot and output below). Report the mean and median for the age distribution, and report the age group (in groups of ten years) which occurs most frequently in the data.
Question 5
As the pandemic has progressed it has become clear that males are at increased risk of death compared with females. But before investigating if that is true in our data, we firstly need to see if age is distributed differently by gender (i.e., is one group older than the other?). To answer this question, make a box plot of age by gender and interpret the plot with regards to the distribution of age by gender.
Question 6
Your manager also wants to know the COVID death rate (percentage of patients who died) in your sample. Do this calculation now and report the percentage who died (hint: you will find the unique() function to be helpful).
Part 2 - Prediction models
In this section, we will start a new and different activity. The code below firstly removes any data from R, then loads the three Synthea data sets, and loads the rpart() and rpart.plot() packages. NOTE: if you have closed R since the last activity you will need to set your working directory again.
rm(list=ls())
patient.data<-read.csv("Synthea_patient_covid.csv", na.strings = "")
obs.data<-read.csv("Synthea_observations_covid.csv", na.strings = "")
con.data<-read.csv("Synthea_conditions_covid.csv", na.strings = "")
library(rpart)
library(rpart.plot)
Now, you are to build a prediction model of death among patients with COVID 19. This will help clinicians understand which COVID patients are at greatest risk of death. The code below reformats our three data sets into one new data set called pred.data which is ready to use for prediction modelling (use this code to prepare your data):
full.data<-Reduce(merge, list(patient.data,obs.data,con.data))
full.data<-full.data[which(full.data$covid_status==1),]
pred.data<-full.data[,c(5,8, 10:11,30,35,43,45,62,64:66,73,75:76,98,108,
115,117,122,134,146,148,173,176,125,194:195,209,217,218,
221:222,225:227,235,236)]
pred.data$GENDER<-as.factor(pred.data$GENDER)
pred.data$MARITAL<-as.factor(pred.data$MARITAL)
pred.data$ETHNICITY<-as.factor(pred.data$ETHNICITY)
Question 7
Explain what each line of code in the above code does (there are six lines, line three takes up three lines) and why we took these steps to produce our pred.data. Explain it line-by-line (e.g., Line 1 does A, B and C. Line 2 does X, Y and Z, and so on…). In your own words.
This data is very similar to the data we prepared in our week 4 exercise, with a few important differences. Here we are predicting death among the COVID+ sample. We have also removed the RACE variable to save time (i.e., we don’t have to dummy code). Before we can run our model we have one last step, we must randomly separate our data into a training and testing set which is done for you in the code below (run the code in R).
set.seed(36457)
pred.data$random<-runif(nrow(pred.data))
train.data<-pred.data[which(pred.data$random<=.7),]
test.data<-pred.data[which(pred.data$random>.7),]
train.data<-train.data[,-39]
test.data<-test.data[,-39]
Next we can construct a very simple classification tree by running the rpart() function with the default setting as below, saving the model in a new object called tree.mod.
set.seed(73525)
tree.mod<-rpart(formula=dead~.,data=train.data)
Question 8
Plot the tree diagram using the package rpart.plot, and underneath the plot write one statement describing which patients are at greatest risk of dying from COVID 19 and their probability of death? In addition, describe which patients are at the least risk of dying from COVID 19 and their risk of death?
Next, we have to obtain the predicted probability of death for each row in our data using the predict() function. After which we must apply a cut-off value to this probability for which we say patients who are above this value are predicted to die from COVID 19. You decide to use a probability cut-off value of 0.5, so that all patients with a predicted probability above 0.5 will be classified as predicted to die. Run the code below which does these steps for you:
train.data$pred.prob<-predict(tree.mod, train.data)
train.data$pred.out<-0
train.data$pred.out[train.data$pred.prob>=0.5]<-1
Lastly, as you will recall from the week 4 notes, we can calculate the sensitivity and specificity of our model from the cross-classification table (see the week 4 notes and lecture if you have forgotten how to interpret this).
table(predicted=train.data$pred.out,actual=train.data$dead)
## actual
## predicted 0 1
## 0 8301 68
## 1 32 142
Remember the sensitivity is simply the number of deceased patients who were correctly predicted to die (142) divided by the total number who actually did die (142+68). Likewise, the specificity is the number of patients who were correctly predicted to not die (8301) divided by the total number who did not die (8301+32).
This gives us: Sensitivity = 142/(142+68) = 0.68 Specificity = 8301/(8301+32) = 0.99
Question 9
Repeat the above process but this time in the test.data, to obtain the sensitivity and specificity in the test data. By comparing the sensitivity and specificity between the train and test set, what do you conclude with regards to ‘model overfit’? Your first line of code is given below, it calculates the predicted probability of the outcome according to the model in the test set (you don’t need to run the classification tree model again, start by generating the predicted probabilities in the test.data which is done in the line of code below).
test.data$pred.prob<-predict(tree.mod, test.data)
Question 10
You manager is not happy with your model’s low sensitivity. S/he tells you to increase the threshold of the probability, so that only patients with a predicted probability above 0.9 will be classified as predicted to die. Why is this a bad idea in your model? (Note: this is an interpretation question - you don’t need to run/present any R code - but experimenting with this cut-off may help you discover the answer).
Question 11
Your manager accepts your answer above, and now suggests that you instead make a more complex model to improve the sensitivity. As discussed in the week 4 exercise, the cp argument (short for complexity parameter) controls the number of splits in our classification tree, by controlling how big of a reduction in model error is required by each split to accept that split. Large values of cp mean that we require larger reductions in error to accept a new split and will lead to a less complex model. Conversely, small values of cp mean we require smaller reductions in error to accept a new split and will lead to a more complex model. The default value for cp is 0.01.
Run a new classification tree in which you reduce the value of cp by one half the default value and interpret the plot (do not change the default values for minsplit and minbucket as we did in the exercise). Using the same probability thresholds as above (0.5), calculate the sensitivity and specificity for this new model in the test and train data. Explain whether the model is better or worse by considering both changes in model performance (i.e., sensitivity and specificity) and over-fit.
IMPORTANT: firstly you will need to remove the pred.prob and pred.out variables from the train.data and test.data. The plot should work OK if you increase the plot panel window - but if your plot is a mess just include it anyway and I will know from your code if you did it correctly.
Part 3 - Imputation of missing data
In this section we’re going to combine the ideas and skills covered in week 3 and 4 to do something new (i.e., not done in class). We will use rpart to impute missing data, which just means we will use rpart to predict the missing values. Importantly however, here we will use trees to predict a continuous outcome (thus we will use a regression tree not a classification tree). Let’s start by loading and preparing our data, we’ll use the same data as in part two, except we’ll keep the full sample of COVID positive and negative patients (obtain the data by running all the code given below):
rm(list=ls())
library(rpart)
library(rpart.plot)
patient.data<-read.csv("Synthea_patient_covid.csv", na.strings = "")
obs.data<-read.csv("Synthea_observations_covid.csv", na.strings = "")
con.data<-read.csv("Synthea_conditions_covid.csv", na.strings = "")
full.data<-Reduce(merge, list(patient.data,obs.data,con.data))
full.data<-full.data[,c(4,5,8, 10:11,30,35,43,45,62,64:66,73,75:76,98,108,
115,117,122,134,146,148,173,176,125,194:195,209,217,218,
221:222,225:227,235,236)]
full.data$GENDER<-as.factor(full.data$GENDER)
full.data$MARITAL<-as.factor(full.data$MARITAL)
full.data$ETHNICITY<-as.factor(full.data$ETHNICITY)
Now, our goal is to provide an estimate of the increased heart rate in COVID patients compared with non-COVID patients. However, we note that our outcome variable Heart.rate has quite a lot of missing data [confirm this: table(is.na(full.data$Heart.rate))]. Considering how much other patient information we have that is related to Heart.rate we believe we can predict the values for those with missing data (also known as imputation). Furthermore, we can do this using our supervised machine learning technique, regression trees, which may build a good model for us.
Just as you did for part two, we now use rpart to build a prediction model (save it again as tree.mod), but this time you need to predict the continuous variable Heart.rate (also include cp=0.001 in the rpart function to make a more complex model). Further, just like before use predict() to get the predicted Heart.rate according to the model, and save it in a new variable called pred.heart.rate. The good news is that all of this is done for you in the code below (just run it):
set.seed(73525)
tree.mod<-rpart(formula=Heart.rate~.,data=full.data, cp=0.001)
full.data$pred.heart.rate<-predict(tree.mod, full.data)
Question 13
Plot the tree using rpart.plot and explain which group had the lowest heart rate, what was the mean heart rate in this group, and what proportion of the sample did they represent? Remember the outcome here is continuous, not binary, so think carefully about what the numbers in the nodes at the bottom of the plot represent. You may need to make your plot window larger or save the image as a PDF to see it properly.
Question 14
Produce a scatter plot comparing Heart.rate (x-axis) with pred.heart.rate (y-axis) and comment on what you see. Why does the plot look the way it does? Do you think the method we chose to create the imputations (regression trees) was a good one? Why or why not?
Question 15
Last of all, you need to make a final variable impute.heart.rate which includes the actual value of Heart.rate for those without missing values and the predicted value of pred.heart.rate for those with missing values (thus imputing the missing values). After you’ve done this, estimate the increased heart rate among COVID patients using the t.test, once using the original variable Heart.rate (i.e., using listwise deletion) and once using the imputed variable impute.heart.rate.
Comment on any difference between the comparison of heart rate between listwise deletion and regression tree imputation. Why do you think the difference occurred?
This PUBH2005 - Data Anaysis has been solved by our PhD Experts at My Uni Paper. Our Assignment Writing Experts are efficient to provide a fresh solution to this question. We are serving more than 10000+ Student in Australia, UK and US by helping them to score HD in their academics. Our Experts are well trained to follow all marking rubrics and referencing style.
© Copyright 2026 My Uni Papers – Student Hustle Made Hassle Free. All rights reserved.