Survival analysis in R

I wrote a tutorial on survival analysis using R, which is located on my RPubs page. The R Markdown code is located on my GitHub site.

I provide an introduction to survivor and hazard functions, Kaplan-Meier curves, and Cox proportional hazards model.

R tutorials on confounding/interaction and linear regression model - Updates

Last year, I created several tutorials on how to use R for identifying confounding/interaction and visualizing linear regression models. I updated these tutorials recently to address some errors and mistakes. They have received new hyperlinks:

R tutorial on confounding and interactions using the epitool and epiR packages is located on my RPubs page here. The R Markdown code is located on my GitHub page here.

R tutorial on linear regression model is located on my RPubs page here. The R Markdown code is located on my GitHub page here.

Communicating data effectively with data visualization: Part 39 (Heatmaps of COVID-19 deaths)

INTRODUCTION

I wanted to incorporate a heatmap that illustrated the death rates (per 100,000 population) across time in the United States. But I also wanted to show when the coronavirus pandemic 2019 (COVID-19) vaccine was introduced and how it impacted death rates. I thought that a heatmap would do a nice job of illustrating this.

The data visualization by Tynan DeBold and Dov Friedman from the Wall Street Journal has a great visualization on the impact of vaccines for various disease from the measles to smallpox on death rates (See figure below). This heatmap shows the number of measles cases per 100,000 population between 1920 and 2000. Each row represents a state or territory of the United States (U.S.). In this tutorial, we’ll create a similar heatmap for COVID-19 deaths.

Source: Tynan DeBold and Dov Friedman, Wall Street Journal (link)

I set out to create my own heat map with COVID-related death rates using data from the Centers for Disease Prevention and Protection (CDC). The CDC provides a dashboard to visualize the trends in death rates by states and U.S. territories (Compare Trends in COVID-19 Cases and Deaths in the US). However, the data was not compiled in an easy manner. You can only visualize 6 territories at a time. I was able to download all the data and compile this into a single file for this tutorial, which you can download here. Use the file with the *.xlsm extension, which supports macros.

TUTORIAL

Step 1. Download the Excel file with the data. Use the data from the “data” tab. Inspect the data. The columns represent the weekly death rate (7-day average number of deaths per 100,000 population). The rows represents the states and U.S. territories.

Step 2. Use the VBA macro. In a previous article, I explained how to create a heatmap with different gradient levels. We will use a modified version of the macro for this exercise.

This is the VBA macro that we’ll use (link). Don’t be intimidated by this. I’ll go over how to use this code

I start by determining the number of gradient levels for the heatmap. The average death rate was 0.35 per 100,000 population, so I generated 20 levels of gradient (0 to 0.999, 1.0 to 1.999, 2.0 to 2.999, etc). I wanted a “blue” shade for this heatmap, so I had to figure out the RGB scheme for each level. I identified the RGB color scheme using a gradient generator by ColorDesigner. RGB code uses three values to represents the main color on the spectrum (red, green, blue).

Once you have the RGB codes for the gradient levels, you can edit the VBA macro.

In the Developer tab, click on “Visual Basic.” Make sure that the Developer tab is viewable on the Ribbon. If it is not, then you can activate this by going to the File > Options > Customize Ribbon and activate it by entering a check by the Ribbon box.

The Visual Basic interface is a separate window that pops up.

In the “Sub ChangeCellColor()” macro, we’re going to include 20 gradient levels. It’s important to make sure the Range() includes the data that we’re interested in modifying. Since the first cell is in A1 and the last cell is in CZ61, the range is Range("A1:CZ61").

Then we include the 20 gradient levels by changing the Case statement with the corresponding RGB codes. As you modify each Case statement, make sure to change the value ranges for each statement. For example, if you want to apply the RGB code for (18, 123, 141), the death rate range is 1.8 to 1.8999999. You can do this for all the gradient levels.

Here is an example:

    Case 1.8 To 1.8999999
         oCell.Interior.Color = RGB(18, 23, 141)
         oCell.Font.Bold = True
         oCell.Font.Color = RGB(18, 23, 141)
         oCell.Font.Name = "Times New Roman"
         oCell.HorizontalAlignment = xlCenter

Case 1.8 to 1.8999999 denotes the range of the values in each cell (7-day average deaths per 100,000 population).

oCell.Interior.Color = RGB(18, 23, 141) denotes the RGB color scheme for our gradient

oCell.Font.Bold = True denotes that the font is bolded

oCell.Font.Color = RGB(18, 23, 141) denotes that the font color matches the cell color

oCell.Font.Name = "Times New Roman" denotes that the font is Times New Roman

