Section 11.3: Assessing trend in cancer rates in the US

University of Southampton
S.K.Sahu@soton.ac.uk

Abstract

This file contains all the code for reproducing the figures and tables for assessing trend in US cancer rates as discussed in Section 11.3.

## Note the start time
start.time<-proc.time()[3]
# figurepath <- "YourPath"   # Folder to save or read the figures from
# filepath <- "YourFilePath" # Folder to save or read the tables from 
figurepath <- "figures/chap11figures/usa/"
filepath <- "txttables/"
dpath <- "datafiles/" # Folder containing the data files 
longrun <- FALSE  # Should we run the model to check the results. 
# If FALSE we use the results, e.g. tables, figures and model fits from previous runs  
if (!longrun) load(file="mcmc_output_files/uscancer/m422.RData")
load(file=paste0(dpath, "C11us48cancerdata.RData"))
colpalette <- c("dodgerblue4",  "dodgerblue2",  "firebrick2",   "firebrick4",   "purple")     

library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(doBy)
library(usmap)

1 Code to reproduce Figures in Section 11.4

head(us48cancertot)
##          Name state fips  totcount   totpop  totrate   AIAN.HL  AIAN.NHL
## 1     Alabama    AL    1  25611.33  4736261 540.1333 0.1272515 0.5892589
## 2     Arizona    AZ    4  29450.87  6363790 462.0467 1.3033103 4.2843893
## 3    Arkansas    AR    5  15938.33  2894066 550.0467 0.1995402 0.7781714
## 4  California    CA    6 165372.00 37325590 442.7733 1.2801805 0.5671331
## 5    Colorado    CO    8  22118.60  5047750 437.6800 0.9988665 0.7593700
## 6 Connecticut    CT    9  21617.73  3555387 607.9467 0.2868551 0.2620113
##       API.HL   API.NHL    BAA.HL   BAA.NHL  White.HL White.NHL
## 1 0.08855458  1.264443 0.2878583 26.567691  3.073277  68.00167
## 2 0.37435044  3.157346 0.6966057  4.198725 27.147177  58.83810
## 3 0.08576009  1.550466 0.2291451 15.805733  5.679722  75.67146
## 4 0.81832412 14.032092 0.8891200  6.373004 34.193313  41.84683
## 5 0.25679949  3.199021 0.5706648  4.324408 18.605895  71.28497
## 6 0.18256169  4.045786 1.5652350 10.071039 11.344880  72.24163
plot_usmap(data = us48cancertot, values = "totrate", color = "red", exclude=c("AK", "HI")) +
  # scale_fill_continuous(name = "Cancer rate", label = scales::comma) +
  scale_fill_gradientn(colours=colpalette,na.value="black",limits=range(us48cancertot$totrate),
  name = "Cancer rate") +
  theme(legend.position = "right")

ggsave(paste0(figurepath, "uscancer_rate.png"),  width=8.27, height=3.44, dpi=300)


us48cancertot$white <- us48cancertot$White.HL + us48cancertot$White.NHL
com1 <- c("gray15","gray17","gray19","gray21",  "white")#colour palette

plot_usmap(data = us48cancertot, values = "White.NHL", color = "red", exclude=c("AK", "HI")) +
  # scale_fill_continuous(name = "Cancer rate", label = scales::comma) +
  scale_fill_gradientn(colours=com1,na.value="black",limits=range(us48cancertot$White.NHL),
                       name = "Percentage of White") +
  theme(legend.position = "right")

ggsave(paste0(figurepath, "percentage_of_white.png"),  width=8.27, height=3.44, dpi=300)

cor(us48cancertot$totrate, us48cancertot$white)
## [1] 0.1255233

cor(us48cancertot$totrate, us48cancertot$White.NHL)
## [1] 0.4569091

plot_usmap(data = us48cancertot, values = "BAA.NHL", color = "red", exclude=c("AK", "HI")) +
  # scale_fill_continuous(name = "Cancer rate", label = scales::comma) +
  scale_fill_gradientn(colours=rev(com1),na.value="black",limits=range(us48cancertot$BAA.NHL),
                       name = "Percentage") +
  theme(legend.position = "right")

ggsave(paste0(figurepath, "baa_nhl.png"),  width=8.27, height=3.44, dpi=300)


## Now the plots
# Plot the observed smr for all the years


