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
i <- 2014
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2014 <- 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)
smr2014
i <- 2015
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2015 <- 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)
smr2015
i <- 2016
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2016 <- 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)
smr2016
i <- 2017
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrobs")]
#head(dp)
smr2017 <- 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)
smr2017
library(ggpubr)
# a <- get_legend(smr2003, position = "bottom")
a <- get_legend(smr2017, position = c(0.25, 0.1))
ggarrange(
smr2003, smr2004, smr2005, smr2006, smr2007,
smr2008, smr2009, smr2010, smr2011, smr2012,
smr2013, smr2014, smr2015, smr2016, smr2017, as_ggplot(a),
common.legend = TRUE, legend = "none", nrow = 4, ncol = 4
)
ggsave(filename=paste0(figurepath, "uscancer_obs_smr.png"), device = "png")
## Saving 7 x 5 in image
### Raw SMR's
p2 <- ggplot(data=us48cancer0317, aes(x=Year, y=unemployment, color=state, group=state)) +
# geom_boxplot() + facet_wrap(~year, ncol=5) +
geom_line() +
labs(x="Year", y="Annual unemployment") +
scale_colour_discrete(guide="none")
# geom_dl(data=dp2, aes(x=Year, y=unemployment, label =state), method = list(dl.combine("last.bumpup")), cex = 0.8)
p2
ggsave(filename=paste0(figurepath, "unemployment.png"))
## Saving 7 x 5 in image
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 exptot
## 1 0.2035754 26.22443 1.975512 70.00944 6.0 23779.00
## 2 0.2141819 26.25379 2.147623 69.71302 5.7 23922.82
## 3 0.2288719 26.26904 2.345592 69.40907 4.5 24108.62
## 4 0.2426236 26.34064 2.557496 69.03839 4.0 24441.60
## 5 0.2551981 26.36388 2.771569 68.71885 4.0 24673.18
## 6 0.2730275 26.41459 2.978441 68.37359 5.7 24912.72
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
us48cancer0317$sids <- rep(1:48, each=15)
mus48cancer0317 <- us48cancer0317
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" "sids"
for (i in 13:21) {
mus48cancer0317[, i] <- (mus48cancer0317[, i]- mean(mus48cancer0317[, i]))/sd(mus48cancer0317[, i])
}
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
## 1 10244 2324069 440.8 -0.8951115 -0.3829786 -0.6807173
## 2 10927 2337857 467.4 -0.8663559 -0.3807052 -0.6205997
## 3 11170 2354514 474.4 -0.8367739 -0.3786847 -0.5731817
## 4 11315 2385480 474.3 -0.8056172 -0.3762056 -0.5147688
## 5 11850 2407275 492.3 -0.7766837 -0.3735007 -0.4584002
## 6 12385 2430257 509.6 -0.7385912 -0.3712892 -0.4154168
## API.NHL BAA.HL BAA.NHL White.HL White.NHL unemployment exptot
## 1 -0.8989287 -0.6490985 1.603981 -0.7959856 -0.2662471 0.03851875 23779.00
## 2 -0.8764698 -0.6312760 1.607064 -0.7776644 -0.2875463 -0.10821935 23922.82
## 3 -0.8563441 -0.6065917 1.608667 -0.7565908 -0.3093866 -0.69517173 24108.62
## 4 -0.8384781 -0.5834842 1.616187 -0.7340337 -0.3360211 -0.93973522 24441.60
## 5 -0.8218118 -0.5623547 1.618628 -0.7112457 -0.3589813 -0.93973522 24673.18
## 6 -0.8052010 -0.5323952 1.623955 -0.6892243 -0.3837895 -0.10821935 24912.72
## sids
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
yt <- 2003:2017
mus48cancer0317$tyear <- (mus48cancer0317$Year - mean(yt))/sd(yt)
us48cancer0317$tyear <- (mus48cancer0317$Year - mean(yt))/sd(yt)
f0 <- totcount ~ offset(log(exptot)) + tyear + unemployment
f22 <- update(f0, . ~ . + AIAN.HL + API.HL + BAA.HL)
if (longrun) {
burn.in <- 5000
N <- 205000
thin <- 20
m422 <- Bcartime(formula = f22, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol="sids", tcol="Year", N=N,
burn.in = burn.in, thin=thin, verbose=FALSE)
summary(m422)
table11.2 <- m422$params[, 1:3]
round(table11.2, 3)
dput(table11.2, file=paste0(filepath, "table11.2.txt"))
save(m422, file="uscancer/m422.RData")
} else {
table11.2 <- dget(file=paste0(filepath, "table11.2.txt"))
}
table11.2
## Median 2.5% 97.5%
## (Intercept) 0.0167 0.0162 0.0172
## tyear 0.0404 0.0335 0.0507
## unemployment 0.0082 0.0014 0.0191
## AIAN.HL -0.0643 -0.0754 -0.0558
## API.HL -0.0327 -0.0414 -0.0221
## BAA.HL 0.0267 0.0156 0.0448
## tau2.S 0.0418 0.0291 0.0637
## tau2.T 0.1135 0.0640 0.2324
## tau2.I 0.0092 0.0082 0.0103
## rho.S 0.8998 0.6787 0.9901
## rho.T 0.5078 0.0620 0.9131
# Inference with m422 instead
# load(file="uscancer/m422.RData")
modfit <- m422$fit
names(modfit)
## [1] "summary.results" "samples" "fitted.values"
## [4] "residuals" "modelfit" "accept"
## [7] "localised.structure" "formula" "model"
## [10] "X"
u <- modfit$fitted.values
length(u)
## [1] 720
mus48cancer0317$fits <- modfit$fitted.values
mus48cancer0317$smrfit <- mus48cancer0317$fits/mus48cancer0317$exptot
#mus48cancer0317$obs_smr <- mus48cancer0317$totcount/mus48cancer0317$exptot
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
## 1 10244 2324069 440.8 -0.8951115 -0.3829786 -0.6807173
## 2 10927 2337857 467.4 -0.8663559 -0.3807052 -0.6205997
## 3 11170 2354514 474.4 -0.8367739 -0.3786847 -0.5731817
## 4 11315 2385480 474.3 -0.8056172 -0.3762056 -0.5147688
## 5 11850 2407275 492.3 -0.7766837 -0.3735007 -0.4584002
## 6 12385 2430257 509.6 -0.7385912 -0.3712892 -0.4154168
## API.NHL BAA.HL BAA.NHL White.HL White.NHL unemployment exptot
## 1 -0.8989287 -0.6490985 1.603981 -0.7959856 -0.2662471 0.03851875 23779.00
## 2 -0.8764698 -0.6312760 1.607064 -0.7776644 -0.2875463 -0.10821935 23922.82
## 3 -0.8563441 -0.6065917 1.608667 -0.7565908 -0.3093866 -0.69517173 24108.62
## 4 -0.8384781 -0.5834842 1.616187 -0.7340337 -0.3360211 -0.93973522 24441.60
## 5 -0.8218118 -0.5623547 1.618628 -0.7112457 -0.3589813 -0.93973522 24673.18
## 6 -0.8052010 -0.5323952 1.623955 -0.6892243 -0.3837895 -0.10821935 24912.72
## sids tyear fits smrfit
## 1 1 -1.5652476 21485.16 0.9035351
## 2 1 -1.3416408 22959.94 0.9597506
## 3 1 -1.1180340 23350.05 0.9685354
## 4 1 -0.8944272 24120.83 0.9868760
## 5 1 -0.6708204 25090.72 1.0169226
## 6 1 -0.4472136 26095.15 1.0474628
b <- c(-Inf, 0.8, 1.0, 1.2, Inf)
labss <- c("<0.8", "0.8-1", "1-1.2", ">1.2")
mus48cancer0317$catsmrfit <- factor(cut(mus48cancer0317$smrfit, breaks=b, labels=labss))
table(mus48cancer0317$catsmrfit)
##
## <0.8 0.8-1 1-1.2 >1.2
## 30 237 400 53
levels(mus48cancer0317$catsmrfit)
## [1] "<0.8" "0.8-1" "1-1.2" ">1.2"
com <- c("green4", "green2", "yellow", "red2")
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" "sids" "tyear" "fits"
## [26] "smrfit" "catsmrfit"
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" "sids" "tyear" "fits"
## [26] "smrfit" "catsmrfit"
i <- 2003
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
head(dp)
## fips state catsmrfit
## 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 = "catsmrfit", color = "red", exclude=c("AK", "HI")) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr, name = "smr") +
scale_fill_manual(values =com, guides("SMR"), guide = guide_legend(reverse=TRUE), na.translate = FALSE) +
# theme(legend.position = "right") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
labs(title= i)
smr2003
i <- 2003
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
head(dp)
## fips state catsmrfit
## 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 = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2004 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2005 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2006 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2007 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2008 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2009 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2010 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2011 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2012 <- plot_usmap(data = dp, values = "catsmrfit", 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", "catsmrfit")]
#head(dp)
smr2013 <- plot_usmap(data = dp, values = "catsmrfit", 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
i <- 2014
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
#head(dp)
smr2014 <- plot_usmap(data = dp, values = "catsmrfit", 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)
smr2014
i <- 2015
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
#head(dp)
smr2015 <- plot_usmap(data = dp, values = "catsmrfit", 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)
smr2015
i <- 2016
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
#head(dp)
smr2016 <- plot_usmap(data = dp, values = "catsmrfit", 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)
smr2016
i <- 2017
dp <- mus48cancer0317[mus48cancer0317$Year==i, c("fips", "state", "catsmrfit")]
#head(dp)
smr2017 <- plot_usmap(data = dp, values = "catsmrfit", 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)
smr2017
library(ggpubr)
# a <- get_legend(smr2003, position = "bottom")
a <- get_legend(smr2017, position = c(0.25, 0.1))
ggarrange(
smr2003, smr2004, smr2005, smr2006, smr2007,
smr2008, smr2009, smr2010, smr2011, smr2012,
smr2013, smr2014, smr2015, smr2016, smr2017, as_ggplot(a),
common.legend = TRUE, legend = "none", nrow = 4, ncol = 4
)
ggsave(filename=paste0(figurepath, "uscancerfits.png"), device = "png")
## Saving 7 x 5 in image
Not reported models
head(us48cancertot)
us48cancertot$totcount <- round(us48cancertot$totcount)
a1 <- sum(us48cancertot$totcount)
a2 <- sum(us48cancertot$totpop)
rr <- a1/a2
us48cancertot$exptot <- rr * us48cancertot$totpop
summary(us48cancertot$exptot)
colnames(us48cancertot)
us48cancertot$logexptot <- log(us48cancertot$exptot)
burn.in <- 5000
N <- 25000
thin <- 10
family <- "poisson"
scol <- "sids"
tcol <- "Year"
link <- "log"
us48cancertot$sids <- 1:48
f2 <- totcount ~ offset(log(exptot)) + BAA.HL + BAA.NHL + White.HL + White.NHL
m0 <- Bcartime(formula = f2, data =us48cancertot, family ="poisson",
N=N, burn.in = burn.in, thin=thin)
summary(m0)
# for INLA
f2 <- totcount ~ BAA.HL + BAA.NHL + White.HL + White.NHL
model <- c("bym")
colnames(us48cancertot)
d1 <- Bcartime(data=us48cancertot, formula=f2, W=W, scol =scol,
offsetcol="logexptot", model=model, link=link, family=family, package="inla")
summary(d1)
rownames(d1$fit$summary.hyperpar)
## Analysis of the cancer totals spatio-temporal
us48cancer0317$sids <- rep(1:48, each=15) ## Needs a proper scol argument
scol <- "sids"
colnames(us48cancer0317)
unique(us48cancer0317$fips)
mus48cancer0317 <- us48cancer0317
colnames(mus48cancer0317)
for (i in 13:21) {
mus48cancer0317[, i] <- (mus48cancer0317[, i]- mean(mus48cancer0317[, i]))/sd(mus48cancer0317[, i])
}
head( mus48cancer0317)
yt <- 2003:2017
mus48cancer0317$tyear <- (mus48cancer0317$Year - mean(yt))/sd(yt)
us48cancer0317$tyear <- (mus48cancer0317$Year - mean(yt))/sd(yt)
f1 <- totcount ~ offset(log(exptot)) + tyear + AIAN.NHL + API.NHL + BAA.NHL + White.NHL+ unemployment
m01 <- Bcartime(formula = f1, data =mus48cancer0317, family ="poisson",
N=N, burn.in = burn.in, thin=thin)
summary(m01)
f2 <- totcount ~ offset(log(exptot)) + tyear+ AIAN.HL + API.HL + BAA.HL + White.HL+ unemployment
m02 <- Bcartime(formula = f2, data =mus48cancer0317, family ="poisson",
N=N, burn.in = burn.in, thin=thin)
summary(m02)
round(m01$mchoice, 2)
round(m02$mchoice, 2)
# f2 <- totcount ~ offset(log(exptot)) + BAA.HL + BAA.NHL + White.HL + White.NHL+ unemployment
m01 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family = family,
link= link, N=N, burn.in = burn.in, thin=thin)
summary(m01)
m2 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family = family,
model="linear", package="CARBayesST", link= link, scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m2)
m3 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family = family,
model="sepspatial", package="CARBayesST",
link= link, scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m3)
m4 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m4)
m5 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family = family,
model="ar", package="CARBayesST",
link= link, scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m5)
m51 <- Bcartime(formula = f1, data =mus48cancer0317, W = W, family = family,
model="ar", package="CARBayesST",
link= link, scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m51)
f3 <- totcount ~ offset(log(exptot)) + tyear+ AIAN.HL + API.HL + BAA.HL + BAA.NHL + White.NHL + White.HL+ unemployment
m53 <- Bcartime(formula = f3, data =mus48cancer0317, W = W, family = "poisson",
model="ar", package="CARBayesST",
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m53)
v <- us48cancer0317[, 13:20]
head(v)
a <- apply(v, 1, sum)
library(GGally)
ggpairs(v)
round(cor(v), 2)
v1 <- v[, c(2, 4, 6, 8)]
head(v1)
round(cor(v1), 2)
v2 <- v[, c(1, 3, 5, 7)]
head(v2)
round(cor(v2), 2)
u <- v[, c(1, 3, 5)]
cor(u)
## INLA did not work
mus48cancer0317$logexptot <- log(mus48cancer0317$exptot)
f2 <- totcount ~ BAA.HL + BAA.NHL + White.HL + White.NHL
model <- c("iid", "ar1")
# model <- c("bym", "ar1")
d1 <- Bcartime(data=mus48cancer0317, formula=f2, W=W, scol =scol, tcol=tcol,
offsetcol="logexptot", model=model, link=link, family=family, package="inla")
summary(d1)
rownames(d1$fit$summary.hyperpar)
## Try the Anova model only
f0 <- totcount ~ offset(log(exptot)) + tyear + unemployment
f1 <- update(f0, . ~ . + AIAN.NHL + API.NHL + BAA.NHL + White.NHL)
f2 <- update(f0, . ~ . + AIAN.HL + API.HL + BAA.HL + White.HL)
f3 <- update(f0, . ~ . + BAA.HL + BAA.NHL + White.NHL + White.HL)
f22 <- update(f0, . ~ . + AIAN.HL + API.HL + BAA.HL)
m40 <- Bcartime(formula = f0, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m40)
m41 <- Bcartime(formula = f1, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m41)
m42 <- Bcartime(formula = f2, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m42)
m422 <- Bcartime(formula = f22, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m422)
m43 <- Bcartime(formula = f3, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m43)
m40$mchoice
m41$mchoice
m42$mchoice
m43$mchoice
m42$mchoice
m422$mchoice
f4 <- totcount ~ offset(log(exptot)) + unemployment + AIAN.HL + API.HL + BAA.HL + White.HL
m44 <- Bcartime(formula = f4, data =mus48cancer0317, W = W, family ="poisson",
model="anova", package="CARBayesST", interaction=T,
scol=scol, tcol=tcol, N=N,
burn.in = burn.in, thin=thin)
summary(m44)
# m422 is the chosen model
## Try INLA
f22
u <- m422$params[, 1:3]
u
library(xtable)
xtable(u, digits=3)
f22inla <- totcount ~ tyear + unemployment + AIAN.HL + API.HL + BAA.HL
model <- c("bym", "ar1")
colnames(mus48cancer0317)
mus48cancer0317$logexptot <- log(mus48cancer0317$exptot)
d22 <- Bcartime(data=mus48cancer0317, formula=f22inla, W=W, scol =scol, tcol=tcol,
offsetcol="logexptot", model=model, link="log", family="poisson", package="inla")
summary(d22)
rownames(d22$fit$summary.hyperpar)
## Bursted parameter estimates
# All done
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
## elapsed
## 0.3530667