oCell.HorizontalAlignment = xlCenter denotes that the value is aligned in the center

After you’ve adjusted your code, you can execute the macro. To execute the macro, go to the Ribbon and select “Macros.” The Macro window will appear with three macros. Select the “ChangeCellColor” macro and click “Run.” This should execute the macro, and you will notice that the data will start to change color to the corresponding gradient values.

To sort by the last column, select the “SortColumn” macro and click “Run.”

To create white borders around the cells, select on the “WhiteOutlineCells” macro and click “Run.”

Step 4. Final touches. You can select the columns and change width to 2.

Once you have the correct cell sizes, you can start to add labels to the file. I included a line to delineate when the first vaccine was introduced and a line for when the president announced that COVID-19 was a national emergency. I also added labels to the bottom part of the table to indicate dates along the timeline. The rows represented the states and U.S. territories.

CONCLUSIONS

The number of deaths was high early in the pandemic in a few select places in the U.S. As the vaccine is introduced, the number of deaths reached a zenith around December 2020 before falling to low levels in February 2021. Then, the death rate started to increase around the beginning of July 2021. Based on the heatmap, the vaccine may have resulted in a decrease in deaths. But the death rate increased approximately 6 months later in what appears to be the beginning of a seasonal pattern. It is unclear whether the introduction of new variants causes the increased death rate, but there is speculation that it may be a contributor. This heatmap does not generate any claims to what is actually happening; it only provides a visual of the patterns that are reported across each U.S. state and territory.

REFERENCES

I took inspiration from the data visualization by Tynan DeBold and Dov Friedman from the Wall Street Journal.

Date for this exercise came from the CDC (link).  

A previous article on how to create heatmaps is available here.

I used the Gradient Generator by ColorDesigner to find out the RGB values for my gradient levels.

Sample size estimation and Power analysis in R

I wrote a tutorial on how to perform sample size estimations and power analysis using R “pwr” package. These are simple examples that will hopefully lead to more complicated estimations. The tutorial is available on RPubs (link). The R Markdown code is available on my Github site (link).

Logistic regression in R - Part 2 (Goodness of fit tests)

In a previous tutorial, I discussed how to perform logistic regression using R. I wrote a follow-up tutorial on how to conduct goodness of fit tests for logistic regression models in R and posted it on RPubs. The R Markdown code is available on my Github site.

I’ve learned how to assess model fit using Pearson correlations, deviance, and modified Hosmer-Lemeshow Goodness Of Fit (GOF) tests. I think these are important tools when assessing the fit of a logistic regression model. However, I wanted to focus on the HL GOF tests for this tutorial because there are a lot of nuances that I learned and wanted to share.

Additionally, I added the usefulness of visualizing whether the model over- or under-predicts the actual observed data using the calibration plot in R.

Logistic regression in R

I wrote a tutorial on how to construct logistic regression models in R. This tutorial was published on RPubs (link). I go through the use of the glm() command to perform a crude logistic regression model and a multivariable logistic regression model. The data (diabetes.csv) that I used for the motivating example is located here.

The RMarkdown code I used to create the tutorial is located on my GitHub site (link).

Visualizing linear regression models using R - Part 2

I continue my previous blog post on visualizing linear regression models using R (link). Part 2 focuses on using visualization to assess whether the model’s residuals were associated with the predicted values and whether they are normally distributed.

The R Markdown code that I wrote to create this tutorial is located on my GitHub site (link).

You can find the tutorials on my RPubs site:

  • Part 1 - Visualizing linear regression model using R (link)

  • Part 2 - Visualizing linear regression model using R (link)

(NOTE: on 30 January 2022, I updated these tutorials and they can be found in my RPubs page here. The R Markdown code is saved on my GitHub page here.)

Visualizing linear regression models using R - Part 1

I wrote a tutorial on how to visualize linear regression models using R. In the tutorial I used the lm() command and the predict3d package to generate the models and visualize them using R. You can view the RPubs tutorial here. (NOTE: on 30 January 2022, I updated this tutorial and it can be found in my RPubs page here.) I created this tutorial using R Markdown, and the codes are available on my GitHub site (link).

R tutorial on using the epitools package to assess confounding and interaction

I created an R tutorial on using the epitools and epiR packages to assess confounding and interaction. It is located on my RPubs page here. I used R Markdown to create this tutorial and I uploaded my code on my GitHub page here. Here is a figure summarizing the lessons from the tutorial.

Reproduction number—COVID-19

BACKGROUND