head(us48cancertot)
##          Name state fips  totcount   totpop  totrate   AIAN.HL  AIAN.NHL
## 1     Alabama    AL    1  25611.33  4736261 540.1333 0.1272515 0.5892589
## 2     Arizona    AZ    4  29450.87  6363790 462.0467 1.3033103 4.2843893
## 3    Arkansas    AR    5  15938.33  2894066 550.0467 0.1995402 0.7781714
## 4  California    CA    6 165372.00 37325590 442.7733 1.2801805 0.5671331
## 5    Colorado    CO    8  22118.60  5047750 437.6800 0.9988665 0.7593700
## 6 Connecticut    CT    9  21617.73  3555387 607.9467 0.2868551 0.2620113
##       API.HL   API.NHL    BAA.HL   BAA.NHL  White.HL White.NHL    white
## 1 0.08855458  1.264443 0.2878583 26.567691  3.073277  68.00167 71.07494
## 2 0.37435044  3.157346 0.6966057  4.198725 27.147177  58.83810 85.98527
## 3 0.08576009  1.550466 0.2291451 15.805733  5.679722  75.67146 81.35118
## 4 0.81832412 14.032092 0.8891200  6.373004 34.193313  41.84683 76.04015
## 5 0.25679949  3.199021 0.5706648  4.324408 18.605895  71.28497 89.89087
## 6 0.18256169  4.045786 1.5652350 10.071039 11.344880  72.24163 83.58651
us48cancertot$totcount <- round(us48cancertot$totcount)
a1 <- sum(us48cancertot$totcount)
a2 <- sum(us48cancertot$totpop)
rr <- a1/a2
us48cancertot$exptot <- rr * us48cancertot$totpop
summary(us48cancertot$exptot)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2926   10453   24135   33637   37189  197084
colnames(us48cancertot)
##  [1] "Name"      "state"     "fips"      "totcount"  "totpop"    "totrate"  
##  [7] "AIAN.HL"   "AIAN.NHL"  "API.HL"    "API.NHL"   "BAA.HL"    "BAA.NHL"  
## [13] "White.HL"  "White.NHL" "white"     "exptot"

us48cancertot$logexptot <- log(us48cancertot$exptot)


head(us48cancer0317)
##   fips state Year totcount  totpop totrate malescount malespop malesrate
## 1    1    AL 2003    21472 4503491   476.8      11228  2179422     515.2
## 2    1    AL 2004    22951 4530729   506.6      12024  2192872     548.3
## 3    1    AL 2005    23342 4565917   511.2      12172  2211403     550.4
## 4    1    AL 2006    24113 4628981   520.9      12798  2243501     570.4
## 5    1    AL 2007    25086 4672840   536.8      13236  2265565     584.2
## 6    1    AL 2008    26093 4718206   553.0      13708  2287949     599.1
##   femalescount femalespop femalesrate    AIAN.HL  AIAN.NHL     API.HL   API.NHL
## 1        10244    2324069       440.8 0.06708129 0.5561019 0.05142677 0.9124255
## 2        10927    2337857       467.4 0.07808898 0.5614108 0.06045385 0.9714331
## 3        11170    2354514       474.4 0.08941301 0.5661292 0.06757400 1.0243107
## 4        11315    2385480       474.3 0.10133980 0.5719185 0.07634510 1.0712509
## 5        11850    2407275       492.3 0.11241558 0.5782351 0.08480924 1.1150392
## 6        12385    2430257       509.6 0.12699742 0.5833997 0.09126350 1.1586819
##      BAA.HL  BAA.NHL White.HL White.NHL unemployment
## 1 0.2035754 26.22443 1.975512  70.00944          6.0
## 2 0.2141819 26.25379 2.147623  69.71302          5.7
## 3 0.2288719 26.26904 2.345592  69.40907          4.5
## 4 0.2426236 26.34064 2.557496  69.03839          4.0
## 5 0.2551981 26.36388 2.771569  68.71885          4.0
## 6 0.2730275 26.41459 2.978441  68.37359          5.7
summary(us48cancer0317$totcount)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2370    9755   24241   33637   38868  176983
summary(us48cancer0317$unemployment)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.600   4.500   5.400   5.921   7.000  13.700
a1 <- sum(us48cancer0317$totcount)
a2 <- sum(us48cancer0317$totpop)
rr <- a1/a2
us48cancer0317$exptot <- rr * us48cancer0317$totpop
summary(us48cancer0317$exptot)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2658   10061   23959   33637   38976  208034


mus48cancer0317 <- us48cancer0317
mus48cancer0317$obs_smr <- mus48cancer0317$totcount/mus48cancer0317$exptot

mr <- range(mus48cancer0317$obs_smr)
mr
## [1] 0.6187699 1.3236774
summary(mus48cancer0317$obs_smr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6188  0.9578  1.0312  1.0251  1.1107  1.3237

b <- c(-Inf, 0.8, 1.0, 1.2, Inf)
labss <- c("<0.8", "0.8-1", "1-1.2", ">1.2")

mus48cancer0317$catsmrobs <- factor(cut(mus48cancer0317$obs_smr, breaks=b,  labels=labss))
table(mus48cancer0317$catsmrobs)
## 
##  <0.8 0.8-1 1-1.2  >1.2 
##    30   238   397    55
levels(mus48cancer0317$catsmrfit)
## NULL
table(mus48cancer0317$catsmrfit)
## < table of extent 0 >
head(mus48cancer0317)
##   fips state Year totcount  totpop totrate malescount malespop malesrate
## 1    1    AL 2003    21472 4503491   476.8      11228  2179422     515.2
## 2    1    AL 2004    22951 4530729   506.6      12024  2192872     548.3
## 3    1    AL 2005    23342 4565917   511.2      12172  2211403     550.4
## 4    1    AL 2006    24113 4628981   520.9      12798  2243501     570.4
## 5    1    AL 2007    25086 4672840   536.8      13236  2265565     584.2
## 6    1    AL 2008    26093 4718206   553.0      13708  2287949     599.1
##   femalescount femalespop femalesrate    AIAN.HL  AIAN.NHL     API.HL   API.NHL
## 1        10244    2324069       440.8 0.06708129 0.5561019 0.05142677 0.9124255
## 2        10927    2337857       467.4 0.07808898 0.5614108 0.06045385 0.9714331
## 3        11170    2354514       474.4 0.08941301 0.5661292 0.06757400 1.0243107
## 4        11315    2385480       474.3 0.10133980 0.5719185 0.07634510 1.0712509
## 5        11850    2407275       492.3 0.11241558 0.5782351 0.08480924 1.1150392
## 6        12385    2430257       509.6 0.12699742 0.5833997 0.09126350 1.1586819
##      BAA.HL  BAA.NHL White.HL White.NHL unemployment   exptot   obs_smr
## 1 0.2035754 26.22443 1.975512  70.00944          6.0 23779.00 0.9029816
## 2 0.2141819 26.25379 2.147623  69.71302          5.7 23922.82 0.9593769
## 3 0.2288719 26.26904 2.345592  69.40907          4.5 24108.62 0.9682015
## 4 0.2426236 26.34064 2.557496  69.03839          4.0 24441.60 0.9865556
## 5 0.2551981 26.36388 2.771569  68.71885          4.0 24673.18 1.0167314
## 6 0.2730275 26.41459 2.978441  68.37359          5.7 24912.72 1.0473765
##   catsmrobs
## 1     0.8-1
## 2     0.8-1
## 3     0.8-1
## 4     0.8-1
## 5     1-1.2
## 6     1-1.2
dim(mus48cancer0317)
## [1] 720  24
colnames(mus48cancer0317)
##  [1] "fips"         "state"        "Year"         "totcount"     "totpop"      
##  [6] "totrate"      "malescount"   "malespop"     "malesrate"    "femalescount"
## [11] "femalespop"   "femalesrate"  "AIAN.HL"      "AIAN.NHL"     "API.HL"      
## [16] "API.NHL"      "BAA.HL"       "BAA.NHL"      "White.HL"     "White.NHL"   
## [21] "unemployment" "exptot"       "obs_smr"      "catsmrobs"
com <- c("green4", "green2", "yellow", "red2")

i <- 2003
colnames(mus48cancer0317)
##  [1] "fips"         "state"        "Year"         "totcount"     "totpop"      
##  [6] "totrate"      "malescount"   "malespop"     "malesrate"    "femalescount"
## [11] "femalespop"   "femalesrate"  "AIAN.HL"      "AIAN.NHL"     "API.HL"      
## [16] "API.NHL"      "BAA.HL"       "BAA.NHL"      "White.HL"     "White.NHL"   
## [21] "unemployment" "exptot"       "obs_smr"      "catsmrobs"
dim(mus48cancer0317)
## [1] 720  24
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
head(dp)
##    fips state catsmrobs
## 1     1    AL     0.8-1
## 16   10    DE     1-1.2
## 31   12    FL     1-1.2
## 46   13    GA     0.8-1
## 61   16    ID     0.8-1
## 76   17    IL     0.8-1
smr2003 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  #  scale_fill_gradientn(colours=com,na.value="black",limits=mr, name = "smr") +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2003

i <- 2004
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2004 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2004

i <- 2005
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2005 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2005

i <- 2006
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2006 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2006

i <- 2007
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2007 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2007

i <- 2008
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2008 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2008

i <- 2009
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2009 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2009

i <- 2010
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2010 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2010

i <- 2011
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2011 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2011

i <- 2012
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2012 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2012

i <- 2013
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2013 <- plot_usmap(data = dp, values = "catsmrobs", color = "red", exclude=c("AK", "HI")) +
  scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
  theme(legend.position = "right")  +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title= i)
smr2013