As the COVID-19 pandemic, which began in December 2019, continues into its second year, public health measures have been put into place to mitigate its spread. At the time of writing this article, there have been over 4.5 million deaths and over 216 million cases due to COVID-19.[1] Surveillance of COVID-19 remains an important public health measure of understanding the spread and impact. Daily reports such as the John Hopkins COVID-19 dashboard provide end users with visual and statistical information about the surges in cases and deaths associated with COVID-19. However, one measure that is of great interest is the reproduction number or R0.

 

Reproduction number (R0) and effective reproduction number (Rt)

The reproduction number is the number of new cases that is directly caused by exposure to a single case.[2,3] Figure 1 provides a visual explanation of the basic reproduction number. However, the underlying assumption with R0 is that everyone in the population is susceptible to infection. With the introduction of vaccines, the R0 isn’t a good measure of the reproductive capabilities of COVID-19. Instead, the effective reproduction number (Rt) is used to provide a more realistic reproduction number based on the population being infected, recovered, or vaccinated. The Rt changes over time as the population susceptible to infection changes.

Figure 1. Basic reproduction number.

I wanted to create a figure that would highlight the changes associated with the Rt for each state in the United States. To do this, I downloaded the Rt data from the by Xihong Lin's Group in the Department of Biostatistics at the Harvard T.H. Chan School of Public Health. They have an amazing COVID-19 tracker dashboard that captures the changing patterns of Rt for each state. Then I created a Cleveland plot to show where the Rt was near the beginning of the pandemic and where it is currently (August 2021). (Note: I wrote a tutorial on creating Cleveland plots that you can review here.) Here is the final figure (because of the length of the figure, I cropped it to show the first 30 states or territories):

 

Figure 2. Effective reproduction number (Rt) for U.S. states and territories, April 17, 2020 (past) to August 14, 2021 (recent).

The blue dots denote the most recent effective reproduction number (14 August 2021) and the past dots denote the earliest effective reproduction number (17 April 2020).

It seems that some states have gotten worse in terms of increase effective reproduction number since the beginning of the pandemic. This could be due to lack of good data in the early phases of the pandemic. However, what is of concern is the high effective reproduction numbers in some states (Rt > 2), which indicates that the pandemic is still spreading at an alarming rate.

There were some missing data which are identified by a single dot (blue or red) or an empty field in the recent or past effective reproduction number. Rather than fill these in, I left them empty. There may be data in between the two time periods that I could have used, but I left those out.

One thing to mention is that this Cleveland plot only tells us one dimension of the effective reproduction number story (the difference between the most recent Rt and the earliest Rt). It doesn’t tell us much about how the effective reproduction number changes across time. For that, I direct your attention to the Lin’s Laboratory Group at Harvard, they have a great figure that shows the fluctuation of the effective reproduction number for the U.S. and its states/territories (see example):

Source: Lin’s Laboratory Group at Harvard (link). [last accessed on 30 August 2021].

CONCLUSIONS

The effective reproduction number provides us with some interesting patterns in spread of COVID-19 by states/territories. It seems to have worsened over time, but this could be due to poor data early in the pandemic. There are some issues with the us of effective reproduction number for policy decisions. Reporting delays can impact the estimates for the effective reproduction number. A technique called “nowcasting” is used to estimate the reproduction number.[3] But when I explored some of the work in this area, there appears to be a variety of methods for performing this technique. Despite this limitation, the effective reproduction number may be useful to evaluate public health policy decisions to reduce the spread of the COVID-19 pandemic.[4,5]

 

DATA SOURCE

I provided the link to the COVID-19 Spread Tracker from the Lin Lab at Harvard. You can also download a curated version of the data for this article from my Dropbox folder. The data are current as of 17 August 2021. If you’re interested in recreating this Cleveland plot, I recommend downloading the most recent data to see how much the effective reproduction number has changed.

REFERENCES

  1. Worldometeres.info. COVID Live Update: 217,770,381 Cases and 4,521,936 Deaths from the Coronavirus - Worldometer. Accessed August 30, 2021. https://www.worldometers.info/coronavirus/

  2. Lim J-S, Cho S-I, Ryu S, Pak S-I. Interpretation of the Basic and Effective Reproduction Number. J Prev Med Pub Health. 2020;53(6):405-408. doi:10.3961/jpmph.20.288

  3. Adam D. A guide to R — the pandemic’s misunderstood metric. Nature. 2020;583(7816):346-348. doi:10.1038/d41586-020-02009-w

  4. Inglesby TV. Public Health Measures and the Reproduction Number of SARS-CoV-2. JAMA. 2020;323(21):2186-2187. doi:10.1001/jama.2020.7878

  5. Pan A, Liu L, Wang C, et al. Association of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan, China. JAMA. 2020;323(19):1915-1923. doi:10.1001/jama.2020.6130